make test similar to python

unit_test/burn_cell_metal_chem/_parameters

@@ -17,7 +17,8 @@ tfirst        real       0.0
 # number of steps for the single zone burn
 nsteps        int        1000
 
-# initial temperature
+# initial number density and temperature
+ninit         real       1e-1
 temperature   real       1e2
 
 #list of species and their number densities used in the network (39 if including deuterium)

unit_test/burn_cell_metal_chem/burn_cell.H

@@ -23,6 +23,8 @@ auto burn_cell_c() -> int {
   int x;
   int nn;
 
+  amrex::Real metal = network_rp::metallicity;
+
   for (n = 1; n <= NumSpec; ++n) {
     switch (n) {
 
@@ -75,7 +77,7 @@ auto burn_cell_c() -> int {
       numdens[n - 1] = primary_species_16;
       break;
     case 17:
-      numdens[n - 1] = primary_species_17;
+      numdens[n - 1] = primary_species_17*metal;
       break;
     case 18:
       numdens[n - 1] = primary_species_18;
@@ -93,7 +95,7 @@ auto burn_cell_c() -> int {
       numdens[n - 1] = primary_species_22;
       break;
     case 23:
-      numdens[n - 1] = primary_species_23;
+      numdens[n - 1] = primary_species_23*metal;
       break;
     case 24:
       numdens[n - 1] = primary_species_24;
@@ -131,10 +133,25 @@ auto burn_cell_c() -> int {
     }
   }
 
+  //scale number densities by initial ninit
+  for (n = 0; n < NumSpec; ++n) {
+    numdens[n] *= ninit;
+  }
+
+  //if metallicity is 0, reset metal number densities to 0
+  if (metal == 0) {
+    for (n = 0; n < NumSpec; ++n) {
+      if ((n < 2) || (n > 15)) {
+        state.xn[n] = 0.0;
+      }
+    }
+  }
+
+
   // Echo initial conditions at burn and fill burn state input
 
   std::cout << "Redshift: " << network_rp::redshift << std::endl;
-  std::cout << "Metallicity: " << network_rp::metallicity << std::endl;
+  std::cout << "Metallicity: " << metal << std::endl;
   std::cout << "Dust2gas Ratio: " << network_rp::dust2gas_ratio << std::endl;
   std::cout << " " << std::endl;
 
@@ -286,7 +303,7 @@ auto burn_cell_c() -> int {
     // get number density
     // make abundance 0 for all metals if metallicity is 0
     for (nn = 0; nn < NumSpec; ++nn) {
-      if (network_rp::metallicity == 0 && ((nn < 2) || (nn > 15))) {
+      if (metal == 0 && ((nn < 2) || (nn > 15))) {
         state.xn[nn] = 0.0;
       }
       sum_numdens += state.xn[nn];

unit_test/burn_cell_metal_chem/inputs_metal_chem

@@ -10,51 +10,52 @@ unit_test.nsteps = 1000
 unit_test.tmax = 7.e20
 # initial temperature
 unit_test.temperature = 3e2
-# initial number densities
-unit_test.primary_species_1 = 1e-41
-unit_test.primary_species_2 = 1e-41
-unit_test.primary_species_3 = 1e-5
-unit_test.primary_species_4 = 1e-5
-unit_test.primary_species_5 = 1e-1
-unit_test.primary_species_6 = 1e-41
-unit_test.primary_species_7 = 1e-41
-unit_test.primary_species_8 = 1e-41
-unit_test.primary_species_9 = 1e-40
-unit_test.primary_species_10 = 1e-41
-unit_test.primary_species_11 = 1e-7
-unit_test.primary_species_12 = 1e-41
-unit_test.primary_species_13 = 1e-41
-unit_test.primary_species_14 = 1e-41
-unit_test.primary_species_15 = 1e-41
-unit_test.primary_species_16 = 0.0775e-1
-unit_test.primary_species_17 = 9.27e-6
-unit_test.primary_species_18 = 1e-41
-unit_test.primary_species_19 = 1e-41
-unit_test.primary_species_20 = 1e-41
-unit_test.primary_species_21 = 1e-41
-unit_test.primary_species_22 = 1e-41
-unit_test.primary_species_23 = 3.568e-5
-unit_test.primary_species_24 = 1e-41
-unit_test.primary_species_25 = 1e-41
-unit_test.primary_species_26 = 1e-41
-unit_test.primary_species_27 = 1e-41
-unit_test.primary_species_28 = 1e-41
-unit_test.primary_species_29 = 1e-41
-unit_test.primary_species_30 = 1e-41
-unit_test.primary_species_31 = 1e-41
-unit_test.primary_species_32 = 1e-41
-unit_test.primary_species_33 = 1e-41
-unit_test.primary_species_34 = 1e-41
+unit_test.ninit = 1e-1
+# initial number densities (will be scaled to metallicity provided below automatically by burn_cell)
+unit_test.primary_species_1 = 1e-40     #co_ice
+unit_test.primary_species_2 = 1e-40     #h2o_ice
+unit_test.primary_species_3 = 1e-4      #e
+unit_test.primary_species_4 = 1e-4      #h+
+unit_test.primary_species_5 = 1e0       #h
+unit_test.primary_species_6 = 1e-40     #h-
+unit_test.primary_species_7 = 1e-40     #d+
+unit_test.primary_species_8 = 1e-5      #d
+unit_test.primary_species_9 = 1e-40     #h2+
+unit_test.primary_species_10 = 1e-40    #d-
+unit_test.primary_species_11 = 1e-6     #h2
+unit_test.primary_species_12 = 1e-40    #hd+
+unit_test.primary_species_13 = 1e-40    #hd
+unit_test.primary_species_14 = 1e-40    #he++
+unit_test.primary_species_15 = 1e-40    #he+
+unit_test.primary_species_16 = 0.0775e0 #he
+unit_test.primary_species_17 = 9.27e-5  #c+
+unit_test.primary_species_18 = 1e-40    #c
+unit_test.primary_species_19 = 1e-40    #ch
+unit_test.primary_species_20 = 1e-40    #ch2
+unit_test.primary_species_21 = 1e-40    #ch3
+unit_test.primary_species_22 = 1e-40    #o+
+unit_test.primary_species_23 = 3.568e-4 #o
+unit_test.primary_species_24 = 1e-40    #ch4
+unit_test.primary_species_25 = 1e-40    #oh+
+unit_test.primary_species_26 = 1e-40    #oh
+unit_test.primary_species_27 = 1e-40    #h2o+
+unit_test.primary_species_28 = 1e-40    #h2o
+unit_test.primary_species_29 = 1e-40    #h3o+
+unit_test.primary_species_30 = 1e-40    #co+
+unit_test.primary_species_31 = 1e-40    #co
+unit_test.primary_species_32 = 1e-40    #o2+
+unit_test.primary_species_33 = 1e-40    #o2
+unit_test.primary_species_34 = 1e-40    #co2
 
 # integrator runtime parameters
-# are we using primordial chemistry? then we use number densities
+# are we using metal chemistry? then we use number densities
 integrator.use_number_densities = 1
 # we do not want to subtract the internal energy
 integrator.subtract_internal_energy = 0
 # we do not want to clip species between 0 and 1
 integrator.do_species_clip = 0
 # minimum positive value of number densities 
-integrator.SMALL_X_SAFE = 1e-100
+integrator.SMALL_X_SAFE = 0 #1e-100
 integrator.burner_verbose = 0
 # do you want to use the jacobian calculated in a previous step?
 integrator.use_jacobian_caching = 0
@@ -69,10 +70,12 @@ integrator.jacobian = 2
 integrator.renormalize_abundances = 0
 # tolerances
 integrator.rtol_spec = 1.0e-4
-integrator.atol_spec = 1.0e-4
+integrator.atol_spec = 1.0e-40
 
 #assumed redshift for Pop III star formation
 network.redshift = 0.0
+network.metallicity = 0.0
+network.dust2gas_ratio = 0.0
 
 # amrex runtime parameters
 # these params help debug the code