Added density dependent UVB heating option

Exec/LyA/inputs

@@ -35,8 +35,8 @@ amr.data_log = runlog
 gravity.gravity_type = PoissonGrav
 gravity.no_sync      = 1
 gravity.no_composite = 1
-gravity.solve_with_cpp = 1
-gravity.solve_with_hpgmg = 0
+gravity.solve_with_cpp = 0
+gravity.solve_with_hpgmg = 1
 
 mg.bottom_solver = 4
 
@@ -71,9 +71,13 @@ nyx.comoving_OmM = 0.275
 nyx.comoving_OmB = 0.046
 nyx.comoving_h   = 0.702d0
 
-# Reionization
+# UVB and reionization
 nyx.inhomo_reion     = 0
+nyx.inhomo_zhi_file  = "zhi.bin"
+nyx.inhomo_grid      = 512
 nyx.uvb_rates_file   = "TREECOOL_middle"
+nyx.uvb_density_A    = 1.0
+nyx.uvb_density_B    = 0.0
 nyx.reionization_zHI_flash   = -1.0
 nyx.reionization_zHeII_flash = -1.0
 nyx.reionization_T_zHI       = 2.0e4

Exec/LyA_AGN/inputs

@@ -38,8 +38,8 @@ amr.data_log = runlog
 gravity.gravity_type = PoissonGrav
 gravity.no_sync      = 1
 gravity.no_composite = 1
-gravity.solve_with_cpp = 1
-gravity.solve_with_hpgmg = 0
+gravity.solve_with_cpp = 0
+gravity.solve_with_hpgmg = 1
 
 mg.bottom_solver = 4
 
@@ -74,9 +74,13 @@ nyx.comoving_OmM = 0.275
 nyx.comoving_OmB = 0.046
 nyx.comoving_h   = 0.702d0
 
-# Reionization
+# UVB and reionization
 nyx.inhomo_reion     = 0
+nyx.inhomo_zhi_file  = "zhi.bin"
+nyx.inhomo_grid      = 512
 nyx.uvb_rates_file   = "TREECOOL_middle"
+nyx.uvb_density_A    = 1.0
+nyx.uvb_density_B    = 0.0
 nyx.reionization_zHI_flash   = -1.0
 nyx.reionization_zHeII_flash = -1.0
 nyx.reionization_T_zHI       = 2.0e4

Exec/Scaling/inputs

@@ -35,8 +35,8 @@ amr.data_log = runlog
 gravity.gravity_type = PoissonGrav
 gravity.no_sync      = 1
 gravity.no_composite = 1
-gravity.solve_with_cpp = 1
-gravity.solve_with_hpgmg = 0
+gravity.solve_with_cpp = 0
+gravity.solve_with_hpgmg = 1
 
 mg.bottom_solver = 4
 
@@ -73,9 +73,13 @@ nyx.comoving_OmM = 0.275
 nyx.comoving_OmB = 0.046
 nyx.comoving_h   = 0.702d0
 
-# Reionization
+# UVB and reionization
 nyx.inhomo_reion     = 0
+nyx.inhomo_zhi_file  = "zhi.bin"
+nyx.inhomo_grid      = 512
 nyx.uvb_rates_file   = "TREECOOL_middle"
+nyx.uvb_density_A    = 1.0
+nyx.uvb_density_B    = 0.0
 nyx.reionization_zHI_flash   = -1.0
 nyx.reionization_zHeII_flash = -1.0
 nyx.reionization_T_zHI       = 2.0e4

Source/EOS/atomic_rates.f90

@@ -40,6 +40,8 @@ module atomic_rates_module
 
   real(rt), parameter, public :: MPROTON = 1.6726231d-24, BOLTZMANN = 1.3806e-16
 
+  real(rt), public, save :: uvb_density_A = 1.0d0, uvb_density_B = 0.0d0, mean_rhob
+
   ! Note that XHYDROGEN can be set by a call to set_xhydrogen which now
   ! lives in set_method_params.
   real(rt), public :: XHYDROGEN = 0.76d0
@@ -50,6 +52,9 @@ module atomic_rates_module
       subroutine fort_tabulate_rates() bind(C, name='fort_tabulate_rates')
       use parallel, only: parallel_ioprocessor
       use amrex_parmparse_module
+      use bl_constants_module, only: M_PI
+      use fundamental_constants_module, only: Gconst
+      use comoving_module, only: comoving_h,comoving_OmB
       use reion_aux_module, only: zhi_flash, zheii_flash, T_zhi, T_zheii, &
                                   flash_h, flash_he, inhomogeneous_on
 
@@ -70,6 +75,8 @@ subroutine fort_tabulate_rates() bind(C, name='fort_tabulate_rates')
          call amrex_parmparse_build(pp, "nyx")
          call pp%query("inhomo_reion"             , inhomo_reion)
          call pp%query("uvb_rates_file"           , file_in)
+         call pp%query("uvb_density_A"            , uvb_density_A)
+         call pp%query("uvb_density_B"            , uvb_density_B)
          call pp%query("reionization_zHI_flash"   , zhi_flash)
          call pp%query("reionization_zHeII_flash" , zheii_flash)
          call pp%query("reionization_T_zHI"       , T_zhi)
@@ -82,7 +89,15 @@ subroutine fort_tabulate_rates() bind(C, name='fort_tabulate_rates')
             print*, '    reionization_zHeII_flash   = ', zheii_flash
             print*, '    reionization_T_zHI         = ', T_zhi
             print*, '    reionization_T_zHeII       = ', T_zheii
-         endif
+
+            print*, 'TABULATE_RATES: rho-dependent heating parameters are:'
+            print*, '    A       = ', uvb_density_A
+            print*, '    B       = ', uvb_density_B
+            print*, '    UVB heating rates will be multiplied by A*(rho/rho_mean)**B'
+        endif
+
+        ! Save mean density (in code units) for density-dependent heating
+        mean_rhob = comoving_OmB * 3.d0*(comoving_h*100.d0)**2 / (8.d0*M_PI*Gconst)
 
          ! Set options in reion_aux_module
          !   Hydrogen reionization

Source/HeatCool/f_rhs.f90

@@ -7,6 +7,7 @@ subroutine f_rhs(num_eq, time, e_in, energy, rpar, ipar)
       use eos_module, only: iterate_ne
       use atomic_rates_module, ONLY: TCOOLMIN, TCOOLMAX, NCOOLTAB, deltaT, &
                                      MPROTON, XHYDROGEN, &
+                                     uvb_density_A, uvb_density_B, mean_rhob, &
                                      BetaH0, BetaHe0, BetaHep, Betaff1, Betaff4, &
                                      RecHp, RecHep, RecHepp, &
                                      eh0, ehe0, ehep
@@ -26,7 +27,7 @@ subroutine f_rhs(num_eq, time, e_in, energy, rpar, ipar)
       real(rt) :: ahp, ahep, ahepp, ad, geh0, gehe0, gehep
       real(rt) :: bh0, bhe0, bhep, bff1, bff4, rhp, rhep, rhepp
       real(rt) :: lambda_c, lambda_ff, lambda, heat
-      real(rt) :: rho, U, a
+      real(rt) :: rho, U, a, rho_heat
       real(rt) :: nh, nh0, nhp, nhe0, nhep, nhepp
       integer :: j
 
@@ -98,6 +99,8 @@ subroutine f_rhs(num_eq, time, e_in, energy, rpar, ipar)
 
       ! Heating terms
       heat = JH_vode*nh0*eh0 + JH_vode*nhe0*ehe0 + JHe_vode*nhep*ehep
+      rho_heat = uvb_density_A * (rho_vode/mean_rhob)**uvb_density_B
+      heat = rho_heat*heat
 
       ! Convert back to code units
       ne_vode     = ne_vode / nh