1	/*//////////////////////////////////////////////////////////////////////////////
     2	// <file
     3	//   repository_path  = "$Source: /pub/cvsprojects/nucleo/web/home/modules/libnucnet/src/examples/example36.c,v $"
     4	//   revision         = "$Revision: 1.6 $"
     5	//   date             = "$Date: 2008/07/30 22:07:46 $"
     6	//   tag              = "$Name: libnucnet-0-4 $"
     7	//   template_version = "cp_template.0.12"
     8	//   type             = "public"
     9	// >
    10	//   <license>
    11	//      See the README.txt file in this directory for copyright and license
    12	//      information.
    13	//   </license>
    14	//   <description>
    15	//     <abstract>
    16	//       Example to demonstrate how to use libnucnet routines to create
    17	//       a new Libnucnet structure from an input xml file, evolve a multi-zone
    18	//       nuclear network system, and clear the structure and free the
    19	//       allocated memory.
    20	//     </abstract>
    21	//   </description>
    22	// </file>
    23	//////////////////////////////////////////////////////////////////////////////*/
    24	
    25	/*##############################################################################
    26	// Includes.
    27	//############################################################################*/
    28	
    29	#include <Libnucnet.h>
    30	#include "remove_duplicate.h"
    31	#include "screen.h"
    32	#include "coul_corr.h"
    33	
    34	/*##############################################################################
    35	// Define some parameters.
    36	//############################################################################*/
    37	
    38	#define D_MIN          1.e-4   /* Newton-Raphson convergence criterion */
    39	#define D_Y_MIN        1.e-30  /* Smallest y for convergence */
    40	#define D_Y_MIN_PRINT  1.e-30  /* Smallest y to print out */
    41	#define I_ITMAX        4       /* Maximum number of Newton-Raphson iterations */
    42	
    43	#define D_REGT         0.15    /* Max fractional change in dt */
    44	#define D_REGY         0.15    /* Max fractional change in y */
    45	#define D_YMIN         1.e-10  /* Smallest y for change in dt */
    46	
    47	/*##############################################################################
    48	// Screening.  "no" = no screening, "yes" = screening.
    49	//############################################################################*/
    50	
    51	#define USE_SCREENING  "yes"
    52	
    53	/*##############################################################################
    54	// User structures.
    55	//############################################################################*/
    56	
    57	typedef struct {
    58	  Libnucnet *pLibnucnet;
    59	  gsl_vector *pYOld;
    60	  gsl_vector *pRhs;
    61	  gsl_vector *pSolution;
    62	  WnMatrix *pMatrix;
    63	  double *aT9;
    64	  double *aRho;
    65	  double dMixTime;
    66	  double dDt;
    67	  double dCheck;
    68	} user_data;
    69	
    70	  user_screening_data my_screening_data;
    71	  user_coul_corr_data my_corr_data;
    72	
    73	/*##############################################################################
    74	// Prototypes.
    75	//############################################################################*/
    76	
    77	int evolve_system( Libnucnet *, double *, double *, double, double );
    78	int print_zone( Libnucnet__Zone *, void * );
    79	void print_abunds( Libnucnet__Zone * );
    80	int print_abund( Libnucnet__Species *, Libnucnet__Zone * );
    81	int get_next_timestep( Libnucnet__Zone *, double * );
    82	int set_zone_abundances( Libnucnet__Zone *, gsl_vector * );
    83	int set_zone( Libnucnet__Zone *, user_data * );
    84	int set_mixing_terms( Libnucnet__Zone *, user_data * );
    85	int update_abundances( Libnucnet__Zone *, user_data * );
    86	int update_abundance_changes( Libnucnet__Zone *, gsl_vector * );
    87	int initialize_zone( Libnucnet__Zone *, void * );
    88	
    89	/*##############################################################################
    90	// main()
    91	//############################################################################*/
    92	
    93	int main( int argc, char * argv[] ) {
    94	
    95	  double *a_t9, *a_rho, d_tend, d_t, d_dt;
    96	  double d_t9, d_rho, d_mix_time;
    97	  size_t i_index;
    98	  int i_step;
    99	  Libnucnet *p_my_nucnet;
   100	  FILE *p_in;
   101	
   102	  /*============================================================================
   103	  // Check input.
   104	  //==========================================================================*/
   105	
   106	  if ( argc < 5 || argc > 7 ) {
   107	    fprintf(
   108	      stderr,
   109	      "\nUsage: %s network_file param_file mix_time t_end xpath_nuc xpath_reac\n\n",
   110	      argv[0]
   111	    );
   112	    fprintf(
   113	      stderr, "  network_file = input nuclear data xml filename\n\n"
   114	    );
   115	    fprintf(
   116	      stderr, "  param_file = input t9 and rho data filename\n\n"
   117	    );
   118	    fprintf(
   119	      stderr, "  mix_time = mixing time between zones (0 = no mixing)\n\n"
   120	    );
   121	    fprintf(
   122	      stderr, "  t_end = time to evolve system\n\n"
   123	    );
   124	    fprintf(
   125	      stderr, "  xpath_nuc = nuclear xpath expression (optional--required if xpath_reac specified\n\n"
   126	    );
   127	    fprintf(
   128	      stderr, "  xpath_reac = reaction xpath expression (optional)\n\n"
   129	    );
   130	    return EXIT_FAILURE;
   131	  }
   132	
   133	  /*============================================================================
   134	  // Validate input xml.
   135	  //==========================================================================*/
   136	
   137	  if( !Libnucnet__is_valid_input_xml( argv[1] ) ) {
   138	    fprintf( stderr, "Not valid libnucnet input xml!\n" );
   139	    return EXIT_FAILURE;
   140	  }
   141	
   142	  /*============================================================================
   143	  // Read and store network data.
   144	  //==========================================================================*/
   145	
   146	  if( argc == 5 ) {
   147	    p_my_nucnet = Libnucnet__new_from_xml( argv[1], NULL, NULL );
   148	  } else if( argc == 6 ) {
   149	    p_my_nucnet = Libnucnet__new_from_xml( argv[1], argv[5], NULL );
   150	  } else {
   151	    p_my_nucnet = Libnucnet__new_from_xml( argv[1], argv[5], argv[6] );
   152	  }
   153	 
   154	  /*============================================================================
   155	  // Open input file containing zone data.
   156	  //==========================================================================*/
   157	
   158	  if ( ( p_in = fopen( argv[2], "r" ) ) == NULL ) {
   159	
   160	    fprintf( stderr, "\nCannot open input param file!\n" );
   161	    return EXIT_FAILURE;
   162	
   163	  }
   164	
   165	  /*============================================================================
   166	  // Allocate memory for t9 and rho arrays.
   167	  //==========================================================================*/
   168	
   169	  a_t9 = ( double * ) malloc( 
   170	    sizeof( double ) * Libnucnet__getNumberOfZones( p_my_nucnet ) 
   171	  );
   172	
   173	  if( !a_t9 ) {
   174	    fprintf( stderr, "Couldn't allocate memory for t9 array!\n" ); 
   175	    return EXIT_FAILURE;
   176	  }
   177	
   178	  a_rho = ( double * ) malloc( 
   179	    sizeof( double ) * Libnucnet__getNumberOfZones( p_my_nucnet ) 
   180	  );
   181	
   182	  if( !a_rho ) {
   183	    fprintf( stderr, "Couldn't allocate memory for rho array!\n" );  
   184	    return EXIT_FAILURE;
   185	  }
   186	
   187	  /*============================================================================
   188	  // Read in zone data.
   189	  //==========================================================================*/
   190	
   191	  i_index = 0;
   192	
   193	  while ( !feof( p_in ) ) {
   194	
   195	    if( i_index == Libnucnet__getNumberOfZones( p_my_nucnet ) ) {
   196	       fprintf( stderr, "Too many zones in input file!\n" );
   197	       return EXIT_FAILURE;
   198	    }
   199	
   200	    fscanf( p_in, "%lf %lf\n", &d_t9, &d_rho );
   201	    a_t9[i_index] = d_t9;
   202	    a_rho[i_index] = d_rho;
   203	    ++i_index;
   204	
   205	  }
   206	
   207	  /*============================================================================
   208	  // Close input file.
   209	  //==========================================================================*/
   210	
   211	  fclose( p_in );
   212	
   213	  /*============================================================================
   214	  // Check number of zones.
   215	  //==========================================================================*/
   216	
   217	  if ( Libnucnet__getNumberOfZones( p_my_nucnet ) != i_index ) {
   218	      printf( "Number of zones wrong in param.dat!\n" );
   219	      return EXIT_FAILURE;
   220	  }
   221	
   222	  /*============================================================================
   223	  // Iterate over zones to remove duplicates and set screening.
   224	  //==========================================================================*/
   225	
   226	  Libnucnet__iterateZones(
   227	    p_my_nucnet,
   228	    (Libnucnet__Zone__iterateFunction) initialize_zone,
   229	    NULL
   230	  );
   231	
   232	  /*============================================================================
   233	  // Store mixing time and final time.
   234	  //==========================================================================*/
   235	
   236	  d_mix_time = atof( argv[3] );
   237	  d_tend = atof( argv[4] );
   238	
   239	  /*============================================================================
   240	  // Initialize time.
   241	  //==========================================================================*/
   242	
   243	  d_dt = 1.e-10;
   244	  d_t = 0.0;
   245	
   246	  /*============================================================================
   247	  // Evolve network while t < final t.
   248	  //==========================================================================*/
   249	
   250	  i_step = 0;
   251	
   252	  while ( d_t < d_tend ) {
   253	
   254	    d_t += d_dt;
   255	
   256	    if ( evolve_system( p_my_nucnet, a_t9, a_rho, d_mix_time, d_dt ) != 0 ) {
   257	
   258	      fprintf( stderr, "Error evolving system!\n" );
   259	      return EXIT_FAILURE;
   260	
   261	    }
   262	
   263	  /*============================================================================
   264	  // Print out abundances every 25th time step.
   265	  //==========================================================================*/
   266	
   267	    if( i_step % 25 == 0 ) {
   268	
   269	      printf( "t = %10.4e, dt = %10.4e\n\n", d_t, d_dt );
   270	
   271	      Libnucnet__iterateZones(
   272	        p_my_nucnet,
   273	        (Libnucnet__Zone__iterateFunction) print_zone,
   274	        NULL
   275	      );
   276	
   277	    }
   278	
   279	  /*============================================================================
   280	  // Update timestep.
   281	  //==========================================================================*/
   282	
   283	    Libnucnet__iterateZones(
   284	      p_my_nucnet,
   285	      (Libnucnet__Zone__iterateFunction) get_next_timestep,
   286	      &d_dt
   287	    ); 
   288	
   289	    if ( d_t + d_dt > d_tend ) {
   290	
   291	      d_dt = d_tend - d_t;
   292	
   293	    }
   294	
   295	    i_step++;
   296	
   297	  }  
   298	
   299	  /*============================================================================
   300	  // Print out final abundances.
   301	  //==========================================================================*/
   302	
   303	  printf( "t = %10.4e, dt = %10.4e\n\n", d_t, d_dt );
   304	
   305	  Libnucnet__iterateZones(
   306	    p_my_nucnet,
   307	    (Libnucnet__Zone__iterateFunction) print_zone,
   308	    NULL
   309	  );
   310	
   311	  /*============================================================================
   312	  // Clean up and exit.
   313	  //==========================================================================*/
   314	
   315	  free( a_t9 );
   316	  free( a_rho );
   317	  Libnucnet__free( p_my_nucnet );
   318	  return EXIT_SUCCESS;
   319	
   320	}
   321	
   322	/*##############################################################################
   323	// evolve_system()
   324	//############################################################################*/
   325	
   326	int
   327	evolve_system( 
   328	  Libnucnet *self, double *a_t9, double *a_rho, double d_mix_time, double d_dt 
   329	) {
   330	
   331	  size_t i_species, i_zones, i_dim, i_iter;
   332	  gsl_vector *p_abundance_changes;
   333	  user_data *p_user_data;
   334	  
   335	  /*============================================================================
   336	  // Allocate Memory.
   337	  //==========================================================================*/
   338	
   339	  i_species =
   340	    Libnucnet__Nuc__getNumberOfSpecies(
   341	      Libnucnet__Net__getNuc(
   342	        Libnucnet__getNet( self )
   343	      )
   344	    );
   345	
   346	  i_zones = Libnucnet__getNumberOfZones( self );
   347	
   348	  i_dim = i_species * i_zones;
   349	
   350	  p_abundance_changes = gsl_vector_calloc( i_dim );
   351	
   352	  p_user_data = ( user_data * ) malloc( sizeof( user_data ) );
   353	
   354	  p_user_data->pRhs = gsl_vector_calloc( i_dim );
   355	
   356	  p_user_data->aT9 = a_t9;
   357	  p_user_data->aRho = a_rho;
   358	  p_user_data->dMixTime = d_mix_time;
   359	  p_user_data->dDt = d_dt;
   360	  p_user_data->pYOld = gsl_vector_calloc( i_dim );
   361	
   362	  /*============================================================================
   363	  // Store initial abundances.
   364	  //==========================================================================*/
   365	
   366	  p_user_data->pLibnucnet = self;
   367	
   368	  Libnucnet__iterateZones(
   369	    self,
   370	    (Libnucnet__Zone__iterateFunction) set_zone_abundances,
   371	    p_user_data->pYOld
   372	  );
   373	
   374	  /*============================================================================
   375	  // Newton-Raphson iterations.
   376	  //==========================================================================*/
   377	
   378	  for( i_iter = 0; i_iter < I_ITMAX; i_iter++ ) {
   379	
   380	    p_user_data->pMatrix = WnMatrix__new( i_dim, i_dim );
   381	
   382	    /*==========================================================================
   383	    // Loop Over Zones.
   384	    //========================================================================*/
   385	
   386	    Libnucnet__iterateZones(
   387	      self,
   388	      (Libnucnet__Zone__iterateFunction) set_zone,
   389	      p_user_data
   390	    ); 
   391	
   392	
   393	    /*==========================================================================
   394	    // Assign mixing terms.
   395	    //========================================================================*/
   396	
   397	    if( d_mix_time < 0. ) {
   398	       fprintf( stderr, "Invalid mixing time" );
   399	       exit( EXIT_FAILURE );
   400	    }
   401	
   402	    if( d_mix_time != 0. )
   403	      Libnucnet__iterateZones(
   404	        self,
   405	        (Libnucnet__Zone__iterateFunction) set_mixing_terms,
   406	        p_user_data
   407	      );
   408	
   409	    /*========================================================================
   410	    // Add 1/dt to diagonals.
   411	    //========================================================================*/
   412	
   413	    WnMatrix__addValueToDiagonals( p_user_data->pMatrix, 1.0 / d_dt );
   414	
   415	    /*========================================================================
   416	    // Solve matrix.
   417	    //========================================================================*/
   418	
   419	    p_user_data->pSolution =
   420	      WnMatrix__solve( p_user_data->pMatrix, p_user_data->pRhs );
   421	  
   422	    /*==========================================================================
   423	    // Update abundances.
   424	    //========================================================================*/
   425	
   426	    p_user_data->dCheck = 0.;
   427	
   428	    Libnucnet__iterateZones(
   429	      self,
   430	      (Libnucnet__Zone__iterateFunction) update_abundances,
   431	      p_user_data
   432	    );
   433	
   434	    /*==========================================================================
   435	    // Update abundance changes.
   436	    //========================================================================*/
   437	
   438	    gsl_vector_add( p_abundance_changes, p_user_data->pSolution );
   439	
   440	    /*==========================================================================
   441	    // Free matrix and solution vector.
   442	    //========================================================================*/
   443	
   444	    gsl_vector_free( p_user_data->pSolution );
   445	    WnMatrix__free( p_user_data->pMatrix );
   446	
   447	    /*==========================================================================
   448	    // Break if solution converged.
   449	    //========================================================================*/
   450	
   451	    if( p_user_data->dCheck < D_MIN ) break;
   452	
   453	  }
   454	
   455	  /*============================================================================
   456	  // Update abundance changes.
   457	  //==========================================================================*/
   458	
   459	  Libnucnet__iterateZones(
   460	    self,
   461	    (Libnucnet__Zone__iterateFunction) update_abundance_changes,
   462	    p_abundance_changes
   463	  );
   464	
   465	
   466	  /*==========================================================================
   467	  // Free allocated memory.
   468	  //========================================================================*/
   469	
   470	  gsl_vector_free( p_abundance_changes );
   471	
   472	  gsl_vector_free( p_user_data->pRhs );
   473	  gsl_vector_free( p_user_data->pYOld );
   474	
   475	  free( p_user_data );
   476	
   477	  return 0;
   478	
   479	}
   480	
   481	/*##############################################################################
   482	// print_zone()
   483	//############################################################################*/
   484	
   485	int
   486	print_zone( Libnucnet__Zone *p_zone, void *p_data )
   487	{
   488	
   489	  if( p_data ) {
   490	    fprintf( stderr, "Routine should have no extra data!\n" );
   491	    return 0;
   492	  }
   493	
   494	  fprintf(
   495	    stdout,
   496	    "Zone: %s %s %s\n\n",
   497	    Libnucnet__Zone__getLabel( p_zone, 1 ),
   498	    Libnucnet__Zone__getLabel( p_zone, 2 ),
   499	    Libnucnet__Zone__getLabel( p_zone, 3 )
   500	  );
   501	
   502	  print_abunds( p_zone );
   503	
   504	  return 1;
   505	
   506	}
   507	
   508	/*##############################################################################
   509	// print_abunds()
   510	//############################################################################*/
   511	
   512	void
   513	print_abunds( Libnucnet__Zone *p_zone ) {
   514	
   515	  Libnucnet__Nuc__iterateSpecies(
   516	    Libnucnet__Net__getNuc(
   517	      Libnucnet__Zone__getNet( p_zone )
   518	    ),
   519	    (Libnucnet__Species__iterateFunction) print_abund,
   520	    p_zone
   521	  );
   522	
   523	  fprintf(
   524	    stdout,
   525	    "\n1 - xsum = %g\n",
   526	    1. - Libnucnet__Zone__computeAMoment( p_zone, 1 )
   527	  );
   528	
   529	}
   530	
   531	/*##############################################################################
   532	// print_abund()
   533	//############################################################################*/
   534	
   535	int
   536	print_abund( Libnucnet__Species *p_species, Libnucnet__Zone *p_zone )
   537	{
   538	
   539	  double d_abund;
   540	
   541	  d_abund = Libnucnet__Zone__getSpeciesAbundance( p_zone, p_species );
   542	
   543	  if ( d_abund > D_Y_MIN_PRINT ) {
   544	
   545	    printf(
   546	      "%5d%5d%14.4e%14.4e\n",
   547	      Libnucnet__Species__getZ( p_species ),
   548	      Libnucnet__Species__getA( p_species ),
   549	      d_abund,
   550	      Libnucnet__Zone__getSpeciesAbundanceChange( p_zone, p_species )
   551	    );
   552	
   553	  }
   554	
   555	  return 1;
   556	
   557	}
   558	
   559	/*##############################################################################
   560	// get_next_timestep().
   561	//############################################################################*/
   562	
   563	int
   564	get_next_timestep( Libnucnet__Zone *p_zone, double *p_dt )
   565	{
   566	
   567	  Libnucnet__Zone__updateTimeStep(
   568	    p_zone,
   569	    p_dt,
   570	    D_REGT,
   571	    D_REGY,
   572	    D_YMIN
   573	  );
   574	
   575	  return 1;
   576	
   577	}
   578	
   579	/*##############################################################################
   580	// set_zone_abundances().
   581	//############################################################################*/
   582	
   583	int
   584	set_zone_abundances( Libnucnet__Zone *p_zone, gsl_vector *p_old )
   585	{
   586	
   587	  gsl_vector *p_abundances;
   588	  gsl_vector_view v_old;
   589	  size_t i_zone;
   590	
   591	  i_zone = atoi( Libnucnet__Zone__getLabel( p_zone, 1 ) );
   592	
   593	  p_abundances =
   594	    Libnucnet__Zone__getAbundances( p_zone );
   595	
   596	  v_old =
   597	    gsl_vector_subvector(
   598	      p_old,
   599	      i_zone * WnMatrix__get_gsl_vector_size( p_abundances ),
   600	      WnMatrix__get_gsl_vector_size( p_abundances )
   601	    );
   602	
   603	  gsl_vector_memcpy( &v_old.vector, p_abundances );
   604	
   605	  gsl_vector_free( p_abundances );
   606	
   607	  return 1;
   608	
   609	}
   610	
   611	/*##############################################################################
   612	// set_zone()
   613	//############################################################################*/
   614	
   615	int
   616	set_zone( Libnucnet__Zone *p_zone, user_data *p_user_data )
   617	{
   618	
   619	  size_t i_zone;
   620	  gsl_vector *p_abund, *p_rhs;
   621	  gsl_vector_view v_old, v_rhs;
   622	  WnMatrix *p_matrix;
   623	
   624	  /*========================================================================
   625	  // Get zone.
   626	  //======================================================================*/
   627	
   628	  i_zone = atoi( Libnucnet__Zone__getLabel( p_zone, 1 ) );
   629	
   630	  /*========================================================================
   631	  // If T9 = 0, return.
   632	  //======================================================================*/
   633	
   634	  if( p_user_data->aT9[i_zone] == 0 ) return 1;
   635	
   636	  /*========================================================================
   637	  // Compute rates.
   638	  //======================================================================*/
   639	
   640	  Libnucnet__Zone__computeRates(
   641	    p_zone,
   642	    p_user_data->aT9[i_zone],
   643	    p_user_data->aRho[i_zone]
   644	  );
   645	
   646	  /*======================================================================
   647	  // Get right hand side vector.
   648	  //======================================================================*/
   649	
   650	  p_rhs = Libnucnet__Zone__computeFlowVector( p_zone );
   651	
   652	  p_abund =
   653	    Libnucnet__Zone__getAbundances( p_zone );
   654	
   655	  v_old =
   656	    gsl_vector_subvector(
   657	      p_user_data->pYOld,
   658	      i_zone * WnMatrix__get_gsl_vector_size( p_abund ),
   659	      WnMatrix__get_gsl_vector_size( p_abund )
   660	    );
   661	
   662	  gsl_vector_sub(
   663	    p_abund, &v_old.vector
   664	  );
   665	
   666	  gsl_vector_scale( p_abund, 1. / p_user_data->dDt );
   667	
   668	  gsl_vector_sub( p_rhs, p_abund );
   669	
   670	  /*======================================================================
   671	  // Insert into large vector.
   672	  //======================================================================*/
   673	
   674	  v_rhs =
   675	    gsl_vector_subvector(
   676	      p_user_data->pRhs,
   677	      i_zone * WnMatrix__get_gsl_vector_size( p_rhs ),
   678	      WnMatrix__get_gsl_vector_size( p_rhs )
   679	    );
   680	
   681	  gsl_vector_memcpy( &v_rhs.vector, p_rhs );
   682	
   683	  gsl_vector_free( p_rhs );
   684	
   685	  /*========================================================================
   686	  // Get the Jacobian Matrix.
   687	  //======================================================================*/
   688	
   689	  p_matrix = Libnucnet__Zone__computeJacobianMatrix( p_zone );
   690	
   691	  /*========================================================================
   692	  // Insert single-zone matrix into full matrix.
   693	  //======================================================================*/
   694	
   695	  WnMatrix__insertMatrix(
   696	    p_user_data->pMatrix,
   697	    p_matrix,
   698	    i_zone * WnMatrix__get_gsl_vector_size( p_abund ) + 1,
   699	    i_zone * WnMatrix__get_gsl_vector_size( p_abund ) + 1
   700	  );
   701	
   702	  /*========================================================================
   703	  // Free memory allocated for single-zone matrix and abundances. 
   704	  //======================================================================*/
   705	
   706	  gsl_vector_free( p_abund );
   707	  WnMatrix__free( p_matrix );
   708	
   709	  return 1;
   710	
   711	}
   712	
   713	/*##############################################################################
   714	// set_mixing_terms().
   715	//############################################################################*/
   716	
   717	int
   718	set_mixing_terms( Libnucnet__Zone *p_zone, user_data *p_user_data )
   719	{
   720	
   721	  Libnucnet__Zone *p_previous, *p_next;
   722	  gsl_vector *p_abundances;
   723	  gsl_vector_view v_view;
   724	  size_t i, i_species, i_zone;
   725	  char s_label[256];
   726	
   727	  i_species =
   728	    Libnucnet__Nuc__getNumberOfSpecies(
   729	      Libnucnet__Net__getNuc(
   730	        Libnucnet__Zone__getNet( p_zone )
   731	      )
   732	    );
   733	
   734	  i_zone = atoi( Libnucnet__Zone__getLabel( p_zone, 1 ) );
   735	
   736	  /*============================================================================
   737	  // Add mixing to/from previous zone.
   738	  //==========================================================================*/
   739	
   740	  if( i_zone != 0 ) {
   741	
   742	    sprintf( s_label, "%d", i_zone - 1 );
   743	
   744	    v_view =
   745	      gsl_vector_subvector(
   746	        p_user_data->pRhs,
   747	        i_zone * i_species,
   748	        i_species
   749	      );
   750	  
   751	    p_previous =
   752	      Libnucnet__getZoneByLabels(
   753	        p_user_data->pLibnucnet,
   754	        s_label,
   755	        "0",
   756	        "0"
   757	      );
   758	
   759	    for( i = 0; i < i_species; i++ )
   760	      WnMatrix__assignElement(
   761	        p_user_data->pMatrix,
   762	        i_zone * i_species + i + 1,
   763	        ( i_zone - 1 ) * i_species + i + 1,
   764	        -1. / p_user_data->dMixTime
   765	      );
   766	
   767	    p_abundances = Libnucnet__Zone__getAbundances( p_previous );
   768	
   769	    gsl_vector_scale( p_abundances, 1. / p_user_data->dMixTime );
   770	
   771	    gsl_vector_add( &v_view.vector, p_abundances );
   772	
   773	    gsl_vector_free( p_abundances );
   774	
   775	    for( i = 0; i < i_species; i++ )
   776	      WnMatrix__assignElement(
   777	        p_user_data->pMatrix,
   778	        i_zone * i_species + i + 1,
   779	        i_zone * i_species + i + 1,
   780	        1. / p_user_data->dMixTime
   781	      );
   782	
   783	    p_abundances = Libnucnet__Zone__getAbundances( p_zone );
   784	
   785	    gsl_vector_scale( p_abundances, -1. / p_user_data->dMixTime );
   786	
   787	    gsl_vector_add( &v_view.vector, p_abundances );
   788	
   789	    gsl_vector_free( p_abundances );
   790	
   791	  }
   792	
   793	  /*============================================================================
   794	  // Add mixing to/from next zone.
   795	  //==========================================================================*/
   796	
   797	  if( i_zone < Libnucnet__getNumberOfZones( p_user_data->pLibnucnet ) - 1 ) {
   798	
   799	    sprintf( s_label, "%d", i_zone + 1 );
   800	
   801	    v_view =
   802	      gsl_vector_subvector(
   803	        p_user_data->pRhs,
   804	        i_zone * i_species,
   805	        i_species
   806	      );
   807	  
   808	    p_next =
   809	      Libnucnet__getZoneByLabels(
   810	        p_user_data->pLibnucnet,
   811	        s_label,
   812	        "0",
   813	        "0"
   814	      );
   815	
   816	    for( i = 0; i < i_species; i++ )
   817	      WnMatrix__assignElement(
   818	        p_user_data->pMatrix,
   819	        i_zone * i_species + i + 1,
   820	        ( i_zone + 1 ) * i_species + i + 1,
   821	        -1. / p_user_data->dMixTime
   822	      );
   823	
   824	    p_abundances = Libnucnet__Zone__getAbundances( p_next );
   825	
   826	    gsl_vector_scale( p_abundances, 1. / p_user_data->dMixTime );
   827	
   828	    gsl_vector_add( &v_view.vector, p_abundances );
   829	
   830	    gsl_vector_free( p_abundances );
   831	
   832	    for( i = 0; i < i_species; i++ )
   833	      WnMatrix__assignElement(
   834	        p_user_data->pMatrix,
   835	        i_zone * i_species + i + 1,
   836	        i_zone * i_species + i + 1,
   837	        1. / p_user_data->dMixTime
   838	      );
   839	
   840	    p_abundances = Libnucnet__Zone__getAbundances( p_zone );
   841	
   842	    gsl_vector_scale( p_abundances, -1. / p_user_data->dMixTime );
   843	
   844	    gsl_vector_add( &v_view.vector, p_abundances );
   845	
   846	    gsl_vector_free( p_abundances );
   847	
   848	  }
   849	
   850	  return 1;
   851	
   852	}
   853	
   854	/*##############################################################################
   855	// update_abundances().
   856	//############################################################################*/
   857	
   858	int
   859	update_abundances( Libnucnet__Zone *p_zone, user_data *p_user_data )
   860	{
   861	
   862	  size_t i, i_zone, i_species;
   863	  double d_y, d_checkT;
   864	  gsl_vector *p_new_abundances;
   865	  gsl_vector_view v_view;
   866	
   867	  i_zone = atoi( Libnucnet__Zone__getLabel( p_zone, 1 ) );
   868	
   869	  i_species =
   870	    Libnucnet__Nuc__getNumberOfSpecies(
   871	      Libnucnet__Net__getNuc(
   872	        Libnucnet__Zone__getNet( p_zone )
   873	      )
   874	    );
   875	
   876	  p_new_abundances = Libnucnet__Zone__getAbundances( p_zone );
   877	
   878	  v_view =
   879	    gsl_vector_subvector(
   880	       p_user_data->pSolution,
   881	       i_species * i_zone,
   882	       i_species
   883	    );
   884	
   885	  gsl_vector_add( p_new_abundances, &v_view.vector );
   886	
   887	  Libnucnet__Zone__updateAbundances( p_zone, p_new_abundances );
   888	
   889	  for( i = 0; i < WnMatrix__get_gsl_vector_size( p_new_abundances ); i++ ) {
   890	    d_y = gsl_vector_get( p_new_abundances, i );
   891	    if( d_y > D_Y_MIN ) {
   892	      d_checkT = fabs( gsl_vector_get( &v_view.vector, i ) / d_y );
   893	      if( d_checkT > p_user_data->dCheck )
   894	        p_user_data->dCheck = d_checkT;
   895	    }
   896	
   897	  }
   898	
   899	  gsl_vector_free( p_new_abundances );
   900	
   901	  return 1;
   902	
   903	}
   904	
   905	/*##############################################################################
   906	// update_abundance_changes().
   907	//############################################################################*/
   908	
   909	int
   910	update_abundance_changes( Libnucnet__Zone *p_zone, gsl_vector *p_changes )
   911	{
   912	
   913	  size_t i_zone, i_species;
   914	  gsl_vector_view v_view;
   915	
   916	  i_zone = atoi( Libnucnet__Zone__getLabel( p_zone, 1 ) );
   917	
   918	  i_species =
   919	    Libnucnet__Nuc__getNumberOfSpecies(
   920	      Libnucnet__Net__getNuc(
   921	        Libnucnet__Zone__getNet( p_zone )
   922	      )
   923	    );
   924	
   925	  v_view =
   926	    gsl_vector_subvector(
   927	       p_changes,
   928	       i_zone * i_species,
   929	       i_species
   930	    );
   931	
   932	  Libnucnet__Zone__updateAbundanceChanges( p_zone, &v_view.vector );
   933	
   934	  return 1;
   935	
   936	}
   937	
   938	/*##############################################################################
   939	// initialize_zone().
   940	//############################################################################*/
   941	
   942	int
   943	initialize_zone( Libnucnet__Zone *p_zone, void *p_data )
   944	{
   945	
   946	  Libnucnet__Reac *p_my_duplicates;
   947	
   948	  /*============================================================================
   949	  // Check input.
   950	  //==========================================================================*/
   951	
   952	  if( p_data ) {
   953	    fprintf( stderr, "Routine takes no extra data!\n" );
   954	    exit( EXIT_FAILURE );
   955	  }
   956	
   957	  /*============================================================================
   958	  // Remove duplicate reactions.
   959	  //==========================================================================*/
   960	
   961	  p_my_duplicates =
   962	    Libnucnet__Reac__getDuplicateReactions(
   963	      Libnucnet__Net__getReac(
   964	        Libnucnet__Zone__getNet( p_zone )
   965	      )
   966	    );
   967	
   968	  Libnucnet__Reac__iterateReactions(
   969	    p_my_duplicates,
   970	    (Libnucnet__Reaction__iterateFunction) remove_duplicate,
   971	    Libnucnet__Zone__getNet( p_zone )
   972	  );
   973	
   974	  Libnucnet__Reac__free( p_my_duplicates );
   975	
   976	  /*============================================================================
   977	  // Set screening, if desired.
   978	  //==========================================================================*/
   979	
   980	  if( strcmp( USE_SCREENING, "yes" ) == 0 ) {
   981	
   982	    /*==========================================================================
   983	    // Set screening data and function.
   984	    //========================================================================*/
   985	
   986	    my_screening_data.dYe2 =
   987	      Libnucnet__Zone__computeZMoment( p_zone, 2 );
   988	
   989	    Libnucnet__Zone__setScreeningFunction(
   990	      p_zone,
   991	      (Libnucnet__Net__screeningFunction) my_screening_function,
   992	      &my_screening_data
   993	    );
   994	
   995	    /*==========================================================================
   996	    // Set correction function.
   997	    //========================================================================*/
   998	
   999	    my_corr_data.dA1 = -0.9052;
  1000	    my_corr_data.dA2 = 0.6322;
  1001	
  1002	    Libnucnet__Zone__setNseCorrectionFactorFunction(
  1003	      p_zone,
  1004	      (Libnucnet__Species__nseCorrectionFactorFunction) my_coulomb_correction,
  1005	      &my_corr_data
  1006	    );
  1007	
  1008	  }
  1009	
  1010	  return 1;
  1011	
  1012	}
  1013