Michael Kuhlen avatar Michael Kuhlen committed cf7a1a5 Merge

Merged in my UVbackground mods.

Comments (0)

Files changed (8)

src/clib/Make.mach.linux-gnu

 # Machine-dependent defines
 #-----------------------------------------------------------------------
 
-MACH_DEFINES   = -DLINUX -DH5_USE_16_API 
+MACH_DEFINES   = -DLINUX -DH5_USE_16_API -fPIC
 
 #-----------------------------------------------------------------------
 # Compiler flag settings
 #-----------------------------------------------------------------------
 
 
-MACH_CPPFLAGS = -P -traditional 
+MACH_CPPFLAGS = -P -traditional
 MACH_CFLAGS   = 
-MACH_CXXFLAGS =
+MACH_CXXFLAGS = 
 MACH_FFLAGS   = -fno-second-underscore -ffixed-line-length-132
 MACH_F90FLAGS = -fno-second-underscore
 MACH_LDFLAGS  = 

src/clib/freefall.C

     return FAIL;
   }
 
-  // my_chemistry.UVbackground = 1;
-  // my_chemistry.UVbackground_file = (char*) "UVB_rates_HM2012.hdf5";
-  // if (initialize_UVbackground_data(my_chemistry, my_units, a_value) == FAIL) {
-  //   fprintf(stderr, "Error in initialize_UVbackground_data.\n");
-  //   return FAIL;
-  // }
+  my_chemistry.UVbackground = 1;
+  my_chemistry.UVbackground_file = (char*) "UVB_rates_HM2012.hdf5";
+  if (initialize_UVbackground_data(my_chemistry) == FAIL) {
+    fprintf(stderr, "Error in initialize_UVbackground_data.\n");
+    return FAIL;
+  }
 
-  // update_UVbackground_rates(my_chemistry, my_units, a_value);
+  update_UVbackground_rates(my_chemistry, my_units, a_value);
 
   gr_float *density, *energy, *x_velocity, *y_velocity, *z_velocity;
   gr_float *HI_density, *HII_density, *HM_density,

src/clib/grackle.h

 int initialize_chemistry_data(chemistry_data &my_chemistry,
                               code_units &my_units, gr_float a_value);
 
-int initialize_UVbackground_data(chemistry_data &my_chemistry,
-				 code_units &my_units, gr_float a_value);
+int initialize_UVbackground_data(chemistry_data &my_chemistry);
 int update_UVbackground_rates(chemistry_data &my_chemistry,
 			      code_units &my_units, float a_value);
 

src/clib/initialize_UVbackground_data.C

 
 
 // Initialize UV Background data
-int initialize_UVbackground_data(chemistry_data &my_chemistry,
-				 code_units &my_units, gr_float a_value)
+int initialize_UVbackground_data(chemistry_data &my_chemistry)
 {
   gr_int Nz, i;
 
   H5Fclose(file_id);
 
 
-  // Now convert the rates to code units.
-
-  /* Get conversion units. */
-  double tbase1 = my_units.time_units;
-  double xbase1 = my_units.length_units/(a_value * my_units.a_units);
-  double dbase1 = my_units.density_units*POW(a_value * my_units.a_units, 3);
-  double mh     = 1.67e-24;
-  double CoolingUnits = (POW(my_units.a_units, 5) * xbase1*xbase1 * mh*mh) /
-    (POW(tbase1, 3) * dbase1);
-  
-  for(i=0;i<Nz;i++) {
-    my_chemistry.UVbackground_table.k24[i] *= my_units.time_units;
-    my_chemistry.UVbackground_table.k25[i] *= my_units.time_units;
-    my_chemistry.UVbackground_table.k26[i] *= my_units.time_units;
-
-    if (my_chemistry.primordial_chemistry > 1) {
-      my_chemistry.UVbackground_table.k27[i] *= my_units.time_units;
-      my_chemistry.UVbackground_table.k28[i] *= my_units.time_units;
-      my_chemistry.UVbackground_table.k29[i] *= my_units.time_units;
-      my_chemistry.UVbackground_table.k30[i] *= my_units.time_units;
-      my_chemistry.UVbackground_table.k31[i] *= my_units.time_units;
-    }
-
-    my_chemistry.UVbackground_table.piHI[i] /= CoolingUnits;
-    my_chemistry.UVbackground_table.piHeII[i] /= CoolingUnits;
-    my_chemistry.UVbackground_table.piHeI[i] /= CoolingUnits;
-  }
-
-
   // Get min/max of redshift vector
   my_chemistry.UVbackground_table.zmin = my_chemistry.UVbackground_table.z[0];
   my_chemistry.UVbackground_table.zmax = my_chemistry.UVbackground_table.z[Nz-1];

src/clib/update_UVbackground_rates.C

 
 #include <string.h>
 #include <stdio.h>
+#include <stdlib.h>
 #include <math.h>
 
 #include "ErrorExceptions.h"
        (Redshift > my_chemistry.UVbackground_table.zmax) )
     return SUCCESS;
 
-
-  /* Set units. */
-
   if (!my_units.comoving_coordinates) {
     my_chemistry.UVbackground_redshift_on = Redshift+0.2;
     my_chemistry.UVbackground_redshift_off = 0.0;
     my_chemistry.UVbackground_redshift_fullon = Redshift+0.1;
     my_chemistry.UVbackground_redshift_drop = 0.0;
   }
-    
-  double tbase1 = my_units.time_units;
-  double xbase1 = my_units.length_units/(a_value * my_units.a_units);
-  double dbase1 = my_units.density_units*POW(a_value * my_units.a_units, 3);
-  double mh     = 1.67e-24;
-  double CoolingUnits = (POW(my_units.a_units, 5) * xbase1*xbase1 * mh*mh) /
-                        (POW(tbase1, 3) * dbase1);
 
   /* ------------------------------------------------------------------ */
   /* First, calculate the ramp value, a number between 0 and 1 which
   my_chemistry.piHeI = (Redshift - zvec[index-1]) * slope + my_chemistry.UVbackground_table.piHeI[index-1];
 
 
+  // Now convert the rates to code units.
+  
+  /* Get conversion units. */
+    
+  double tbase1 = my_units.time_units;
+  double xbase1 = my_units.length_units/(a_value * my_units.a_units);
+  double dbase1 = my_units.density_units*POW(a_value * my_units.a_units, 3);
+  double mh     = 1.67262171e-24;
+  double ev2erg = 1.60217653e-12;
+  double CoolingUnits = (POW(my_units.a_units, 5) * xbase1*xbase1 * mh*mh) / (POW(tbase1, 3) * dbase1) / ev2erg;  // compared to Enzo source, there's an additional factor of 1/ev2erg here, because the heating rates are stored as eV/s.
+
+
+  my_chemistry.k24 *= my_units.time_units;
+  my_chemistry.k25 *= my_units.time_units;
+  my_chemistry.k26 *= my_units.time_units;
+  
+  if (my_chemistry.primordial_chemistry > 1) {
+    my_chemistry.k27 *= my_units.time_units;
+    my_chemistry.k28 *= my_units.time_units;
+    my_chemistry.k29 *= my_units.time_units;
+    my_chemistry.k30 *= my_units.time_units;
+    my_chemistry.k31 *= my_units.time_units;
+  }
+    
+  my_chemistry.piHI /= CoolingUnits;
+  my_chemistry.piHeII /= CoolingUnits;
+  my_chemistry.piHeI /= CoolingUnits;
+  
+
   // printf("%e %e %e\n",my_chemistry.k24,my_chemistry.k25,my_chemistry.k26);
-
+  // printf("%e %e %e\n",my_chemistry.piHI,my_chemistry.piHeII,my_chemistry.piHeI);
 
   // Now apply the Ramp factor
   my_chemistry.k24 *= Ramp;
        is the average photon energy in keV, corrected for relativistic
        effects.  Eq.(4) and Eq.(11) of Madau & Efstathiou (1999) */
     
-    my_chemistry.comp_xray = 6.65e-25 * 3.0e10 * 
+    my_chemistry.comp_xray = 4.15e-13 * 3.0e10 * 
       (31.8*POW(1.0+Redshift, 0.3333)/511.0) * 
       (6.3e-5 * 1.6e-12) * 
       POW(1.0 + Redshift, 4) * 

src/python/examples/freefall.py

 
 my_chemistry.initialize(a_value)
 
+
+my_chemistry.UVbackground = 1;
+my_chemistry.UVbackground_file = "UVB_rates_HM2012.hdf5";
+my_chemistry.initialize_UVbackground()
+my_chemistry.update_UVbackground(a_value)
+
+
 fc = FluidContainer(my_chemistry, 1)
 fc["density"][:] = 1.0
 fc["HI"][:] = 0.76 * fc["density"]

src/python/pygrackle/grackle_defs.pxd

         gr_int cie_cooling
         gr_int h2_optical_depth_approximation
         gr_int photoelectric_heating
+        gr_int UVbackground
+        char *UVbackground_file
         # Most of the rest are not user-settable
 
 cdef extern from "code_units.h":
     gr_int initialize_chemistry_data(c_chemistry_data &my_chemistry,
                                   c_code_units &my_units, gr_float a_value)
 
+    gr_int initialize_UVbackground_data(c_chemistry_data &my_chemistry)
+
+    gr_int update_UVbackground_rates(c_chemistry_data &my_chemistry,
+                                  c_code_units &my_units, gr_float a_value)
+
     gr_int c_solve_chemistry "solve_chemistry"(
                 c_chemistry_data &my_chemistry,
                 c_code_units &my_units,

src/python/pygrackle/grackle_wrapper.pyx

     def initialize(self, a_value):
         initialize_chemistry_data(self.data, self.units, a_value)
 
+    def initialize_UVbackground(self):
+        initialize_UVbackground_data(self.data)
+
+    def update_UVbackground(self, a_value):
+        update_UVbackground_rates(self.data, self.units, a_value)
+
     property Gamma:
         def __get__(self):
             return self.data.Gamma
         def __set__(self, val):
             self.data.cloudy_table_file = val
 
+    property UVbackground:
+        def __get__(self):
+            return self.data.UVbackground
+        def __set__(self, val):
+            self.data.UVbackground = val
+
+    property UVbackground_file:
+        def __get__(self):
+            return self.data.UVbackground_file
+        def __set__(self, val):
+            self.data.UVbackground_file = val
+
     property comoving_coordinates:
         def __get__(self):
             return self.units.comoving_coordinates
Tip: Filter by directory path e.g. /media app.js to search for public/media/app.js.
Tip: Use camelCasing e.g. ProjME to search for ProjectModifiedEvent.java.
Tip: Filter by extension type e.g. /repo .js to search for all .js files in the /repo directory.
Tip: Separate your search with spaces e.g. /ssh pom.xml to search for src/ssh/pom.xml.
Tip: Use ↑ and ↓ arrow keys to navigate and return to view the file.
Tip: You can also navigate files with Ctrl+j (next) and Ctrl+k (previous) and view the file with Ctrl+o.
Tip: You can also navigate files with Alt+j (next) and Alt+k (previous) and view the file with Alt+o.