Commits

Britton Smith  committed 610cff8

Cloudy cooling now only supports the latest tables that are either 2D
(density, temperature) or 3D (density, redshift, temperature).
Additionally, the final redshift slice of the 3D data is now the
collisional 2D data, so simulations don't have to be started with the
2D data, run up to the point where the background turns on, and
restarted with the 3D data. Simulations can be run straight from the
beginning with the 3D data.

  • Participants
  • Parent commits be6e382

Comments (0)

Files changed (4)

File src/clib/chemistry_data.h

 #ifndef __CHEMISTRY_DATA_H__
 #define __CHEMISTRY_DATA_H__
-#define CLOUDY_COOLING_MAX_DIMENSION 5
+#define CLOUDY_COOLING_MAX_DIMENSION 3
 
 typedef struct 
 {

File src/clib/cool1d_cloudy.F

 
 !  Locals
 
-      integer i, q, zindex, zmidpt, zhighpt
+      integer i, q, zindex, zmidpt, zhighpt, get_heat
       real dclPar(clGridRank), inv_log10, log10_tCMB
+      logical end_int
 
 !  Slice locals
 
 !\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\/////////////////////////////////
 !=======================================================================
 
+      end_int = .false.
+      get_heat = iClHeat
+
       inv_log10 = 1.d0 / log(10.d0)
       log10_tCMB = log10(comp2)
 
                if (zr .le. clPar2(1)) then
                   zindex = 1
                else if (zr .ge. clPar2(clGridDim(2)-1)) then
-                  zindex = clGridDim(2) - 1
+                  zindex = clGridDim(2)
+                  end_int = .true.
+                  get_heat = 0
+               else if (zr .ge. clPar2(clGridDim(2)-2)) then
+                  zindex = clGridDim(2) - 2
                else
                   zindex = 1
-                  zhighpt = clGridDim(2)
+                  zhighpt = clGridDim(2) - 2
                   do while ((zhighpt - zindex) .gt. 1)
                      zmidpt = int((zhighpt + zindex) / 2)
                      if (zr .ge. clPar2(zmidpt)) then
                   edot_met(i) = edot_met(i) + 10.d0**log_cool_cmb(i)
                endif
 
-               if (iClHeat .eq. 1) then
+               if (get_heat .eq. 1) then
                   call interpolate_1D(log10tem(i), clGridDim, clPar1, 
      &                 dclPar(1), clDataSize, clHeating, 
      &                 log_heat(i))
                   edot_met(i) = edot_met(i) + 10.d0**log_cool_cmb(i)
                endif
 
-               if (iClHeat .eq. 1) then
+               if (get_heat .eq. 1) then
                call interpolate_2D(log_n_h(i), log10tem(i), clGridDim,
      &              clPar1, dclPar(1), clPar2, dclPar(2),
      &              clDataSize, clHeating, log_heat(i))
      &              clPar1, dclPar(1), 
      &              clPar2, zindex,
      &              clPar3, dclPar(3),
-     &              clDataSize, clCooling, log_cool(i))
+     &              clDataSize, clCooling, 
+     &              end_int, log_cool(i))
                edot_met(i) = -10.d0**log_cool(i)
 
 !     Ignore CMB term if T >> T_CMB
      &                 clPar1, dclPar(1),
      &                 clPar2, zindex,
      &                 clPar3, dclPar(3),
-     &                 clDataSize, clCooling, log_cool_cmb(i))
+     &                 clDataSize, clCooling, 
+     &                 end_int, log_cool_cmb(i))
                   edot_met(i) = edot_met(i) + 10.d0**log_cool_cmb(i)
                endif
 
-               if (iClHeat .eq. 1) then
+               if (get_heat .eq. 1) then
                   call interpolate_3Dz(log_n_h(i), zr, log10tem(i),
      &                 clGridDim,
      &                 clPar1, dclPar(1),
      &                 clPar2, zindex,
      &                 clPar3, dclPar(3),
-     &                 clDataSize, clHeating, log_heat(i))
+     &                 clDataSize, clHeating, 
+     &                 end_int, log_heat(i))
                   edot_met(i) = edot_met(i) + 10.d0**log_heat(i)
                endif
 
       end
 !=======================================================================
 !////////////////////  SUBROUTINE INTERPOLATE_3Dz  \\\\\\\\\\\\\\\\\\\\\
-
+!
+!     Similar to interpolate_3D except index2 is calculated
+!     ahead of time because it is the redshift and will not 
+!     change for the entire grid.
+!
       subroutine interpolate_3Dz(input1, input2, input3, gridDim,
      &     gridPar1, dgridPar1,
      &     gridPar2, index2,
      &     gridPar3, dgridPar3,
-     &     dataSize, dataField, value)
+     &     dataSize, dataField, 
+     &     end_int, value)
 
 !  General Arguments
 
 
 !  Locals
 
+      logical end_int
       integer index1, index3, int_index, q, w
       real slope, value3(2), value2(2)
 
 !\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\/////////////////////////////////
 !=======================================================================
 
+      if (end_int) then
+         call interpolate_2Df3D(input1, input2, 
+     &        input3, gridDim,
+     &        gridPar1, dgridPar1,
+     &        gridPar2, index2,
+     &        gridPar3, dgridPar3,
+     &        dataSize, dataField, 
+     &        value)
+         return
+      endif
+
 !     Calculate interpolation indices
 
       index1 = min(gridDim(1)-1,
       end
 
 !=======================================================================
-!////////////////////  SUBROUTINE INTERPOLATE_4D  \\\\\\\\\\\\\\\\\\\\\\
-
-      subroutine interpolate_4D(input1, input2, input3, input4, 
-     &     gridDim,
+!///////////////////  SUBROUTINE INTERPOLATE_2Df3D  \\\\\\\\\\\\\\\\\\\\
+!
+!     Interpolation in 2 dimensions but with a 3D grid.
+!     This is used for interpolating from just the last 
+!     slice in the datacube before the redshift where 
+!     the UV background turns on.
+!
+      subroutine interpolate_2Df3D(input1, input2, input3, gridDim,
      &     gridPar1, dgridPar1,
-     &     gridPar2, dgridPar2,
+     &     gridPar2, index2,
      &     gridPar3, dgridPar3,
-     &     gridPar4, dgridPar4,
-     &     dataSize, dataField, value)
+     &     dataSize, dataField, 
+     &     value)
 
 !  General Arguments
 
-      integer dataSize
-      integer gridDim(4)
-      real input1, input2, input3, input4, value
+      integer dataSize, index2
+      integer gridDim(3)
+      real input1, input2, input3, value
       real gridPar1(gridDim(1)), dgridPar1,
-     &     gridPar2(gridDim(2)), dgridPar2,
-     &     gridPar3(gridDim(3)), dgridPar3,
-     &     gridPar4(gridDim(4)), dgridPar4
+     &     gridPar2(gridDim(2)),
+     &     gridPar3(gridDim(3)), dgridPar3
       real dataField(dataSize)
 
 !  Locals
 
-      integer index1, index2, index3, index4, int_index, q, w, e
-      real slope, value4(2), value3(2), value2(2)
+      logical end_int
+      integer index1, index3, int_index, q
+      real slope, value3(2), value2(2)
 
 !\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\/////////////////////////////////
 !=======================================================================
 
       index1 = min(gridDim(1)-1,
      &     max(1,int((input1-gridPar1(1))/dgridPar1)+1))
-      index2 = min(gridDim(2)-1,
-     &     max(1,int((input2-gridPar2(1))/dgridPar2)+1))
       index3 = min(gridDim(3)-1,
      &     max(1,int((input3-gridPar3(1))/dgridPar3)+1))
-      index4 = min(gridDim(4)-1,
-     &     max(1,int((input4-gridPar4(1))/dgridPar4)+1))
 
       do q=1, 2
 
-         do w=1, 2
-
-            do e=1, 2
-
-!     interpolate over parameter 4
-
-               int_index = (((q+index1-2) * gridDim(2) + (w+index2-2)) * 
-     &              gridDim(3) + (e+index3-2)) * gridDim(4) + index4
-
-               slope = (dataField(int_index+1) - dataField(int_index)) /
-     &              (gridPar4(index4+1) - gridPar4(index4))
-
-               value4(e) = (input4 - gridPar4(index4)) * slope + 
-     &              dataField(int_index)
-
-            enddo
-
 !     interpolate over parameter 3
 
-            slope = (value4(2) - value4(1)) /
+            int_index = ((q+index1-2) * gridDim(2) + (index2-1)) * 
+     &           gridDim(3) + index3
+
+            slope = (dataField(int_index+1) - dataField(int_index)) /
      &           (gridPar3(index3+1) - gridPar3(index3))
 
-            value3(w) = (input3 - gridPar3(index3)) * slope +
-     &           value4(1)
-
-         enddo
-
-!     interpolate over parameter 2
-
-         slope = (value3(2) - value3(1)) /
-     &        (gridPar2(index2+1) - gridPar2(index2))
-
-         value2(q) = (input2 - gridPar2(index2)) * slope + value3(1)
+            value3(q) = (input3 - gridPar3(index3)) * slope +
+     &           dataField(int_index)
 
       enddo
 
 !     interpolate over parameter 1
 
-      slope = (value2(2) - value2(1)) /
+      slope = (value3(2) - value3(1)) /
      &     (gridPar1(index1+1) - gridPar1(index1))
 
-      value = (input1 - gridPar1(index1)) * slope + value2(1)
+      value = (input1 - gridPar1(index1)) * slope + value3(1)
 
       return
       end
-
-!=======================================================================
-!////////////////////  SUBROUTINE INTERPOLATE_5D  \\\\\\\\\\\\\\\\\\\\\\
-
-      subroutine interpolate_5D(input1, input2, input3, input4, input5,
-     &     gridDim,
-     &     gridPar1, dgridPar1,
-     &     gridPar2, dgridPar2,
-     &     gridPar3, dgridPar3,
-     &     gridPar4, dgridPar4,
-     &     gridPar5, dgridPar5,
-     &     dataSize, dataField, value)
-
-!  General Arguments
-
-      integer dataSize
-      integer gridDim(4)
-      real input1, input2, input3, input4, input5, value
-      real gridPar1(gridDim(1)), dgridPar1,
-     &     gridPar2(gridDim(2)), dgridPar2,
-     &     gridPar3(gridDim(3)), dgridPar3,
-     &     gridPar4(gridDim(4)), dgridPar4,
-     &     gridPar5(gridDim(5)), dgridPar5
-      real dataField(dataSize)
-
-!  Locals
-
-      integer index1, index2, index3, index4, index5, 
-     &     int_index, q, w, e, r, midPt, highPt
-      real slope, value5(2), value4(2), value3(2), value2(2)
-
-!\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\/////////////////////////////////
-!=======================================================================
-
-!     Calculate interpolation indices
-
-      index1 = min(gridDim(1)-1,
-     &     max(1,int((input1-gridPar1(1))/dgridPar1)+1))
-      index2 = min(gridDim(2)-1,
-     &     max(1,int((input2-gridPar2(1))/dgridPar2)+1))
-      index3 = min(gridDim(3)-1,
-     &     max(1,int((input3-gridPar3(1))/dgridPar3)+1))
-#define INDEX_4_BISECTION
-#ifdef INDEX_4_BISECTION
-!     get index 4 with bisection, since not evenly spaced
-      if (input4 .le. gridPar4(1)) then
-         index4 = 1
-      else if (input4 .ge. gridPar4(gridDim(4)-1)) then
-         index4 = gridDim(4) - 1
-      else
-         index4 = 1
-         highPt = gridDim(4)
-         do while ((highPt - index4) .gt. 1)
-            midPt = int((highPt + index4) / 2)
-            if (input4 .ge. gridPar4(midPt)) then
-               index4 = midPt
-            else
-               highPt = midPt
-            endif
-         enddo
-      endif
-#else
-      index4 = min(gridDim(4)-1,
-     &     max(1,int((input4-gridPar4(1))/dgridPar4)+1))
-#endif /* INDEX_4_BISECTION */
-      index5 = min(gridDim(5)-1,
-     &     max(1,int((input5-gridPar5(1))/dgridPar5)+1))
-
-      do q=1, 2
-
-         do w=1, 2
-
-            do e=1, 2
-
-               do r=1, 2
-
-!     interpolate over parameter 5
-
-                  int_index = ((((q+index1-2) * gridDim(2) + 
-     &                 (w+index2-2)) * gridDim(3) + (e+index3-2)) * 
-     &                 gridDim(4) + (r+index4-2)) * gridDim(5) +
-     &                 index5
-
-                  slope = (dataField(int_index+1) - 
-     &                 dataField(int_index)) /
-     &                 (gridPar5(index5+1) - gridPar5(index5))
-
-                  value5(r) = (input5 - gridPar5(index5)) * slope +
-     &                 dataField(int_index)
-
-               enddo
-
-!     interpolate over parameter 4
-
-               slope = (value5(2) - value5(1)) /
-     &              (gridPar4(index4+1) - gridPar4(index4))
-
-               value4(e) = (input4 - gridPar4(index4)) * slope +
-     &              value5(1)
-
-            enddo
-
-!     interpolate over parameter 3
-
-            slope = (value4(2) - value4(1)) /
-     &           (gridPar3(index3+1) - gridPar3(index3))
-
-            value3(w) = (input3 - gridPar3(index3)) * slope +
-     &           value4(1)
-
-         enddo
-
-!     interpolate over parameter 2
-
-         slope = (value3(2) - value3(1)) /
-     &        (gridPar2(index2+1) - gridPar2(index2))
-
-         value2(q) = (input2 - gridPar2(index2)) * slope +
-     &        value3(1)
-
-      enddo
-
-!     interpolate over parameter 1
-
-      slope = (value2(2) - value2(1)) /
-     &     (gridPar1(index1+1) - gridPar1(index1))
-
-      value = (input1 - gridPar1(index1)) * slope + value2(1)
-
-      return
-      end

File src/clib/initialize_cloudy_data.C

 /            from file.
 /  Rank = 1: interpolate over temperature.
 /  Rank = 2: interpolate over density and temperature.
-/  Rank = 3: interpolate over density, metallicity, and temperature.
-/  Rank = 4: interpolate over density, metallicity, electron fraction, 
-/            and temperature.
-/  Rank = 5: interpolate over density, metallicity, electron fraction, 
-/            redshift (for Haardt Madau background), and temperature.
+/  Rank = 3: interpolate over density, redshift, and temperature.
 /
 /  RETURNS:
 /    SUCCESS or FAIL

File src/python/examples/cooling_rate.py

 my_chemistry.with_radiative_cooling = 0
 my_chemistry.primordial_chemistry = 3
 my_chemistry.metal_cooling = 1
-my_chemistry.cloudy_table_file = "solar_cie_c10.h5";
+my_chemistry.cloudy_table_file = "hm_2011_plus.h5";
+my_chemistry.include_metal_heating = 1
 
 my_chemistry.comoving_coordinates = 0
-my_chemistry.density_units = 1.67e-24
+my_chemistry.density_units = mass_h
 my_chemistry.length_units = 1.0e16
 my_chemistry.time_units = yr_to_s * 1e6
 my_chemistry.a_units = 1.0
                   "H2I", "H2II", "DI", "DII", "HDI", "de"]:
         fc_last[field] = np.copy(fc[field])
     solve_chemistry(fc, a_value, dt)
-    fc["energy"] = np.copy(initial_energy)    
+    fc["energy"] = np.copy(initial_energy)   
     if check_convergence(fc, fc_last):
         break
     my_time += dt
 
 calculate_temperature(fc)
 calculate_cooling_time(fc, a_value)
+  
+xbase1 = my_chemistry.length_units/(a_value*my_chemistry.a_units)
+dbase1   = my_chemistry.density_units*(a_value*my_chemistry.a_units)**3
+cool_unit = (my_chemistry.a_units**5 * xbase1**2 * mass_h**2) / \
+  (my_chemistry.time_units**3 * dbase1)
 
-cooling_rate = fc["density"] * my_chemistry.density_units * \
-  fc["energy"] * energy_units / \
-  (fc["cooling_time"] * my_chemistry.time_units)
-
+cooling_rate = cool_unit * fc["density"] * fc["energy"] / \
+  fc["cooling_time"]
+  
 from matplotlib import pyplot
 pyplot.loglog(fc['temperature'], cooling_rate)
 pyplot.xlabel('T [K]')
 pyplot.ylabel('$\\Lambda$ [erg s$^{-1}$ cm$^{-3}$]')
-pyplot.ylim(1e-30, 1e-21)
+#pyplot.ylim(1e-30, 1e-21)
 output_file = 'cooling_rate.png'
 print "Writing %s." % output_file
 pyplot.savefig(output_file)