From a178c299f98c43a0993f33095d0426fb5a73dfc5 Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Mon, 9 Apr 2018 21:29:56 -0700 Subject: [PATCH 01/12] We have to redefine grav_n_grow to allow for the fact that ghost particles were created on the coarser level which means they may be (iteration-1) cells further out than ghost_width which means we need to grow grav_n_grow to reflect that. --- Source/Nyx_advance.cpp | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/Source/Nyx_advance.cpp b/Source/Nyx_advance.cpp index 3a83c768..15595c16 100644 --- a/Source/Nyx_advance.cpp +++ b/Source/Nyx_advance.cpp @@ -108,14 +108,16 @@ Nyx::advance_hydro_plus_particles (Real time, // ghost_width + (1-iteration) - 1: // the minus 1 arises because this occurs *after* the move - int where_width = ghost_width + (1-iteration) - 1; + int where_width = ghost_width + (1-iteration) - 1; // *** grav_n_grow *** is used // *) to determine how many ghost cells we need to fill in the MultiFab from // which the particle interpolates its acceleration // *) to set how many cells the Where call in moveKickDrift tests = (grav.nGrow()-2). + // *) the (1-iteration) arises because the ghost particles are created on the coarser + // level which means in iteration 2 the ghost particles may have moved 1 additional cell along - int grav_n_grow = ghost_width + (1-iteration) + + int grav_n_grow = ghost_width + (1-iteration) + (iteration-1) + stencil_interpolation_width ; BL_PROFILE_REGION_START("R::Nyx::advance_hydro_plus_particles"); @@ -250,7 +252,7 @@ Nyx::advance_hydro_plus_particles (Real time, const Real a_half = 0.5 * (a_old + a_new); if (particle_verbose && ParallelDescriptor::IOProcessor()) - std::cout << "moveKickDrift ... updating particle positions and velocity\n"; + std::cout << "moveKickDrift ... updating particle positions and velocity\n"; for (int lev = level; lev <= finest_level_to_advance; lev++) { From 71c5efd3c89e0d89c506372d8b80976d40f2b5bd Mon Sep 17 00:00:00 2001 From: Weiqun Zhang Date: Wed, 11 Apr 2018 12:59:49 -0700 Subject: [PATCH 02/12] USE_CVODE --> AMREX_USE_CVODE --- Exec/HydroTests/TurbForce/Nyx_setup.cpp | 2 +- Source/Initialization/Nyx_setup.cpp | 2 +- Source/Nyx.cpp | 4 ++-- Source/nyx_main.cpp | 4 ++-- 4 files changed, 6 insertions(+), 6 deletions(-) diff --git a/Exec/HydroTests/TurbForce/Nyx_setup.cpp b/Exec/HydroTests/TurbForce/Nyx_setup.cpp index 288adece..32b2a771 100644 --- a/Exec/HydroTests/TurbForce/Nyx_setup.cpp +++ b/Exec/HydroTests/TurbForce/Nyx_setup.cpp @@ -731,7 +731,7 @@ Nyx::no_hydro_setup() } #endif -#ifdef USE_CVODE +#ifdef AMREX_USE_CVODE void Nyx::set_simd_width(const int simd_width) { diff --git a/Source/Initialization/Nyx_setup.cpp b/Source/Initialization/Nyx_setup.cpp index 898a08b7..14561cd4 100644 --- a/Source/Initialization/Nyx_setup.cpp +++ b/Source/Initialization/Nyx_setup.cpp @@ -739,7 +739,7 @@ Nyx::no_hydro_setup() } #endif -#ifdef USE_CVODE +#ifdef AMREX_USE_CVODE void Nyx::set_simd_width(const int simd_width) { diff --git a/Source/Nyx.cpp b/Source/Nyx.cpp index ce0d60ba..fa2c28d8 100644 --- a/Source/Nyx.cpp +++ b/Source/Nyx.cpp @@ -262,7 +262,7 @@ Nyx::read_params () pp.query("strict_subcycling",strict_subcycling); -#ifdef USE_CVODE +#ifdef AMREX_USE_CVODE pp.query("simd_width", simd_width); if (simd_width < 1) amrex::Abort("simd_width must be a positive integer"); set_simd_width(simd_width); @@ -426,7 +426,7 @@ Nyx::read_params () std::cout << std::endl; } -#ifndef USE_CVODE +#ifndef AMREX_USE_CVODE if (heat_cool_type == 5 || heat_cool_type == 7) amrex::Error("Nyx:: cannot set heat_cool_type = 5 or 7 unless USE_CVODE=TRUE"); #endif diff --git a/Source/nyx_main.cpp b/Source/nyx_main.cpp index c6352cfa..a8d44512 100644 --- a/Source/nyx_main.cpp +++ b/Source/nyx_main.cpp @@ -112,7 +112,7 @@ nyx_main (int argc, char* argv[]) const Real time_before_main_loop = ParallelDescriptor::second(); -#ifdef USE_CVODE +#ifdef AMREX_USE_CVODE Nyx::alloc_simd_vec(); #endif @@ -145,7 +145,7 @@ nyx_main (int argc, char* argv[]) } // ---- end while( ! finished) -#ifdef USE_CVODE +#ifdef AMREX_USE_CVODE Nyx::dealloc_simd_vec(); #endif From 3e69c2c6403c00bfa8512df413c532d40305168e Mon Sep 17 00:00:00 2001 From: Weiqun Zhang Date: Wed, 18 Apr 2018 10:31:39 -0700 Subject: [PATCH 03/12] mempool_module --> amrex_mempool_module --- Source/HydroFortran/Nyx_advection_3d.f90 | 6 +++--- Source/HydroFortran/ppm_3d.f90 | 2 +- Source/HydroFortran/trans_3d.f90 | 2 +- 3 files changed, 5 insertions(+), 5 deletions(-) diff --git a/Source/HydroFortran/Nyx_advection_3d.f90 b/Source/HydroFortran/Nyx_advection_3d.f90 index 757500da..1f7d6d58 100644 --- a/Source/HydroFortran/Nyx_advection_3d.f90 +++ b/Source/HydroFortran/Nyx_advection_3d.f90 @@ -19,7 +19,7 @@ subroutine fort_advance_gas(time,lo,hi,& bind(C, name="fort_advance_gas") use amrex_fort_module, only : rt => amrex_real - use mempool_module, only : bl_allocate, bl_deallocate + use amrex_mempool_module, only : bl_allocate, bl_deallocate use meth_params_module, only : QVAR, NVAR, NHYP, normalize_species use enforce_module, only : enforce_nonnegative_species use bl_constants_module @@ -205,7 +205,7 @@ subroutine umeth3d(q, c, csml, flatn, qd_l1, qd_l2, qd_l3, qd_h1, qd_h2, qd_h3, pdivu,a_old,a_new,print_fortran_warnings) use amrex_fort_module, only : rt => amrex_real - use mempool_module, only : bl_allocate, bl_deallocate + use amrex_mempool_module, only : bl_allocate, bl_deallocate use bl_constants_module use meth_params_module, only : QVAR, NVAR, QU, ppm_type, & use_colglaz, corner_coupling, & @@ -1362,7 +1362,7 @@ subroutine cmpflx(qm,qp,qpd_l1,qpd_l2,qpd_l3,qpd_h1,qpd_h2,qpd_h3, & idir,ilo,ihi,jlo,jhi,kc,kflux,k3d,print_fortran_warnings) use amrex_fort_module, only : rt => amrex_real - use mempool_module, only : bl_allocate, bl_deallocate + use amrex_mempool_module, only : bl_allocate, bl_deallocate use bl_constants_module use meth_params_module, only : QVAR, NVAR diff --git a/Source/HydroFortran/ppm_3d.f90 b/Source/HydroFortran/ppm_3d.f90 index d0667256..722ce8f6 100644 --- a/Source/HydroFortran/ppm_3d.f90 +++ b/Source/HydroFortran/ppm_3d.f90 @@ -74,7 +74,7 @@ subroutine ppm_type1(s,s_l1,s_l2,s_l3,s_h1,s_h2,s_h3, & flatn,f_l1,f_l2,f_l3,f_h1,f_h2,f_h3, & Ip,Im,ilo1,ilo2,ihi1,ihi2,dx,dy,dz,dt_over_a,k3d,kc) - use mempool_module, only: bl_allocate, bl_deallocate + use amrex_mempool_module, only: bl_allocate, bl_deallocate use amrex_fort_module, only : rt => amrex_real use meth_params_module, only : ppm_type, ppm_flatten_before_integrals use bl_constants_module diff --git a/Source/HydroFortran/trans_3d.f90 b/Source/HydroFortran/trans_3d.f90 index ee2007ff..0a4a2b81 100644 --- a/Source/HydroFortran/trans_3d.f90 +++ b/Source/HydroFortran/trans_3d.f90 @@ -770,7 +770,7 @@ subroutine transz(qxm,qxmo,qxp,qxpo, & ugdnvz,pgdnvz,pgdz_l1,pgdz_l2,pgdz_l3,pgdz_h1,pgdz_h2,pgdz_h3, & cdtdz,ilo,ihi,jlo,jhi,km,kc,k3d) - use mempool_module, only: bl_allocate, bl_deallocate + use amrex_mempool_module, only: bl_allocate, bl_deallocate use meth_params_module, only : QVAR, NVAR, QRHO, QU, QV, QW, & QPRES, QREINT, & URHO, UMX, UMY, UMZ, UEDEN, & From 61e579b10f29af971f3a355e4f673261310a1f5b Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Fri, 20 Apr 2018 14:24:41 -0700 Subject: [PATCH 04/12] 1) Remove calls to fill_ec_grow 2) Clean up some stuff so we can remove F_BaseLib dependency in the long run --- Source/Constants/constants_cgs.f90 | 1 - Source/EOS/eos_stuff.f90 | 2 -- Source/Forcing/forcing_spect.f90 | 8 +++----- Source/Gravity/Gravity.cpp | 24 ++++++++++++++++++------ Source/Network/network.f90 | 1 - Source/Network/nyx_burner.f90 | 3 --- 6 files changed, 21 insertions(+), 18 deletions(-) diff --git a/Source/Constants/constants_cgs.f90 b/Source/Constants/constants_cgs.f90 index 3eefdb5c..0b7a5708 100644 --- a/Source/Constants/constants_cgs.f90 +++ b/Source/Constants/constants_cgs.f90 @@ -3,7 +3,6 @@ module fundamental_constants_module use amrex_fort_module, only : rt => amrex_real - use bl_types implicit none diff --git a/Source/EOS/eos_stuff.f90 b/Source/EOS/eos_stuff.f90 index 5eab4309..166c0a50 100644 --- a/Source/EOS/eos_stuff.f90 +++ b/Source/EOS/eos_stuff.f90 @@ -26,8 +26,6 @@ module eos_module - use bl_types - use bl_space use bl_constants_module, only: M_PI, ONE use network, only: nspec, aion, zion use atomic_rates_module, only: XHYDROGEN diff --git a/Source/Forcing/forcing_spect.f90 b/Source/Forcing/forcing_spect.f90 index f03d9fa7..2245e8d9 100644 --- a/Source/Forcing/forcing_spect.f90 +++ b/Source/Forcing/forcing_spect.f90 @@ -3,8 +3,6 @@ module forcing_spect_module use amrex_fort_module, only : bl_spacedim, rt => amrex_real - use bl_types - use bl_constants_module, only : ZERO implicit none @@ -46,9 +44,9 @@ subroutine fort_alloc_spect(length) & call bl_abort('number of forcing modes must be positive') end if - modes_even(:,:) = ZERO - modes_even(:,:) = ZERO - wavevectors(:,:) = ZERO + modes_even(:,:) = 0.d0 + modes_even(:,:) = 0.d0 + wavevectors(:,:) = 0.d0 end subroutine fort_alloc_spect diff --git a/Source/Gravity/Gravity.cpp b/Source/Gravity/Gravity.cpp index b76cb1c4..7425c750 100644 --- a/Source/Gravity/Gravity.cpp +++ b/Source/Gravity/Gravity.cpp @@ -958,13 +958,11 @@ Gravity::get_old_grav_vector (int level, // Fill grow cells in grad_phi, will need to compute grad_phi_cc in 1 grow cell const Geometry& geom = parent->Geom(level); +#if 0 if (level == 0) { for (int i = 0; i < BL_SPACEDIM ; i++) - { - grad_phi_prev[level][i]->setBndry(0); - grad_phi_prev[level][i]->FillBoundary(geom.periodicity()); - } + grad_phi_prev[level][i]->setBndry(0.); } else { @@ -973,11 +971,16 @@ Gravity::get_old_grav_vector (int level, fill_ec_grow(level, amrex::GetVecOfPtrs(grad_phi_prev[level]), amrex::GetVecOfPtrs(crse_grad_phi)); } +#endif + + // Fill boundary values at the current level + for (int i = 0; i < BL_SPACEDIM ; i++) + grad_phi_prev[level][i]->FillBoundary(geom.periodicity()); // Average edge-centered gradients to cell centers. amrex::average_face_to_cellcenter(grav_vector, - amrex::GetVecOfConstPtrs(grad_phi_prev[level]), - geom); + amrex::GetVecOfConstPtrs(grad_phi_prev[level]), + geom); #ifdef CGRAV if (gravity_type == "CompositeGrav") @@ -999,6 +1002,8 @@ Gravity::get_old_grav_vector (int level, // This is a hack-y way to fill the ghost cell values of grav_vector // before returning it + // Note that this fills ghost cells over the coarse grid from interpolation, + // not from the ghost cell values previously filled after the fill_ec_grow stuff. AmrLevel* amrlev = &parent->getLevel(level); int ng = grav_vector.nGrow(); AmrLevel::FillPatch(*amrlev,grav_vector,ng,time,Gravity_Type,0,BL_SPACEDIM); @@ -1022,6 +1027,7 @@ Gravity::get_new_grav_vector (int level, // Fill grow cells in `grad_phi`, will need to compute `grad_phi_cc` in // 1 grow cell const Geometry& geom = parent->Geom(level); +#if 0 if (level == 0) { for (int i = 0; i < BL_SPACEDIM ; i++) @@ -1037,6 +1043,10 @@ Gravity::get_new_grav_vector (int level, fill_ec_grow(level, amrex::GetVecOfPtrs(grad_phi_curr[level]), amrex::GetVecOfPtrs(crse_grad_phi)); } +#endif + + for (int i = 0; i < BL_SPACEDIM ; i++) + grad_phi_curr[level][i]->FillBoundary(geom.periodicity()); // Average edge-centered gradients to cell centers, excluding grow cells amrex::average_face_to_cellcenter(grav_vector, @@ -1068,6 +1078,8 @@ Gravity::get_new_grav_vector (int level, // This is a hack-y way to fill the ghost cell values of grav_vector // before returning it + // Note that this fills ghost cells over the coarse grid from interpolation, + // not from the ghost cell values previously filled after the fill_ec_grow stuff. AmrLevel* amrlev = &parent->getLevel(level) ; int ng = grav_vector.nGrow(); AmrLevel::FillPatch(*amrlev,grav_vector,ng,time,Gravity_Type,0,BL_SPACEDIM); diff --git a/Source/Network/network.f90 b/Source/Network/network.f90 index 4be6655c..b2a4bc41 100644 --- a/Source/Network/network.f90 +++ b/Source/Network/network.f90 @@ -22,7 +22,6 @@ module network use amrex_fort_module, only : rt => amrex_real - use bl_types implicit none diff --git a/Source/Network/nyx_burner.f90 b/Source/Network/nyx_burner.f90 index cb7dc07c..4d9ba796 100644 --- a/Source/Network/nyx_burner.f90 +++ b/Source/Network/nyx_burner.f90 @@ -1,9 +1,6 @@ module nyx_burner_module use amrex_fort_module, only : rt => amrex_real - use bl_types - use bl_constants_module - use bl_error_module use eos_module use network From 7d42ce6470f37e66a4e7b531b1e3d04035e19322 Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Fri, 20 Apr 2018 14:56:05 -0700 Subject: [PATCH 05/12] Use the Fortran version of multilevel AssignDensity. --- Source/NyxParticleContainer.H | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/Source/NyxParticleContainer.H b/Source/NyxParticleContainer.H index c5977e4e..d778c155 100644 --- a/Source/NyxParticleContainer.H +++ b/Source/NyxParticleContainer.H @@ -63,11 +63,7 @@ public: virtual void AssignDensity (amrex::Vector >& mf, int lev_min = 0, int ncomp = 1, int finest_level = -1, int ngrow = 1) const override { - if (finestLevel() == 0) { - amrex::NeighborParticleContainer::AssignDensityFort(0, mf, lev_min, ncomp, finest_level); - } else { - amrex::NeighborParticleContainer::AssignDensity(0, sub_cycle, mf, lev_min, ncomp, finest_level, ngrow); - } + amrex::NeighborParticleContainer::AssignDensityFort(0, mf, lev_min, ncomp, finest_level); } void MultiplyParticleMass (int lev, amrex::Real mult); From d0db7c65c81e5a552961a4d63aaabd2ecb2d9d9b Mon Sep 17 00:00:00 2001 From: Weiqun Zhang Date: Wed, 25 Apr 2018 13:36:56 -0700 Subject: [PATCH 06/12] use a different version of complementIn --- Source/NeutrinoParticleContainer.cpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/Source/NeutrinoParticleContainer.cpp b/Source/NeutrinoParticleContainer.cpp index 7775197f..61c1ea7e 100644 --- a/Source/NeutrinoParticleContainer.cpp +++ b/Source/NeutrinoParticleContainer.cpp @@ -253,7 +253,9 @@ NeutrinoParticleContainer::AssignRelativisticDensity (Vector Date: Wed, 2 May 2018 10:26:06 -0700 Subject: [PATCH 07/12] These changes should not change the answer at all, but they are structural changes in prep for merging in the SDC stuff. --- Source/HeatCool/integrate_state_3d.f90 | 8 +- Source/HydroFortran/Make.package | 1 + Source/HydroFortran/Nyx_advection_3d.f90 | 290 +++++++++++++++--- Source/HydroFortran/add_grav_source_3d.f90 | 6 +- Source/HydroFortran/make_hydro_sources_3d.f90 | 143 +++++++++ Source/Make.package | 11 +- Source/Nyx.H | 42 ++- Source/Nyx_F.H | 42 ++- Source/Nyx_advance.cpp | 24 +- Source/compute_hydro_sources.cpp | 116 +++++++ Source/sdc_hydro.cpp | 132 ++++++++ Source/sdc_reactions.cpp | 96 ++++++ Source/{Nyx_hydro.cpp => strang_hydro.cpp} | 18 +- ...ang_splitting.cpp => strang_reactions.cpp} | 4 +- Source/update_state_with_sources.cpp | 41 +++ 15 files changed, 902 insertions(+), 72 deletions(-) create mode 100644 Source/HydroFortran/make_hydro_sources_3d.f90 create mode 100644 Source/compute_hydro_sources.cpp create mode 100644 Source/sdc_hydro.cpp create mode 100644 Source/sdc_reactions.cpp rename Source/{Nyx_hydro.cpp => strang_hydro.cpp} (94%) rename Source/{strang_splitting.cpp => strang_reactions.cpp} (94%) create mode 100644 Source/update_state_with_sources.cpp diff --git a/Source/HeatCool/integrate_state_3d.f90 b/Source/HeatCool/integrate_state_3d.f90 index d463080c..f17f5078 100644 --- a/Source/HeatCool/integrate_state_3d.f90 +++ b/Source/HeatCool/integrate_state_3d.f90 @@ -1,7 +1,7 @@ subroutine integrate_state(lo, hi, & state , s_l1, s_l2, s_l3, s_h1, s_h2, s_h3, & diag_eos, d_l1, d_l2, d_l3, d_h1, d_h2, d_h3, & - dx, time, a, half_dt, min_iter, max_iter) & + a, half_dt, min_iter, max_iter) & bind(C, name="integrate_state") ! @@ -17,10 +17,6 @@ subroutine integrate_state(lo, hi, & ! The state vars ! diag_eos* : double arrays ! Temp and Ne -! src_* : double arrays -! The source terms to be added to state (iterative approx.) -! double array (3) -! The low corner of the entire domain ! a : double ! The current a ! half_dt : double @@ -43,7 +39,7 @@ subroutine integrate_state(lo, hi, & integer , intent(in ) :: d_l1, d_l2, d_l3, d_h1, d_h2, d_h3 real(rt), intent(inout) :: state(s_l1:s_h1, s_l2:s_h2,s_l3:s_h3, NVAR) real(rt), intent(inout) :: diag_eos(d_l1:d_h1, d_l2:d_h2,d_l3:d_h3, NDIAG) - real(rt), intent(in ) :: dx(3), time, a, half_dt + real(rt), intent(in ) :: a, half_dt integer , intent(inout) :: min_iter, max_iter if (heat_cool_type .eq. 1) then diff --git a/Source/HydroFortran/Make.package b/Source/HydroFortran/Make.package index 16f55a30..2932557f 100644 --- a/Source/HydroFortran/Make.package +++ b/Source/HydroFortran/Make.package @@ -2,6 +2,7 @@ f90EXE_sources += add_grav_source_3d.f90 f90EXE_sources += analriem.f90 f90EXE_sources += enforce_minimum_density_3d.f90 f90EXE_sources += flatten_3d.f90 +f90EXE_sources += make_hydro_sources_3d.f90 f90EXE_sources += normalize_species_3d.f90 f90EXE_sources += Nyx_advection_3d.f90 f90EXE_sources += ppm_3d.f90 diff --git a/Source/HydroFortran/Nyx_advection_3d.f90 b/Source/HydroFortran/Nyx_advection_3d.f90 index 1f7d6d58..bd320066 100644 --- a/Source/HydroFortran/Nyx_advection_3d.f90 +++ b/Source/HydroFortran/Nyx_advection_3d.f90 @@ -56,8 +56,8 @@ subroutine fort_advance_gas(time,lo,hi,& real(rt), pointer :: flatn(:,:,:) real(rt), pointer :: c(:,:,:) real(rt), pointer :: csml(:,:,:) - real(rt), pointer :: div(:,:,:) - real(rt), pointer :: pdivu(:,:,:) + real(rt), pointer :: divu_nd(:,:,:) + real(rt), pointer :: divu_cc(:,:,:) real(rt), pointer :: srcQ(:,:,:,:) real(rt) dx,dy,dz @@ -89,8 +89,8 @@ subroutine fort_advance_gas(time,lo,hi,& call bl_allocate( srcQ, lo-1, hi+1, QVAR) - call bl_allocate( div, lo, hi+1) - call bl_allocate( pdivu, lo, hi) + call bl_allocate(divu_nd, lo, hi+1) + call bl_allocate(divu_cc, lo, hi) dx = delta(1) dy = delta(2) @@ -119,11 +119,11 @@ subroutine fort_advance_gas(time,lo,hi,& ugdnvx_out,ugdnvx_l1,ugdnvx_l2,ugdnvx_l3,ugdnvx_h1,ugdnvx_h2,ugdnvx_h3, & ugdnvy_out,ugdnvy_l1,ugdnvy_l2,ugdnvy_l3,ugdnvy_h1,ugdnvy_h2,ugdnvy_h3, & ugdnvz_out,ugdnvz_l1,ugdnvz_l2,ugdnvz_l3,ugdnvz_h1,ugdnvz_h2,ugdnvz_h3, & - pdivu,a_old,a_new,print_fortran_warnings) + divu_cc,a_old,a_new,print_fortran_warnings) ! Compute divergence of velocity field (on surroundingNodes(lo,hi)) - call divu(lo,hi,q,q_l1,q_l2,q_l3,q_h1,q_h2,q_h3, & - dx,dy,dz,div,lo(1),lo(2),lo(3),hi(1)+1,hi(2)+1,hi(3)+1) + call make_divu_nd(lo,hi,q,q_l1,q_l2,q_l3,q_h1,q_h2,q_h3, & + dx,dy,dz,divu_nd,lo(1),lo(2),lo(3),hi(1)+1,hi(2)+1,hi(3)+1) ! Conservative update call consup(uin,uin_l1,uin_l2,uin_l3,uin_h1,uin_h2,uin_h3, & @@ -132,16 +132,16 @@ subroutine fort_advance_gas(time,lo,hi,& flux1,flux1_l1,flux1_l2,flux1_l3,flux1_h1,flux1_h2,flux1_h3, & flux2,flux2_l1,flux2_l2,flux2_l3,flux2_h1,flux2_h2,flux2_h3, & flux3,flux3_l1,flux3_l2,flux3_l3,flux3_h1,flux3_h2,flux3_h3, & - div,pdivu,lo,hi,dx,dy,dz,dt,a_old,a_new) + divu_nd,divu_cc,lo,hi,dx,dy,dz,dt,a_old,a_new) ! We are done with these here so can go ahead and free up the space. call bl_deallocate(q) call bl_deallocate(flatn) call bl_deallocate(c) call bl_deallocate(csml) - call bl_deallocate(div) + call bl_deallocate(divu_nd) call bl_deallocate(srcQ) - call bl_deallocate(pdivu) + call bl_deallocate(divu_cc) ! Enforce the density >= small_dens. Make sure we do this immediately after consup. call enforce_minimum_density(uin, uin_l1, uin_l2, uin_l3, uin_h1, uin_h2, uin_h3, & @@ -202,7 +202,7 @@ subroutine umeth3d(q, c, csml, flatn, qd_l1, qd_l2, qd_l3, qd_h1, qd_h2, qd_h3, ugdnvy_h1,ugdnvy_h2,ugdnvy_h3, & ugdnvz_out,ugdnvz_l1,ugdnvz_l2,ugdnvz_l3, & ugdnvz_h1,ugdnvz_h2,ugdnvz_h3, & - pdivu,a_old,a_new,print_fortran_warnings) + divu_cc,a_old,a_new,print_fortran_warnings) use amrex_fort_module, only : rt => amrex_real use amrex_mempool_module, only : bl_allocate, bl_deallocate @@ -244,7 +244,7 @@ subroutine umeth3d(q, c, csml, flatn, qd_l1, qd_l2, qd_l3, qd_h1, qd_h2, qd_h3, real(rt) ugdnvx_out(ugdnvx_l1:ugdnvx_h1,ugdnvx_l2:ugdnvx_h2,ugdnvx_l3:ugdnvx_h3) real(rt) ugdnvy_out(ugdnvy_l1:ugdnvy_h1,ugdnvy_l2:ugdnvy_h2,ugdnvy_l3:ugdnvy_h3) real(rt) ugdnvz_out(ugdnvz_l1:ugdnvz_h1,ugdnvz_l2:ugdnvz_h2,ugdnvz_l3:ugdnvz_h3) - real(rt) pdivu(ilo1:ihi1,ilo2:ihi2,ilo3:ihi3) + real(rt) divu_cc(ilo1:ihi1,ilo2:ihi2,ilo3:ihi3) real(rt) dx, dy, dz, dt real(rt) dtdx, dtdy, dtdz, hdt real(rt) cdtdx, cdtdy, cdtdz @@ -415,8 +415,8 @@ subroutine umeth3d(q, c, csml, flatn, qd_l1, qd_l2, qd_l3, qd_h1, qd_h2, qd_h3, cdtdy = THIRD*dtdy/a_half cdtdz = THIRD*dtdz/a_half - ! Initialize pdivu to zero - pdivu(:,:,:) = ZERO + ! Initialize divu_cc to zero + divu_cc(:,:,:) = ZERO ! Initialize kc (current k-level) and km (previous k-level) kc = 1 @@ -646,11 +646,7 @@ subroutine umeth3d(q, c, csml, flatn, qd_l1, qd_l2, qd_l3, qd_h1, qd_h2, qd_h3, if (k3d .ge. ilo3+1 .and. k3d .le. ihi3+1) then do j = ilo2,ihi2 do i = ilo1,ihi1 -! pdivu(i,j,k3d-1) = pdivu(i,j,k3d-1) + & -! HALF*(pgdnvzf(i,j,kc)+pgdnvzf(i,j,km)) * & -! (ugdnvzf(i,j,kc)-ugdnvzf(i,j,km))/dz - pdivu(i,j,k3d-1) = pdivu(i,j,k3d-1) + & - (ugdnvzf(i,j,kc)-ugdnvzf(i,j,km))/dz + divu_cc(i,j,k3d-1) = divu_cc(i,j,k3d-1) + (ugdnvzf(i,j,kc)-ugdnvzf(i,j,km))/dz end do end do end if @@ -772,12 +768,7 @@ subroutine umeth3d(q, c, csml, flatn, qd_l1, qd_l2, qd_l3, qd_h1, qd_h2, qd_h3, do j = ilo2,ihi2 do i = ilo1,ihi1 -! pdivu(i,j,k3d-1) = pdivu(i,j,k3d-1) + & -! HALF*(pgdnvxf(i+1,j,km) + pgdnvxf(i,j,km)) * & -! (ugdnvxf(i+1,j,km)-ugdnvxf(i,j,km))/dx + & -! HALF*(pgdnvyf(i,j+1,km) + pgdnvyf(i,j,km)) * & -! (ugdnvyf(i,j+1,km)-ugdnvyf(i,j,km))/dy - pdivu(i,j,k3d-1) = pdivu(i,j,k3d-1) + & + divu_cc(i,j,k3d-1) = divu_cc(i,j,k3d-1) + & (ugdnvxf(i+1,j,km)-ugdnvxf(i,j,km))/dx + & (ugdnvyf(i,j+1,km)-ugdnvyf(i,j,km))/dy end do @@ -1165,7 +1156,7 @@ subroutine consup(uin,uin_l1,uin_l2,uin_l3,uin_h1,uin_h2,uin_h3, & flux1,flux1_l1,flux1_l2,flux1_l3,flux1_h1,flux1_h2,flux1_h3, & flux2,flux2_l1,flux2_l2,flux2_l3,flux2_h1,flux2_h2,flux2_h3, & flux3,flux3_l1,flux3_l2,flux3_l3,flux3_h1,flux3_h2,flux3_h3, & - div,pdivu,lo,hi,dx,dy,dz,dt,a_old,a_new) + divu_nd,divu_cc,lo,hi,dx,dy,dz,dt,a_old,a_new) use amrex_fort_module, only : rt => amrex_real use bl_constants_module @@ -1188,8 +1179,8 @@ subroutine consup(uin,uin_l1,uin_l2,uin_l3,uin_h1,uin_h2,uin_h3, & real(rt) flux1(flux1_l1:flux1_h1,flux1_l2:flux1_h2,flux1_l3:flux1_h3,NVAR) real(rt) flux2(flux2_l1:flux2_h1,flux2_l2:flux2_h2,flux2_l3:flux2_h3,NVAR) real(rt) flux3(flux3_l1:flux3_h1,flux3_l2:flux3_h2,flux3_l3:flux3_h3,NVAR) - real(rt) div(lo(1):hi(1)+1,lo(2):hi(2)+1,lo(3):hi(3)+1) - real(rt) pdivu(lo(1):hi(1),lo(2):hi(2),lo(3):hi(3)) + real(rt) divu_nd(lo(1):hi(1)+1,lo(2):hi(2)+1,lo(3):hi(3)+1) + real(rt) divu_cc(lo(1):hi(1),lo(2):hi(2),lo(3):hi(3)) real(rt) dx, dy, dz, dt real(rt) a_old, a_new @@ -1214,7 +1205,7 @@ subroutine consup(uin,uin_l1,uin_l2,uin_l3,uin_h1,uin_h2,uin_h3, & do k = lo(3),hi(3) do j = lo(2),hi(2) do i = lo(1),hi(1)+1 - div1 = FOURTH*(div(i,j,k) + div(i,j+1,k) + div(i,j,k+1) + div(i,j+1,k+1)) + div1 = FOURTH*(divu_nd(i,j,k) + divu_nd(i,j+1,k) + divu_nd(i,j,k+1) + divu_nd(i,j+1,k+1)) div1 = difmag*min(ZERO,div1) flux1(i,j,k,n) = flux1(i,j,k,n) + dx*div1*(uin(i,j,k,n)-uin(i-1,j,k,n)) flux1(i,j,k,n) = flux1(i,j,k,n) * area1 * dt @@ -1224,7 +1215,7 @@ subroutine consup(uin,uin_l1,uin_l2,uin_l3,uin_h1,uin_h2,uin_h3, & do k = lo(3),hi(3) do j = lo(2),hi(2)+1 do i = lo(1),hi(1) - div1 = FOURTH*(div(i,j,k) + div(i+1,j,k) + div(i,j,k+1) + div(i+1,j,k+1)) + div1 = FOURTH*(divu_nd(i,j,k) + divu_nd(i+1,j,k) + divu_nd(i,j,k+1) + divu_nd(i+1,j,k+1)) div1 = difmag*min(ZERO,div1) flux2(i,j,k,n) = flux2(i,j,k,n) + dy*div1*(uin(i,j,k,n)-uin(i,j-1,k,n)) flux2(i,j,k,n) = flux2(i,j,k,n) * area2 * dt @@ -1234,7 +1225,7 @@ subroutine consup(uin,uin_l1,uin_l2,uin_l3,uin_h1,uin_h2,uin_h3, & do k = lo(3),hi(3)+1 do j = lo(2),hi(2) do i = lo(1),hi(1) - div1 = FOURTH*(div(i,j,k) + div(i+1,j,k) + div(i,j+1,k) + div(i+1,j+1,k)) + div1 = FOURTH*(divu_nd(i,j,k) + divu_nd(i+1,j,k) + divu_nd(i,j+1,k) + divu_nd(i+1,j+1,k)) div1 = difmag*min(ZERO,div1) flux3(i,j,k,n) = flux3(i,j,k,n) + dz*div1*(uin(i,j,k,n)-uin(i,j,k-1,n)) flux3(i,j,k,n) = flux3(i,j,k,n) * area3 * dt @@ -1301,16 +1292,15 @@ subroutine consup(uin,uin_l1,uin_l2,uin_l3,uin_h1,uin_h2,uin_h3, & + a_half * dt * src(i,j,k,n) ! ********************************************************************************* - ! This is the version where "pdivu" is actually just divu uout(i,j,k,n) = uout(i,j,k,n) & - - a_half * dt * (HALF * gamma_minus_1 * uin(i,j,k,n)) * pdivu(i,j,k) + - a_half * dt * (HALF * gamma_minus_1 * uin(i,j,k,n)) * divu_cc(i,j,k) uout(i,j,k,n) = uout(i,j,k,n) / & - ( ONE + a_half * dt * (HALF * gamma_minus_1 * pdivu(i,j,k)) * a_newsq_inv ) + ( ONE + a_half * dt * (HALF * gamma_minus_1 * divu_cc(i,j,k)) * a_newsq_inv ) ! ********************************************************************************* ! This is the original version - ! uout(i,j,k,n) = uout(i,j,k,n) - a_half * dt * pdivu(i,j,k) + ! uout(i,j,k,n) = uout(i,j,k,n) - a_half * dt * divu_cc(i,j,k) ! ********************************************************************************* uout(i,j,k,n) = uout(i,j,k,n) * a_newsq_inv @@ -1351,6 +1341,119 @@ subroutine consup(uin,uin_l1,uin_l2,uin_l3,uin_h1,uin_h2,uin_h3, & end subroutine consup +! ::: +! ::: ------------------------------------------------------------------ +! ::: + + subroutine update_state(lo,hi, & + uin,uin_l1,uin_l2,uin_l3,uin_h1,uin_h2,uin_h3, & + uout,uout_l1,uout_l2,uout_l3,uout_h1,uout_h2,uout_h3, & + src ,src_l1,src_l2,src_l3,src_h1,src_h2,src_h3, & + hydro_src ,hsrc_l1,hsrc_l2,hsrc_l3,hsrc_h1,hsrc_h2,hsrc_h3, & + divu_cc,d_l1,d_l2,d_l3,d_h1,d_h2,d_h3, & + dt,a_old,a_new,print_fortran_warnings) & + bind(C, name="fort_update_state") + + use amrex_fort_module, only : rt => amrex_real + use bl_constants_module + use meth_params_module, only : NVAR, URHO, UMX, UMZ, UEDEN, UEINT, UFS, & + gamma_minus_1, normalize_species + use enforce_module, only : enforce_nonnegative_species + + implicit none + + integer, intent(in) :: lo(3), hi(3) + integer, intent(in) ::print_fortran_warnings + integer, intent(in) :: uin_l1, uin_l2, uin_l3, uin_h1, uin_h2, uin_h3 + integer, intent(in) :: uout_l1, uout_l2, uout_l3, uout_h1, uout_h2, uout_h3 + integer, intent(in) :: src_l1, src_l2, src_l3, src_h1, src_h2, src_h3 + integer, intent(in) :: hsrc_l1, hsrc_l2, hsrc_l3, hsrc_h1, hsrc_h2, hsrc_h3 + integer, intent(in) :: d_l1,d_l2,d_l3,d_h1,d_h2,d_h3 + + real(rt), intent(in) :: uin( uin_l1: uin_h1, uin_l2: uin_h2, uin_l3: uin_h3,NVAR) + real(rt), intent(out) :: uout(uout_l1:uout_h1,uout_l2:uout_h2,uout_l3:uout_h3,NVAR) + real(rt), intent(in) :: src( src_l1: src_h1, src_l2: src_h2, src_l3: src_h3,NVAR) + real(rt), intent(in) :: hydro_src(hsrc_l1:hsrc_h1,hsrc_l2:hsrc_h2,hsrc_l3:hsrc_h3,NVAR) + real(rt), intent(in) :: divu_cc( d_l1: d_h1, d_l2: d_h2, d_l3: d_h3) + real(rt), intent(in) :: dt, a_old, a_new + + real(rt) :: a_half, a_oldsq, a_newsq + real(rt) :: a_new_inv, a_newsq_inv, a_half_inv, dt_a_new + integer :: i, j, k, n + + a_half = HALF * (a_old + a_new) + a_oldsq = a_old * a_old + a_newsq = a_new * a_new + + a_half_inv = ONE / a_half + a_new_inv = ONE / a_new + a_newsq_inv = ONE / a_newsq + + do n = 1, NVAR + + ! Actually do the update here + do k = lo(3),hi(3) + do j = lo(2),hi(2) + do i = lo(1),hi(1) + + ! Density + if (n .eq. URHO) then + uout(i,j,k,n) = uin(i,j,k,n) + dt * hydro_src(i,j,k,n) & + + dt * src(i,j,k,n) * a_half_inv + + ! Momentum + else if (n .ge. UMX .and. n .le. UMZ) then + uout(i,j,k,n) = a_old * uin(i,j,k,n) + dt * hydro_src(i,j,k,n) & + + dt * src(i,j,k,n) + uout(i,j,k,n) = uout(i,j,k,n) * a_new_inv + + ! (rho E) + else if (n .eq. UEDEN) then + uout(i,j,k,n) = a_oldsq * uin(i,j,k,n) + dt * hydro_src(i,j,k,n) & + + a_half * dt * src(i,j,k,n) + uout(i,j,k,n) = uout(i,j,k,n) * a_newsq_inv + + ! (rho e) + else if (n .eq. UEINT) then + + uout(i,j,k,n) = a_oldsq*uin(i,j,k,n) + dt * hydro_src(i,j,k,n) & + + a_half * dt * src(i,j,k,n) + + uout(i,j,k,n) = uout(i,j,k,n) * a_newsq_inv + +! We don't do this here because we are adding all of this term explicitly in hydro_src +! uout(i,j,k,n) = uout(i,j,k,n) * a_newsq_inv / & +! ( ONE + a_half * dt * (HALF * gamma_minus_1 * divu_cc(i,j,k)) * a_newsq_inv ) + + ! (rho X_i) and (rho adv_i) and (rho aux_i) + else + uout(i,j,k,n) = uin(i,j,k,n) + dt * hydro_src(i,j,k,n) & + + dt * src(i,j,k,n) * a_half_inv + + endif + + enddo + enddo + enddo + enddo + + ! Enforce the density >= small_dens. Make sure we do this immediately after consup. + call enforce_minimum_density(uin, uin_l1, uin_l2, uin_l3, uin_h1, uin_h2, uin_h3, & + uout,uout_l1,uout_l2,uout_l3,uout_h1,uout_h2,uout_h3, & + lo,hi,print_fortran_warnings) + + ! Enforce species >= 0 + call enforce_nonnegative_species(uout,uout_l1,uout_l2,uout_l3, & + uout_h1,uout_h2,uout_h3,lo,hi,0) + + ! Re-normalize the species + if (normalize_species .eq. 1) then + call normalize_new_species(uout,uout_l1,uout_l2,uout_l3,uout_h1,uout_h2,uout_h3, & + lo,hi) + end if + + end subroutine update_state + ! ::: ! ::: ------------------------------------------------------------------ ! ::: @@ -1827,8 +1930,8 @@ end subroutine riemannus ! ::: ------------------------------------------------------------------ ! ::: - subroutine divu(lo,hi,q,q_l1,q_l2,q_l3,q_h1,q_h2,q_h3,dx,dy,dz, & - div,div_l1,div_l2,div_l3,div_h1,div_h2,div_h3) + subroutine make_divu_nd(lo,hi,q,q_l1,q_l2,q_l3,q_h1,q_h2,q_h3,dx,dy,dz, & + divu_nd,div_l1,div_l2,div_l3,div_h1,div_h2,div_h3) use amrex_fort_module, only : rt => amrex_real use bl_constants_module @@ -1840,7 +1943,7 @@ subroutine divu(lo,hi,q,q_l1,q_l2,q_l3,q_h1,q_h2,q_h3,dx,dy,dz, & integer :: q_l1,q_l2,q_l3,q_h1,q_h2,q_h3 integer :: div_l1,div_l2,div_l3,div_h1,div_h2,div_h3 real(rt) :: dx, dy, dz - real(rt) :: div(div_l1:div_h1,div_l2:div_h2,div_l3:div_h3) + real(rt) :: divu_nd(div_l1:div_h1,div_l2:div_h2,div_l3:div_h3) real(rt) :: q(q_l1:q_h1,q_l2:q_h2,q_l3:q_h3,*) integer :: i, j, k @@ -1869,11 +1972,118 @@ subroutine divu(lo,hi,q,q_l1,q_l2,q_l3,q_h1,q_h2,q_h3,dx,dy,dz, & + q(i-1,j ,k ,QW) - q(i-1,j ,k-1,QW) & + q(i-1,j-1,k ,QW) - q(i-1,j-1,k-1,QW) ) * dzinv - div(i,j,k) = FOURTH*(ux + vy + wz) + divu_nd(i,j,k) = FOURTH*(ux + vy + wz) + + enddo + enddo + enddo + + end subroutine make_divu_nd + +! ::: +! ::: ------------------------------------------------------------------ +! ::: + + subroutine add_grav_source(uin,uin_l1,uin_l2,uin_l3,uin_h1,uin_h2,uin_h3, & + uout,uout_l1,uout_l2,uout_l3,uout_h1,uout_h2,uout_h3, & + grav, gv_l1, gv_l2, gv_l3, gv_h1, gv_h2, gv_h3, & + lo,hi,dx,dy,dz,dt,a_old,a_new,e_added,ke_added) + + use amrex_fort_module, only : rt => amrex_real + use eos_module + use meth_params_module, only : NVAR, URHO, UMX, UMY, UMZ, & + UEDEN, grav_source_type + + implicit none + + integer lo(3), hi(3) + integer uin_l1,uin_l2,uin_l3,uin_h1,uin_h2,uin_h3 + integer uout_l1, uout_l2, uout_l3, uout_h1, uout_h2, uout_h3 + integer gv_l1, gv_l2, gv_l3, gv_h1, gv_h2, gv_h3 + + real(rt) uin( uin_l1: uin_h1, uin_l2: uin_h2, uin_l3: uin_h3,NVAR) + real(rt) uout(uout_l1:uout_h1,uout_l2:uout_h2,uout_l3:uout_h3,NVAR) + real(rt) grav( gv_l1: gv_h1, gv_l2: gv_h2, gv_l3: gv_h3,3) + real(rt) dx, dy, dz, dt + real(rt) a_old, a_new + real(rt) e_added,ke_added + + real(rt) :: a_half, a_oldsq, a_newsq, a_newsq_inv + real(rt) :: rho + real(rt) :: SrU, SrV, SrW, SrE + real(rt) :: rhoInv, dt_a_new + real(rt) :: old_rhoeint, new_rhoeint, old_ke, new_ke + integer :: i, j, k + + a_half = 0.5d0 * (a_old + a_new) + a_oldsq = a_old * a_old + a_newsq = a_new * a_new + a_newsq_inv = 1.d0 / a_newsq + + dt_a_new = dt / a_new + + ! Gravitational source options for how to add the work to (rho E): + ! grav_source_type = + ! 1: Original version ("does work") + ! 3: Puts all gravitational work into KE, not (rho e) + + ! Add gravitational source terms + do k = lo(3),hi(3) + do j = lo(2),hi(2) + do i = lo(1),hi(1) + + ! **** Start Diagnostics **** + old_ke = 0.5d0 * (uout(i,j,k,UMX)**2 + uout(i,j,k,UMY)**2 + uout(i,j,k,UMZ)**2) / & + uout(i,j,k,URHO) + old_rhoeint = uout(i,j,k,UEDEN) - old_ke + ! **** End Diagnostics **** + + rho = uin(i,j,k,URHO) + rhoInv = 1.0d0 / rho + + SrU = rho * grav(i,j,k,1) + SrV = rho * grav(i,j,k,2) + SrW = rho * grav(i,j,k,3) + + ! We use a_new here because we think of d/dt(a rho u) = ... + (rho g) + uout(i,j,k,UMX) = uout(i,j,k,UMX) + SrU * dt_a_new + uout(i,j,k,UMY) = uout(i,j,k,UMY) + SrV * dt_a_new + uout(i,j,k,UMZ) = uout(i,j,k,UMZ) + SrW * dt_a_new + + if (grav_source_type .eq. 1) then + + ! This does work (in 1-d) + ! Src = rho u dot g, evaluated with all quantities at t^n + SrE = uin(i,j,k,UMX) * grav(i,j,k,1) + & + uin(i,j,k,UMY) * grav(i,j,k,2) + & + uin(i,j,k,UMZ) * grav(i,j,k,3) + uout(i,j,k,UEDEN) = (a_newsq*uout(i,j,k,UEDEN) + SrE * (dt*a_half)) * a_newsq_inv + + else if (grav_source_type .eq. 3) then + + new_ke = 0.5d0 * (uout(i,j,k,UMX)**2 + uout(i,j,k,UMY)**2 + uout(i,j,k,UMZ)**2) / & + uout(i,j,k,URHO) + uout(i,j,k,UEDEN) = old_rhoeint + new_ke + + else + call bl_error("Error:: Nyx_advection_3d.f90 :: bogus grav_source_type") + end if + + ! **** Start Diagnostics **** + ! This is the new (rho e) as stored in (rho E) after the gravitational work is added + new_ke = 0.5d0 * (uout(i,j,k,UMX)**2 + uout(i,j,k,UMY)**2 + uout(i,j,k,UMZ)**2) / & + uout(i,j,k,URHO) + new_rhoeint = uout(i,j,k,UEDEN) - new_ke + + e_added = e_added + (new_rhoeint - old_rhoeint) + ke_added = ke_added + (new_ke - old_ke ) + ! **** End Diagnostics **** enddo enddo enddo - end subroutine divu + ! print *,' EADDED ',lo(1),lo(2),lo(3), e_added + ! print *,'KEADDED ',lo(1),lo(2),lo(3),ke_added + end subroutine add_grav_source diff --git a/Source/HydroFortran/add_grav_source_3d.f90 b/Source/HydroFortran/add_grav_source_3d.f90 index 7f037889..71201f48 100644 --- a/Source/HydroFortran/add_grav_source_3d.f90 +++ b/Source/HydroFortran/add_grav_source_3d.f90 @@ -2,13 +2,11 @@ ! ::: ------------------------------------------------------------------ ! ::: - !=========================================================================== - ! This is called from within threaded loops in advance_gas_tile so *no* OMP here ... - !=========================================================================== subroutine add_grav_source(uin,uin_l1,uin_l2,uin_l3,uin_h1,uin_h2,uin_h3, & uout,uout_l1,uout_l2,uout_l3,uout_h1,uout_h2,uout_h3, & grav, gv_l1, gv_l2, gv_l3, gv_h1, gv_h2, gv_h3, & - lo,hi,dx,dy,dz,dt,a_old,a_new,e_added,ke_added) + lo,hi,dx,dy,dz,dt,a_old,a_new,e_added,ke_added) & + bind(C, name="fort_add_grav_source") use amrex_fort_module, only : rt => amrex_real use eos_module diff --git a/Source/HydroFortran/make_hydro_sources_3d.f90 b/Source/HydroFortran/make_hydro_sources_3d.f90 new file mode 100644 index 00000000..07e23484 --- /dev/null +++ b/Source/HydroFortran/make_hydro_sources_3d.f90 @@ -0,0 +1,143 @@ + +! ::: +! ::: ---------------------------------------------------------------- +! ::: + + subroutine fort_make_hydro_sources(time,lo,hi,& + uin,uin_l1,uin_l2,uin_l3,uin_h1,uin_h2,uin_h3, & + ugdnvx_out,ugdnvx_l1,ugdnvx_l2,ugdnvx_l3,ugdnvx_h1,ugdnvx_h2,ugdnvx_h3, & + ugdnvy_out,ugdnvy_l1,ugdnvy_l2,ugdnvy_l3,ugdnvy_h1,ugdnvy_h2,ugdnvy_h3, & + ugdnvz_out,ugdnvz_l1,ugdnvz_l2,ugdnvz_l3,ugdnvz_h1,ugdnvz_h2,ugdnvz_h3, & + src ,src_l1,src_l2,src_l3,src_h1,src_h2,src_h3, & + hydro_src ,hsrc_l1,hsrc_l2,hsrc_l3,hsrc_h1,hsrc_h2,hsrc_h3, & + divu_cc,d_l1,d_l2,d_l3,d_h1,d_h2,d_h3, & + grav,gv_l1,gv_l2,gv_l3,gv_h1,gv_h2,gv_h3, & + delta,dt, & + flux1,flux1_l1,flux1_l2,flux1_l3,flux1_h1,flux1_h2,flux1_h3, & + flux2,flux2_l1,flux2_l2,flux2_l3,flux2_h1,flux2_h2,flux2_h3, & + flux3,flux3_l1,flux3_l2,flux3_l3,flux3_h1,flux3_h2,flux3_h3, & + courno,a_old,a_new,print_fortran_warnings,do_grav) & + bind(C, name="fort_make_hydro_sources") + + use amrex_fort_module, only : rt => amrex_real + use amrex_mempool_module, only : bl_allocate, bl_deallocate + use meth_params_module, only : QVAR, NVAR, NHYP, normalize_species + use bl_constants_module + + implicit none + + integer lo(3),hi(3),print_fortran_warnings,do_grav + integer uin_l1,uin_l2,uin_l3,uin_h1,uin_h2,uin_h3 + integer ugdnvx_l1,ugdnvx_l2,ugdnvx_l3,ugdnvx_h1,ugdnvx_h2,ugdnvx_h3 + integer ugdnvy_l1,ugdnvy_l2,ugdnvy_l3,ugdnvy_h1,ugdnvy_h2,ugdnvy_h3 + integer ugdnvz_l1,ugdnvz_l2,ugdnvz_l3,ugdnvz_h1,ugdnvz_h2,ugdnvz_h3 + integer flux1_l1,flux1_l2,flux1_l3,flux1_h1,flux1_h2,flux1_h3 + integer flux2_l1,flux2_l2,flux2_l3,flux2_h1,flux2_h2,flux2_h3 + integer flux3_l1,flux3_l2,flux3_l3,flux3_h1,flux3_h2,flux3_h3 + integer src_l1,src_l2,src_l3,src_h1,src_h2,src_h3 + integer hsrc_l1,hsrc_l2,hsrc_l3,hsrc_h1,hsrc_h2,hsrc_h3 + integer d_l1,d_l2,d_l3,d_h1,d_h2,d_h3 + integer gv_l1,gv_l2,gv_l3,gv_h1,gv_h2,gv_h3 + real(rt) uin( uin_l1:uin_h1, uin_l2:uin_h2, uin_l3:uin_h3, NVAR) + real(rt) ugdnvx_out(ugdnvx_l1:ugdnvx_h1,ugdnvx_l2:ugdnvx_h2,ugdnvx_l3:ugdnvx_h3) + real(rt) ugdnvy_out(ugdnvy_l1:ugdnvy_h1,ugdnvy_l2:ugdnvy_h2,ugdnvy_l3:ugdnvy_h3) + real(rt) ugdnvz_out(ugdnvz_l1:ugdnvz_h1,ugdnvz_l2:ugdnvz_h2,ugdnvz_l3:ugdnvz_h3) + real(rt) src( src_l1:src_h1, src_l2:src_h2, src_l3:src_h3, NVAR) + real(rt) hydro_src(hsrc_l1:hsrc_h1,hsrc_l2:hsrc_h2,hsrc_l3:hsrc_h3,NVAR) + real(rt) divu_cc(d_l1:d_h1,d_l2:d_h2,d_l3:d_h3) + real(rt) grav( gv_l1:gv_h1, gv_l2:gv_h2, gv_l3:gv_h3, 3) + real(rt) flux1(flux1_l1:flux1_h1,flux1_l2:flux1_h2, flux1_l3:flux1_h3,NVAR) + real(rt) flux2(flux2_l1:flux2_h1,flux2_l2:flux2_h2, flux2_l3:flux2_h3,NVAR) + real(rt) flux3(flux3_l1:flux3_h1,flux3_l2:flux3_h2, flux3_l3:flux3_h3,NVAR) + real(rt) delta(3),dt,time,courno + real(rt) a_old, a_new + + ! Automatic arrays for workspace + real(rt), pointer :: q(:,:,:,:) + real(rt), pointer :: flatn(:,:,:) + real(rt), pointer :: c(:,:,:) + real(rt), pointer :: csml(:,:,:) + real(rt), pointer :: divu_nd(:,:,:) + real(rt), pointer :: srcQ(:,:,:,:) + + real(rt) dx,dy,dz + integer ngq,ngf + integer q_l1, q_l2, q_l3, q_h1, q_h2, q_h3 + integer srcq_l1, srcq_l2, srcq_l3, srcq_h1, srcq_h2, srcq_h3 + + ngq = NHYP + ngf = 1 + + q_l1 = lo(1)-NHYP + q_l2 = lo(2)-NHYP + q_l3 = lo(3)-NHYP + q_h1 = hi(1)+NHYP + q_h2 = hi(2)+NHYP + q_h3 = hi(3)+NHYP + + srcq_l1 = lo(1)-1 + srcq_l2 = lo(2)-1 + srcq_l3 = lo(3)-1 + srcq_h1 = hi(1)+1 + srcq_h2 = hi(2)+1 + srcq_h3 = hi(3)+1 + + call bl_allocate( q, lo-NHYP, hi+NHYP, QVAR) + call bl_allocate( flatn, lo-NHYP, hi+NHYP ) + call bl_allocate( c, lo-NHYP, hi+NHYP ) + call bl_allocate( csml, lo-NHYP, hi+NHYP ) + + call bl_allocate( srcQ, lo-1, hi+1, QVAR) + call bl_allocate(divu_nd, lo , hi+1) + + dx = delta(1) + dy = delta(2) + dz = delta(3) + + ! 1) Translate conserved variables (u) to primitive variables (q). + ! 2) Compute sound speeds (c) + ! Note that (q,c,csml,flatn) are all dimensioned the same + ! and set to correspond to coordinates of (lo:hi) + ! 3) Translate source terms + call ctoprim(lo,hi,uin,uin_l1,uin_l2,uin_l3,uin_h1,uin_h2,uin_h3, & + q,c,csml,flatn,q_l1,q_l2,q_l3,q_h1,q_h2,q_h3, & + src , src_l1, src_l2, src_l3, src_h1, src_h2, src_h3, & + srcQ,srcq_l1,srcq_l2,srcq_l3,srcq_h1,srcq_h2,srcq_h3, & + grav,gv_l1,gv_l2,gv_l3,gv_h1,gv_h2,gv_h3, & + courno,dx,dy,dz,dt,ngq,ngf,a_old,a_new) + + ! Compute hyperbolic fluxes using unsplit Godunov + call umeth3d(q,c,csml,flatn,q_l1,q_l2,q_l3,q_h1,q_h2,q_h3, & + srcQ,srcq_l1,srcq_l2,srcq_l3,srcq_h1,srcq_h2,srcq_h3, & + lo(1),lo(2),lo(3),hi(1),hi(2),hi(3),dx,dy,dz,dt, & + flux1,flux1_l1,flux1_l2,flux1_l3,flux1_h1,flux1_h2,flux1_h3, & + flux2,flux2_l1,flux2_l2,flux2_l3,flux2_h1,flux2_h2,flux2_h3, & + flux3,flux3_l1,flux3_l2,flux3_l3,flux3_h1,flux3_h2,flux3_h3, & + grav,gv_l1,gv_l2,gv_l3,gv_h1,gv_h2,gv_h3, & + ugdnvx_out,ugdnvx_l1,ugdnvx_l2,ugdnvx_l3,ugdnvx_h1,ugdnvx_h2,ugdnvx_h3, & + ugdnvy_out,ugdnvy_l1,ugdnvy_l2,ugdnvy_l3,ugdnvy_h1,ugdnvy_h2,ugdnvy_h3, & + ugdnvz_out,ugdnvz_l1,ugdnvz_l2,ugdnvz_l3,ugdnvz_h1,ugdnvz_h2,ugdnvz_h3, & + divu_cc,a_old,a_new,print_fortran_warnings) + + ! Compute divergence of velocity field (on surroundingNodes(lo,hi)) + call make_divu_nd(lo,hi,q,q_l1,q_l2,q_l3,q_h1,q_h2,q_h3, & + dx,dy,dz,divu_nd,lo(1),lo(2),lo(3),hi(1)+1,hi(2)+1,hi(3)+1) + + ! Conservative update + call consup(uin,uin_l1,uin_l2,uin_l3,uin_h1,uin_h2,uin_h3, & + hydro_src , hsrc_l1, hsrc_l2, hsrc_l3, hsrc_h1, hsrc_h2, hsrc_h3, & + flux1,flux1_l1,flux1_l2,flux1_l3,flux1_h1,flux1_h2,flux1_h3, & + flux2,flux2_l1,flux2_l2,flux2_l3,flux2_h1,flux2_h2,flux2_h3, & + flux3,flux3_l1,flux3_l2,flux3_l3,flux3_h1,flux3_h2,flux3_h3, & + divu_nd,divu_cc,lo,hi,dx,dy,dz,dt,a_old,a_new) + + ! We are done with these here so can go ahead and free up the space. + call bl_deallocate(q) + call bl_deallocate(flatn) + call bl_deallocate(c) + call bl_deallocate(csml) + call bl_deallocate(divu_nd) + call bl_deallocate(srcQ) + + end subroutine fort_make_hydro_sources + diff --git a/Source/Make.package b/Source/Make.package index 6ec5c214..b6f7c562 100644 --- a/Source/Make.package +++ b/Source/Make.package @@ -20,10 +20,17 @@ CEXE_sources += write_info.cpp CEXE_sources += bl_random_c.cpp ifneq ($(NO_HYDRO), TRUE) -CEXE_sources += Nyx_hydro.cpp +CEXE_sources += compute_hydro_sources.cpp +CEXE_sources += update_state_with_sources.cpp +CEXE_sources += strang_hydro.cpp +CEXE_sources += strang_reactions.cpp CEXE_sources += sum_integrated_quantities.cpp CEXE_sources += sum_utils.cpp -CEXE_sources += strang_splitting.cpp +endif + +ifeq ($(USE_SDC), TRUE) +CEXE_sources += sdc_hydro.cpp +CEXE_sources += sdc_reactions.cpp endif CEXE_headers += Nyx.H diff --git a/Source/Nyx.H b/Source/Nyx.H index 4c015a95..9e530c1e 100644 --- a/Source/Nyx.H +++ b/Source/Nyx.H @@ -38,6 +38,9 @@ enum StateType { #ifdef GRAVITY PhiGrav_Type, Gravity_Type, +#endif +#ifdef SDC + SDC_IR_Type, #endif NUM_STATE_TYPE }; @@ -379,10 +382,33 @@ public: amrex::Real advance_no_hydro(amrex::Real time, amrex::Real dt, int iteration, int ncycle); amrex::Real advance_hydro_plus_particles(amrex::Real time, amrex::Real dt, int iteration, int ncycle); - void just_the_hydro(amrex::Real time, amrex::Real dt, amrex::Real a_old, amrex::Real a_new); + + void strang_hydro(amrex::Real time, amrex::Real dt, amrex::Real a_old, amrex::Real a_new); +#ifdef SDC + void sdc_hydro(amrex::Real time, amrex::Real dt, amrex::Real a_old, amrex::Real a_new); +#endif + + void compute_hydro_sources(amrex::Real time, amrex::Real dt, amrex::Real a_old, amrex::Real a_new, + amrex::MultiFab& S_border, amrex::MultiFab& D_border, + amrex::MultiFab& ext_src_old, amrex::MultiFab& hydro_src, + amrex::MultiFab& grav, amrex::MultiFab& divu_cc, + bool init_flux_register, bool add_to_flux_register); + + void update_state_with_sources( amrex::MultiFab& S_old , amrex::MultiFab& S_new, + amrex::MultiFab& ext_src_old, amrex::MultiFab& hydro_src, + amrex::MultiFab& grav , amrex::MultiFab& divu_cc, + amrex::Real dt, amrex::Real a_old, amrex::Real a_new); + void strang_first_step (amrex::Real time, amrex::Real dt, amrex::MultiFab& state, amrex::MultiFab& dstate); void strang_second_step (amrex::Real time, amrex::Real dt, amrex::MultiFab& state, amrex::MultiFab& dstate); +#ifdef SDC + void sdc_reactions ( amrex::MultiFab& state_old, amrex::MultiFab& state_new, amrex::MultiFab& dstate, + amrex::MultiFab& hydro_src, amrex::MultiFab& IR, + amrex::Real dt, amrex::Real a_old, amrex::Real a_new, + int sdc_iter); +#endif + amrex::Real advance_particles_only (amrex::Real time, amrex::Real dt, int iteration, int ncycle); void moveKickDriftExact(amrex::Real dt); @@ -700,9 +726,13 @@ protected: static int do_forcing; // if true , incorporate the source term through Strang-splitting - // if false, incorporate the source term through predictor-corrector methodology static int strang_split; +#ifdef SDC + // if true , use SDC to couple hydro and reactions + static int sdc_split; +#endif + #ifdef GRAVITY // There can be only one Gravity object, it covers all levels: static class Gravity *gravity; @@ -765,9 +795,9 @@ Nyx::num_grow() inline Nyx& -Nyx::get_level(int level) +Nyx::get_level(int my_level) { - return *(Nyx *) &parent->getLevel(level); + return *(Nyx *) &parent->getLevel(my_level); } #ifndef NO_HYDRO @@ -781,9 +811,9 @@ Nyx::get_flux_reg() inline amrex::FluxRegister& -Nyx::get_flux_reg(int level) +Nyx::get_flux_reg(int my_level) { - return get_level(level).get_flux_reg(); + return get_level(my_level).get_flux_reg(); } #endif // NO_HYDRO diff --git a/Source/Nyx_F.H b/Source/Nyx_F.H index 730eb794..47035359 100644 --- a/Source/Nyx_F.H +++ b/Source/Nyx_F.H @@ -138,6 +138,47 @@ extern "C" const amrex::Real dx[], const amrex::Real glo[], const amrex::Real* time, const int bc[]); + void fort_make_hydro_sources + (const amrex::Real* time, + const int lo[], const int hi[], + const BL_FORT_FAB_ARG(state), + const BL_FORT_FAB_ARG(ugdnvx), + const BL_FORT_FAB_ARG(ugdnvy), + const BL_FORT_FAB_ARG(ugdnvz), + const BL_FORT_FAB_ARG(src), + BL_FORT_FAB_ARG(hydro_src), + BL_FORT_FAB_ARG(divu_cc), + const BL_FORT_FAB_ARG(grav), + const amrex::Real dx[], const amrex::Real* dt, + D_DECL(const BL_FORT_FAB_ARG(xflux), + const BL_FORT_FAB_ARG(yflux), + const BL_FORT_FAB_ARG(zflux)), + const amrex::Real* cflLoc, + const amrex::Real* a_old, const amrex::Real* a_new, + const int* print_fortran_warnings, + const int* do_grav); + + void fort_update_state + ( const int lo[], const int hi[], + const BL_FORT_FAB_ARG(u_in), + BL_FORT_FAB_ARG(u_out), + const BL_FORT_FAB_ARG(src), + const BL_FORT_FAB_ARG(hydro_src), + const BL_FORT_FAB_ARG(divu_cc), + const amrex::Real* dt, + const amrex::Real* a_old, + const amrex::Real* a_new, + const int* print_fortran_warnings); + + void fort_add_grav_source + ( const int lo[], const int hi[], + const BL_FORT_FAB_ARG(u_in), + BL_FORT_FAB_ARG(u_out), + const BL_FORT_FAB_ARG(grav), + const amrex::Real* dt, + const amrex::Real* a_old, + const amrex::Real* a_new); + void fort_advance_gas (const amrex::Real* time, const int lo[], const int hi[], @@ -248,7 +289,6 @@ extern "C" (const int* lo, const int* hi, BL_FORT_FAB_ARG(state), BL_FORT_FAB_ARG(diag_eos), - const amrex::Real* dx, const amrex::Real* time, const amrex::Real* z, const amrex::Real* dt, const int* min_iter, const int* max_iter); diff --git a/Source/Nyx_advance.cpp b/Source/Nyx_advance.cpp index 15595c16..449aa479 100644 --- a/Source/Nyx_advance.cpp +++ b/Source/Nyx_advance.cpp @@ -285,7 +285,16 @@ Nyx::advance_hydro_plus_particles (Real time, BL_PROFILE_VAR("just_the_hydro", just_the_hydro); for (int lev = level; lev <= finest_level_to_advance; lev++) { - get_level(lev).just_the_hydro(time, dt, a_old, a_new); +#ifdef SDC + if (sdc_split > 0) + { + get_level(lev).sdc_hydro(time, dt, a_old, a_new); + } else { + get_level(lev).strang_hydro(time, dt, a_old, a_new); + } +#else + get_level(lev).strang_hydro(time, dt, a_old, a_new); +#endif } BL_PROFILE_VAR_STOP(just_the_hydro); @@ -527,7 +536,18 @@ Nyx::advance_hydro (Real time, #endif // Call the hydro advance itself - just_the_hydro(time, dt, a_old, a_new); + BL_PROFILE_VAR("just_the_hydro", just_the_hydro); +#ifdef SDC + if (sdc_split > 0) + { + sdc_hydro(time, dt, a_old, a_new); + } else { + strang_hydro(time, dt, a_old, a_new); + } +#else + strang_hydro(time, dt, a_old, a_new); +#endif + BL_PROFILE_VAR_STOP(just_the_hydro); MultiFab& S_new = get_new_data(State_Type); MultiFab& D_new = get_new_data(DiagEOS_Type); diff --git a/Source/compute_hydro_sources.cpp b/Source/compute_hydro_sources.cpp new file mode 100644 index 00000000..3189a3a7 --- /dev/null +++ b/Source/compute_hydro_sources.cpp @@ -0,0 +1,116 @@ +#include "Nyx.H" +#include "Nyx_F.H" + +using namespace amrex; + +using std::string; + +void +Nyx::compute_hydro_sources(amrex::Real time, amrex::Real dt, amrex::Real a_old, amrex::Real a_new, + MultiFab& S_border, MultiFab& D_border, + MultiFab& ext_src_old, MultiFab& hydro_src, + MultiFab& grav_vector, MultiFab& divu_cc, + bool init_flux_register, bool add_to_flux_register) +{ + amrex::Print() << "Computing the hydro sources ... " << std::endl; + const int finest_level = parent->finestLevel(); + const Real* dx = geom.CellSize(); + Real courno = -1.0e+200; + + MultiFab fluxes[BL_SPACEDIM]; + for (int j = 0; j < BL_SPACEDIM; j++) + { + fluxes[j].define(getEdgeBoxArray(j), dmap, NUM_STATE, 0); + fluxes[j].setVal(0.0); + } + + // + // Get pointers to Flux registers, or set pointer to zero if not there. + // + FluxRegister* fine = 0; + FluxRegister* current = 0; + + if (do_reflux) + { + if (level < finest_level) + { + fine = &get_flux_reg(level+1); + if (init_flux_register) + fine->setVal(0); + + } else if (level > 0) { + current = &get_flux_reg(level); + } + } + +#ifdef _OPENMP +#pragma omp parallel reduction(max:courno) +#endif + { + FArrayBox flux[BL_SPACEDIM], u_gdnv[BL_SPACEDIM]; + Real cflLoc = -1.e+200; + + for (MFIter mfi(S_border,true); mfi.isValid(); ++mfi) + { + const Box& bx = mfi.tilebox(); + + FArrayBox& state = S_border[mfi]; + FArrayBox& dstate = D_border[mfi]; + + // Allocate fabs for fluxes. + for (int i = 0; i < BL_SPACEDIM ; i++) { + const Box &bxtmp = amrex::surroundingNodes(bx, i); + flux[i].resize(bxtmp, NUM_STATE); + u_gdnv[i].resize(amrex::grow(bxtmp, 1), 1); + u_gdnv[i].setVal(1.e200); + } + fort_make_hydro_sources + (&time, bx.loVect(), bx.hiVect(), + BL_TO_FORTRAN(state), + BL_TO_FORTRAN(u_gdnv[0]), + BL_TO_FORTRAN(u_gdnv[1]), + BL_TO_FORTRAN(u_gdnv[2]), + BL_TO_FORTRAN(ext_src_old[mfi]), + BL_TO_FORTRAN(hydro_src[mfi]), + BL_TO_FORTRAN(divu_cc[mfi]), + BL_TO_FORTRAN(grav_vector[mfi]), + dx, &dt, + BL_TO_FORTRAN(flux[0]), + BL_TO_FORTRAN(flux[1]), + BL_TO_FORTRAN(flux[2]), + &cflLoc, &a_old, &a_new, + &print_fortran_warnings, &do_grav); + + for (int i = 0; i < BL_SPACEDIM; ++i) + fluxes[i][mfi].copy(flux[i], mfi.nodaltilebox(i)); + + } // end of MFIter loop + + courno = std::max(courno, cflLoc); + + } // end of parallel + + if (add_to_flux_register) + { + if (do_reflux) { + if (current) { + for (int i = 0; i < BL_SPACEDIM ; i++) { + current->FineAdd(fluxes[i], i, 0, 0, NUM_STATE, 1); + } + } + if (fine) { + for (int i = 0; i < BL_SPACEDIM ; i++) { + fine->CrseInit(fluxes[i],i,0,0,NUM_STATE,-1.,FluxRegister::ADD); + } + } + } + } + + if (courno > 1.0) + { + if (ParallelDescriptor::IOProcessor()) + std::cout << "OOPS -- EFFECTIVE CFL AT THIS LEVEL " << level + << " IS " << courno << '\n'; + amrex::Abort("CFL is too high at this level -- go back and restart with lower cfl number"); + } +} diff --git a/Source/sdc_hydro.cpp b/Source/sdc_hydro.cpp new file mode 100644 index 00000000..db4a6f78 --- /dev/null +++ b/Source/sdc_hydro.cpp @@ -0,0 +1,132 @@ + +#include "Nyx.H" +#include "Nyx_F.H" +#include + +#ifdef GRAVITY +#include "Gravity.H" +#endif + +using namespace amrex; + +using std::string; + +#ifndef NO_HYDRO +void +Nyx::sdc_hydro (Real time, + Real dt, + Real a_old, + Real a_new) +{ + BL_PROFILE("Nyx::sdc_hydro()"); + + BL_ASSERT(NUM_GROW == 4); + + int sdc_iter; + Real IR_fac; + const Real prev_time = state[State_Type].prevTime(); + const Real cur_time = state[State_Type].curTime(); + + MultiFab& S_old = get_old_data(State_Type); + MultiFab& S_new = get_new_data(State_Type); + + MultiFab& D_old = get_old_data(DiagEOS_Type); + MultiFab& D_new = get_new_data(DiagEOS_Type); + + MultiFab& IR_old = get_old_data(SDC_IR_Type); + MultiFab& IR_new = get_new_data(SDC_IR_Type); + +#if 0 + VisMF::Write(IR_old,"IR_OLD"); +#endif + + // It's possible for interpolation to create very small negative values for + // species so we make sure here that all species are non-negative after this + // point + enforce_nonnegative_species(S_old); + + MultiFab ext_src_old(grids, dmap, NUM_STATE, 3); + ext_src_old.setVal(0.); + + // Define the gravity vector + MultiFab grav_vector(grids, dmap, BL_SPACEDIM, 3); + grav_vector.setVal(0.); + +#ifdef GRAVITY + gravity->get_old_grav_vector(level, grav_vector, time); + grav_vector.FillBoundary(geom.periodicity()); +#endif + + // Create FAB for extended grid values (including boundaries) and fill. + MultiFab S_old_tmp(S_old.boxArray(), S_old.DistributionMap(), NUM_STATE, NUM_GROW); + FillPatch(*this, S_old_tmp, NUM_GROW, time, State_Type, 0, NUM_STATE); + + MultiFab D_old_tmp(D_old.boxArray(), D_old.DistributionMap(), D_old.nComp(), NUM_GROW); + FillPatch(*this, D_old_tmp, NUM_GROW, time, DiagEOS_Type, 0, D_old.nComp()); + + MultiFab hydro_src(grids, dmap, NUM_STATE, 0); + hydro_src.setVal(0.); + + MultiFab divu_cc(grids, dmap, 1, 0); + divu_cc.setVal(0.); + + //Begin loop over SDC iterations + int sdc_iter_max = 2; + + for (sdc_iter = 0; sdc_iter < sdc_iter_max; sdc_iter++) + { + amrex::Print() << "STARTING SDC_ITER LOOP " << sdc_iter << std::endl; + + // Here we bundle all the source terms into ext_src_old + // FillPatch the IR term into ext_src_old for both Eint and Eden + MultiFab IR_tmp(grids, dmap, 1, NUM_GROW); + FillPatch(*this, IR_tmp, NUM_GROW, time, SDC_IR_Type, 0, 1); + + MultiFab::Add(ext_src_old,IR_tmp,0,Eden,1,0); + MultiFab::Add(ext_src_old,IR_tmp,0,Eint,1,0); + + ext_src_old.FillBoundary(geom.periodicity()); + + bool init_flux_register = (sdc_iter == 0); + bool add_to_flux_register = (sdc_iter == sdc_iter_max-1); + compute_hydro_sources(time,dt,a_old,a_new,S_old_tmp,D_old_tmp, + ext_src_old,hydro_src,grav_vector,divu_cc, + init_flux_register, add_to_flux_register); + + // We subtract IR_tmp before we add the rest of the source term to (rho e) and (rho E) + MultiFab::Subtract(ext_src_old,IR_tmp,0,Eden,1,0); + MultiFab::Subtract(ext_src_old,IR_tmp,0,Eint,1,0); + update_state_with_sources(S_old_tmp,S_new, + ext_src_old,hydro_src,grav_vector,divu_cc, + dt,a_old,a_new); + + // We copy old Temp and Ne to new Temp and Ne so that they can be used + // as guesses when we next need them. + MultiFab::Copy(D_new,D_old,0,0,D_old.nComp(),0); + + // This step needs to do the update of (rho), (rho e) and (rho E) + // AND needs to return an updated value of I_R in the old SDC_IR statedata. + BL_PROFILE_VAR("sdc_reactions", sdc_reactions); + sdc_reactions(S_old_tmp, S_new, D_new, hydro_src, IR_old, dt, a_old, a_new, sdc_iter); + BL_PROFILE_VAR_STOP(sdc_reactions); + +#if 0 + if (sdc_iter == 0) + VisMF::Write(IR_old,"IR_AFTER_SDC1"); + if (sdc_iter == 1) + VisMF::Write(IR_old,"IR_AFTER_SDC2"); +#endif + + // Internal to sdc_reactions, we have added a_half * IR to asq*(rho e) and asq*(rho E) + // I_R satisfies the equation anewsq * (rho_out e_out ) = + // aoldsq * (rho_orig e_orig) + dt * a_half * I_R + dt * H_{rho e} + + } //End loop over SDC iterations + + // Copy IR_old (the current IR) into IR_new here so that when the pointer swap occurs + // we can use IR in the next timestep + MultiFab::Copy(IR_new,IR_old,0,0,1,0); + + grav_vector.clear(); +} +#endif diff --git a/Source/sdc_reactions.cpp b/Source/sdc_reactions.cpp new file mode 100644 index 00000000..e016075b --- /dev/null +++ b/Source/sdc_reactions.cpp @@ -0,0 +1,96 @@ + +#include "Nyx.H" +#include "Nyx_F.H" + +using namespace amrex; +using std::string; + +void +Nyx::sdc_reactions (MultiFab& S_old, MultiFab& S_new, MultiFab& D_new, + MultiFab& hydro_src, MultiFab& IR, + Real delta_time, Real a_old, Real a_new, int sdc_iter) +{ + BL_PROFILE("Nyx::sdc_reactions()"); + + const Real* dx = geom.CellSize(); + + MultiFab reset_src(grids, dmap, 1, NUM_GROW); + reset_src.setVal(0.0); + compute_new_temp(S_new,D_new,reset_src); + // MultiFab reset_src(grids, dmap, 1, NUM_GROW); + // reset_src.setVal(0.0); + +#ifndef FORCING + { + const Real z = 1.0/a_old - 1.0; + fort_interp_to_this_z(&z); + } +#endif + + int min_iter = 100000; + int max_iter = 0; + + int min_iter_grid, max_iter_grid; + + /////////////////////Consider adding ifdefs for whether CVODE is compiled in for these statements + if(heat_cool_type == 3) + { +#ifdef _OPENMP +#pragma omp parallel +#endif + for (MFIter mfi(S_old,true); mfi.isValid(); ++mfi) + { + // Note that this "bx" is only the valid region (unlike for Strang) + const Box& bx = mfi.tilebox(); + + min_iter_grid = 100000; + max_iter_grid = 0; + + integrate_state_with_source + (bx.loVect(), bx.hiVect(), + BL_TO_FORTRAN(S_old[mfi]), + BL_TO_FORTRAN(S_new[mfi]), + BL_TO_FORTRAN(D_new[mfi]), + BL_TO_FORTRAN(hydro_src[mfi]), + BL_TO_FORTRAN(reset_src[mfi]), + BL_TO_FORTRAN(IR[mfi]), + &a_old, &delta_time, &min_iter_grid, &max_iter_grid); + + min_iter = std::min(min_iter,min_iter_grid); + max_iter = std::max(max_iter,max_iter_grid); + } + } + else if(heat_cool_type == 5) + { +#ifdef _OPENMP +#pragma omp parallel +#endif + for (MFIter mfi(S_old,true); mfi.isValid(); ++mfi) + { + // Note that this "bx" is only the valid region (unlike for Strang) + const Box& bx = mfi.tilebox(); + + min_iter_grid = 100000; + max_iter_grid = 0; + + integrate_state_fcvode_with_source + (bx.loVect(), bx.hiVect(), + BL_TO_FORTRAN(S_old[mfi]), + BL_TO_FORTRAN(S_new[mfi]), + BL_TO_FORTRAN(D_new[mfi]), + BL_TO_FORTRAN(hydro_src[mfi]), + BL_TO_FORTRAN(reset_src[mfi]), + BL_TO_FORTRAN(IR[mfi]), + &a_old, &delta_time, &min_iter_grid, &max_iter_grid); + + min_iter = std::min(min_iter,min_iter_grid); + max_iter = std::max(max_iter,max_iter_grid); + } + + } + + ParallelDescriptor::ReduceIntMax(max_iter); + ParallelDescriptor::ReduceIntMin(min_iter); + + amrex::Print() << "Min/Max Number of Iterations in SDC: " << min_iter << " " << max_iter << std::endl; +} diff --git a/Source/Nyx_hydro.cpp b/Source/strang_hydro.cpp similarity index 94% rename from Source/Nyx_hydro.cpp rename to Source/strang_hydro.cpp index c50b9825..56925ada 100644 --- a/Source/Nyx_hydro.cpp +++ b/Source/strang_hydro.cpp @@ -13,12 +13,12 @@ using std::string; #ifndef NO_HYDRO void -Nyx::just_the_hydro (Real time, - Real dt, - Real a_old, - Real a_new) +Nyx::strang_hydro (Real time, + Real dt, + Real a_old, + Real a_new) { - BL_PROFILE("Nyx::just_the_hydro()"); + BL_PROFILE("Nyx::strang_hydro()"); const Real prev_time = state[State_Type].prevTime(); const Real cur_time = state[State_Type].curTime(); @@ -32,10 +32,10 @@ Nyx::just_the_hydro (Real time, { if (ParallelDescriptor::IOProcessor()) { - std::cout << "just_the_hydro: prev_time = " << prev_time << std::endl; - std::cout << "just_the_hydro: time = " << time << std::endl; + std::cout << "strang_hydro: prev_time = " << prev_time << std::endl; + std::cout << "strang_hydro: time = " << time << std::endl; } - amrex::Abort("time should equal prev_time in just_the_hydro!"); + amrex::Abort("time should equal prev_time in strang_hydro!"); } #ifndef NDEBUG @@ -44,7 +44,7 @@ Nyx::just_the_hydro (Real time, for (int i = 0; i < S_old.nComp(); i++) { if (ParallelDescriptor::IOProcessor()) - std::cout << "just_the_hydro: testing component " << i << " for NaNs" << std::endl; + std::cout << "strang_hydro: testing component " << i << " for NaNs" << std::endl; if (S_old.contains_nan(Density+i,1,0)) amrex::Abort("S_old has NaNs in this component"); } diff --git a/Source/strang_splitting.cpp b/Source/strang_reactions.cpp similarity index 94% rename from Source/strang_splitting.cpp rename to Source/strang_reactions.cpp index 4d831fd9..432b2f85 100644 --- a/Source/strang_splitting.cpp +++ b/Source/strang_reactions.cpp @@ -36,7 +36,7 @@ Nyx::strang_first_step (Real time, Real dt, MultiFab& S_old, MultiFab& D_old) (bx.loVect(), bx.hiVect(), BL_TO_FORTRAN(S_old[mfi]), BL_TO_FORTRAN(D_old[mfi]), - dx, &time, &a, &half_dt, &min_iter, &max_iter); + &a, &half_dt, &min_iter, &max_iter); #ifndef NDEBUG if (S_old[mfi].contains_nan()) @@ -85,7 +85,7 @@ Nyx::strang_second_step (Real time, Real dt, MultiFab& S_new, MultiFab& D_new) (bx.loVect(), bx.hiVect(), BL_TO_FORTRAN(S_new[mfi]), BL_TO_FORTRAN(D_new[mfi]), - dx, &time, &a, &half_dt, &min_iter_grid, &max_iter_grid); + &a, &half_dt, &min_iter_grid, &max_iter_grid); if (S_new[mfi].contains_nan(bx,0,S_new.nComp())) { diff --git a/Source/update_state_with_sources.cpp b/Source/update_state_with_sources.cpp new file mode 100644 index 00000000..9cbceec1 --- /dev/null +++ b/Source/update_state_with_sources.cpp @@ -0,0 +1,41 @@ +#include "Nyx.H" +#include "Nyx_F.H" + +using namespace amrex; + +using std::string; + +void +Nyx::update_state_with_sources( MultiFab& S_old, MultiFab& S_new, + MultiFab& ext_src_old, MultiFab& hydro_src, + MultiFab& grav, MultiFab& divu_cc, + amrex::Real dt, amrex::Real a_old, amrex::Real a_new) +{ + amrex::Print() << "Updating state with the hydro sources ... " << std::endl; + +#ifdef _OPENMP +#pragma omp parallel +#endif + for (MFIter mfi(S_old,true); mfi.isValid(); ++mfi) + { + const Box& bx = mfi.tilebox(); + fort_update_state ( + bx.loVect(), bx.hiVect(), + BL_TO_FORTRAN(S_old[mfi]), + BL_TO_FORTRAN(S_new[mfi]), + BL_TO_FORTRAN(ext_src_old[mfi]), + BL_TO_FORTRAN(hydro_src[mfi]), + BL_TO_FORTRAN(divu_cc[mfi]), + &dt, &a_old, &a_new, &print_fortran_warnings); + + // Note this increments S_new, it doesn't add source to S_old + // However we create the source term using rho_old + if (do_grav) + fort_add_grav_source ( + bx.loVect(), bx.hiVect(), + BL_TO_FORTRAN(S_old[mfi]), + BL_TO_FORTRAN(S_new[mfi]), + BL_TO_FORTRAN(grav[mfi]), + &dt, &a_old, &a_new); + } +} From b53170510fff6b9e31416012fcc9ec0c50e7a829 Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Wed, 2 May 2018 12:08:40 -0700 Subject: [PATCH 08/12] More sync-ing of current code with SDC version. Specifically we separate the call to reset_internal_energy out of the call to compute_new_temp. In the call to reset_internal_energy we pass reset_e_src which holds the amount by which we have changed e in the process of resetting it. --- Source/HydroFortran/Nyx_advection_3d.f90 | 21 +--- Source/Initialization/Nyx_initdata.cpp | 25 ++++- Source/Nyx.H | 15 +-- Source/Nyx.cpp | 68 ++++++++----- Source/Nyx_F.H | 4 +- Source/Nyx_advance.cpp | 113 +++++---------------- Source/Nyx_nd.f90 | 3 + Source/SourceTerms/Nyx_grav_sources_3d.f90 | 14 +-- Source/sdc_reactions.cpp | 15 +-- Source/strang_hydro.cpp | 83 +-------------- Source/strang_reactions.cpp | 5 +- Source/sum_integrated_quantities.cpp | 10 +- Source/write_info.cpp | 8 +- 13 files changed, 142 insertions(+), 242 deletions(-) diff --git a/Source/HydroFortran/Nyx_advection_3d.f90 b/Source/HydroFortran/Nyx_advection_3d.f90 index bd320066..e55a290f 100644 --- a/Source/HydroFortran/Nyx_advection_3d.f90 +++ b/Source/HydroFortran/Nyx_advection_3d.f90 @@ -15,7 +15,7 @@ subroutine fort_advance_gas(time,lo,hi,& flux1,flux1_l1,flux1_l2,flux1_l3,flux1_h1,flux1_h2,flux1_h3, & flux2,flux2_l1,flux2_l2,flux2_l3,flux2_h1,flux2_h2,flux2_h3, & flux3,flux3_l1,flux3_l2,flux3_l3,flux3_h1,flux3_h2,flux3_h3, & - courno,a_old,a_new,e_added,ke_added,print_fortran_warnings,do_grav) & + courno,a_old,a_new,print_fortran_warnings,do_grav) & bind(C, name="fort_advance_gas") use amrex_fort_module, only : rt => amrex_real @@ -49,7 +49,6 @@ subroutine fort_advance_gas(time,lo,hi,& real(rt) flux3(flux3_l1:flux3_h1,flux3_l2:flux3_h2, flux3_l3:flux3_h3,NVAR) real(rt) delta(3),dt,time,courno real(rt) a_old, a_new - real(rt) e_added,ke_added ! Automatic arrays for workspace real(rt), pointer :: q(:,:,:,:) @@ -152,7 +151,7 @@ subroutine fort_advance_gas(time,lo,hi,& call add_grav_source(uin,uin_l1,uin_l2,uin_l3,uin_h1,uin_h2,uin_h3, & uout,uout_l1,uout_l2,uout_l3,uout_h1,uout_h2,uout_h3, & grav, gv_l1, gv_l2, gv_l3, gv_h1, gv_h2, gv_h3, & - lo,hi,dx,dy,dz,dt,a_old,a_new,e_added,ke_added) + lo,hi,dx,dy,dz,dt,a_old,a_new) ! Enforce species >= 0 call enforce_nonnegative_species(uout,uout_l1,uout_l2,uout_l3, & @@ -1987,7 +1986,7 @@ end subroutine make_divu_nd subroutine add_grav_source(uin,uin_l1,uin_l2,uin_l3,uin_h1,uin_h2,uin_h3, & uout,uout_l1,uout_l2,uout_l3,uout_h1,uout_h2,uout_h3, & grav, gv_l1, gv_l2, gv_l3, gv_h1, gv_h2, gv_h3, & - lo,hi,dx,dy,dz,dt,a_old,a_new,e_added,ke_added) + lo,hi,dx,dy,dz,dt,a_old,a_new) use amrex_fort_module, only : rt => amrex_real use eos_module @@ -2006,7 +2005,6 @@ subroutine add_grav_source(uin,uin_l1,uin_l2,uin_l3,uin_h1,uin_h2,uin_h3, & real(rt) grav( gv_l1: gv_h1, gv_l2: gv_h2, gv_l3: gv_h3,3) real(rt) dx, dy, dz, dt real(rt) a_old, a_new - real(rt) e_added,ke_added real(rt) :: a_half, a_oldsq, a_newsq, a_newsq_inv real(rt) :: rho @@ -2069,21 +2067,8 @@ subroutine add_grav_source(uin,uin_l1,uin_l2,uin_l3,uin_h1,uin_h2,uin_h3, & call bl_error("Error:: Nyx_advection_3d.f90 :: bogus grav_source_type") end if - ! **** Start Diagnostics **** - ! This is the new (rho e) as stored in (rho E) after the gravitational work is added - new_ke = 0.5d0 * (uout(i,j,k,UMX)**2 + uout(i,j,k,UMY)**2 + uout(i,j,k,UMZ)**2) / & - uout(i,j,k,URHO) - new_rhoeint = uout(i,j,k,UEDEN) - new_ke - - e_added = e_added + (new_rhoeint - old_rhoeint) - ke_added = ke_added + (new_ke - old_ke ) - ! **** End Diagnostics **** - enddo enddo enddo - ! print *,' EADDED ',lo(1),lo(2),lo(3), e_added - ! print *,'KEADDED ',lo(1),lo(2),lo(3),ke_added - end subroutine add_grav_source diff --git a/Source/Initialization/Nyx_initdata.cpp b/Source/Initialization/Nyx_initdata.cpp index 76ceec0f..3f1056e9 100644 --- a/Source/Initialization/Nyx_initdata.cpp +++ b/Source/Initialization/Nyx_initdata.cpp @@ -180,9 +180,10 @@ Nyx::initData () return; } + MultiFab& S_new = get_new_data(State_Type); + #ifndef NO_HYDRO // We need this because otherwise we might operate on uninitialized data. - MultiFab& S_new = get_new_data(State_Type); S_new.setVal(0.0); #endif @@ -226,7 +227,12 @@ Nyx::initData () if (inhomo_reion) init_zhi(); - compute_new_temp(); + // First reset internal energy before call to compute_temp + MultiFab reset_e_src(grids, dmap, 1, NUM_GROW); + reset_e_src.setVal(0.0); + + reset_internal_energy(S_new,D_new,reset_e_src); + compute_new_temp (S_new,D_new); enforce_consistent_e(S_new); } else @@ -268,6 +274,14 @@ Nyx::initData () #endif +#ifdef SDC + // + // Initialize this to zero before we use it in advance + // + MultiFab& IR_new = get_new_data(SDC_IR_Type); + IR_new.setVal(0.0); +#endif + #ifndef NO_HYDRO // // Read in initial conditions from a file. @@ -283,10 +297,10 @@ Nyx::initData () VisMF::Read(mf, mfDirName.c_str()); - MultiFab& S_new = get_level(0).get_new_data(State_Type); + MultiFab& S_new_crse = get_level(0).get_new_data(State_Type); - S_new.copy(mf, 0, 0, 6); - S_new.copy(mf, 0, FirstSpec, 1); + S_new_crse.copy(mf, 0, 0, 6); + S_new_crse.copy(mf, 0, FirstSpec, 1); if (do_hydro == 1) { @@ -430,4 +444,5 @@ Nyx::init_from_plotfile () std::cout << "Done initializing the particles from the plotfile " << std::endl; std::cout << " " << std::endl; } + } diff --git a/Source/Nyx.H b/Source/Nyx.H index 9e530c1e..0940bb40 100644 --- a/Source/Nyx.H +++ b/Source/Nyx.H @@ -389,12 +389,12 @@ public: #endif void compute_hydro_sources(amrex::Real time, amrex::Real dt, amrex::Real a_old, amrex::Real a_new, - amrex::MultiFab& S_border, amrex::MultiFab& D_border, - amrex::MultiFab& ext_src_old, amrex::MultiFab& hydro_src, + amrex::MultiFab& S_border, amrex::MultiFab& D_border, + amrex::MultiFab& ext_src_old, amrex::MultiFab& hydro_src, amrex::MultiFab& grav, amrex::MultiFab& divu_cc, bool init_flux_register, bool add_to_flux_register); - void update_state_with_sources( amrex::MultiFab& S_old , amrex::MultiFab& S_new, + void update_state_with_sources( amrex::MultiFab& S_old , amrex::MultiFab& S_new, amrex::MultiFab& ext_src_old, amrex::MultiFab& hydro_src, amrex::MultiFab& grav , amrex::MultiFab& divu_cc, amrex::Real dt, amrex::Real a_old, amrex::Real a_new); @@ -404,8 +404,8 @@ public: #ifdef SDC void sdc_reactions ( amrex::MultiFab& state_old, amrex::MultiFab& state_new, amrex::MultiFab& dstate, - amrex::MultiFab& hydro_src, amrex::MultiFab& IR, - amrex::Real dt, amrex::Real a_old, amrex::Real a_new, + amrex::MultiFab& hydro_src, amrex::MultiFab& IR, + amrex::Real dt, amrex::Real a_old, amrex::Real a_new, int sdc_iter); #endif @@ -506,9 +506,10 @@ public: static int num_grow(); // Synchronize (rho e) and (rho E) so they are consistent with each other - void reset_internal_energy(amrex::MultiFab& State, amrex::MultiFab& DiagEOS); + void reset_internal_energy(amrex::MultiFab& State, amrex::MultiFab& DiagEOS, amrex::MultiFab& reset_e_src); - void compute_new_temp(); + // Note: this no longer includes the call to reset_internal_energy + void compute_new_temp(amrex::MultiFab& S_new, amrex::MultiFab& D_new); void compute_rho_temp(amrex::Real& rho_T_avg, amrex::Real& T_avg, amrex::Real& Tinv_avg, amrex::Real& T_meanrho); void compute_gas_fractions(amrex::Real T_cut, amrex::Real rho_cut, diff --git a/Source/Nyx.cpp b/Source/Nyx.cpp index fa2c28d8..60ba2beb 100644 --- a/Source/Nyx.cpp +++ b/Source/Nyx.cpp @@ -123,6 +123,9 @@ int Nyx::do_hydro = -1; int Nyx::add_ext_src = 0; int Nyx::heat_cool_type = 0; int Nyx::strang_split = 0; +#ifdef SDC +int Nyx::sdc_split = 0; +#endif Real Nyx::average_gas_density = 0; Real Nyx::average_dm_density = 0; @@ -358,6 +361,18 @@ Nyx::read_params () pp.query("add_ext_src", add_ext_src); pp.query("strang_split", strang_split); +#ifdef SDC + pp.query("sdc_split", sdc_split); + if (sdc_split == 1 && strang_split == 1) + amrex::Error("Cant have strang_split == 1 and sdc_split == 1"); + if (sdc_split == 0 && strang_split == 0) + amrex::Error("Cant have strang_split == 0 and sdc_split == 0"); + if (sdc_split != 1 && strang_split != 1) + amrex::Error("Cant have strang_split != 1 and sdc_split != 1"); +#else + if (strang_split != 1) + amrex::Error("Cant have strang_split != 1 with USE_SDC != TRUE"); +#endif #ifdef FORCING pp.get("do_forcing", do_forcing); @@ -402,17 +417,14 @@ Nyx::read_params () #ifdef HEATCOOL if (heat_cool_type > 0 && add_ext_src == 0) amrex::Error("Nyx::must set add_ext_src to 1 if heat_cool_type > 0"); - if (heat_cool_type != 1 && heat_cool_type != 3 && heat_cool_type != 5 && heat_cool_type != 7) - amrex::Error("Nyx:: nonzero heat_cool_type must equal 1 or 3 or 5 or 7"); + if (heat_cool_type != 3 && heat_cool_type != 5 && heat_cool_type != 7) + amrex::Error("Nyx:: nonzero heat_cool_type must equal 3 or 5 or 7"); if (heat_cool_type == 0) amrex::Error("Nyx::contradiction -- HEATCOOL is defined but heat_cool_type == 0"); if (ParallelDescriptor::IOProcessor()) { std::cout << "Integrating heating/cooling method with the following method: "; switch (heat_cool_type) { - case 1: - std::cout << "HC"; - break; case 3: std::cout << "VODE"; break; @@ -429,6 +441,9 @@ Nyx::read_params () #ifndef AMREX_USE_CVODE if (heat_cool_type == 5 || heat_cool_type == 7) amrex::Error("Nyx:: cannot set heat_cool_type = 5 or 7 unless USE_CVODE=TRUE"); +#else + if (heat_cool_type == 7 && sdc_split == 1) + amrex::Error("Nyx:: cannot set heat_cool_type = 7 with sdc_split = 1"); #endif #else @@ -617,7 +632,7 @@ Nyx::Nyx (Amr& papa, #ifdef HEATCOOL // Initialize "this_z" in the atomic_rates_module - if (heat_cool_type == 1 || heat_cool_type == 3 || heat_cool_type == 5 || heat_cool_type == 7) + if (heat_cool_type == 3 || heat_cool_type == 5 || heat_cool_type == 7) fort_interp_to_this_z(&initial_z); #endif @@ -760,9 +775,14 @@ Nyx::init (AmrLevel& old) } #endif - // Set E in terms of e + kinetic energy - // if (do_hydro) - // enforce_consistent_e(S_new); +#ifdef SDC + MultiFab& IR_new = get_new_data(SDC_IR_Type); + for (FillPatchIterator fpi(old, IR_new, 0, cur_time, SDC_IR_Type, 0, 1); + fpi.isValid(); ++fpi) + { + IR_new[fpi].copy(fpi()); + } +#endif } // @@ -800,10 +820,6 @@ Nyx::init () // We set dt to be large for this new level to avoid screwing up // computeNewDt. parent->setDtLevel(1.e100, level); - - // Set E in terms of e + kinetic energy - // if (do_hydro) - // enforce_consistent_e(S_new); } Real @@ -1439,8 +1455,16 @@ Nyx::post_timestep (int iteration) #ifndef NO_HYDRO if (do_hydro) { + MultiFab& S_new = get_new_data(State_Type); + MultiFab& D_new = get_new_data(DiagEOS_Type); + + // First reset internal energy before call to compute_temp + MultiFab reset_e_src(grids, dmap, 1, NUM_GROW); + reset_e_src.setVal(0.0); + reset_internal_energy(S_new,D_new,reset_e_src); + // Re-compute temperature after all the other updates. - compute_new_temp(); + compute_new_temp(S_new,D_new); } #endif } @@ -2180,7 +2204,7 @@ Nyx::network_init () #ifndef NO_HYDRO void -Nyx::reset_internal_energy (MultiFab& S_new, MultiFab& D_new) +Nyx::reset_internal_energy (MultiFab& S_new, MultiFab& D_new, MultiFab& reset_e_src) { BL_PROFILE("Nyx::reset_internal_energy()"); // Synchronize (rho e) and (rho E) so they are consistent with each other @@ -2204,6 +2228,7 @@ Nyx::reset_internal_energy (MultiFab& S_new, MultiFab& D_new) Real se = 0; reset_internal_e (BL_TO_FORTRAN(S_new[mfi]), BL_TO_FORTRAN(D_new[mfi]), + BL_TO_FORTRAN(reset_e_src[mfi]), bx.loVect(), bx.hiVect(), &print_fortran_warnings, &a, &s, &se); sum_energy_added += s; @@ -2237,20 +2262,15 @@ Nyx::reset_internal_energy (MultiFab& S_new, MultiFab& D_new) #ifndef NO_HYDRO void -Nyx::compute_new_temp () +Nyx::compute_new_temp (MultiFab& S_new, MultiFab& D_new) { BL_PROFILE("Nyx::compute_new_temp()"); - MultiFab& S_new = get_new_data(State_Type); - MultiFab& D_new = get_new_data(DiagEOS_Type); - - Real cur_time = state[State_Type].curTime(); - - reset_internal_energy(S_new,D_new); - Real a = get_comoving_a(cur_time); + Real cur_time = state[State_Type].curTime(); + Real a = get_comoving_a(cur_time); #ifdef HEATCOOL - if (heat_cool_type == 1 || heat_cool_type == 3 || heat_cool_type == 5 || heat_cool_type == 7) { + if (heat_cool_type == 3 || heat_cool_type == 5 || heat_cool_type == 7) { const Real z = 1.0/a - 1.0; fort_interp_to_this_z(&z); } diff --git a/Source/Nyx_F.H b/Source/Nyx_F.H index 47035359..9cbd525f 100644 --- a/Source/Nyx_F.H +++ b/Source/Nyx_F.H @@ -115,6 +115,7 @@ extern "C" void reset_internal_e (BL_FORT_FAB_ARG(S_new), BL_FORT_FAB_ARG(D_new), + BL_FORT_FAB_ARG(reset_e_src), const int lo[], const int hi[], const int* print_fortran_warnings, amrex::Real* comoving_a, amrex::Real* sum_energy_added, @@ -194,7 +195,6 @@ extern "C" const BL_FORT_FAB_ARG(zflux)), const amrex::Real* cflLoc, const amrex::Real* a_old, const amrex::Real* a_new, - const amrex::Real* e_added, const amrex::Real* ke_added, const int* print_fortran_warnings, const int* do_gas); @@ -217,7 +217,7 @@ extern "C" const BL_FORT_FAB_ARG(grad_phi_cc), const BL_FORT_FAB_ARG(S_old), BL_FORT_FAB_ARG(S_new), const amrex::Real* a_old, const amrex::Real* a_new, - const amrex::Real* dt, const amrex::Real* e_added, const amrex::Real* ke_added); + const amrex::Real* dt); void fort_syncgsrc (const int lo[], const int hi[], diff --git a/Source/Nyx_advance.cpp b/Source/Nyx_advance.cpp index 449aa479..ac5684d9 100644 --- a/Source/Nyx_advance.cpp +++ b/Source/Nyx_advance.cpp @@ -286,12 +286,12 @@ Nyx::advance_hydro_plus_particles (Real time, for (int lev = level; lev <= finest_level_to_advance; lev++) { #ifdef SDC - if (sdc_split > 0) - { + if (sdc_split > 0) + { get_level(lev).sdc_hydro(time, dt, a_old, a_new); - } else { + } else { get_level(lev).strang_hydro(time, dt, a_old, a_new); - } + } #else get_level(lev).strang_hydro(time, dt, a_old, a_new); #endif @@ -365,6 +365,9 @@ Nyx::advance_hydro_plus_particles (Real time, { MultiFab& S_old = get_level(lev).get_old_data(State_Type); MultiFab& S_new = get_level(lev).get_new_data(State_Type); + MultiFab& D_new = get_level(lev).get_new_data(DiagEOS_Type); + MultiFab reset_e_src(grids, dmap, 1, NUM_GROW); + reset_e_src.setVal(0.0); const auto& ba = get_level(lev).get_new_data(State_Type).boxArray(); const auto& dm = get_level(lev).get_new_data(State_Type).DistributionMap(); @@ -377,12 +380,10 @@ Nyx::advance_hydro_plus_particles (Real time, get_level(lev).gravity->get_old_grav_vector(lev, grav_vec_old, time); get_level(lev).gravity->get_new_grav_vector(lev, grav_vec_new, cur_time); - Real e_added = 0; - Real ke_added = 0; const Real* dx = get_level(lev).Geom().CellSize(); #ifdef _OPENMP -#pragma omp parallel reduction(+:e_added,ke_added) +#pragma omp parallel #endif for (MFIter mfi(S_new,true); mfi.isValid(); ++mfi) { @@ -394,40 +395,12 @@ Nyx::advance_hydro_plus_particles (Real time, fort_correct_gsrc (bx.loVect(), bx.hiVect(), BL_TO_FORTRAN(grav_vec_old[mfi]), BL_TO_FORTRAN(grav_vec_new[mfi]), BL_TO_FORTRAN(S_old[mfi]), - BL_TO_FORTRAN(S_new[mfi]), &a_old, &a_new, &dt, &se, &ske); - - e_added += se; - ke_added += ske; - } - - if (verbose > 2) - { - Real added[2] = {e_added,ke_added}; - - const int IOProc = ParallelDescriptor::IOProcessorNumber(); - - ParallelDescriptor::ReduceRealSum(added, 2, IOProc); - - if (ParallelDescriptor::IOProcessor()) - { - const Real vol = D_TERM(dx[0],*dx[1],*dx[2]); - - e_added = vol*added[0]; - ke_added = vol*added[1]; - - const Real sum_added = e_added + ke_added; - - std::cout << "Grav. correct work at level " - << lev - << " is " - << e_added/sum_added*100 - << " % into (rho e) and " - << ke_added/sum_added*100 - << " % into (KE) " << '\n'; - } + BL_TO_FORTRAN(S_new[mfi]), &a_old, &a_new, &dt); } - get_level(lev).compute_new_temp(); + // First reset internal energy before call to compute_temp + get_level(lev).reset_internal_energy(S_new,D_new,reset_e_src); + get_level(lev).compute_new_temp(S_new,D_new); } // Must average down again after doing the gravity correction; @@ -476,7 +449,11 @@ Nyx::advance_hydro_plus_particles (Real time, { MultiFab& S_new = get_level(lev).get_new_data(State_Type); MultiFab& D_new = get_level(lev).get_new_data(DiagEOS_Type); - get_level(lev).reset_internal_energy(S_new,D_new); + MultiFab reset_e_src(grids, dmap, 1, NUM_GROW); + reset_e_src.setVal(0.0); + + get_level(lev).reset_internal_energy(S_new,D_new,reset_e_src); + } BL_PROFILE_REGION_STOP("R::Nyx::advance_hydro_plus_particles"); @@ -538,12 +515,12 @@ Nyx::advance_hydro (Real time, // Call the hydro advance itself BL_PROFILE_VAR("just_the_hydro", just_the_hydro); #ifdef SDC - if (sdc_split > 0) - { + if (sdc_split > 0) + { sdc_hydro(time, dt, a_old, a_new); - } else { + } else { strang_hydro(time, dt, a_old, a_new); - } + } #else strang_hydro(time, dt, a_old, a_new); #endif @@ -551,6 +528,8 @@ Nyx::advance_hydro (Real time, MultiFab& S_new = get_new_data(State_Type); MultiFab& D_new = get_new_data(DiagEOS_Type); + MultiFab reset_e_src(grids, dmap, 1, NUM_GROW); + reset_e_src.setVal(0.0); #ifdef GRAVITY if (verbose && ParallelDescriptor::IOProcessor()) @@ -576,61 +555,25 @@ Nyx::advance_hydro (Real time, MultiFab grav_vec_new(grids,dmap,BL_SPACEDIM,0); gravity->get_new_grav_vector(level,grav_vec_new,cur_time); - Real e_added = 0; - Real ke_added = 0; - const Real* dx = geom.CellSize(); - // Now do corrector part of source term update #ifdef _OPENMP -#pragma omp parallel reduction(+:e_added,ke_added) +#pragma omp parallel #endif for (MFIter mfi(S_new,true); mfi.isValid(); ++mfi) { const Box& bx = mfi.tilebox(); - Real se = 0; - Real ske = 0; - fort_correct_gsrc (bx.loVect(), bx.hiVect(), BL_TO_FORTRAN(grav_vec_old[mfi]), BL_TO_FORTRAN(grav_vec_new[mfi]), BL_TO_FORTRAN(S_old[mfi]), - BL_TO_FORTRAN(S_new[mfi]), &a_old, &a_new, &dt, &se, &ske); - - e_added += se; - ke_added += ske; - } - - if (verbose > 2) - { - Real added[2] = {e_added,ke_added}; - - const int IOProc = ParallelDescriptor::IOProcessorNumber(); - - ParallelDescriptor::ReduceRealSum(added, 2, IOProc); - - if (ParallelDescriptor::IOProcessor()) - { - const Real vol = D_TERM(dx[0],*dx[1],*dx[2]); - - e_added = vol*added[0]; - ke_added = vol*added[1]; - - const Real sum_added = e_added + ke_added; - - std::cout << "Grav. correct work at level " - << level - << " is " - << e_added/sum_added*100 - << " % into (rho e) and " - << ke_added/sum_added*100 - << " % into (KE) " << '\n'; - } + BL_TO_FORTRAN(S_new[mfi]), &a_old, &a_new, &dt); } - compute_new_temp(); #endif /*GRAVITY*/ - reset_internal_energy(S_new,D_new); + // First reset internal energy before call to compute_temp + reset_internal_energy(S_new,D_new,reset_e_src); + compute_new_temp(S_new,D_new); return dt; } diff --git a/Source/Nyx_nd.f90 b/Source/Nyx_nd.f90 index d1f07c2a..2244b56d 100644 --- a/Source/Nyx_nd.f90 +++ b/Source/Nyx_nd.f90 @@ -234,6 +234,9 @@ subroutine fort_set_method_params( & difmag = 0.1d0 grav_source_type = 1 + ! We may want to default to 3 when using SDC because then the gravity updates + ! do not change the internal energy (rho e) -- but this needs further investigation + ! grav_source_type = 3 comoving_type = 1 diff --git a/Source/SourceTerms/Nyx_grav_sources_3d.f90 b/Source/SourceTerms/Nyx_grav_sources_3d.f90 index 9dc3878c..5be0ec25 100644 --- a/Source/SourceTerms/Nyx_grav_sources_3d.f90 +++ b/Source/SourceTerms/Nyx_grav_sources_3d.f90 @@ -11,7 +11,7 @@ subroutine fort_correct_gsrc(lo,hi, & gnew,gnew_l1,gnew_l2,gnew_l3,gnew_h1,gnew_h2,gnew_h3, & uold,uold_l1,uold_l2,uold_l3,uold_h1,uold_h2,uold_h3, & unew,unew_l1,unew_l2,unew_l3,unew_h1,unew_h2,unew_h3, & - a_old,a_new,dt,e_added,ke_added) & + a_old,a_new,dt) & bind(C, name="fort_correct_gsrc") use amrex_fort_module, only : rt => amrex_real @@ -28,7 +28,7 @@ subroutine fort_correct_gsrc(lo,hi, & real(rt) gnew(gnew_l1:gnew_h1,gnew_l2:gnew_h2,gnew_l3:gnew_h3,3) real(rt) uold(uold_l1:uold_h1,uold_l2:uold_h2,uold_l3:uold_h3,NVAR) real(rt) unew(unew_l1:unew_h1,unew_l2:unew_h2,unew_l3:unew_h3,NVAR) - real(rt) a_old,a_new,dt,e_added,ke_added + real(rt) a_old,a_new,dt integer i,j,k real(rt) SrU_old, SrV_old, SrW_old @@ -109,16 +109,6 @@ subroutine fort_correct_gsrc(lo,hi, & call bl_error("Error:: Nyx_advection_3d.f90 :: bogus grav_source_type") end if - ! **** Start Diagnostics **** - ! This is the new (rho e) as stored in (rho E) after the gravitational work is added - new_ke = 0.5d0 * (unew(i,j,k,UMX)**2 + unew(i,j,k,UMY)**2 + unew(i,j,k,UMZ)**2) / & - unew(i,j,k,URHO) - new_rhoeint = unew(i,j,k,UEDEN) - new_ke - - e_added = e_added + (new_rhoeint - old_rhoeint) - ke_added = ke_added + (new_ke - old_ke ) - ! **** End Diagnostics **** - enddo enddo enddo diff --git a/Source/sdc_reactions.cpp b/Source/sdc_reactions.cpp index e016075b..abbb9a31 100644 --- a/Source/sdc_reactions.cpp +++ b/Source/sdc_reactions.cpp @@ -14,11 +14,12 @@ Nyx::sdc_reactions (MultiFab& S_old, MultiFab& S_new, MultiFab& D_new, const Real* dx = geom.CellSize(); - MultiFab reset_src(grids, dmap, 1, NUM_GROW); - reset_src.setVal(0.0); - compute_new_temp(S_new,D_new,reset_src); - // MultiFab reset_src(grids, dmap, 1, NUM_GROW); - // reset_src.setVal(0.0); + // First reset internal energy before call to compute_temp + MultiFab reset_e_src(grids, dmap, 1, NUM_GROW); + reset_e_src.setVal(0.0); + + reset_internal_energy(S_new,D_new,reset_e_src); + compute_new_temp (S_new,D_new); #ifndef FORCING { @@ -52,7 +53,7 @@ Nyx::sdc_reactions (MultiFab& S_old, MultiFab& S_new, MultiFab& D_new, BL_TO_FORTRAN(S_new[mfi]), BL_TO_FORTRAN(D_new[mfi]), BL_TO_FORTRAN(hydro_src[mfi]), - BL_TO_FORTRAN(reset_src[mfi]), + BL_TO_FORTRAN(reset_e_src[mfi]), BL_TO_FORTRAN(IR[mfi]), &a_old, &delta_time, &min_iter_grid, &max_iter_grid); @@ -79,7 +80,7 @@ Nyx::sdc_reactions (MultiFab& S_old, MultiFab& S_new, MultiFab& D_new, BL_TO_FORTRAN(S_new[mfi]), BL_TO_FORTRAN(D_new[mfi]), BL_TO_FORTRAN(hydro_src[mfi]), - BL_TO_FORTRAN(reset_src[mfi]), + BL_TO_FORTRAN(reset_e_src[mfi]), BL_TO_FORTRAN(IR[mfi]), &a_old, &delta_time, &min_iter_grid, &max_iter_grid); diff --git a/Source/strang_hydro.cpp b/Source/strang_hydro.cpp index 56925ada..a3efa5fb 100644 --- a/Source/strang_hydro.cpp +++ b/Source/strang_hydro.cpp @@ -80,24 +80,6 @@ Nyx::strang_hydro (Real time, MultiFab ext_src_old(grids, dmap, NUM_STATE, 3); ext_src_old.setVal(0); - if (add_ext_src && ParallelDescriptor::IOProcessor()) - { - if (strang_split) - { - std::cout << "Source terms are handled with strang splitting" << std::endl; - } else { - std::cout << "Source terms are handled with predictor/corrector" << std::endl; - } - } - - if (add_ext_src && !strang_split) - { -#ifndef NO_OLD_SRC - get_old_source(prev_time, dt, ext_src_old); -#endif //#ifndef NO_OLD_SRC - ext_src_old.FillBoundary(); - } - // Define the gravity vector so we can pass this to ca_umdrv. MultiFab grav_vector(grids, dmap, BL_SPACEDIM, 3); grav_vector.setVal(0.); @@ -116,9 +98,6 @@ Nyx::strang_hydro (Real time, BL_ASSERT(NUM_GROW == 4); - Real e_added = 0; - Real ke_added = 0; - // Create FAB for extended grid values (including boundaries) and fill. MultiFab S_old_tmp(S_old.boxArray(), S_old.DistributionMap(), NUM_STATE, NUM_GROW); FillPatch(*this, S_old_tmp, NUM_GROW, time, State_Type, 0, NUM_STATE); @@ -126,11 +105,10 @@ Nyx::strang_hydro (Real time, MultiFab D_old_tmp(D_old.boxArray(), D_old.DistributionMap(), D_old.nComp(), NUM_GROW); FillPatch(*this, D_old_tmp, NUM_GROW, time, DiagEOS_Type, 0, D_old.nComp()); - if (add_ext_src && strang_split) - strang_first_step(time,dt,S_old_tmp,D_old_tmp); + strang_first_step(time,dt,S_old_tmp,D_old_tmp); #ifdef _OPENMP -#pragma omp parallel reduction(max:courno) reduction(+:e_added,ke_added) +#pragma omp parallel reduction(max:courno) #endif { FArrayBox flux[BL_SPACEDIM], u_gdnv[BL_SPACEDIM]; @@ -173,14 +151,12 @@ Nyx::strang_hydro (Real time, BL_TO_FORTRAN(flux[0]), BL_TO_FORTRAN(flux[1]), BL_TO_FORTRAN(flux[2]), - &cflLoc, &a_old, &a_new, &se, &ske, &print_fortran_warnings, &do_grav); + &cflLoc, &a_old, &a_new, &print_fortran_warnings, &do_grav); for (int i = 0; i < BL_SPACEDIM; ++i) { fluxes[i][mfi].copy(flux[i], mfi.nodaltilebox(i)); } - e_added += se; - ke_added += ske; } // end of MFIter loop courno = std::max(courno, cflLoc); @@ -204,38 +180,6 @@ Nyx::strang_hydro (Real time, } } -#ifdef GRAVITY - if (verbose > 1) - { - Real added[2] = {e_added,ke_added}; - - const int IOProc = ParallelDescriptor::IOProcessorNumber(); - - ParallelDescriptor::ReduceRealSum(added, 2, IOProc); - - if (ParallelDescriptor::IOProcessor()) - { - const Real vol = D_TERM(dx[0],*dx[1],*dx[2]); - - e_added = vol*added[0]; - ke_added = vol*added[1]; - - const Real sum_added = std::abs(e_added) + std::abs(ke_added); - - if (sum_added > 0) - { - std::cout << "Gravitational work at level " - << level - << " is " - << std::abs( e_added)/sum_added*100 - << " % into (rho e) and " - << std::abs(ke_added)/sum_added*100 - << " % into (KE) " << '\n'; - } - } - } -#endif /*GRAVITY*/ - grav_vector.clear(); ParallelDescriptor::ReduceRealMax(courno); @@ -265,32 +209,13 @@ Nyx::strang_hydro (Real time, } #endif - if (add_ext_src && !strang_split) - { - get_old_source(prev_time, dt, ext_src_old); - // Must compute new temperature in case it is needed in the source term - // evaluation - compute_new_temp(); - - // Compute source at new time (no ghost cells needed) - MultiFab ext_src_new(grids, dmap, NUM_STATE, 0); - ext_src_new.setVal(0); - - get_new_source(prev_time, cur_time, dt, ext_src_new); - - time_center_source_terms(S_new, ext_src_old, ext_src_new, dt); - - compute_new_temp(); - } // end if (add_ext_src && !strang_split) - #ifndef NDEBUG if (S_new.contains_nan(Density, S_new.nComp(), 0)) amrex::Abort("S_new has NaNs before the second strang call"); #endif // This returns updated (rho e), (rho E), and Temperature - if (add_ext_src && strang_split) - strang_second_step(cur_time,dt,S_new,D_new); + strang_second_step(cur_time,dt,S_new,D_new); #ifndef NDEBUG if (S_new.contains_nan(Density, S_new.nComp(), 0)) diff --git a/Source/strang_reactions.cpp b/Source/strang_reactions.cpp index 432b2f85..d226705e 100644 --- a/Source/strang_reactions.cpp +++ b/Source/strang_reactions.cpp @@ -61,7 +61,10 @@ Nyx::strang_second_step (Real time, Real dt, MultiFab& S_new, MultiFab& D_new) const Real a = get_comoving_a(time-half_dt); const Real* dx = geom.CellSize(); - compute_new_temp(); + MultiFab reset_e_src(grids, dmap, 1, NUM_GROW); + reset_e_src.setVal(0.0); + reset_internal_energy(S_new,D_new,reset_e_src); + compute_new_temp (S_new,D_new); #ifndef FORCING { diff --git a/Source/sum_integrated_quantities.cpp b/Source/sum_integrated_quantities.cpp index bc17c456..056131ba 100644 --- a/Source/sum_integrated_quantities.cpp +++ b/Source/sum_integrated_quantities.cpp @@ -275,7 +275,15 @@ Nyx::compute_average_temperature (Real& average_temperature) for (int lev = 0; lev <= finest_level; lev++) { Nyx& nyx_lev = get_level(lev); - nyx_lev.compute_new_temp(); + MultiFab& S_new = nyx_lev.get_new_data(State_Type); + MultiFab& D_new = nyx_lev.get_new_data(DiagEOS_Type); + + MultiFab reset_e_src(grids, dmap, 1, NUM_GROW); + reset_e_src.setVal(0.0); + nyx_lev.reset_internal_energy(S_new,D_new,reset_e_src); + + nyx_lev.compute_new_temp (S_new,D_new); + average_temperature += nyx_lev.vol_weight_sum("Temp",time,true); } diff --git a/Source/write_info.cpp b/Source/write_info.cpp index 184c510b..340baad6 100644 --- a/Source/write_info.cpp +++ b/Source/write_info.cpp @@ -11,14 +11,20 @@ Nyx::write_info () if (ndatalogs > 0) { #ifndef NO_HYDRO + MultiFab& S_new = get_new_data(State_Type); MultiFab& D_new = get_new_data(DiagEOS_Type); Real max_t = 0; Real rho_T_avg=0.0, T_avg=0.0, Tinv_avg=0.0, T_meanrho=0.0; Real whim_mass_frac, whim_vol_frac, hh_mass_frac, hh_vol_frac, igm_mass_frac, igm_vol_frac; + if (do_hydro) { - compute_new_temp(); + // First reset internal energy before call to compute_temp + MultiFab reset_e_src(grids, dmap, 1, NUM_GROW); + reset_e_src.setVal(0.0); + reset_internal_energy(S_new,D_new,reset_e_src); + compute_new_temp (S_new,D_new); max_t = D_new.norm0(Temp_comp); compute_rho_temp(rho_T_avg, T_avg, Tinv_avg, T_meanrho); compute_gas_fractions(1.0e5, 120.0, whim_mass_frac, whim_vol_frac, From f2e24e89e9fc17e8cfdb932f0a25ae1ec3a977ab Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Wed, 2 May 2018 13:29:43 -0700 Subject: [PATCH 09/12] 1) Some cleanup 2) We need to use the boxArray and DistributionMap of the right level when constructing reset_e_src so we use S_new.boxArray and S_new.DistributionMap rather than grids and dmap. --- Source/Initialization/Nyx_initdata.cpp | 15 +- Source/Nyx.cpp | 178 +++++++++------------ Source/NyxParticles.cpp | 6 +- Source/Nyx_F.H | 7 +- Source/Nyx_advance.cpp | 6 +- Source/Src_3d/reset_internal_energy_3d.f90 | 18 +-- Source/sdc_reactions.cpp | 2 +- Source/strang_hydro.cpp | 4 + Source/strang_reactions.cpp | 4 +- Source/sum_integrated_quantities.cpp | 22 +-- Source/write_info.cpp | 2 +- 11 files changed, 109 insertions(+), 155 deletions(-) diff --git a/Source/Initialization/Nyx_initdata.cpp b/Source/Initialization/Nyx_initdata.cpp index 3f1056e9..106b4b82 100644 --- a/Source/Initialization/Nyx_initdata.cpp +++ b/Source/Initialization/Nyx_initdata.cpp @@ -143,13 +143,13 @@ Nyx::init_zhi () MultiFab& D_new = get_new_data(DiagEOS_Type); int nd = D_new.nComp(); - const BoxArray& ba = D_new.boxArray(); - const DistributionMapping& dmap = D_new.DistributionMap(); + const BoxArray& my_ba = D_new.boxArray(); + const DistributionMapping& my_dmap = D_new.DistributionMap(); - BL_ASSERT(ba.coarsenable(ratio)); - BoxArray coarse_ba = ba; + BL_ASSERT(my_ba.coarsenable(ratio)); + BoxArray coarse_ba = my_ba; coarse_ba.coarsen(ratio); - MultiFab zhi(coarse_ba, dmap, 1, 0); + MultiFab zhi(coarse_ba, my_dmap, 1, 0); MultiFab zhi_from_file; VisMF::Read(zhi_from_file, inhomo_zhi_file); @@ -228,7 +228,7 @@ Nyx::initData () if (inhomo_reion) init_zhi(); // First reset internal energy before call to compute_temp - MultiFab reset_e_src(grids, dmap, 1, NUM_GROW); + MultiFab reset_e_src(S_new.boxArray(), S_new.DistributionMap(), 1, NUM_GROW); reset_e_src.setVal(0.0); reset_internal_energy(S_new,D_new,reset_e_src); @@ -355,10 +355,7 @@ Nyx::initData () // Need to compute this in case we want to use overdensity for regridding. // if (level == 0) - { compute_average_density(); -// compute_level_averages(); - } #endif if (verbose && ParallelDescriptor::IOProcessor()) diff --git a/Source/Nyx.cpp b/Source/Nyx.cpp index 60ba2beb..0df9199b 100644 --- a/Source/Nyx.cpp +++ b/Source/Nyx.cpp @@ -122,7 +122,7 @@ Real Nyx::comoving_h; int Nyx::do_hydro = -1; int Nyx::add_ext_src = 0; int Nyx::heat_cool_type = 0; -int Nyx::strang_split = 0; +int Nyx::strang_split = 1; #ifdef SDC int Nyx::sdc_split = 0; #endif @@ -244,29 +244,29 @@ Nyx::read_params () done = true; // ? - ParmParse pp("nyx"); + ParmParse pp_nyx("nyx"); - pp.query("v", verbose); - pp.get("init_shrink", init_shrink); - pp.get("cfl", cfl); - pp.query("change_max", change_max); - pp.query("fixed_dt", fixed_dt); - pp.query("initial_dt", initial_dt); - pp.query("sum_interval", sum_interval); - pp.query("do_reflux", do_reflux); + pp_nyx.query("v", verbose); + pp_nyx.get("init_shrink", init_shrink); + pp_nyx.get("cfl", cfl); + pp_nyx.query("change_max", change_max); + pp_nyx.query("fixed_dt", fixed_dt); + pp_nyx.query("initial_dt", initial_dt); + pp_nyx.query("sum_interval", sum_interval); + pp_nyx.query("do_reflux", do_reflux); do_reflux = (do_reflux ? 1 : 0); - pp.get("dt_cutoff", dt_cutoff); + pp_nyx.get("dt_cutoff", dt_cutoff); - pp.query("dump_old", dump_old); + pp_nyx.query("dump_old", dump_old); - pp.query("small_dens", small_dens); - pp.query("small_temp", small_temp); - pp.query("gamma", gamma); + pp_nyx.query("small_dens", small_dens); + pp_nyx.query("small_temp", small_temp); + pp_nyx.query("gamma", gamma); - pp.query("strict_subcycling",strict_subcycling); + pp_nyx.query("strict_subcycling",strict_subcycling); #ifdef AMREX_USE_CVODE - pp.query("simd_width", simd_width); + pp_nyx.query("simd_width", simd_width); if (simd_width < 1) amrex::Abort("simd_width must be a positive integer"); set_simd_width(simd_width); @@ -277,8 +277,8 @@ Nyx::read_params () // Get boundary conditions Vector lo_bc(BL_SPACEDIM), hi_bc(BL_SPACEDIM); - pp.getarr("lo_bc", lo_bc, 0, BL_SPACEDIM); - pp.getarr("hi_bc", hi_bc, 0, BL_SPACEDIM); + pp_nyx.getarr("lo_bc", lo_bc, 0, BL_SPACEDIM); + pp_nyx.getarr("hi_bc", hi_bc, 0, BL_SPACEDIM); for (int i = 0; i < BL_SPACEDIM; i++) { phys_bc.setLo(i, lo_bc[i]); @@ -339,15 +339,15 @@ Nyx::read_params () } } - pp.get("comoving_OmB", comoving_OmB); - pp.get("comoving_OmM", comoving_OmM); - pp.get("comoving_h", comoving_h); + pp_nyx.get("comoving_OmB", comoving_OmB); + pp_nyx.get("comoving_OmM", comoving_OmM); + pp_nyx.get("comoving_h", comoving_h); fort_set_omb(comoving_OmB); fort_set_omm(comoving_OmM); fort_set_hubble(comoving_h); - pp.get("do_hydro", do_hydro); + pp_nyx.get("do_hydro", do_hydro); #ifdef NO_HYDRO if (do_hydro == 1) amrex::Error("Cant have do_hydro == 1 when NO_HYDRO is true"); @@ -359,10 +359,10 @@ Nyx::read_params () #endif #endif - pp.query("add_ext_src", add_ext_src); - pp.query("strang_split", strang_split); + pp_nyx.query("add_ext_src", add_ext_src); + pp_nyx.query("strang_split", strang_split); #ifdef SDC - pp.query("sdc_split", sdc_split); + pp_nyx.query("sdc_split", sdc_split); if (sdc_split == 1 && strang_split == 1) amrex::Error("Cant have strang_split == 1 and sdc_split == 1"); if (sdc_split == 0 && strang_split == 0) @@ -375,7 +375,7 @@ Nyx::read_params () #endif #ifdef FORCING - pp.get("do_forcing", do_forcing); + pp_nyx.get("do_forcing", do_forcing); #ifdef NO_HYDRO if (do_forcing == 1) amrex::Error("Cant have do_forcing == 1 when NO_HYDRO is true "); @@ -387,7 +387,7 @@ Nyx::read_params () amrex::Error("Nyx::you set do_forcing = 1 but forgot to set USE_FORCING = TRUE "); #endif - pp.query("heat_cool_type", heat_cool_type); + pp_nyx.query("heat_cool_type", heat_cool_type); if (heat_cool_type == 7) { amrex::Print() << "----- WARNING WARNING WARNING WARNING WARNING -----" << std::endl; @@ -397,21 +397,22 @@ Nyx::read_params () amrex::Print() << " " << std::endl; amrex::Print() << "----- WARNING WARNING WARNING WARNING WARNING -----" << std::endl; Vector n_cell(BL_SPACEDIM); - ParmParse pp("amr"); - pp.getarr("n_cell", n_cell, 0, BL_SPACEDIM); + + ParmParse pp_amr("amr"); + pp_amr.getarr("n_cell", n_cell, 0, BL_SPACEDIM); if (n_cell[0] % simd_width) { const std::string errmsg = "Currently the SIMD CVODE solver requires that n_cell[0] % simd_width = 0"; amrex::Abort(errmsg); } } - pp.query("use_exact_gravity", use_exact_gravity); + pp_nyx.query("use_exact_gravity", use_exact_gravity); - pp.query("inhomo_reion", inhomo_reion); + pp_nyx.query("inhomo_reion", inhomo_reion); if (inhomo_reion) { - pp.get("inhomo_zhi_file", inhomo_zhi_file); - pp.get("inhomo_grid", inhomo_grid); + pp_nyx.get("inhomo_zhi_file", inhomo_zhi_file); + pp_nyx.get("inhomo_grid", inhomo_grid); } #ifdef HEATCOOL @@ -453,23 +454,23 @@ Nyx::read_params () amrex::Error("Nyx::you set inhomo_reion > 0 but forgot to set USE_HEATCOOL = TRUE"); #endif - pp.query("allow_untagging", allow_untagging); - pp.query("use_const_species", use_const_species); - pp.query("normalize_species", normalize_species); - pp.query("ppm_type", ppm_type); - pp.query("ppm_reference", ppm_reference); - pp.query("ppm_flatten_before_integrals", ppm_flatten_before_integrals); - pp.query("use_flattening", use_flattening); - pp.query("use_colglaz", use_colglaz); - pp.query("version_2", version_2); - pp.query("corner_coupling", corner_coupling); + pp_nyx.query("allow_untagging", allow_untagging); + pp_nyx.query("use_const_species", use_const_species); + pp_nyx.query("normalize_species", normalize_species); + pp_nyx.query("ppm_type", ppm_type); + pp_nyx.query("ppm_reference", ppm_reference); + pp_nyx.query("ppm_flatten_before_integrals", ppm_flatten_before_integrals); + pp_nyx.query("use_flattening", use_flattening); + pp_nyx.query("use_colglaz", use_colglaz); + pp_nyx.query("version_2", version_2); + pp_nyx.query("corner_coupling", corner_coupling); if (do_hydro == 1) { if (do_hydro == 1 && use_const_species == 1) { - pp.get("h_species" , h_species); - pp.get("he_species", he_species); + pp_nyx.get("h_species" , h_species); + pp_nyx.get("he_species", he_species); fort_set_xhydrogen(h_species); if (ParallelDescriptor::IOProcessor()) { @@ -519,35 +520,35 @@ Nyx::read_params () read_init_params(); - pp.query("write_parameter_file",write_parameters_in_plotfile); - pp.query("print_fortran_warnings",print_fortran_warnings); + pp_nyx.query("write_parameter_file",write_parameters_in_plotfile); + pp_nyx.query("print_fortran_warnings",print_fortran_warnings); read_comoving_params(); - if (pp.contains("plot_z_values")) + if (pp_nyx.contains("plot_z_values")) { - int num_z_values = pp.countval("plot_z_values"); + int num_z_values = pp_nyx.countval("plot_z_values"); plot_z_values.resize(num_z_values); - pp.queryarr("plot_z_values",plot_z_values,0,num_z_values); + pp_nyx.queryarr("plot_z_values",plot_z_values,0,num_z_values); } - if (pp.contains("analysis_z_values")) + if (pp_nyx.contains("analysis_z_values")) { - int num_z_values = pp.countval("analysis_z_values"); + int num_z_values = pp_nyx.countval("analysis_z_values"); analysis_z_values.resize(num_z_values); - pp.queryarr("analysis_z_values",analysis_z_values,0,num_z_values); + pp_nyx.queryarr("analysis_z_values",analysis_z_values,0,num_z_values); } // How often do we want to write x,y,z 2-d slices of S_new - pp.query("slice_int", slice_int); - pp.query("slice_file", slice_file); - pp.query("slice_nfiles", slice_nfiles); + pp_nyx.query("slice_int", slice_int); + pp_nyx.query("slice_file", slice_file); + pp_nyx.query("slice_nfiles", slice_nfiles); - pp.query("gimlet_int", gimlet_int); + pp_nyx.query("gimlet_int", gimlet_int); #ifdef AGN - pp.query("mass_halo_min", mass_halo_min); - pp.query("mass_seed", mass_seed); + pp_nyx.query("mass_halo_min", mass_halo_min); + pp_nyx.query("mass_seed", mass_seed); #endif } @@ -1302,13 +1303,9 @@ Nyx::post_timestep (int iteration) if ((iteration < ncycle and level < finest_level) || level == 0) { for (int i = 0; i < theActiveParticles().size(); i++) - { - int ngrow = (level == 0) ? 0 : iteration; - theActiveParticles()[i]->Redistribute(level, theActiveParticles()[i]->finestLevel(), iteration); - } } #ifndef NO_HYDRO @@ -1459,7 +1456,7 @@ Nyx::post_timestep (int iteration) MultiFab& D_new = get_new_data(DiagEOS_Type); // First reset internal energy before call to compute_temp - MultiFab reset_e_src(grids, dmap, 1, NUM_GROW); + MultiFab reset_e_src(S_new.boxArray(), S_new.DistributionMap(), 1, NUM_GROW); reset_e_src.setVal(0.0); reset_internal_energy(S_new,D_new,reset_e_src); @@ -1601,8 +1598,6 @@ Nyx::postCoarseTimeStep (Real cumtime) AmrLevel::postCoarseTimeStep(cumtime); - const Real cur_time = state[State_Type].curTime(); - #ifdef AGN halo_find(parent->dtLevel(level)); #endif @@ -1754,10 +1749,12 @@ Nyx::post_regrid (int lbase, particle_redistribute(lbase, false); } +#ifdef GRAVITY + int which_level_being_advanced = parent->level_being_advanced(); -#ifdef GRAVITY bool do_grav_solve_here; + if (which_level_being_advanced >= 0) { do_grav_solve_here = (level == which_level_being_advanced) && (lbase == which_level_being_advanced); @@ -2211,51 +2208,19 @@ Nyx::reset_internal_energy (MultiFab& S_new, MultiFab& D_new, MultiFab& reset_e_ const Real cur_time = state[State_Type].curTime(); Real a = get_comoving_a(cur_time); - const Real* dx = geom.CellSize(); - const Real vol = D_TERM(dx[0],*dx[1],*dx[2]); - - Real sum_energy_added = 0; - Real sum_energy_total = 0; #ifdef _OPENMP -#pragma omp parallel reduction(+:sum_energy_added,sum_energy_total) +#pragma omp parallel #endif for (MFIter mfi(S_new,true); mfi.isValid(); ++mfi) { const Box& bx = mfi.tilebox(); - Real s = 0; - Real se = 0; reset_internal_e - (BL_TO_FORTRAN(S_new[mfi]), BL_TO_FORTRAN(D_new[mfi]), + (bx.loVect(), bx.hiVect(), + BL_TO_FORTRAN(S_new[mfi]), BL_TO_FORTRAN(D_new[mfi]), BL_TO_FORTRAN(reset_e_src[mfi]), - bx.loVect(), bx.hiVect(), - &print_fortran_warnings, &a, &s, &se); - sum_energy_added += s; - sum_energy_total += se; - } - - if (verbose > 1) - { - Real sums[2] = {sum_energy_added,sum_energy_total}; - - const int IOProc = ParallelDescriptor::IOProcessorNumber(); - - ParallelDescriptor::ReduceRealSum(sums,2,IOProc); - - if (ParallelDescriptor::IOProcessor()) - { - sum_energy_added = vol*sums[0]; - sum_energy_total = vol*sums[1]; - - if (sum_energy_added > (1.e-12)*sum_energy_total) - { - std::cout << "Adding to (rho E) " - << sum_energy_added - << " out of total (rho E) " - << sum_energy_total << '\n'; - } - } + &print_fortran_warnings, &a); } } #endif @@ -2269,8 +2234,9 @@ Nyx::compute_new_temp (MultiFab& S_new, MultiFab& D_new) Real cur_time = state[State_Type].curTime(); Real a = get_comoving_a(cur_time); -#ifdef HEATCOOL - if (heat_cool_type == 3 || heat_cool_type == 5 || heat_cool_type == 7) { +#ifdef HEATCOOL + if (heat_cool_type == 3 || heat_cool_type == 5 || heat_cool_type == 7) + { const Real z = 1.0/a - 1.0; fort_interp_to_this_z(&z); } diff --git a/Source/NyxParticles.cpp b/Source/NyxParticles.cpp index f658391d..7c4f4db5 100644 --- a/Source/NyxParticles.cpp +++ b/Source/NyxParticles.cpp @@ -907,16 +907,16 @@ Nyx::particle_est_time_step (Real& est_dt) #endif void -Nyx::particle_redistribute (int lbase, bool init) +Nyx::particle_redistribute (int lbase, bool my_init) { BL_PROFILE("Nyx::particle_redistribute()"); if (DMPC) { // - // If we are calling with init = true, then we want to force the redistribute + // If we are calling with my_init = true, then we want to force the redistribute // without checking whether the grids have changed. // - if (init) + if (my_init) { DMPC->Redistribute(lbase); return; diff --git a/Source/Nyx_F.H b/Source/Nyx_F.H index 9cbd525f..d6df42a7 100644 --- a/Source/Nyx_F.H +++ b/Source/Nyx_F.H @@ -113,13 +113,12 @@ extern "C" const int* print_fortran_warnings); void reset_internal_e - (BL_FORT_FAB_ARG(S_new), + (const int lo[], const int hi[], + BL_FORT_FAB_ARG(S_new), BL_FORT_FAB_ARG(D_new), BL_FORT_FAB_ARG(reset_e_src), - const int lo[], const int hi[], const int* print_fortran_warnings, - amrex::Real* comoving_a, amrex::Real* sum_energy_added, - amrex::Real* sum_energy_total); + amrex::Real* comoving_a); void hypfill (BL_FORT_FAB_ARG(state), diff --git a/Source/Nyx_advance.cpp b/Source/Nyx_advance.cpp index ac5684d9..56fec9f8 100644 --- a/Source/Nyx_advance.cpp +++ b/Source/Nyx_advance.cpp @@ -366,7 +366,7 @@ Nyx::advance_hydro_plus_particles (Real time, MultiFab& S_old = get_level(lev).get_old_data(State_Type); MultiFab& S_new = get_level(lev).get_new_data(State_Type); MultiFab& D_new = get_level(lev).get_new_data(DiagEOS_Type); - MultiFab reset_e_src(grids, dmap, 1, NUM_GROW); + MultiFab reset_e_src(S_new.boxArray(), S_new.DistributionMap(), 1, NUM_GROW); reset_e_src.setVal(0.0); const auto& ba = get_level(lev).get_new_data(State_Type).boxArray(); @@ -449,7 +449,7 @@ Nyx::advance_hydro_plus_particles (Real time, { MultiFab& S_new = get_level(lev).get_new_data(State_Type); MultiFab& D_new = get_level(lev).get_new_data(DiagEOS_Type); - MultiFab reset_e_src(grids, dmap, 1, NUM_GROW); + MultiFab reset_e_src(S_new.boxArray(), S_new.DistributionMap(), 1, NUM_GROW); reset_e_src.setVal(0.0); get_level(lev).reset_internal_energy(S_new,D_new,reset_e_src); @@ -528,7 +528,7 @@ Nyx::advance_hydro (Real time, MultiFab& S_new = get_new_data(State_Type); MultiFab& D_new = get_new_data(DiagEOS_Type); - MultiFab reset_e_src(grids, dmap, 1, NUM_GROW); + MultiFab reset_e_src(S_new.boxArray(), S_new.DistributionMap(), 1, NUM_GROW); reset_e_src.setVal(0.0); #ifdef GRAVITY diff --git a/Source/Src_3d/reset_internal_energy_3d.f90 b/Source/Src_3d/reset_internal_energy_3d.f90 index 5218057b..4453c43a 100644 --- a/Source/Src_3d/reset_internal_energy_3d.f90 +++ b/Source/Src_3d/reset_internal_energy_3d.f90 @@ -1,7 +1,7 @@ - subroutine reset_internal_e(u,u_l1,u_l2,u_l3,u_h1,u_h2,u_h3, & - d,d_l1,d_l2,d_l3,d_h1,d_h2,d_h3,lo,hi, & - print_fortran_warnings,& - comoving_a,sum_energy_added,sum_energy_total) & + subroutine reset_internal_e(lo,hi,& + u,u_l1,u_l2,u_l3,u_h1,u_h2,u_h3, & + d,d_l1,d_l2,d_l3,d_h1,d_h2,d_h3,& + print_fortran_warnings, comoving_a) & bind(C, name="reset_internal_e") use amrex_fort_module, only : rt => amrex_real @@ -20,8 +20,6 @@ subroutine reset_internal_e(u,u_l1,u_l2,u_l3,u_h1,u_h2,u_h3, & real(rt) :: d(d_l1:d_h1,d_l2:d_h2,d_l3:d_h3,NDIAG) real(rt), intent(in ) :: comoving_a - real(rt), intent(inout) :: sum_energy_added - real(rt), intent(inout) :: sum_energy_total ! Local variables integer :: i,j,k @@ -49,9 +47,6 @@ subroutine reset_internal_e(u,u_l1,u_l2,u_l3,u_h1,u_h2,u_h3, & ! If (e from E) < 0 or (e from E) < .0001*E but (e from e) > 0. else if (u(i,j,k,UEINT) .gt. 0.d0) then - ! Keep track of how much energy we are adding to (rho E) - sum_energy_added = sum_energy_added + (u(i,j,k,UEINT) + ke - u(i,j,k,UEDEN)) - u(i,j,k,UEDEN) = u(i,j,k,UEINT) + ke ! If not resetting and little e is negative ... @@ -70,15 +65,10 @@ subroutine reset_internal_e(u,u_l1,u_l2,u_l3,u_h1,u_h2,u_h3, & u(i,j,k,UEINT) = u(i,j,k,URHO) * eint_new - ! Keep track of how much energy we are adding to (rho E) - sum_energy_added = sum_energy_added + (u(i,j,k,UEINT) + ke - u(i,j,k,UEDEN)) - u(i,j,k,UEDEN) = u(i,j,k,UEINT) + ke end if - sum_energy_total = sum_energy_total + u(i,j,k,UEDEN) - enddo enddo enddo diff --git a/Source/sdc_reactions.cpp b/Source/sdc_reactions.cpp index abbb9a31..f8212020 100644 --- a/Source/sdc_reactions.cpp +++ b/Source/sdc_reactions.cpp @@ -15,7 +15,7 @@ Nyx::sdc_reactions (MultiFab& S_old, MultiFab& S_new, MultiFab& D_new, const Real* dx = geom.CellSize(); // First reset internal energy before call to compute_temp - MultiFab reset_e_src(grids, dmap, 1, NUM_GROW); + MultiFab reset_e_src(S_new.boxArray(), S_new.DistributionMap(), 1, NUM_GROW); reset_e_src.setVal(0.0); reset_internal_energy(S_new,D_new,reset_e_src); diff --git a/Source/strang_hydro.cpp b/Source/strang_hydro.cpp index a3efa5fb..7220ac94 100644 --- a/Source/strang_hydro.cpp +++ b/Source/strang_hydro.cpp @@ -105,7 +105,9 @@ Nyx::strang_hydro (Real time, MultiFab D_old_tmp(D_old.boxArray(), D_old.DistributionMap(), D_old.nComp(), NUM_GROW); FillPatch(*this, D_old_tmp, NUM_GROW, time, DiagEOS_Type, 0, D_old.nComp()); +#ifdef HEATCOOL strang_first_step(time,dt,S_old_tmp,D_old_tmp); +#endif #ifdef _OPENMP #pragma omp parallel reduction(max:courno) @@ -214,8 +216,10 @@ Nyx::strang_hydro (Real time, amrex::Abort("S_new has NaNs before the second strang call"); #endif +#ifdef HEATCOOL // This returns updated (rho e), (rho E), and Temperature strang_second_step(cur_time,dt,S_new,D_new); +#endif #ifndef NDEBUG if (S_new.contains_nan(Density, S_new.nComp(), 0)) diff --git a/Source/strang_reactions.cpp b/Source/strang_reactions.cpp index d226705e..3c25c66c 100644 --- a/Source/strang_reactions.cpp +++ b/Source/strang_reactions.cpp @@ -12,7 +12,6 @@ Nyx::strang_first_step (Real time, Real dt, MultiFab& S_old, MultiFab& D_old) Real half_dt = 0.5*dt; const Real a = get_comoving_a(time); - const Real* dx = geom.CellSize(); #ifndef FORCING { @@ -59,9 +58,8 @@ Nyx::strang_second_step (Real time, Real dt, MultiFab& S_new, MultiFab& D_new) // Set a at the half of the time step in the second strang const Real a = get_comoving_a(time-half_dt); - const Real* dx = geom.CellSize(); - MultiFab reset_e_src(grids, dmap, 1, NUM_GROW); + MultiFab reset_e_src(S_new.boxArray(), S_new.DistributionMap(), 1, NUM_GROW); reset_e_src.setVal(0.0); reset_internal_energy(S_new,D_new,reset_e_src); compute_new_temp (S_new,D_new); diff --git a/Source/sum_integrated_quantities.cpp b/Source/sum_integrated_quantities.cpp index 056131ba..93d0fb76 100644 --- a/Source/sum_integrated_quantities.cpp +++ b/Source/sum_integrated_quantities.cpp @@ -178,7 +178,7 @@ Nyx::compute_average_density () { int finest_level = parent->finestLevel(); Real time = state[State_Type].curTime(); - const Geometry& geom = parent->Geom(0); + const Geometry& crse_geom = parent->Geom(0); // This is a static in the Nyx class average_gas_density = 0; @@ -235,9 +235,9 @@ Nyx::compute_average_density () // Divide by physical volume of domain. if (do_hydro == 1) - average_gas_density /= geom.ProbSize(); - average_dm_density /= geom.ProbSize(); - average_neutr_density /= geom.ProbSize(); + average_gas_density /= crse_geom.ProbSize(); + average_dm_density /= crse_geom.ProbSize(); + average_neutr_density /= crse_geom.ProbSize(); // Define the total density = gas density + dark matter density if (do_hydro == 1) @@ -268,7 +268,7 @@ Nyx::compute_average_temperature (Real& average_temperature) { int finest_level = parent->finestLevel(); Real time = state[State_Type].curTime(); - const Geometry& geom = parent->Geom(0); + const Geometry& crse_geom = parent->Geom(0); // Add up the temperature -- this is just volume-weighted, not mass-weighted average_temperature = 0; @@ -278,17 +278,17 @@ Nyx::compute_average_temperature (Real& average_temperature) MultiFab& S_new = nyx_lev.get_new_data(State_Type); MultiFab& D_new = nyx_lev.get_new_data(DiagEOS_Type); - MultiFab reset_e_src(grids, dmap, 1, NUM_GROW); + MultiFab reset_e_src(S_new.boxArray(), S_new.DistributionMap(), 1, NUM_GROW); reset_e_src.setVal(0.0); - nyx_lev.reset_internal_energy(S_new,D_new,reset_e_src); + nyx_lev.reset_internal_energy(S_new,D_new,reset_e_src); nyx_lev.compute_new_temp (S_new,D_new); average_temperature += nyx_lev.vol_weight_sum("Temp",time,true); } // Divide by physical volume of domain. - average_temperature = average_temperature / geom.ProbSize(); + average_temperature = average_temperature / crse_geom.ProbSize(); if (verbose > 0 && ParallelDescriptor::IOProcessor()) { std::cout << "Average temperature " << average_temperature << '\n'; @@ -312,8 +312,8 @@ Nyx::compute_average_species (int nspec, } else { - Real time = state[State_Type].curTime(); - const Geometry& geom = parent->Geom(0); + Real time = state[State_Type].curTime(); + const Geometry& crse_geom = parent->Geom(0); // // Get the species names from the network model. // @@ -361,7 +361,7 @@ Nyx::compute_average_species (int nspec, delete [] name; // Divide by physical volume of domain. - average_species[i] = average_species[i] / geom.ProbSize(); + average_species[i] = average_species[i] / crse_geom.ProbSize(); if (verbose > 0 && ParallelDescriptor::IOProcessor()) std::cout << "Average species " << i << ": " << average_species[i] << '\n'; diff --git a/Source/write_info.cpp b/Source/write_info.cpp index 340baad6..c449206d 100644 --- a/Source/write_info.cpp +++ b/Source/write_info.cpp @@ -21,7 +21,7 @@ Nyx::write_info () if (do_hydro) { // First reset internal energy before call to compute_temp - MultiFab reset_e_src(grids, dmap, 1, NUM_GROW); + MultiFab reset_e_src(S_new.boxArray(), S_new.DistributionMap(), 1, NUM_GROW); reset_e_src.setVal(0.0); reset_internal_energy(S_new,D_new,reset_e_src); compute_new_temp (S_new,D_new); From 186c09019584e9a2bd53bcb97b1cb5f6299c86fe Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Wed, 2 May 2018 13:44:27 -0700 Subject: [PATCH 10/12] Oops -- missed one of the pp --> pp_nyx --- Source/Nyx.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Source/Nyx.cpp b/Source/Nyx.cpp index 0df9199b..5ee007b0 100644 --- a/Source/Nyx.cpp +++ b/Source/Nyx.cpp @@ -513,7 +513,7 @@ Nyx::read_params () if (do_hydro == 0) do_reflux = 0; #ifdef GRAVITY - pp.get("do_grav", do_grav); + pp_nyx.get("do_grav", do_grav); #endif read_particle_params(); From 43aaa244f7e924417c215bbd0ffe1450d90bbc81 Mon Sep 17 00:00:00 2001 From: Zarija Lukic Date: Mon, 7 May 2018 13:28:55 -0700 Subject: [PATCH 11/12] updating setups --- Exec/AMR-density/GNUmakefile | 4 ++-- Exec/LyA/inputs | 6 ++---- 2 files changed, 4 insertions(+), 6 deletions(-) diff --git a/Exec/AMR-density/GNUmakefile b/Exec/AMR-density/GNUmakefile index 53a57329..2fbfecbb 100644 --- a/Exec/AMR-density/GNUmakefile +++ b/Exec/AMR-density/GNUmakefile @@ -9,8 +9,8 @@ TOP = ../.. # compilation options COMP = intel # gnu -USE_MPI = FALSE -USE_OMP = FALSE +USE_MPI = TRUE +USE_OMP = TRUE PROFILE = TRUE TRACE_PROFILE = FALSE diff --git a/Exec/LyA/inputs b/Exec/LyA/inputs index c3b653db..bd18bcd7 100644 --- a/Exec/LyA/inputs +++ b/Exec/LyA/inputs @@ -1,7 +1,5 @@ # ------------------ INPUTS TO MAIN PROGRAM ------------------- max_step = 10000000 -max_step = 0 -max_step = 4 amr.mffile_nstreams = 4 amr.precreateDirectories = 1 @@ -54,8 +52,8 @@ amr.data_log = runlog gravity.gravity_type = PoissonGrav gravity.no_sync = 1 gravity.no_composite = 1 -gravity.solve_with_cpp = 0 -gravity.solve_with_hpgmg = 1 +gravity.solve_with_mlmg = 1 +gravity.solve_with_hpgmg = 0 mg.bottom_solver = 4 From 6dfe7172f1adff63467e72b94f84ce01e597ff7d Mon Sep 17 00:00:00 2001 From: Zarija Lukic Date: Mon, 7 May 2018 13:29:49 -0700 Subject: [PATCH 12/12] license info in README --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index 462468df..c518128b 100644 --- a/README.md +++ b/README.md @@ -89,7 +89,7 @@ be run either "in situ" or "in-transit", or with a combination of both. ## License -See the [license.txt](license.txt) file for details. +Nyx is released under the LBL's modified BSD license, see the [license.txt](license.txt) file for details. ## Contact