From 22c9d3b7d0d2a5eb545836fc8d23d053aa09437f Mon Sep 17 00:00:00 2001 From: Claire Kopenhafer Date: Thu, 29 Jun 2023 17:31:18 -0400 Subject: [PATCH] Add correct, newer gen eq table script! --- input/gen_equilibrium_table.py | 59 ++++++++++++++++++--------------- src/enzo/ReadEquilibriumTable.C | 21 +----------- 2 files changed, 33 insertions(+), 47 deletions(-) diff --git a/input/gen_equilibrium_table.py b/input/gen_equilibrium_table.py index 4417afaff..994ef272a 100644 --- a/input/gen_equilibrium_table.py +++ b/input/gen_equilibrium_table.py @@ -8,7 +8,7 @@ # As is, this script makes a table spanning 1e-29 to 1e-23 g/cm^3 # and 1e4 to 1e6 K using 50 cells on a side. The constant metallicity is # 0.27 Zsun and the HM2012 UV background is used. The table is calculated for -# redshift 0. The table size and are included in the resulting HDF5 filename. +# redshift 0. The table size and metallicity are included in the resulting HDF5 filename. # ############################################################################### @@ -19,18 +19,22 @@ chemistry_data, \ setup_fluid_container from pygrackle.utilities.physical_constants import \ + cm_per_mpc as Mpc, \ + amu_cgs as m_u, \ + mass_sun_cgs as Msun, \ + keV_per_K as kboltz, \ sec_per_Myr, \ cm_per_kpc -size_1d = 50 # square table; length of each side -z = 0.27 # solar units +# 1e-29 to 1e-23 g/cm**3 +# 6.5e4 to 1.5e6 K -# Table's density & temperature grid is in log space -dens = np.logspace(-29, -23, size_1d) # cgs -temp = np.logspace(4, 6, size_1d) # K -d, t = np.meshgrid(dens, temp) +size_1d = 60 # square table; length of each side +Z = 0.3 # solar units -current_redshift = 0. +dens = np.logspace(np.log10(1e-31), np.log10(1e-23), size_1d) # cgs +temp = np.logspace(np.log10(3e3), np.log10(3e6), size_1d) # K +d, t = np.meshgrid(dens, temp) # Set solver parameters chem = chemistry_data() @@ -40,8 +44,7 @@ chem.metal_cooling = 1 chem.UVbackground = 1 chem.cmb_temperature_floor = 1 -chem.h2_on_dust = 1 -chem.grackle_data_file = "/PATH/TO/GRACKLE/input/CloudyData_UVB=HM2012.h5" +chem.grackle_data_file = b"/PATH/TO/GRACKLE/input/CloudyData_UVB=HM2012.h5" chem.use_specific_heating_rate = 0 chem.use_volumetric_heating_rate = 0 @@ -49,7 +52,7 @@ # Set units chem.comoving_coordinates = 0 # proper units chem.a_units = 1.0 -chem.a_value = 1.0 / (1.0 + current_redshift) / chem.a_units +chem.a_value = 1.0 chem.density_units = 1.67e-27 chem.length_units = 1638.4 * cm_per_kpc chem.time_units = sec_per_Myr @@ -60,22 +63,24 @@ # Call convenience function for setting up a fluid container. # This container holds the solver parameters, units, and fields. fc = setup_fluid_container(chem, - density=d.flatten(), - temperature=t.flatten(), - metal_mass_fraction=0.02014*z, - converge=True) + density=d.ravel(), + temperature=t.ravel(), + metal_mass_fraction=0.01295*Z, + converge=True, + tolerance=1e-7, # default: 0.01 + max_iterations=10000000000) -# Save in cgs -dataset = {'HI' :fc['HI'] *fc.chemistry_data.density_units, - 'HII' :fc['HII'] *fc.chemistry_data.density_units, - 'HeI' :fc['HeI'] *fc.chemistry_data.density_units, - 'HeII' :fc['HeII'] *fc.chemistry_data.density_units, - 'HeIII':fc['HeIII']*fc.chemistry_data.density_units, - 'HM' :fc['HM'] *fc.chemistry_data.density_units, - 'H2I' :fc['H2I'] *fc.chemistry_data.density_units, - 'H2II' :fc['H2II'] *fc.chemistry_data.density_units, - 'de' :fc['de'] *fc.chemistry_data.density_units, - 'metal':fc['metal']*fc.chemistry_data.density_units, +# Save as fraction +dataset = {'HI' :fc['HI'] /fc['density'], + 'HII' :fc['HII'] /fc['density'], + 'HeI' :fc['HeI'] /fc['density'], + 'HeII' :fc['HeII'] /fc['density'], + 'HeIII':fc['HeIII']/fc['density'], + 'HM' :fc['HM'] /fc['density'], + 'H2I' :fc['H2I'] /fc['density'], + 'H2II' :fc['H2II'] /fc['density'], + 'de' :fc['de'] /fc['density'], + 'metal':fc['metal']/fc['density'], 'density': dens, 'temperature': temp} @@ -92,5 +97,5 @@ 'density': 'indexer', 'temperature': 'indexer'} -yt.save_as_dataset(None, f'equilibrium_table_{size_1d}_0{z*100:.0f}-Zsun.h5', dataset, +yt.save_as_dataset(None, f'equilibrium_table_{size_1d}_0{Z*100:.0f}-Zsun.h5', dataset, field_types = field_types) diff --git a/src/enzo/ReadEquilibriumTable.C b/src/enzo/ReadEquilibriumTable.C index e29296076..acffdb9da 100644 --- a/src/enzo/ReadEquilibriumTable.C +++ b/src/enzo/ReadEquilibriumTable.C @@ -422,26 +422,7 @@ int ReadEquilibriumTable(char* name, FLOAT Time) EquilibriumTable.density[i] /= DensityUnits; //EquilibriumTable.temperature[i] /= TemperatureUnits; // keep in K } - // for (int i=0; i 1) { - // EquilibriumTable.HM[i] /= DensityUnits; - // EquilibriumTable.H2I[i] /= DensityUnits; - // EquilibriumTable.H2II[i] /= DensityUnits; - // } - // if (MultiSpecies > 2) { - // EquilibriumTable.DI[i] /= DensityUnits; - // EquilibriumTable.DII[i] /= DensityUnits; - // EquilibriumTable.HDI[i] /= DensityUnits; - // } - // } - // } + return SUCCESS; }