diff --git a/Exec/AMR-density/GNUmakefile b/Exec/AMR-density/GNUmakefile index 2fbfecbb..fc2ee684 100644 --- a/Exec/AMR-density/GNUmakefile +++ b/Exec/AMR-density/GNUmakefile @@ -1,7 +1,6 @@ # AMREX_HOME defines the directory in which we will find all the AMReX code AMREX_HOME ?= ../../../amrex -HPGMG_DIR ?= ../../Util/hpgmg/finite-volume CVODE_LIB_DIR ?= $(CVODE_LIB) # TOP defines the directory in which we will find Source, Exec, etc @@ -24,8 +23,6 @@ DEBUG = FALSE GIMLET = FALSE REEBER = FALSE -USE_HPGMG = FALSE - # physics DIM = 3 USE_GRAV = TRUE diff --git a/Exec/AMR-density/inputs b/Exec/AMR-density/inputs index eb15f41f..cc86f475 100644 --- a/Exec/AMR-density/inputs +++ b/Exec/AMR-density/inputs @@ -155,7 +155,6 @@ amr.plot_file = plt amr.plot_int = -1 amr.plot_nfiles = 64 nyx.plot_z_values = 7.0 6.0 5.0 4.0 3.0 2.0 -particles.write_in_plotfile = 1 amr.plot_vars = density xmom ymom zmom rho_e Temp phi_grav amr.derive_plot_vars = particle_mass_density diff --git a/Exec/AMR-density/inputs.rt b/Exec/AMR-density/inputs.rt index a5e7cc9f..2247a1fe 100644 --- a/Exec/AMR-density/inputs.rt +++ b/Exec/AMR-density/inputs.rt @@ -136,7 +136,6 @@ amr.plot_files_output = 1 amr.plot_file = plt amr.plot_int = 10 amr.plot_nfiles = 64 -particles.write_in_plotfile = 1 amr.plot_vars = density xmom ymom zmom rho_e Temp phi_grav amr.derive_plot_vars = particle_mass_density diff --git a/Exec/AMR-zoom/GNUmakefile b/Exec/AMR-zoom/GNUmakefile index 53a57329..e92b29c1 100644 --- a/Exec/AMR-zoom/GNUmakefile +++ b/Exec/AMR-zoom/GNUmakefile @@ -1,7 +1,6 @@ # AMREX_HOME defines the directory in which we will find all the AMReX code AMREX_HOME ?= ../../../amrex -HPGMG_DIR ?= ../../Util/hpgmg/finite-volume CVODE_LIB_DIR ?= $(CVODE_LIB) # TOP defines the directory in which we will find Source, Exec, etc @@ -24,8 +23,6 @@ DEBUG = FALSE GIMLET = FALSE REEBER = FALSE -USE_HPGMG = FALSE - # physics DIM = 3 USE_GRAV = TRUE diff --git a/Exec/AMR-zoom/inputs b/Exec/AMR-zoom/inputs index 577d7325..4c56540e 100644 --- a/Exec/AMR-zoom/inputs +++ b/Exec/AMR-zoom/inputs @@ -155,7 +155,6 @@ amr.plot_file = plt amr.plot_int = -1 amr.plot_nfiles = 64 nyx.plot_z_values = 7.0 6.0 5.0 4.0 3.0 2.0 -particles.write_in_plotfile = 1 amr.plot_vars = density xmom ymom zmom rho_e Temp phi_grav amr.derive_plot_vars = particle_mass_density diff --git a/Exec/GravityTests/MacLaurin/Tagging_3d.f90 b/Exec/GravityTests/MacLaurin/Tagging_3d.f90 index 0a4d787e..263be699 100644 --- a/Exec/GravityTests/MacLaurin/Tagging_3d.f90 +++ b/Exec/GravityTests/MacLaurin/Tagging_3d.f90 @@ -18,10 +18,11 @@ ! ::: time => problem evolution time ! ::: level => refinement level of this array ! ::: ----------------------------------------------------------- + subroutine tag_overdensity(tag,tagl1,tagl2,tagl3,tagh1,tagh2,tagh3, & set,clear, & den,denl1,denl2,denl3,denh1,denh2,denh3, & - lo,hi,nc,delta,level,avg_den) + lo,hi,nc,domlo,domhi,delta,level,avg_den) use amrex_fort_module, only : rt => amrex_real use probdata_module @@ -30,7 +31,7 @@ subroutine tag_overdensity(tag,tagl1,tagl2,tagl3,tagh1,tagh2,tagh3, & integer set, clear, nc, level integer tagl1,tagl2,tagl3,tagh1,tagh2,tagh3 integer denl1,denl2,denl3,denh1,denh2,denh3 - integer lo(3), hi(3) + integer lo(3), hi(3), domlo(3), domhi(3) integer tag(tagl1:tagh1,tagl2:tagh2,tagl3:tagh3) real(rt) den(denl1:denh1,denl2:denh2,denl3:denh3,nc) real(rt) delta(3), avg_den @@ -74,21 +75,17 @@ end subroutine tag_overdensity ! ::: level => refinement level of this array ! ::: ----------------------------------------------------------- subroutine tag_region(tag,tagl1,tagl2,tagl3,tagh1,tagh2,tagh3, & - set,clear, & - var,varl1,varl2,varl3,varh1,varh2,varh3, & - lo,hi,nd,domlo,domhi, & - delta,xlo,problo,time,level) + set,lo,hi,domlo,domhi, & + delta,xlo,problo,level) use amrex_fort_module, only : rt => amrex_real use probdata_module implicit none - integer :: set, clear, nd, level + integer :: set,level integer :: tagl1,tagl2,tagl3,tagh1,tagh2,tagh3 - integer :: varl1,varl2,varl3,varh1,varh2,varh3 integer :: lo(3), hi(3), domlo(3), domhi(3) integer :: tag(tagl1:tagh1,tagl2:tagh2,tagl3:tagh3) - real(rt) :: var(varl1:varh1,varl2:varh2,varl3:varh3,nd) - real(rt) :: delta(3), xlo(3), problo(3), time + real(rt) :: delta(3), xlo(3), problo(3) integer :: ilo,jlo,klo,ihi,jhi,khi integer :: i,j,k diff --git a/Exec/HDF5_example/GNUmakefile b/Exec/HDF5_example/GNUmakefile index 63167c9e..ff7627d4 100644 --- a/Exec/HDF5_example/GNUmakefile +++ b/Exec/HDF5_example/GNUmakefile @@ -1,7 +1,6 @@ # AMREX_HOME defines the directory in which we will find all the AMReX code AMREX_HOME ?= ../../../amrex -HPGMG_DIR ?= ../../Util/hpgmg/finite-volume CVODE_LIB_DIR ?= $(CVODE_LIB) # TOP defines the directory in which we will find Source, Exec, etc @@ -24,8 +23,6 @@ DEBUG = FALSE GIMLET = FALSE REEBER = FALSE -USE_HPGMG = TRUE - # physics DIM = 3 USE_GRAV = TRUE diff --git a/Exec/HDF5_example/inputs b/Exec/HDF5_example/inputs index bd18bcd7..3eeaba98 100644 --- a/Exec/HDF5_example/inputs +++ b/Exec/HDF5_example/inputs @@ -155,7 +155,6 @@ amr.plot_file = plt amr.plot_int = -1 amr.plot_nfiles = 64 nyx.plot_z_values = 7.0 6.0 5.0 4.0 3.0 2.0 -particles.write_in_plotfile = 1 amr.plot_vars = density xmom ymom zmom rho_e Temp phi_grav amr.derive_plot_vars = particle_mass_density diff --git a/Exec/LyA/GNUmakefile b/Exec/LyA/GNUmakefile index 57922262..77419fdc 100644 --- a/Exec/LyA/GNUmakefile +++ b/Exec/LyA/GNUmakefile @@ -1,7 +1,6 @@ # AMREX_HOME defines the directory in which we will find all the AMReX code AMREX_HOME ?= ../../../amrex -HPGMG_DIR ?= ../../Util/hpgmg/finite-volume CVODE_LIB_DIR ?= $(CVODE_LIB) # TOP defines the directory in which we will find Source, Exec, etc @@ -26,8 +25,6 @@ DEBUG = FALSE GIMLET = FALSE REEBER = FALSE -USE_HPGMG = TRUE - # physics DIM = 3 USE_GRAV = TRUE diff --git a/Exec/LyA/inputs b/Exec/LyA/inputs index 370e2c87..52d9b32b 100644 --- a/Exec/LyA/inputs +++ b/Exec/LyA/inputs @@ -155,7 +155,6 @@ amr.plot_file = plt amr.plot_int = -1 amr.plot_nfiles = 64 nyx.plot_z_values = 7.0 6.0 5.0 4.0 3.0 2.0 -particles.write_in_plotfile = 1 amr.plot_vars = density xmom ymom zmom rho_e Temp phi_grav amr.derive_plot_vars = particle_mass_density diff --git a/Exec/LyA/inputs_gimlet_in_transit.dsc b/Exec/LyA/inputs_gimlet_in_transit.dsc index 561287cf..befbbd73 100644 --- a/Exec/LyA/inputs_gimlet_in_transit.dsc +++ b/Exec/LyA/inputs_gimlet_in_transit.dsc @@ -127,7 +127,6 @@ amr.plot_file = plt amr.plot_int = 20 amr.plot_nfiles = 86 #nyx.plot_z_values = 4.0 3.5 3.0 2.5 2.25 2.0 -#particles.write_in_plotfile = 1 amr.plot_vars = density #amr.derive_plot_vars = particle_mass_density diff --git a/Exec/LyA_AGN/GNUmakefile b/Exec/LyA_AGN/GNUmakefile index 17c80a87..e94354fb 100644 --- a/Exec/LyA_AGN/GNUmakefile +++ b/Exec/LyA_AGN/GNUmakefile @@ -1,7 +1,6 @@ # AMREX_HOME defines the directory in which we will find all the AMReX code AMREX_HOME ?= ../../../amrex -HPGMG_DIR ?= ../../Util/hpgmg/finite-volume CVODE_LIB_DIR ?= $(CVODE_LIB) # TOP defines the directory in which we will find Source, Exec, etc @@ -10,7 +9,7 @@ TOP = ../.. # compilation options COMP = intel # gnu USE_MPI = TRUE -USE_OMP = TRUE +USE_OMP = FALSE PROFILE = TRUE TRACE_PROFILE = FALSE @@ -24,8 +23,6 @@ DEBUG = FALSE GIMLET = FALSE REEBER = TRUE -USE_HPGMG = TRUE - # physics DIM = 3 USE_GRAV = TRUE diff --git a/Exec/LyA_AGN/inputs b/Exec/LyA_AGN/inputs index 25a9ed53..e18160a4 100644 --- a/Exec/LyA_AGN/inputs +++ b/Exec/LyA_AGN/inputs @@ -144,7 +144,6 @@ amr.plot_file = plt amr.plot_int = -1 amr.plot_nfiles = 64 nyx.plot_z_values = 7.0 6.0 5.0 4.0 3.0 2.0 -particles.write_in_plotfile = 1 amr.plot_vars = density xmom ymom zmom rho_e Temp phi_grav amr.derive_plot_vars = particle_mass_density diff --git a/Exec/MiniSB/GNUmakefile b/Exec/MiniSB/GNUmakefile index 488bdf5a..050a66c3 100644 --- a/Exec/MiniSB/GNUmakefile +++ b/Exec/MiniSB/GNUmakefile @@ -1,8 +1,6 @@ # AMREX_HOME defines the directory in which we will find all the BoxLib code AMREX_HOME ?= ../../../amrex -HPGMG_DIR ?= ../../Util/hpgmg/finite-volume - # TOP defines the directory in which we will find Source, Exec, etc TOP = ../.. @@ -15,8 +13,6 @@ USE_OMP = FALSE # Analysis REEBER = FALSE -USE_HPGMG = FALSE - PRECISION = DOUBLE DEBUG = FALSE diff --git a/Exec/RegressionTest/GNUmakefile b/Exec/RegressionTest/GNUmakefile index 9a35b23f..208e8c58 100644 --- a/Exec/RegressionTest/GNUmakefile +++ b/Exec/RegressionTest/GNUmakefile @@ -12,6 +12,7 @@ USE_OMP = FALSE PRECISION = DOUBLE DEBUG = FALSE +DEBUG = TRUE # physics DIM = 3 diff --git a/Exec/RegressionTest/Prob_3d.f90 b/Exec/RegressionTest/Prob_3d.f90 index ee749c1c..f88e2c7a 100644 --- a/Exec/RegressionTest/Prob_3d.f90 +++ b/Exec/RegressionTest/Prob_3d.f90 @@ -81,24 +81,33 @@ subroutine fort_initdata(level,time,lo,hi, & integer i,j,k - do k = lo(3), hi(3) - do j = lo(2), hi(2) - do i = lo(1), hi(1) - state(i,j,k,URHO) = 1.d0 - state(i,j,k,UMX:UMZ) = 0.0d0 - - state(i,j,k,UEINT) = 10.d0 - state(i,j,k,UEDEN) = state(i,j,k,UEINT) + 0.5d0 / state(i,j,k,URHO) * & - (state(i,j,k,UMX)**2 + state(i,j,k,UMY)**2 + state(i,j,k,UMZ)**2) - - if (UFS .gt. -1) then - state(i,j,k,UFS ) = XHYDROGEN - state(i,j,k,UFS+1) = (1.d0 - XHYDROGEN) - end if - - diag_eos(i,j,k,TEMP_COMP) = 1000.d0 - enddo - enddo - enddo + if (ns.gt.1) then + + do k = lo(3), hi(3) + do j = lo(2), hi(2) + do i = lo(1), hi(1) + state(i,j,k,URHO) = 1.d0 + state(i,j,k,UMX:UMZ) = 0.0d0 + + state(i,j,k,UEINT) = 10.d0 + state(i,j,k,UEDEN) = state(i,j,k,UEINT) + 0.5d0 / state(i,j,k,URHO) * & + (state(i,j,k,UMX)**2 + state(i,j,k,UMY)**2 + state(i,j,k,UMZ)**2) + + if (UFS .gt. -1) then + state(i,j,k,UFS ) = XHYDROGEN + state(i,j,k,UFS+1) = (1.d0 - XHYDROGEN) + end if + + diag_eos(i,j,k,TEMP_COMP) = 1000.d0 + enddo + enddo + enddo + + else + + state(lo(1):hi(1),lo(2):hi(2),lo(3):hi(3),URHO) = 1.d0 + + endif + end subroutine fort_initdata diff --git a/Exec/RegressionTest/inputs.mass.nosub b/Exec/RegressionTest/inputs.mass.nosub index bdc97fa9..3d074266 100644 --- a/Exec/RegressionTest/inputs.mass.nosub +++ b/Exec/RegressionTest/inputs.mass.nosub @@ -69,8 +69,6 @@ nyx.particle_init_type = AsciiFile nyx.ascii_particle_file = particle_file.mass particles.particle_output_file = final_particles -particles.write_in_plotfile = 1 - # >>>>>>>>>>>>> PARTICLE AGGREGATION OPTIONS <<<<<<<<<<<<<<<< # "None" "Cell" "Flow" # >>>>>>>>>>>>> PARTICLE AGGREGATION OPTIONS <<<<<<<<<<<<<<<< diff --git a/Exec/RegressionTest/inputs.mass.sub b/Exec/RegressionTest/inputs.mass.sub index b82351ca..aa0def2f 100644 --- a/Exec/RegressionTest/inputs.mass.sub +++ b/Exec/RegressionTest/inputs.mass.sub @@ -71,8 +71,6 @@ nyx.particle_init_type = AsciiFile nyx.ascii_particle_file = particle_file.mass particles.particle_output_file = final_particles -particles.write_in_plotfile = 1 - # >>>>>>>>>>>>> PARTICLE AGGREGATION OPTIONS <<<<<<<<<<<<<<<< # "None" "Cell" "Flow" # >>>>>>>>>>>>> PARTICLE AGGREGATION OPTIONS <<<<<<<<<<<<<<<< diff --git a/Exec/SanityCheck/GNUmakefile b/Exec/SanityCheck/GNUmakefile new file mode 100644 index 00000000..37108945 --- /dev/null +++ b/Exec/SanityCheck/GNUmakefile @@ -0,0 +1,39 @@ +# AMREX_HOME defines the directory in which we will find all the AMReX code +AMREX_HOME ?= ../../../amrex + +CVODE_LIB_DIR ?= $(CVODE_LIB) + +# TOP defines the directory in which we will find Source, Exec, etc +TOP = ../.. + +# compilation options +COMP = gnu +USE_MPI = TRUE +USE_OMP = FALSE + +PROFILE = FALSE +TRACE_PROFILE = FALSE +COMM_PROFILE = FALSE +TINY_PROFILE = FALSE + +PRECISION = DOUBLE +USE_SINGLE_PRECISION_PARTICLES = TRUE +DEBUG = FALSE + +GIMLET = FALSE +REEBER = FALSE + +USE_SDC = TRUE + +# physics +DIM = 3 +USE_GRAV = TRUE +USE_HEATCOOL = TRUE +USE_AGN = FALSE +USE_CVODE = FALSE + +Bpack := ./Make.package +Blocs := . + +include $(TOP)/Exec/Make.Nyx + diff --git a/Exec/SantaBarbara/GNUmakefile b/Exec/SantaBarbara/GNUmakefile index 19dd050e..d6ad140a 100644 --- a/Exec/SantaBarbara/GNUmakefile +++ b/Exec/SantaBarbara/GNUmakefile @@ -1,8 +1,6 @@ # AMREX_HOME defines the directory in which we will find all the BoxLib code AMREX_HOME ?= ../../../amrex -HPGMG_DIR ?= ../../Util/hpgmg/finite-volume - # TOP defines the directory in which we will find Source, Exec, etc TOP = ../.. @@ -15,8 +13,6 @@ PRECISION = DOUBLE USE_SINGLE_PRECISION_PARTICLES = FALSE DEBUG = FALSE -USE_HPGMG = FALSE - # physics DIM = 3 USE_GRAV = TRUE diff --git a/Exec/Scaling/GNUmakefile b/Exec/Scaling/GNUmakefile index 3b67f54c..261db364 100644 --- a/Exec/Scaling/GNUmakefile +++ b/Exec/Scaling/GNUmakefile @@ -1,16 +1,15 @@ # AMREX_HOME defines the directory in which we will find all the AMReX code AMREX_HOME ?= ../../../amrex -HPGMG_DIR ?= ../../Util/hpgmg/finite-volume CVODE_LIB_DIR ?= $(CVODE_LIB) # TOP defines the directory in which we will find Source, Exec, etc TOP = ../.. # compilation options -COMP = intel # gnu -USE_MPI = TRUE -USE_OMP = TRUE +COMP = gnu +USE_MPI = FALSE +USE_OMP = FALSE PROFILE = TRUE TRACE_PROFILE = FALSE @@ -24,8 +23,6 @@ DEBUG = FALSE GIMLET = FALSE REEBER = FALSE -USE_HPGMG = TRUE - # physics DIM = 3 USE_GRAV = TRUE diff --git a/Exec/Scaling/inputs b/Exec/Scaling/inputs index 7fa6a0b8..ed41baf2 100644 --- a/Exec/Scaling/inputs +++ b/Exec/Scaling/inputs @@ -139,7 +139,6 @@ amr.plot_file = plt amr.plot_int = -1 amr.plot_nfiles = 64 nyx.plot_z_values = 7.0 6.0 5.0 4.0 3.0 2.0 -particles.write_in_plotfile = 1 amr.plot_vars = density xmom ymom zmom rho_e Temp phi_grav amr.derive_plot_vars = particle_mass_density diff --git a/README.md b/README.md index c518128b..58b1a437 100644 --- a/README.md +++ b/README.md @@ -66,11 +66,13 @@ For the reaction and thermal rates of the primordial chemical composition gas Lukic, Stark, Nugent, White, Meiksin & Almgren (2015), MNRAS, 446, 3697: http://adsabs.harvard.edu/abs/2015MNRAS.446.3697L -For considerations regarding the synthesis model of the UV background, +For considerations regarding the spatially uniform synthesis model of the UV background, which provides the photo-ionization and photo-heating rates, see Onorbe, Hennawi & Lukic (2017), ApJ, 837, 106: http://adsabs.harvard.edu/abs/2017ApJ...837..106O +We have also implemented non-radiative transfer methods to model inhomogeneous reionization, +the paper is in preparation. ## Output diff --git a/Source/Initialization/Nyx_setup.cpp b/Source/Initialization/Nyx_setup.cpp index 6e77b1f9..6527624f 100644 --- a/Source/Initialization/Nyx_setup.cpp +++ b/Source/Initialization/Nyx_setup.cpp @@ -646,7 +646,9 @@ Nyx::no_hydro_setup() use_const_species, gamma, normalize_species, heat_cool_type, inhomo_reion); +#ifdef HEATCOOL fort_tabulate_rates(); +#endif int coord_type = Geometry::Coord(); fort_set_problem_params(dm, phys_bc.lo(), phys_bc.hi(), Outflow, Symmetry, coord_type); diff --git a/Source/Nyx.H b/Source/Nyx.H index 58284579..e2026439 100644 --- a/Source/Nyx.H +++ b/Source/Nyx.H @@ -179,6 +179,16 @@ public: // static amrex::Real relative_max_change_a; + // + // Absolute change in a allowed in one timestep for fixed delta_a + // + static amrex::Real absolute_max_change_a; + + // + // Positive number means use powers of 2 binning for relative dt + // + static amrex::Real dt_binpow; + // // Old and new times at which "old_a" and "new_a" are defined. // diff --git a/Source/Nyx_F.H b/Source/Nyx_F.H index e8c02546..e87e3c4c 100644 --- a/Source/Nyx_F.H +++ b/Source/Nyx_F.H @@ -18,6 +18,9 @@ extern "C" void fort_integrate_comoving_a_to_z (amrex::Real* old_a, amrex::Real* z_value, amrex::Real* dt); + void fort_integrate_comoving_a_to_a + (amrex::Real* old_a, amrex::Real* a_value, amrex::Real* dt); + void set_simd(const int *simd_width); // void fort_get_omm (amrex::Real* omm); @@ -35,7 +38,7 @@ extern "C" void fort_estdt_comoving_a (amrex::Real* old_a, amrex::Real* new_dummy_a, amrex::Real* dt, amrex::Real* change_allowed, - amrex::Real* final_a, int* dt_modified); + amrex::Real* fixed_da_interval, amrex::Real* final_a, int* dt_modified); void fort_get_method_params(int* HYP_GROW); diff --git a/Source/Nyx_nd.f90 b/Source/Nyx_nd.f90 index be4a5887..e57a5727 100644 --- a/Source/Nyx_nd.f90 +++ b/Source/Nyx_nd.f90 @@ -147,7 +147,6 @@ subroutine fort_set_small_values(average_dens, average_temp, a, & use amrex_fort_module, only : rt => amrex_real use meth_params_module use eos_module - use network, only : nspec, naux implicit none @@ -372,26 +371,27 @@ subroutine fort_set_method_params( & heat_cool_type = heat_cool_in inhomo_reion = inhomo_reion_in - end if - - ! Easy indexing for the passively advected quantities. - ! This lets us loop over all four groups (advected, species, aux) - ! in a single loop. - allocate(qpass_map(QVAR)) - allocate(upass_map(NVAR)) - npassive = 0 - do iadv = 1, nadv - upass_map(npassive + iadv) = UFA + iadv - 1 - qpass_map(npassive + iadv) = QFA + iadv - 1 - enddo - npassive = npassive + nadv - if(QFS > -1) then - do ispec = 1, nspec+naux - upass_map(npassive + ispec) = UFS + ispec - 1 - qpass_map(npassive + ispec) = QFS + ispec - 1 + ! Easy indexing for the passively advected quantities. + ! This lets us loop over all four groups (advected, species, aux) + ! in a single loop. + allocate(qpass_map(QVAR)) + allocate(upass_map(NVAR)) + npassive = 0 + do iadv = 1, nadv + upass_map(npassive + iadv) = UFA + iadv - 1 + qpass_map(npassive + iadv) = QFA + iadv - 1 enddo - npassive = npassive + nspec + naux - endif + npassive = npassive + nadv + if(QFS > -1) then + print *,'NSPEC NAUX ',nspec, naux + do ispec = 1, nspec+naux + upass_map(npassive + ispec) = UFS + ispec - 1 + qpass_map(npassive + ispec) = QFS + ispec - 1 + enddo + npassive = npassive + nspec + naux + endif + + end if end subroutine fort_set_method_params diff --git a/Source/comoving.cpp b/Source/comoving.cpp index 42f233a7..767d2c9c 100644 --- a/Source/comoving.cpp +++ b/Source/comoving.cpp @@ -8,6 +8,8 @@ Real Nyx::initial_z = -1.0; Real Nyx::final_a = -1.0; Real Nyx::final_z = -1.0; Real Nyx::relative_max_change_a = 0.01; +Real Nyx::absolute_max_change_a = -1.0; +Real Nyx::dt_binpow = -1.0; void Nyx::read_comoving_params () @@ -32,6 +34,9 @@ Nyx::read_comoving_params () } pp.query("relative_max_change_a", relative_max_change_a); + pp.query("absolute_max_change_a", absolute_max_change_a); + pp.query("dt_binpow", dt_binpow); + // for shrinking box tests, initial_z < 0 is ok if (initial_z < 0) @@ -45,6 +50,7 @@ void Nyx::comoving_est_time_step (Real& cur_time, Real& estdt) { Real change_allowed = relative_max_change_a; + Real fixed_da = absolute_max_change_a; Real dt = estdt; Real new_dummy_a; int dt_modified; @@ -56,7 +62,45 @@ Nyx::comoving_est_time_step (Real& cur_time, Real& estdt) // "old_a" and "new_a" -- we can't do that until after we compute dt and then // integrate a forward. fort_estdt_comoving_a - (&new_a, &new_dummy_a, &dt, &change_allowed, &final_a, &dt_modified); + (&new_a, &new_dummy_a, &dt, &change_allowed, &fixed_da, &final_a, &dt_modified); + + if(dt_binpow >= 0) + { + if(estdt>=dt) + estdt=dt; + else if(estdt>.5*dt) + { + estdt=.5*dt; + // std::cout << "Lavel = 1" <.25*dt) + { + estdt=.25*dt; + // std::cout << "Lavel = 2" <.125*dt) + { + estdt=.125*dt; + // std::cout << "Lavel = 3" <.0625*dt) + { + estdt=.0625*dt; + // std::cout << "Lavel = 4" < 4" <= 0) + { + if(estdt>=dt) + estdt=dt; + else if(estdt>.5*dt) + { + estdt=.5*dt; + // std::cout << "Lavel = 1" <.25*dt) + { + estdt=.25*dt; + // std::cout << "Lavel = 2" <.125*dt) + { + estdt=.125*dt; + // std::cout << "Lavel = 3" <.0625*dt) + { + estdt=.0625*dt; + // std::cout << "Lavel = 4" < 4" < 0) then + start_slope = H_0*dsqrt(comoving_OmM / start_a + OmL*start_a**2) + else + start_slope = comoving_h + end if + + ! Compute a provisional value of ln(a) at the new time + end_a = start_a + start_slope * Delta_t + + ! Compute the slope at the new time + if (comoving_type > 0) then + end_slope = H_0*dsqrt(comoving_OmM / end_a + OmL*end_a**2) + else + end_slope = comoving_h + end if + + ! Now recompute a at the new time using the average of the two slopes + end_a = start_a + 0.5d0 * (start_slope + end_slope) * Delta_t + + ! We have crossed from a too small to a too big in this step + if ( (end_a - a_value) * (start_a - a_value) < 0) then + dt = ( ( end_a - a_value) * dble(j ) + & + (a_value - start_a) * dble(j+1) ) / (end_a - start_a) * Delta_t + exit + end if + end do + + end subroutine fort_integrate_comoving_a_to_a +! ::: +! ::: ---------------------------------------------------------------- ! ::: subroutine fort_est_maxdt_comoving_a(old_a,dt) @@ -190,29 +267,84 @@ end subroutine fort_est_maxdt_comoving_a ! ::: ---------------------------------------------------------------- ! ::: - subroutine fort_estdt_comoving_a(old_a,new_a,dt,change_allowed,final_a,dt_modified) & + ! This might only work for t=0=> a=.00625, although constant canceled + subroutine fort_est_lindt_comoving_a(old_a,new_a,dt) + + use fundamental_constants_module, only: Hubble_const + use comoving_module , only: comoving_h, comoving_OmM, comoving_type + + implicit none + + real(rt), intent(in ) :: old_a,new_a + real(rt), intent(inout) :: dt + + real(rt) :: H_0, OmL + real(rt) :: lin_dt + + OmL = 1.d0 - comoving_OmM + + ! This subroutine computes dt based on not changing a by more than 5% + ! if we use forward Euler integration + ! d(ln(a)) / dt = H_0 * sqrt(OmM/a^3 + OmL) + + H_0 = comoving_h * Hubble_const + + ! Could definately be optimized better + if (H_0 .ne. 0.0d0) then + + lin_dt= ((new_a/(.75**(2/3)*(OmL+ comoving_OmM)**(1/3)))**(.75) - & + ((old_a/(.75**(2/3)*(OmL+ comoving_OmM)**(1/3)))**(.75) ) ) /H_0 + dt=lin_dt + + else + + ! dt is unchanged + + end if + + end subroutine fort_est_lindt_comoving_a + +! ::: +! ::: ---------------------------------------------------------------- +! ::: + + subroutine fort_estdt_comoving_a(old_a,new_a,dt,change_allowed,fixed_da,final_a,dt_modified) & bind(C, name="fort_estdt_comoving_a") use comoving_module , only: comoving_h implicit none - real(rt), intent(in ) :: old_a, change_allowed, final_a + real(rt), intent(in ) :: old_a, change_allowed, fixed_da, final_a real(rt), intent(inout) :: dt real(rt), intent( out) :: new_a integer , intent( out) :: dt_modified + real(rt) a_value + real(rt) max_dt + max_dt = dt if (comoving_h .ne. 0.0d0) then - ! First call this to make sure dt that we send to integration routine isnt outrageous - call fort_est_maxdt_comoving_a(old_a,dt) + if( fixed_da .le. 0.0d0) then + ! First call this to make sure dt that we send to integration routine isnt outrageous + call fort_est_maxdt_comoving_a(old_a,dt) + + ! Initial call to see if existing dt will work + call fort_integrate_comoving_a(old_a,new_a,dt) + + ! Make sure a isn't growing too fast + call enforce_percent_change(old_a,new_a,dt,change_allowed) + else + ! First call this to make sure dt that we send to integration routine isnt outrageous + new_a = (old_a + fixed_da); + call fort_est_lindt_comoving_a(old_a,new_a,dt) + call fort_est_maxdt_comoving_a(old_a,dt) - ! Initial call to see if existing dt will work - call fort_integrate_comoving_a(old_a,new_a,dt) + ! Then integrate old_a to a_value using dt as a guess for the maximum dt + ! Output dt is based on a fraction of the input dt + call fort_integrate_comoving_a_to_a(old_a,new_a,dt) + endif - ! Make sure a isn't growing too fast - call enforce_percent_change(old_a,new_a,dt,change_allowed) - ! Make sure we don't go past final_a (if final_a is set) if (final_a > 0.0d0) & call enforce_final_a(old_a,new_a,dt,final_a) diff --git a/Source/sum_utils.cpp b/Source/sum_utils.cpp index df2f0366..94bf5ffc 100644 --- a/Source/sum_utils.cpp +++ b/Source/sum_utils.cpp @@ -32,7 +32,7 @@ Nyx::vol_weight_sum (const std::string& name, } #ifdef _OPENMP -#pragma omp if (!system::regtest_reduction) parallel reduction(+:sum) +#pragma omp parallel if (!system::regtest_reduction) reduction(+:sum) #endif for (MFIter mfi(*mf,true); mfi.isValid(); ++mfi) { @@ -77,7 +77,7 @@ Nyx::vol_weight_sum (MultiFab& mf, bool masked) } #ifdef _OPENMP -#pragma omp if (!system::regtest_reduction) parallel reduction(+:sum) +#pragma omp parallel if (!system::regtest_reduction) reduction(+:sum) #endif for (MFIter mfi(mf,true); mfi.isValid(); ++mfi) { @@ -124,7 +124,7 @@ Nyx::vol_weight_squared_sum_level (const std::string& name, Real lev_vol = parent->boxArray(level).d_numPts() * dx[0] * dx[1] * dx[2]; #ifdef _OPENMP -#pragma omp if (!system::regtest_reduction) parallel reduction(+:sum) +#pragma omp parallel if (!system::regtest_reduction) reduction(+:sum) #endif for (MFIter mfi(*mf,true); mfi.isValid(); ++mfi) { @@ -168,7 +168,7 @@ Nyx::vol_weight_squared_sum (const std::string& name, } #ifdef _OPENMP -#pragma omp if (!system::regtest_reduction) parallel reduction(+:sum) +#pragma omp parallel if (!system::regtest_reduction) reduction(+:sum) #endif for (MFIter mfi(*mf,true); mfi.isValid(); ++mfi) {