Add correct, newer gen eq table script!

clairekope · Jun 29, 2023 · 22c9d3b · 22c9d3b
1 parent ca4a04c
commit 22c9d3b
Show file tree

Hide file tree

Showing 2 changed files with 33 additions and 47 deletions.
diff --git 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,16 +44,15 @@
 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
 
 # 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
@@ -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<EquilibriumTable.dim_size*EquilibriumTable.dim_size; ++i){
-  //   if (MultiSpecies) {
-  //     EquilibriumTable.HI[i] /= DensityUnits;
-  //     EquilibriumTable.HII[i] /= DensityUnits;
-  //     EquilibriumTable.HeI[i] /= DensityUnits;
-  //     EquilibriumTable.HeII[i] /= DensityUnits;
-  //     EquilibriumTable.HeIII[i] /= DensityUnits;
-  //     EquilibriumTable.de[i] /= DensityUnits;
-  //     if (MultiSpecies > 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;
 }