diff --git a/Exec/LyA/inputs b/Exec/LyA/inputs index e16f4b7c..539fa4a7 100644 --- a/Exec/LyA/inputs +++ b/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 diff --git a/Exec/LyA_AGN/inputs b/Exec/LyA_AGN/inputs index 6b9163a1..4f95a3b0 100644 --- a/Exec/LyA_AGN/inputs +++ b/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 diff --git a/Exec/Scaling/inputs b/Exec/Scaling/inputs index 118b807e..48808055 100644 --- a/Exec/Scaling/inputs +++ b/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 diff --git a/Source/EOS/atomic_rates.f90 b/Source/EOS/atomic_rates.f90 index fbd9af7d..425b877c 100644 --- a/Source/EOS/atomic_rates.f90 +++ b/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 diff --git a/Source/HeatCool/f_rhs.f90 b/Source/HeatCool/f_rhs.f90 index 92be3a4e..d1d41c8d 100644 --- a/Source/HeatCool/f_rhs.f90 +++ b/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