From b7875cdb9561d123167e738abb8287d9e6dc76d4 Mon Sep 17 00:00:00 2001 From: Piyush Sharda Date: Fri, 23 Aug 2024 11:58:07 +0200 Subject: [PATCH] make test similar to python --- unit_test/burn_cell_metal_chem/_parameters | 3 +- unit_test/burn_cell_metal_chem/burn_cell.H | 25 +++++- .../burn_cell_metal_chem/inputs_metal_chem | 79 ++++++++++--------- 3 files changed, 64 insertions(+), 43 deletions(-) diff --git a/unit_test/burn_cell_metal_chem/_parameters b/unit_test/burn_cell_metal_chem/_parameters index 42fd7be3b..4d602783f 100644 --- a/unit_test/burn_cell_metal_chem/_parameters +++ b/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) diff --git a/unit_test/burn_cell_metal_chem/burn_cell.H b/unit_test/burn_cell_metal_chem/burn_cell.H index 82cf4055b..c5dc598ed 100644 --- a/unit_test/burn_cell_metal_chem/burn_cell.H +++ b/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]; diff --git a/unit_test/burn_cell_metal_chem/inputs_metal_chem b/unit_test/burn_cell_metal_chem/inputs_metal_chem index 666880ee7..483e55b42 100644 --- a/unit_test/burn_cell_metal_chem/inputs_metal_chem +++ b/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