From 4a4c136ecf04979b14d79edec2b1122124196ed6 Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Tue, 29 Oct 2024 15:31:46 -0700 Subject: [PATCH] Generalize fillpatching (#1916) * generalize fp, etc * fix re-used variable --- CMake/BuildERFExe.cmake | 1 + .../ERF_BoundaryConditions_basestate.cpp | 126 +++++++++- .../ERF_FillCoarsePatch.cpp | 17 +- .../ERF_FillIntermediatePatch.cpp | 168 ++++++++++++-- Source/BoundaryConditions/ERF_FillPatch.cpp | 215 +++++++++++++----- Source/BoundaryConditions/ERF_PhysBCFunct.H | 5 +- Source/BoundaryConditions/ERF_PhysBCFunct.cpp | 11 + Source/ERF.H | 18 ++ Source/ERF.cpp | 140 +++--------- Source/ERF_make_new_arrays.cpp | 10 +- Source/ERF_make_new_level.cpp | 55 ++--- Source/IO/ERF_Checkpoint.cpp | 64 +++--- Source/IO/ERF_NCColumnFile.cpp | 14 +- Source/IO/ERF_Plotfile.cpp | 18 +- Source/Initialization/ERF_init1d.cpp | 93 +------- Source/TimeIntegration/ERF_Advance.cpp | 24 +- Source/Utils/ERF_AverageDown.cpp | 162 +++++++++++++ Source/Utils/Make.package | 1 + 18 files changed, 780 insertions(+), 362 deletions(-) create mode 100644 Source/Utils/ERF_AverageDown.cpp diff --git a/CMake/BuildERFExe.cmake b/CMake/BuildERFExe.cmake index 6c5d6195d..78d438d47 100644 --- a/CMake/BuildERFExe.cmake +++ b/CMake/BuildERFExe.cmake @@ -187,6 +187,7 @@ function(build_erf_lib erf_lib_name) ${SRC_DIR}/TimeIntegration/ERF_fast_rhs_N.cpp ${SRC_DIR}/TimeIntegration/ERF_fast_rhs_T.cpp ${SRC_DIR}/TimeIntegration/ERF_fast_rhs_MT.cpp + ${SRC_DIR}/Utils/ERF_AverageDown.cpp ${SRC_DIR}/Utils/ERF_ChopGrids.cpp ${SRC_DIR}/Utils/ERF_MomentumToVelocity.cpp ${SRC_DIR}/Utils/ERF_TerrainMetrics.cpp diff --git a/Source/BoundaryConditions/ERF_BoundaryConditions_basestate.cpp b/Source/BoundaryConditions/ERF_BoundaryConditions_basestate.cpp index fbd64a53c..d2101ad4c 100644 --- a/Source/BoundaryConditions/ERF_BoundaryConditions_basestate.cpp +++ b/Source/BoundaryConditions/ERF_BoundaryConditions_basestate.cpp @@ -1,5 +1,6 @@ #include "AMReX_PhysBCFunct.H" #include +#include using namespace amrex; @@ -15,9 +16,14 @@ void ERFPhysBCFunct_base::impose_lateral_basestate_bcs (const Array4& dest int ncomp, const IntVect& nghost) { BL_PROFILE_VAR("impose_lateral_base_bcs()",impose_lateral_base_bcs); + // + // Note that the "bx" that comes in here has already been grown in the lateral directions + // but not in the vertical + // const int* bxlo = bx.loVect(); const int* bxhi = bx.hiVect(); + const int* dlo = domain.loVect(); const int* dhi = domain.hiVect(); @@ -63,6 +69,84 @@ void ERFPhysBCFunct_base::impose_lateral_basestate_bcs (const Array4& dest // Populate ghost cells on lo-x and hi-x domain boundaries Box bx_xlo(bx); bx_xlo.setBig (0,dom_lo.x-nghost[0]); Box bx_xhi(bx); bx_xhi.setSmall(0,dom_hi.x+nghost[0]); + + ParallelFor( + bx_xlo, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) + { + int dest_comp = n; + int l_bc_type = bc_ptr[n].lo(0); + int iflip = dom_lo.x - 1 - i; + if (l_bc_type == ERFBCType::foextrap) { + dest_arr(i,j,k,dest_comp) = dest_arr(dom_lo.x,j,k,dest_comp); + } else if (l_bc_type == ERFBCType::open) { + dest_arr(i,j,k,dest_comp) = dest_arr(dom_lo.x,j,k,dest_comp); + } else if (l_bc_type == ERFBCType::reflect_even) { + dest_arr(i,j,k,dest_comp) = dest_arr(iflip,j,k,dest_comp); + } + }, + bx_xhi, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) + { + int dest_comp = n; + int h_bc_type = bc_ptr[n].hi(0); + int iflip = 2*dom_hi.x + 1 - i; + if (h_bc_type == ERFBCType::foextrap) { + dest_arr(i,j,k,dest_comp) = dest_arr(dom_hi.x,j,k,dest_comp); + } else if (h_bc_type == ERFBCType::open) { + dest_arr(i,j,k,dest_comp) = dest_arr(dom_hi.x,j,k,dest_comp); + } else if (h_bc_type == ERFBCType::reflect_even) { + dest_arr(i,j,k,dest_comp) = dest_arr(iflip,j,k,dest_comp); + } + } + ); + } + + if (!is_periodic_in_y) + { + // Populate ghost cells on lo-y and hi-y domain boundaries + Box bx_ylo(bx); bx_ylo.setBig (1,dom_lo.y-nghost[1]); + Box bx_yhi(bx); bx_yhi.setSmall(1,dom_hi.y+nghost[1]); + if (bx_ylo.smallEnd(2) != domain.smallEnd(2)) bx_ylo.growLo(2,nghost[2]); + if (bx_ylo.bigEnd(2) != domain.bigEnd(2)) bx_ylo.growHi(2,nghost[2]); + if (bx_yhi.smallEnd(2) != domain.smallEnd(2)) bx_yhi.growLo(2,nghost[2]); + if (bx_yhi.bigEnd(2) != domain.bigEnd(2)) bx_yhi.growHi(2,nghost[2]); + ParallelFor( + bx_ylo, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) + { + int dest_comp = n; + int l_bc_type = bc_ptr[n].lo(1); + int jflip = dom_lo.y - 1 - j; + if (l_bc_type == ERFBCType::foextrap) { + dest_arr(i,j,k,dest_comp) = dest_arr(i,dom_lo.y,k,dest_comp); + } else if (l_bc_type == ERFBCType::open) { + dest_arr(i,j,k,dest_comp) = dest_arr(i,dom_lo.y,k,dest_comp); + } else if (l_bc_type == ERFBCType::reflect_even) { + dest_arr(i,j,k,dest_comp) = dest_arr(i,jflip,k,dest_comp); + } + + }, + bx_yhi, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) + { + int dest_comp = n; + int h_bc_type = bc_ptr[n].hi(1); + int jflip = 2*dom_hi.y + 1 - j; + if (h_bc_type == ERFBCType::foextrap) { + dest_arr(i,j,k,dest_comp) = dest_arr(i,dom_hi.y,k,dest_comp); + } else if (h_bc_type == ERFBCType::open) { + dest_arr(i,j,k,dest_comp) = dest_arr(i,dom_hi.y,k,dest_comp); + } else if (h_bc_type == ERFBCType::reflect_even) { + dest_arr(i,j,k,dest_comp) = dest_arr(i,jflip,k,dest_comp); + } + } + ); + } + + // Next do ghost cells in x-direction but not reaching out in y + // The corners we miss here will be covered in the y-loop below or by periodicity + if (!is_periodic_in_x) + { + // Populate ghost cells on lo-x and hi-x domain boundaries + Box bx_xlo(bx); bx_xlo.setBig (0,dom_lo.x-1); + Box bx_xhi(bx); bx_xhi.setSmall(0,dom_hi.x+1); if (bx_xlo.smallEnd(2) != domain.smallEnd(2)) bx_xlo.growLo(2,nghost[2]); if (bx_xlo.bigEnd(2) != domain.bigEnd(2)) bx_xlo.growHi(2,nghost[2]); if (bx_xhi.smallEnd(2) != domain.smallEnd(2)) bx_xhi.growLo(2,nghost[2]); @@ -79,6 +163,11 @@ void ERFPhysBCFunct_base::impose_lateral_basestate_bcs (const Array4& dest dest_arr(i,j,k,dest_comp) = dest_arr(dom_lo.x,j,k,dest_comp); } else if (l_bc_type == ERFBCType::reflect_even) { dest_arr(i,j,k,dest_comp) = dest_arr(iflip,j,k,dest_comp); + } else if (l_bc_type == ERFBCType::reflect_odd) { + dest_arr(i,j,k,dest_comp) = -dest_arr(iflip,j,k,dest_comp); + } else if (l_bc_type == ERFBCType::hoextrapcc) { + Real delta_i = (dom_lo.x - i); + dest_arr(i,j,k,dest_comp) = (1.0 + delta_i)*dest_arr(dom_lo.x,j,k,dest_comp) - delta_i*dest_arr(dom_lo.x+1,j,k,dest_comp) ; } }, bx_xhi, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) @@ -92,6 +181,11 @@ void ERFPhysBCFunct_base::impose_lateral_basestate_bcs (const Array4& dest dest_arr(i,j,k,dest_comp) = dest_arr(dom_hi.x,j,k,dest_comp); } else if (h_bc_type == ERFBCType::reflect_even) { dest_arr(i,j,k,dest_comp) = dest_arr(iflip,j,k,dest_comp); + } else if (h_bc_type == ERFBCType::reflect_odd) { + dest_arr(i,j,k,dest_comp) = -dest_arr(iflip,j,k,dest_comp); + } else if (h_bc_type == ERFBCType::hoextrapcc) { + Real delta_i = (i - dom_hi.x); + dest_arr(i,j,k,dest_comp) = (1.0 + delta_i)*dest_arr(dom_hi.x,j,k,dest_comp) - delta_i*dest_arr(dom_hi.x-1,j,k,dest_comp) ; } } ); @@ -100,8 +194,8 @@ void ERFPhysBCFunct_base::impose_lateral_basestate_bcs (const Array4& dest if (!is_periodic_in_y) { // Populate ghost cells on lo-y and hi-y domain boundaries - Box bx_ylo(bx); bx_ylo.setBig (1,dom_lo.y-nghost[1]); - Box bx_yhi(bx); bx_yhi.setSmall(1,dom_hi.y+nghost[1]); + Box bx_ylo(bx); bx_ylo.setBig (1,dom_lo.y-1); + Box bx_yhi(bx); bx_yhi.setSmall(1,dom_hi.y+1); if (bx_ylo.smallEnd(2) != domain.smallEnd(2)) bx_ylo.growLo(2,nghost[2]); if (bx_ylo.bigEnd(2) != domain.bigEnd(2)) bx_ylo.growHi(2,nghost[2]); if (bx_yhi.smallEnd(2) != domain.smallEnd(2)) bx_yhi.growLo(2,nghost[2]); @@ -118,6 +212,11 @@ void ERFPhysBCFunct_base::impose_lateral_basestate_bcs (const Array4& dest dest_arr(i,j,k,dest_comp) = dest_arr(i,dom_lo.y,k,dest_comp); } else if (l_bc_type == ERFBCType::reflect_even) { dest_arr(i,j,k,dest_comp) = dest_arr(i,jflip,k,dest_comp); + } else if (l_bc_type == ERFBCType::reflect_odd) { + dest_arr(i,j,k,dest_comp) = -dest_arr(i,jflip,k,dest_comp); + } else if (l_bc_type == ERFBCType::hoextrapcc) { + Real delta_j = (dom_lo.y - j); + dest_arr(i,j,k,dest_comp) = (1.0 + delta_j)*dest_arr(i,dom_lo.y,k,dest_comp) - delta_j*dest_arr(i,dom_lo.y+1,k,dest_comp) ; } }, @@ -132,6 +231,11 @@ void ERFPhysBCFunct_base::impose_lateral_basestate_bcs (const Array4& dest dest_arr(i,j,k,dest_comp) = dest_arr(i,dom_hi.y,k,dest_comp); } else if (h_bc_type == ERFBCType::reflect_even) { dest_arr(i,j,k,dest_comp) = dest_arr(i,jflip,k,dest_comp); + } else if (h_bc_type == ERFBCType::reflect_odd) { + dest_arr(i,j,k,dest_comp) = -dest_arr(i,jflip,k,dest_comp); + } else if (h_bc_type == ERFBCType::hoextrapcc) { + Real delta_j = (j - dom_hi.y); + dest_arr(i,j,k,dest_comp) = (1.0 + delta_j)*dest_arr(i,dom_hi.y,k,dest_comp) - delta_j*dest_arr(i,dom_hi.y-1,k,dest_comp); } } ); @@ -147,8 +251,22 @@ void ERFPhysBCFunct_base::impose_vertical_basestate_bcs (const Array4& des const auto& dom_lo = lbound(domain); const auto& dom_hi = ubound(domain); + const Real hz = Real(0.5) * m_geom.CellSize(2); + + Box bx_zlo1(bx); bx_zlo1.setBig(2,dom_lo.z-1); if (bx_zlo1.ok()) bx_zlo1.setSmall(2,dom_lo.z-1); + ParallelFor( + bx_zlo1, [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + dest_arr(i,j,k,BaseState::r0_comp) = dest_arr(i,j,dom_lo.z,BaseState::r0_comp); + dest_arr(i,j,k,BaseState::p0_comp) = p_0 - + dest_arr(i,j,k,BaseState::r0_comp) * hz * 9.81; + dest_arr(i,j,k,BaseState::pi0_comp) = dest_arr(i,j,dom_lo.z,BaseState::pi0_comp); + dest_arr(i,j,k,BaseState::th0_comp) = dest_arr(i,j,dom_lo.z,BaseState::th0_comp); + } + ); + Box bx_zlo(bx); bx_zlo.setBig(2,dom_lo.z-2); - Box bx_zhi(bx); bx_zhi.setSmall(2,dom_hi.z+2); + Box bx_zhi(bx); bx_zhi.setSmall(2,dom_hi.z+1); ParallelFor( bx_zlo, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) { @@ -156,7 +274,7 @@ void ERFPhysBCFunct_base::impose_vertical_basestate_bcs (const Array4& des }, bx_zhi, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) { - dest_arr(i,j,k,n) = dest_arr(i,j,dom_hi.z+1,n); + dest_arr(i,j,k,n) = dest_arr(i,j,dom_hi.z,n); } ); } diff --git a/Source/BoundaryConditions/ERF_FillCoarsePatch.cpp b/Source/BoundaryConditions/ERF_FillCoarsePatch.cpp index 9815d050a..fb900ea7d 100644 --- a/Source/BoundaryConditions/ERF_FillCoarsePatch.cpp +++ b/Source/BoundaryConditions/ERF_FillCoarsePatch.cpp @@ -29,11 +29,18 @@ ERF::FillCoarsePatch (int lev, Real time) //**************************************************************************************************************** // bool cons_only = false; - FillPatch(lev-1, time, {&vars_new[lev-1][Vars::cons], &vars_new[lev-1][Vars::xvel], - &vars_new[lev-1][Vars::yvel], &vars_new[lev-1][Vars::zvel]}, - {&vars_new[lev-1][Vars::cons], - &rU_new[lev-1], &rV_new[lev-1], &rW_new[lev-1]}, - false, cons_only); + if (lev == 1) { + FillPatch(lev-1, time, {&vars_new[lev-1][Vars::cons], &vars_new[lev-1][Vars::xvel], + &vars_new[lev-1][Vars::yvel], &vars_new[lev-1][Vars::zvel]}, + cons_only); + } else { + FillPatch(lev-1, time, {&vars_new[lev-1][Vars::cons], &vars_new[lev-1][Vars::xvel], + &vars_new[lev-1][Vars::yvel], &vars_new[lev-1][Vars::zvel]}, + {&vars_new[lev-1][Vars::cons], + &rU_new[lev-1], &rV_new[lev-1], &rW_new[lev-1]}, + base_state[lev-1], base_state[lev-1], + false, cons_only); + } // // ************************************************ diff --git a/Source/BoundaryConditions/ERF_FillIntermediatePatch.cpp b/Source/BoundaryConditions/ERF_FillIntermediatePatch.cpp index 99be97002..6b8c013e2 100644 --- a/Source/BoundaryConditions/ERF_FillIntermediatePatch.cpp +++ b/Source/BoundaryConditions/ERF_FillIntermediatePatch.cpp @@ -61,28 +61,24 @@ ERF::FillIntermediatePatch (int lev, Real time, } } + // amrex::Print() << "CONS ONLY " << cons_only << std::endl; + // amrex::Print() << "ICOMP NCOMP " << icomp_cons << " " << ncomp_cons << " NGHOST " << ng_cons << std::endl; + AMREX_ALWAYS_ASSERT(mfs_mom.size() == IntVars::NumTypes); AMREX_ALWAYS_ASSERT(mfs_vel.size() == Vars::NumTypes); // Enforce no penetration for thin immersed body - if (xflux_imask[lev]) { - ApplyMask(*mfs_mom[IntVars::xmom], *xflux_imask[lev]); - } - if (yflux_imask[lev]) { - ApplyMask(*mfs_mom[IntVars::ymom], *yflux_imask[lev]); - } - if (zflux_imask[lev]) { - ApplyMask(*mfs_mom[IntVars::zmom], *zflux_imask[lev]); - } - - // We always come in to this call with updated momenta but we need to create updated velocity - // in order to impose the rest of the bc's if (!cons_only) { - // This only fills VALID region of velocity - MomentumToVelocity(*mfs_vel[Vars::xvel], *mfs_vel[Vars::yvel], *mfs_vel[Vars::zvel], - *mfs_vel[Vars::cons], - *mfs_mom[IntVars::xmom], *mfs_mom[IntVars::ymom], *mfs_mom[IntVars::zmom], - Geom(lev).Domain(), domain_bcs_type); + // Enforce no penetration for thin immersed body + if (xflux_imask[lev]) { + ApplyMask(*mfs_mom[IntVars::xmom], *xflux_imask[lev]); + } + if (yflux_imask[lev]) { + ApplyMask(*mfs_mom[IntVars::ymom], *yflux_imask[lev]); + } + if (zflux_imask[lev]) { + ApplyMask(*mfs_mom[IntVars::zmom], *zflux_imask[lev]); + } } // @@ -93,19 +89,82 @@ ERF::FillIntermediatePatch (int lev, Real time, // We don't do anything here because we will call the physbcs routines below, // which calls FillBoundary and fills other domain boundary conditions // Physical boundaries will be filled below + + if (!cons_only) + { + // *************************************************************************** + // We always come in to this call with updated momenta but we need to create updated velocity + // in order to impose the rest of the bc's + // *************************************************************************** + // This only fills VALID region of velocity + MomentumToVelocity(*mfs_vel[Vars::xvel], *mfs_vel[Vars::yvel], *mfs_vel[Vars::zvel], + *mfs_vel[Vars::cons], + *mfs_mom[IntVars::xmom], *mfs_mom[IntVars::ymom], *mfs_mom[IntVars::zmom], + Geom(lev).Domain(), domain_bcs_type); + } } else { - MultiFab& mf = *mfs_vel[Vars::cons]; - - Vector fmf = {&mf,&mf}; + // + // We must fill a temporary then copy it back so we don't double add/subtract + // + MultiFab mf(mfs_vel[Vars::cons]->boxArray(),mfs_vel[Vars::cons]->DistributionMap(), + mfs_vel[Vars::cons]->nComp() ,mfs_vel[Vars::cons]->nGrowVect()); + // + // Set all components to 1.789e19, then copy just the density from *mfs_vel[Vars::cons] + // + mf.setVal(1.789e19); + MultiFab::Copy(mf,*mfs_vel[Vars::cons],Rho_comp,Rho_comp,1,mf.nGrowVect()); + + Vector fmf = {mfs_vel[Vars::cons],mfs_vel[Vars::cons]}; Vector cmf = {&vars_old[lev-1][Vars::cons], &vars_new[lev-1][Vars::cons]}; Vector ctime = {t_old[lev-1], t_new[lev-1]}; Vector ftime = {time,time}; // Impose physical bc's on coarse data (note time and 0 are not used) - (*physbcs_cons[lev-1])(vars_old[lev-1][Vars::cons],0,ncomp_cons,IntVect{ng_cons},time,BCVars::cons_bc,true); - (*physbcs_cons[lev-1])(vars_new[lev-1][Vars::cons],0,ncomp_cons,IntVect{ng_cons},time,BCVars::cons_bc,true); + (*physbcs_cons[lev-1])(vars_old[lev-1][Vars::cons],icomp_cons,ncomp_cons,IntVect{ng_cons},time,BCVars::cons_bc,true); + (*physbcs_cons[lev-1])(vars_new[lev-1][Vars::cons],icomp_cons,ncomp_cons,IntVect{ng_cons},time,BCVars::cons_bc,true); + + // Impose physical bc's on fine data (note time and 0 are not used) + (*physbcs_cons[lev])(*mfs_vel[Vars::cons],icomp_cons,ncomp_cons,IntVect{ng_cons},time,BCVars::cons_bc,true); + + if ( (icomp_cons+ncomp_cons > 1) && (interpolation_type == StateInterpType::Perturbational) ) + { + // Divide (rho theta) by rho to get theta + MultiFab::Divide(*mfs_vel[Vars::cons],*mfs_vel[Vars::cons],Rho_comp,RhoTheta_comp,1,IntVect{0}); + + // Subtract theta_0 from theta + MultiFab::Subtract(*mfs_vel[Vars::cons],base_state[lev],BaseState::th0_comp,RhoTheta_comp,1,IntVect{0}); + + if (!amrex::almostEqual(time,ctime[1])) { + MultiFab::Divide(vars_old[lev-1][Vars::cons], vars_old[lev-1][Vars::cons], + Rho_comp,RhoTheta_comp,1,vars_old[lev-1][Vars::cons].nGrowVect()); + MultiFab::Subtract(vars_old[lev-1][Vars::cons], base_state[lev-1], + BaseState::th0_comp,RhoTheta_comp,1,vars_old[lev-1][Vars::cons].nGrowVect()); + } + if (!amrex::almostEqual(time,ctime[0])) { + MultiFab::Divide(vars_new[lev-1][Vars::cons], vars_new[lev-1][Vars::cons], + Rho_comp,RhoTheta_comp,1,vars_new[lev-1][Vars::cons].nGrowVect()); + MultiFab::Subtract(vars_new[lev-1][Vars::cons], base_state[lev-1], + BaseState::th0_comp,RhoTheta_comp,1,vars_new[lev-1][Vars::cons].nGrowVect()); + } + } + + // Subtract rho_0 from rho before we interpolate -- note we only subtract + // on valid region of mf since the ghost cells will be filled below + if (icomp_cons == 0 && (interpolation_type == StateInterpType::Perturbational)) + { + MultiFab::Subtract(*mfs_vel[Vars::cons],base_state[lev],BaseState::r0_comp,Rho_comp,1,IntVect{0}); + + if (!amrex::almostEqual(time,ctime[1])) { + MultiFab::Subtract(vars_old[lev-1][Vars::cons], base_state[lev-1], + BaseState::r0_comp,Rho_comp,1,vars_old[lev-1][Vars::cons].nGrowVect()); + } + if (!amrex::almostEqual(time,ctime[0])) { + MultiFab::Subtract(vars_new[lev-1][Vars::cons], base_state[lev-1], + BaseState::r0_comp,Rho_comp,1,vars_new[lev-1][Vars::cons].nGrowVect()); + } + } // Call FillPatchTwoLevels which ASSUMES that all ghost cells have already been filled mapper = &cell_cons_interp; @@ -115,10 +174,73 @@ ERF::FillIntermediatePatch (int lev, Real time, refRatio(lev-1), mapper, domain_bcs_type, icomp_cons); + if (icomp_cons == 0 && (interpolation_type == StateInterpType::Perturbational)) + { + // Restore the coarse values to what they were + if (!amrex::almostEqual(time,ctime[1])) { + MultiFab::Add(vars_old[lev-1][Vars::cons], base_state[lev-1], + BaseState::r0_comp,Rho_comp,1,vars_old[lev-1][Vars::cons].nGrowVect()); + } + if (!amrex::almostEqual(time,ctime[0])) { + MultiFab::Add(vars_new[lev-1][Vars::cons], base_state[lev-1], + BaseState::r0_comp,Rho_comp,1,vars_new[lev-1][Vars::cons].nGrowVect()); + } + + // Set values in the cells outside the domain boundary so that we can do the Add + // without worrying about uninitialized values outside the domain -- these + // will be filled in the physbcs call + mf.setDomainBndry(1.234e20,Rho_comp,1,geom[lev]); + + // Add rho_0 back to rho after we interpolate -- on all the valid + ghost region + MultiFab::Add(mf, base_state[lev],BaseState::r0_comp,Rho_comp,1,IntVect{ng_cons}); + } + + if ( (icomp_cons+ncomp_cons > 1) && (interpolation_type == StateInterpType::Perturbational) ) + { + // Add theta_0 to theta + if (!amrex::almostEqual(time,ctime[1])) { + MultiFab::Add(vars_old[lev-1][Vars::cons], base_state[lev-1], + BaseState::th0_comp,RhoTheta_comp,1,vars_old[lev-1][Vars::cons].nGrowVect()); + MultiFab::Multiply(vars_old[lev-1][Vars::cons], vars_old[lev-1][Vars::cons], + Rho_comp,RhoTheta_comp,1,vars_old[lev-1][Vars::cons].nGrowVect()); + } + if (!amrex::almostEqual(time,ctime[0])) { + MultiFab::Add(vars_new[lev-1][Vars::cons], base_state[lev-1], + BaseState::th0_comp,RhoTheta_comp,1,vars_new[lev-1][Vars::cons].nGrowVect()); + MultiFab::Multiply(vars_new[lev-1][Vars::cons], vars_new[lev-1][Vars::cons], + Rho_comp,RhoTheta_comp,1,vars_new[lev-1][Vars::cons].nGrowVect()); + } + + // Multiply theta by rho to get (rho theta) + MultiFab::Multiply(*mfs_vel[Vars::cons],*mfs_vel[Vars::cons],Rho_comp,RhoTheta_comp,1,IntVect{0}); + + // Add theta_0 to theta + MultiFab::Add(*mfs_vel[Vars::cons],base_state[lev],BaseState::th0_comp,RhoTheta_comp,1,IntVect{0}); + + // Add theta_0 back to theta + MultiFab::Add(mf,base_state[lev],BaseState::th0_comp,RhoTheta_comp,1,IntVect{ng_cons}); + + // Multiply (theta) by rho to get (rho theta) + MultiFab::Multiply(mf,mf,Rho_comp,RhoTheta_comp,1,IntVect{ng_cons}); + } + + // Make sure to only copy back the components we worked on + MultiFab::Copy(*mfs_vel[Vars::cons],mf,icomp_cons,icomp_cons,ncomp_cons,IntVect{ng_cons}); + // ***************************************************************************************** if (!cons_only) { + // *************************************************************************** + // We always come in to this call with updated momenta but we need to create updated velocity + // in order to impose the rest of the bc's + // *************************************************************************** + // This only fills VALID region of velocity + MomentumToVelocity(*mfs_vel[Vars::xvel], *mfs_vel[Vars::yvel], *mfs_vel[Vars::zvel], + *mfs_vel[Vars::cons], + *mfs_mom[IntVars::xmom], *mfs_mom[IntVars::ymom], *mfs_mom[IntVars::zmom], + Geom(lev).Domain(), domain_bcs_type); + mapper = &face_cons_linear_interp; // @@ -201,7 +323,7 @@ ERF::FillIntermediatePatch (int lev, Real time, #ifdef ERF_USE_NETCDF // We call this here because it is an ERF routine if (use_real_bcs && (lev==0)) { - fill_from_realbdy(mfs_vel,time,cons_only,icomp_cons,ncomp_cons,ngvect_cons, ngvect_vels); + fill_from_realbdy(mfs_vel,time,cons_only,icomp_cons,ncomp_cons,ngvect_cons,ngvect_vels); do_fb = false; } #endif diff --git a/Source/BoundaryConditions/ERF_FillPatch.cpp b/Source/BoundaryConditions/ERF_FillPatch.cpp index 03f272b55..5913f110b 100644 --- a/Source/BoundaryConditions/ERF_FillPatch.cpp +++ b/Source/BoundaryConditions/ERF_FillPatch.cpp @@ -20,11 +20,15 @@ void ERF::FillPatch (int lev, Real time, const Vector& mfs_vel, // This includes cc quantities and VELOCITIES const Vector& mfs_mom, // This includes cc quantities and MOMENTA + const MultiFab& old_base_state, + const MultiFab& new_base_state, bool fillset, bool cons_only) { BL_PROFILE_VAR("ERF::FillPatch()",ERF_FillPatch); Interpolater* mapper = nullptr; + AMREX_ALWAYS_ASSERT(lev > 0); + PhysBCFunctNoOp null_bc; // @@ -38,7 +42,7 @@ ERF::FillPatch (int lev, Real time, // conditions are imposed on velocity, so we convert to momentum here then // convert back. // *************************************************************************** - if (lev>0 && fillset) { + if (fillset) { if (cf_set_width > 0) { FPr_c[lev-1].FillSet(*mfs_vel[Vars::cons], time, null_bc, domain_bcs_type); } @@ -70,48 +74,56 @@ ERF::FillPatch (int lev, Real time, IntVect ngvect_cons = mfs_vel[Vars::cons]->nGrowVect(); IntVect ngvect_vels = mfs_vel[Vars::xvel]->nGrowVect(); - if (lev == 0) { Vector ftime = {t_old[lev], t_new[lev]}; - - // - // Below we call FillPatchSingleLevel which does NOT fill ghost cells outside the domain - // - - Vector fmf = {&vars_old[lev][Vars::cons], &vars_new[lev][Vars::cons]}; - const int ncomp = mfs_vel[Vars::cons]->nComp(); - - FillPatchSingleLevel(*mfs_vel[Vars::cons], ngvect_cons, time, fmf, IntVect(0,0,0), ftime, - 0, 0, ncomp, geom[lev]); - - if (!cons_only) { - fmf = {&vars_old[lev][Vars::xvel], &vars_new[lev][Vars::xvel]}; - FillPatchSingleLevel(*mfs_vel[Vars::xvel], ngvect_vels, time, fmf, - IntVect(0,0,0), ftime, 0, 0, 1, geom[lev]); - - fmf = {&vars_old[lev][Vars::yvel], &vars_new[lev][Vars::yvel]}; - FillPatchSingleLevel(*mfs_vel[Vars::yvel], ngvect_vels, time, fmf, - IntVect(0,0,0), ftime, 0, 0, 1, geom[lev]); - - fmf = {&vars_old[lev][Vars::zvel], &vars_new[lev][Vars::zvel]}; - FillPatchSingleLevel(*mfs_vel[Vars::zvel], ngvect_vels, time, fmf, - IntVect(0,0,0), ftime, 0, 0, 1, geom[lev]); - } // !cons_only - - } else { - - Vector ftime = {t_old[lev], t_new[lev]}; Vector ctime = {t_old[lev-1], t_new[lev-1]}; Vector fmf = {&vars_old[lev ][Vars::cons], &vars_new[lev ][Vars::cons]}; Vector cmf = {&vars_old[lev-1][Vars::cons], &vars_new[lev-1][Vars::cons]}; - MultiFab& mf_c = *mfs_vel[Vars::cons]; + + // We must fill a temporary then copy it back so we don't double add/subtract + MultiFab mf_c(mfs_vel[Vars::cons]->boxArray(),mfs_vel[Vars::cons]->DistributionMap(), + mfs_vel[Vars::cons]->nComp() ,mfs_vel[Vars::cons]->nGrowVect()); + mapper = &cell_cons_interp; - // Impose physical bc's on coarse data (note time and 0 are not used) - // Note that we call FillBoundary inside the physbcs call - // We should not need to call this on old data since that would have been filled before the timestep started - (*physbcs_cons[lev-1])(vars_new[lev-1][Vars::cons],0,mf_c.nComp(),ngvect_cons,time,BCVars::cons_bc,true); + if (interpolation_type == StateInterpType::Perturbational) + { + // Divide (rho theta) by rho to get theta (before we subtract rho0 from rho!) + if (!amrex::almostEqual(time,ctime[1])) { + MultiFab::Divide(vars_old[lev-1][Vars::cons],vars_old[lev-1][Vars::cons], + Rho_comp,RhoTheta_comp,1,ngvect_cons); + MultiFab::Subtract(vars_old[lev-1][Vars::cons],base_state[lev-1], + BaseState::r0_comp,Rho_comp,1,ngvect_cons); + MultiFab::Subtract(vars_old[lev-1][Vars::cons],base_state[lev-1], + BaseState::th0_comp,RhoTheta_comp,1,ngvect_cons); + } + if (!amrex::almostEqual(time,ctime[0])) { + MultiFab::Divide(vars_new[lev-1][Vars::cons],vars_new[lev-1][Vars::cons], + Rho_comp,RhoTheta_comp,1,ngvect_cons); + MultiFab::Subtract(vars_new[lev-1][Vars::cons],base_state[lev-1], + BaseState::r0_comp,Rho_comp,1,ngvect_cons); + MultiFab::Subtract(vars_new[lev-1][Vars::cons],base_state[lev-1], + BaseState::th0_comp,RhoTheta_comp,1,ngvect_cons); + } + + if (!amrex::almostEqual(time,ftime[1])) { + MultiFab::Divide(vars_old[lev ][Vars::cons],vars_old[lev ][Vars::cons], + Rho_comp,RhoTheta_comp,1,IntVect{0}); + MultiFab::Subtract(vars_old[lev ][Vars::cons],old_base_state, + BaseState::r0_comp,Rho_comp,1,IntVect{0}); + MultiFab::Subtract(vars_old[lev ][Vars::cons],old_base_state, + BaseState::th0_comp,RhoTheta_comp,1,IntVect{0}); + } + if (!amrex::almostEqual(time,ftime[0])) { + MultiFab::Divide(vars_new[lev ][Vars::cons],vars_new[lev ][Vars::cons], + Rho_comp,RhoTheta_comp,1,IntVect{0}); + MultiFab::Subtract(vars_new[lev ][Vars::cons],old_base_state, + BaseState::r0_comp,Rho_comp,1,IntVect{0}); + MultiFab::Subtract(vars_new[lev ][Vars::cons],old_base_state, + BaseState::th0_comp,RhoTheta_comp,1,IntVect{0}); + } + } // Call FillPatchTwoLevels which ASSUMES that all ghost cells have already been filled FillPatchTwoLevels(mf_c, ngvect_cons, IntVect(0,0,0), @@ -120,6 +132,56 @@ ERF::FillPatch (int lev, Real time, refRatio(lev-1), mapper, domain_bcs_type, BCVars::cons_bc); + if (interpolation_type == StateInterpType::Perturbational) + { + // Restore the coarse values to what they were + if (!amrex::almostEqual(time,ctime[1])) { + MultiFab::Add(vars_old[lev-1][Vars::cons], base_state[lev-1], + BaseState::r0_comp,Rho_comp,1,ngvect_cons); + MultiFab::Add(vars_old[lev-1][Vars::cons], base_state[lev-1], + BaseState::th0_comp,RhoTheta_comp,1,ngvect_cons); + MultiFab::Multiply(vars_old[lev-1][Vars::cons], vars_old[lev-1][Vars::cons], + Rho_comp,RhoTheta_comp,1,ngvect_cons); + } + if (!amrex::almostEqual(time,ctime[0])) { + MultiFab::Add(vars_new[lev-1][Vars::cons], base_state[lev-1], + BaseState::r0_comp,Rho_comp,1,vars_new[lev-1][Vars::cons].nGrowVect()); + MultiFab::Add(vars_new[lev-1][Vars::cons], base_state[lev-1], + BaseState::th0_comp,RhoTheta_comp,1,vars_new[lev-1][Vars::cons].nGrowVect()); + MultiFab::Multiply(vars_new[lev-1][Vars::cons], vars_new[lev-1][Vars::cons], + Rho_comp,RhoTheta_comp,1,ngvect_cons); + } + + if (!amrex::almostEqual(time,ftime[1])) { + MultiFab::Add(vars_old[lev][Vars::cons],base_state[lev ],BaseState::r0_comp,Rho_comp,1,ngvect_cons); + MultiFab::Add(vars_old[lev][Vars::cons],base_state[lev ],BaseState::th0_comp,RhoTheta_comp,1,ngvect_cons); + MultiFab::Multiply(vars_old[lev][Vars::cons], vars_old[lev][Vars::cons], + Rho_comp,RhoTheta_comp,1,ngvect_cons); + } + if (!amrex::almostEqual(time,ftime[0])) { + MultiFab::Add(vars_new[lev][Vars::cons], base_state[lev],BaseState::r0_comp,Rho_comp,1,ngvect_cons); + MultiFab::Add(vars_new[lev][Vars::cons], base_state[lev],BaseState::th0_comp,RhoTheta_comp,1,ngvect_cons); + MultiFab::Multiply(vars_new[lev][Vars::cons], vars_new[lev][Vars::cons], + Rho_comp,RhoTheta_comp,1,ngvect_cons); + } + + // Set values in the cells outside the domain boundary so that we can do the Add + // without worrying about uninitialized values outside the domain -- these + // will be filled in the physbcs call + mf_c.setDomainBndry(1.234e20,0,2,geom[lev]); // Do both rho and (rho theta) together + + // Add rho_0 back to rho and theta_0 back to theta + MultiFab::Add(mf_c, new_base_state,BaseState::r0_comp,Rho_comp,1,ngvect_cons); + MultiFab::Add(mf_c, new_base_state,BaseState::th0_comp,RhoTheta_comp,1,ngvect_cons); + + // Multiply (theta) by rho to get (rho theta) + MultiFab::Multiply(mf_c,mf_c,Rho_comp,RhoTheta_comp,1,ngvect_cons); + } + + MultiFab::Copy(*mfs_vel[Vars::cons],mf_c,0,0,mf_c.nComp(),mf_c.nGrowVect()); + + // *************************************************************************************** + if (!cons_only) { mapper = &face_cons_linear_interp; @@ -130,13 +192,6 @@ ERF::FillPatch (int lev, Real time, // ********************************************************************** - cmf = {&vars_old[lev-1][Vars::xvel], &vars_new[lev-1][Vars::xvel]}; - - // Impose physical bc's on coarse data (note time and 0 are not used) - // Note that we call FillBoundary inside the physbcs call - // We should not need to call this on old data since that would have been filled before the timestep started - (*physbcs_u[lev-1])(vars_new[lev-1][Vars::xvel],0,1,ngvect_vels,time,BCVars::xvel_bc,true); - fmf = {&vars_old[lev ][Vars::xvel], &vars_new[lev ][Vars::xvel]}; cmf = {&vars_old[lev-1][Vars::xvel], &vars_new[lev-1][Vars::xvel]}; @@ -149,13 +204,6 @@ ERF::FillPatch (int lev, Real time, // ********************************************************************** - cmf = {&vars_old[lev-1][Vars::yvel], &vars_new[lev-1][Vars::yvel]}; - - // Impose physical bc's on coarse data (note time and 0 are not used) - // Note that we call FillBoundary inside the physbcs call - // We should not need to call this on old data since that would have been filled before the timestep started - (*physbcs_v[lev-1])(vars_new[lev-1][Vars::yvel],0,1,ngvect_vels,time,BCVars::yvel_bc,true); - fmf = {&vars_old[lev ][Vars::yvel], &vars_new[lev ][Vars::yvel]}; cmf = {&vars_old[lev-1][Vars::yvel], &vars_new[lev-1][Vars::yvel]}; @@ -168,16 +216,6 @@ ERF::FillPatch (int lev, Real time, // ********************************************************************** - cmf = {&vars_old[lev-1][Vars::zvel], &vars_new[lev-1][Vars::zvel]}; - - // Impose physical bc's on coarse data (note time and 0 are not used) - // Note that we call FillBoundary inside the physbcs call - // We should not need to call this on old data since that would have been filled before the timestep started - (*physbcs_w[lev-1])(vars_new[lev-1][Vars::zvel], - vars_new[lev-1][Vars::xvel], - vars_new[lev-1][Vars::yvel], - ngvect_vels,time,BCVars::zvel_bc,true); - fmf = {&vars_old[lev ][Vars::zvel], &vars_new[lev ][Vars::zvel]}; cmf = {&vars_old[lev-1][Vars::zvel], &vars_new[lev-1][Vars::zvel]}; @@ -198,6 +236,65 @@ ERF::FillPatch (int lev, Real time, bool do_fb = true; + if (m_r2d) fill_from_bndryregs(mfs_vel,time); + + // We call these even if init_type == InitType::Real because these will fill the vertical bcs + // Note that we call FillBoundary inside the physbcs call + (*physbcs_cons[lev])(*mfs_vel[Vars::cons],icomp_cons,ncomp_cons,ngvect_cons,time,BCVars::cons_bc, do_fb); + if (!cons_only) { + (*physbcs_u[lev])(*mfs_vel[Vars::xvel],0,1,ngvect_vels,time,BCVars::xvel_bc, do_fb); + (*physbcs_v[lev])(*mfs_vel[Vars::yvel],0,1,ngvect_vels,time,BCVars::yvel_bc, do_fb); + (*physbcs_w[lev])(*mfs_vel[Vars::zvel],*mfs_vel[Vars::xvel],*mfs_vel[Vars::yvel], + ngvect_vels,time,BCVars::zvel_bc, do_fb); + } +} + +void +ERF::FillPatch (int lev, Real time, + const Vector& mfs_vel, // This includes cc quantities and VELOCITIES + bool cons_only) +{ + BL_PROFILE_VAR("ERF::FillPatch()",ERF_FillPatch); + + AMREX_ALWAYS_ASSERT(lev == 0); + + IntVect ngvect_cons = mfs_vel[Vars::cons]->nGrowVect(); + IntVect ngvect_vels = mfs_vel[Vars::xvel]->nGrowVect(); + + Vector ftime = {t_old[lev], t_new[lev]}; + + // + // Below we call FillPatchSingleLevel which does NOT fill ghost cells outside the domain + // + + Vector fmf = {&vars_old[lev][Vars::cons], &vars_new[lev][Vars::cons]}; + const int ncomp = mfs_vel[Vars::cons]->nComp(); + + FillPatchSingleLevel(*mfs_vel[Vars::cons], ngvect_cons, time, fmf, IntVect(0,0,0), ftime, + 0, 0, ncomp, geom[lev]); + + if (!cons_only) { + fmf = {&vars_old[lev][Vars::xvel], &vars_new[lev][Vars::xvel]}; + FillPatchSingleLevel(*mfs_vel[Vars::xvel], ngvect_vels, time, fmf, + IntVect(0,0,0), ftime, 0, 0, 1, geom[lev]); + + fmf = {&vars_old[lev][Vars::yvel], &vars_new[lev][Vars::yvel]}; + FillPatchSingleLevel(*mfs_vel[Vars::yvel], ngvect_vels, time, fmf, + IntVect(0,0,0), ftime, 0, 0, 1, geom[lev]); + + fmf = {&vars_old[lev][Vars::zvel], &vars_new[lev][Vars::zvel]}; + FillPatchSingleLevel(*mfs_vel[Vars::zvel], ngvect_vels, time, fmf, + IntVect(0,0,0), ftime, 0, 0, 1, geom[lev]); + } // !cons_only + + // *************************************************************************** + // Physical bc's at domain boundary + // *************************************************************************** + int icomp_cons = 0; + int ncomp_cons = mfs_vel[Vars::cons]->nComp(); + + bool do_fb = true; + #ifdef ERF_USE_NETCDF // We call this here because it is an ERF routine if (use_real_bcs && (lev==0)) { diff --git a/Source/BoundaryConditions/ERF_PhysBCFunct.H b/Source/BoundaryConditions/ERF_PhysBCFunct.H index 3666ea06c..756e4725d 100644 --- a/Source/BoundaryConditions/ERF_PhysBCFunct.H +++ b/Source/BoundaryConditions/ERF_PhysBCFunct.H @@ -262,8 +262,10 @@ class ERFPhysBCFunct_base public: ERFPhysBCFunct_base (const int lev, const amrex::Geometry& geom, const amrex::Vector& domain_bcs_type, - const amrex::Gpu::DeviceVector& domain_bcs_type_d) + const amrex::Gpu::DeviceVector& domain_bcs_type_d, + bool moving_terrain) : m_lev(lev), m_geom(geom), + m_moving_terrain(moving_terrain), m_domain_bcs_type(domain_bcs_type), m_domain_bcs_type_d(domain_bcs_type_d) {} @@ -291,6 +293,7 @@ public: private: int m_lev; amrex::Geometry m_geom; + bool m_moving_terrain; amrex::Vector m_domain_bcs_type; amrex::Gpu::DeviceVector m_domain_bcs_type_d; }; diff --git a/Source/BoundaryConditions/ERF_PhysBCFunct.cpp b/Source/BoundaryConditions/ERF_PhysBCFunct.cpp index 3b055a9fa..fb5a18f8f 100644 --- a/Source/BoundaryConditions/ERF_PhysBCFunct.cpp +++ b/Source/BoundaryConditions/ERF_PhysBCFunct.cpp @@ -372,6 +372,11 @@ void ERFPhysBCFunct_base::operator() (MultiFab& mf, int /*icomp*/, int ncomp, In if (m_geom.isAllPeriodic()) return; + if (m_moving_terrain) { + mf.FillBoundary(m_geom.periodicity()); + return; + } + const auto& domain = m_geom.Domain(); // Create a grown domain box containing valid + periodic cells @@ -382,6 +387,12 @@ void ERFPhysBCFunct_base::operator() (MultiFab& mf, int /*icomp*/, int ncomp, In } } + // + // We fill all of the interior and periodic ghost cells first, so we can fill + // those directly inside the lateral and vertical calls. + // + mf.FillBoundary(m_geom.periodicity()); + #ifdef AMREX_USE_OMP #pragma omp parallel if (Gpu::notInLaunchRegion()) #endif diff --git a/Source/ERF.H b/Source/ERF.H index 3c3dda751..2af6c99f9 100644 --- a/Source/ERF.H +++ b/Source/ERF.H @@ -78,6 +78,13 @@ AMREX_ENUM(InitType, None, Input_Sounding, Ideal, Real, Metgrid, Uniform ); +/** + * Enum of possible interpolation types between coarse/fine +*/ +AMREX_ENUM(StateInterpType, + FullState, Perturbational +); + /** * Enum of possible plotfile types */ @@ -86,6 +93,7 @@ AMREX_ENUM(PlotFileType, ); +#if 0 /** * Enum of possible coarse/fine interpolation options */ @@ -102,6 +110,7 @@ namespace InterpType { FaceConserativeLinear }; } +#endif /** * Main class in ERF code, instantiated from main.cpp @@ -516,9 +525,17 @@ private: // // NOTE: FillPatch takes in an empty MF, and returns cell-centered + velocities (not momenta) // + // This one works only at level = 0 (base state does not change) + void FillPatch (int lev, amrex::Real time, + const amrex::Vector& mfs_vel, + bool cons_only=false); + + // This one works only at level > 0 (base state does change) void FillPatch (int lev, amrex::Real time, const amrex::Vector& mfs_vel, const amrex::Vector& mfs_mom, + const amrex::MultiFab& old_base_state, + const amrex::MultiFab& new_base_state, bool fillset=true, bool cons_only=false); // Compute new multifabs by copying data from valid region and filling ghost cells. @@ -978,6 +995,7 @@ private: static PlotFileType plotfile_type_2; static InitType init_type; + static StateInterpType interpolation_type; // sponge_type: "input_sponge" static std::string sponge_type; diff --git a/Source/ERF.cpp b/Source/ERF.cpp index 1c81680fb..a39415562 100644 --- a/Source/ERF.cpp +++ b/Source/ERF.cpp @@ -45,6 +45,7 @@ PlotFileType ERF::plotfile_type_1 = PlotFileType::None; PlotFileType ERF::plotfile_type_2 = PlotFileType::None; InitType ERF::init_type; +StateInterpType ERF::interpolation_type; // use_real_bcs: only true if 1) ( (init_type == InitType::Real) or (init_type == InitGrid::Metgrid) ) // AND 2) we want to use the bc's from the WRF bdy file @@ -874,6 +875,18 @@ ERF::InitData_post () int ncomp = lev_new[Vars::cons].nComp(); + // *************************************************************************** + // Physical bc's at domain boundary + // *************************************************************************** + IntVect ngvect_cons = vars_new[lev][Vars::cons].nGrowVect(); + IntVect ngvect_vels = vars_new[lev][Vars::xvel].nGrowVect(); + + (*physbcs_cons[lev])(lev_new[Vars::cons],0,ncomp,ngvect_cons,t_new[lev],BCVars::cons_bc,true); + ( *physbcs_u[lev])(lev_new[Vars::xvel],0,1 ,ngvect_vels,t_new[lev],BCVars::xvel_bc,true); + ( *physbcs_v[lev])(lev_new[Vars::yvel],0,1 ,ngvect_vels,t_new[lev],BCVars::yvel_bc,true); + ( *physbcs_w[lev])(lev_new[Vars::zvel],lev_new[Vars::xvel],lev_new[Vars::yvel], + ngvect_vels,t_new[lev],BCVars::zvel_bc,true); + MultiFab::Copy(lev_old[Vars::cons],lev_new[Vars::cons],0,0,ncomp,lev_new[Vars::cons].nGrowVect()); MultiFab::Copy(lev_old[Vars::xvel],lev_new[Vars::xvel],0,0, 1,lev_new[Vars::xvel].nGrowVect()); MultiFab::Copy(lev_old[Vars::yvel],lev_new[Vars::yvel],0,0, 1,lev_new[Vars::yvel].nGrowVect()); @@ -905,10 +918,16 @@ ERF::InitData_post () // Fill boundary conditions -- not sure why we need this here // bool fillset = false; - FillPatch(lev, t_new[lev], - {&lev_new[Vars::cons],&lev_new[Vars::xvel],&lev_new[Vars::yvel],&lev_new[Vars::zvel]}, - {&lev_new[Vars::cons],&rU_new[lev],&rV_new[lev],&rW_new[lev]}, - fillset); + if (lev == 0) { + FillPatch(lev, t_new[lev], + {&lev_new[Vars::cons],&lev_new[Vars::xvel],&lev_new[Vars::yvel],&lev_new[Vars::zvel]}); + } else { + FillPatch(lev, t_new[lev], + {&lev_new[Vars::cons],&lev_new[Vars::xvel],&lev_new[Vars::yvel],&lev_new[Vars::zvel]}, + {&lev_new[Vars::cons],&rU_new[lev],&rV_new[lev],&rW_new[lev]}, + base_state[lev], base_state[lev], + fillset); + } // // We do this here to make sure level (lev-1) boundary conditions are filled @@ -1460,6 +1479,10 @@ ERF::ReadParameters () // Flag to trigger initialization from input_sounding like WRF's ideal.exe pp.query("init_sounding_ideal", init_sounding_ideal); + // Set default to FullState for now ... later we will try Perturbation + interpolation_type = StateInterpType::FullState; + pp.query_enum_case_insensitive("interpolation_type" ,interpolation_type); + PlotFileType plotfile_type_temp = PlotFileType::None; pp.query_enum_case_insensitive("plotfile_type" ,plotfile_type_temp); pp.query_enum_case_insensitive("plotfile_type_1",plotfile_type_1); @@ -1814,115 +1837,6 @@ ERF::MakeDiagnosticAverage (Vector& h_havg, MultiFab& S, int n) } } -// Set covered coarse cells to be the average of overlying fine cells for all levels -void -ERF::AverageDown () -{ - AMREX_ALWAYS_ASSERT(solverChoice.coupling_type == CouplingType::TwoWay); - int src_comp = 0; - int num_comp = vars_new[0][Vars::cons].nComp(); - for (int lev = finest_level-1; lev >= 0; --lev) - { - AverageDownTo(lev,src_comp,num_comp); - } -} - -// Set covered coarse cells to be the average of overlying fine cells at level crse_lev -void -ERF::AverageDownTo (int crse_lev, int scomp, int ncomp) // NOLINT -{ - AMREX_ALWAYS_ASSERT(solverChoice.coupling_type == CouplingType::TwoWay); - - // ****************************************************************************************** - // First do cell-centered quantities - // The quantity that is conserved is not (rho S), but rather (rho S / m^2) where - // m is the map scale factor at cell centers - // Here we pre-divide (rho S) by m^2 before average down - // ****************************************************************************************** - for (int lev = crse_lev; lev <= crse_lev+1; lev++) { - for (MFIter mfi(vars_new[lev][Vars::cons], TilingIfNotGPU()); mfi.isValid(); ++mfi) { - const Box& bx = mfi.tilebox(); - const Array4< Real> cons_arr = vars_new[lev][Vars::cons].array(mfi); - const Array4 mapfac_arr = mapfac_m[lev]->const_array(mfi); - if (solverChoice.use_terrain) { - const Array4 detJ_arr = detJ_cc[lev]->const_array(mfi); - ParallelFor(bx, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept - { - cons_arr(i,j,k,scomp+n) *= detJ_arr(i,j,k) / (mapfac_arr(i,j,0)*mapfac_arr(i,j,0)); - }); - } else { - ParallelFor(bx, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept - { - cons_arr(i,j,k,scomp+n) /= (mapfac_arr(i,j,0)*mapfac_arr(i,j,0)); - }); - } - } // mfi - } // lev - - average_down(vars_new[crse_lev+1][Vars::cons], - vars_new[crse_lev ][Vars::cons], - scomp, ncomp, refRatio(crse_lev)); - - // Here we multiply (rho S) by m^2 after average down - for (int lev = crse_lev; lev <= crse_lev+1; lev++) { - for (MFIter mfi(vars_new[lev][Vars::cons], TilingIfNotGPU()); mfi.isValid(); ++mfi) { - const Box& bx = mfi.tilebox(); - const Array4< Real> cons_arr = vars_new[lev][Vars::cons].array(mfi); - const Array4 mapfac_arr = mapfac_m[lev]->const_array(mfi); - if (solverChoice.use_terrain) { - const Array4 detJ_arr = detJ_cc[lev]->const_array(mfi); - ParallelFor(bx, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept - { - cons_arr(i,j,k,scomp+n) *= (mapfac_arr(i,j,0)*mapfac_arr(i,j,0)) / detJ_arr(i,j,k); - }); - } else { - ParallelFor(bx, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept - { - cons_arr(i,j,k,scomp+n) *= (mapfac_arr(i,j,0)*mapfac_arr(i,j,0)); - }); - } - } // mfi - } // lev - - // ****************************************************************************************** - // Now average down momenta. - // Note that vars_new holds velocities not momenta, but we want to do conservative - // averaging so we first convert to momentum, then average down, then convert - // back to velocities -- only on the valid region - // ****************************************************************************************** - for (int lev = crse_lev; lev <= crse_lev+1; lev++) - { - // FillBoundary for density so we can go back and forth between velocity and momentum - vars_new[lev][Vars::cons].FillBoundary(geom[lev].periodicity()); - - VelocityToMomentum(vars_new[lev][Vars::xvel], IntVect(0,0,0), - vars_new[lev][Vars::yvel], IntVect(0,0,0), - vars_new[lev][Vars::zvel], IntVect(0,0,0), - vars_new[lev][Vars::cons], - rU_new[lev], - rV_new[lev], - rW_new[lev], - Geom(lev).Domain(), - domain_bcs_type); - } - - average_down_faces(rU_new[crse_lev+1], rU_new[crse_lev], refRatio(crse_lev), geom[crse_lev]); - average_down_faces(rV_new[crse_lev+1], rV_new[crse_lev], refRatio(crse_lev), geom[crse_lev]); - average_down_faces(rW_new[crse_lev+1], rW_new[crse_lev], refRatio(crse_lev), geom[crse_lev]); - - for (int lev = crse_lev; lev <= crse_lev+1; lev++) { - MomentumToVelocity(vars_new[lev][Vars::xvel], - vars_new[lev][Vars::yvel], - vars_new[lev][Vars::zvel], - vars_new[lev][Vars::cons], - rU_new[lev], - rV_new[lev], - rW_new[lev], - Geom(lev).Domain(), - domain_bcs_type); - } -} - void ERF::Construct_ERFFillPatchers (int lev) { diff --git a/Source/ERF_make_new_arrays.cpp b/Source/ERF_make_new_arrays.cpp index 39c89acf5..a93d83780 100644 --- a/Source/ERF_make_new_arrays.cpp +++ b/Source/ERF_make_new_arrays.cpp @@ -26,13 +26,16 @@ ERF::init_stuff (int lev, const BoxArray& ba, const DistributionMapping& dm, std::unique_ptr& tmp_zphys_nd) { // ******************************************************************************************** - // Base state holds r_0, pres_0, pi_0 (in that order) + // Base state holds r_0, pres_0, pi_0, th_0 (in that order) + // + // Here is where we set 3 ghost cells for the base state! + // // ******************************************************************************************** tmp_base_state.define(ba,dm,BaseState::num_comps,3); tmp_base_state.setVal(0.); if (solverChoice.use_terrain && solverChoice.terrain_type == TerrainType::Moving) { - base_state_new[lev].define(ba,dm,BaseState::num_comps,3); + base_state_new[lev].define(ba,dm,BaseState::num_comps,base_state[lev].nGrowVect()); base_state_new[lev].setVal(0.); } @@ -557,5 +560,6 @@ ERF::make_physbcs (int lev) m_bc_extdir_vals, m_bc_neumann_vals, solverChoice.terrain_type, z_phys_nd[lev], use_real_bcs, zvel_bc_data[lev].data()); - physbcs_base[lev] = std::make_unique (lev, geom[lev], domain_bcs_type, domain_bcs_type_d); + physbcs_base[lev] = std::make_unique (lev, geom[lev], domain_bcs_type, domain_bcs_type_d, + (solverChoice.terrain_type == TerrainType::Moving)); } diff --git a/Source/ERF_make_new_level.cpp b/Source/ERF_make_new_level.cpp index f248edf5b..282114b7d 100644 --- a/Source/ERF_make_new_level.cpp +++ b/Source/ERF_make_new_level.cpp @@ -47,7 +47,7 @@ void ERF::MakeNewLevelFromScratch (int lev, Real time, const BoxArray& ba_in, // Define dmap[lev] to be dm SetDistributionMap(lev, dm); - amrex::Print() <<" BA FROM SCRATCH AT LEVEL " << lev << " " << ba << std::endl; + // amrex::Print() <<" BA FROM SCRATCH AT LEVEL " << lev << " " << ba << std::endl; if (lev == 0) init_bcs(); @@ -283,7 +283,7 @@ ERF::MakeNewLevelFromCoarse (int lev, Real time, const BoxArray& ba, void ERF::RemakeLevel (int lev, Real time, const BoxArray& ba, const DistributionMapping& dm) { - amrex::Print() <<" REMAKING WITH NEW BA AT LEVEL " << lev << " " << ba << std::endl; + // amrex::Print() <<" REMAKING WITH NEW BA AT LEVEL " << lev << " " << ba << std::endl; AMREX_ALWAYS_ASSERT(lev > 0); AMREX_ALWAYS_ASSERT(solverChoice.terrain_type != TerrainType::Moving); @@ -332,39 +332,40 @@ ERF::RemakeLevel (int lev, Real time, const BoxArray& ba, const DistributionMapp // ***************************************************************************************************** make_physbcs(lev); - // ************************************************************************************************* - // This will fill the temporary MultiFabs with data from vars_new - // NOTE: the momenta here are only used as scratch space, the momenta themselves are not fillpatched - // ************************************************************************************************* - FillPatch(lev, time, {&temp_lev_new[Vars::cons],&temp_lev_new[Vars::xvel], - &temp_lev_new[Vars::yvel],&temp_lev_new[Vars::zvel]}, - {&temp_lev_new[Vars::cons],&rU_new[lev],&rV_new[lev],&rW_new[lev]}, - false); - // ******************************************************************************************** // Update the base state at this level by interpolation from coarser level AND copy // from previous (pre-regrid) base_state array // ******************************************************************************************** - if (lev > 0) { - Interpolater* mapper = &cell_cons_interp; + Interpolater* mapper = &cell_cons_interp; - Vector fmf = {&base_state[lev ], &base_state[lev ]}; - Vector cmf = {&base_state[lev-1], &base_state[lev-1]}; - Vector ftime = {time, time}; - Vector ctime = {time, time}; + Vector fmf = {&base_state[lev ], &base_state[lev ]}; + Vector cmf = {&base_state[lev-1], &base_state[lev-1]}; + Vector ftime = {time, time}; + Vector ctime = {time, time}; - // Call FillPatch which ASSUMES that all ghost cells at lev-1 have already been filled - FillPatchTwoLevels(temp_base_state, temp_base_state.nGrowVect(), IntVect(0,0,0), - time, cmf, ctime, fmf, ftime, - 0, 0, temp_base_state.nComp(), geom[lev-1], geom[lev], - refRatio(lev-1), mapper, domain_bcs_type, - BCVars::base_bc); + // Call FillPatch which ASSUMES that all ghost cells at lev-1 have already been filled + FillPatchTwoLevels(temp_base_state, temp_base_state.nGrowVect(), IntVect(0,0,0), + time, cmf, ctime, fmf, ftime, + 0, 0, temp_base_state.nComp(), geom[lev-1], geom[lev], + refRatio(lev-1), mapper, domain_bcs_type, + BCVars::base_bc); - // Impose bc's outside the domain - (*physbcs_base[lev])(temp_base_state,0,temp_base_state.nComp(),base_state[lev].nGrowVect()); + // Impose bc's outside the domain + (*physbcs_base[lev])(temp_base_state,0,temp_base_state.nComp(),base_state[lev].nGrowVect()); - std::swap(temp_base_state, base_state[lev]); - } + // ************************************************************************************************* + // This will fill the temporary MultiFabs with data from vars_new + // NOTE: the momenta here are only used as scratch space, the momenta themselves are not fillpatched + // NOTE: we must create the new base state before calling FillPatch because we will + // interpolate perturbational quantities + // ************************************************************************************************* + FillPatch(lev, time, {&temp_lev_new[Vars::cons],&temp_lev_new[Vars::xvel], + &temp_lev_new[Vars::yvel],&temp_lev_new[Vars::zvel]}, + {&temp_lev_new[Vars::cons],&rU_new[lev],&rV_new[lev],&rW_new[lev]}, + base_state[lev], temp_base_state, false); + + // Now swap the pointers since we needed both old and new in the FillPatch + std::swap(temp_base_state, base_state[lev]); // ******************************************************************************************** // Copy from new into old just in case diff --git a/Source/IO/ERF_Checkpoint.cpp b/Source/IO/ERF_Checkpoint.cpp index af211bf22..9137bace2 100644 --- a/Source/IO/ERF_Checkpoint.cpp +++ b/Source/IO/ERF_Checkpoint.cpp @@ -132,22 +132,22 @@ ERF::WriteCheckpointFile () const // Note that we write the ghost cells of the base state (unlike above) // For backward compatibility we only write the first components and 1 ghost cell - IntVect ng; int ncomp; + IntVect ng_base; int ncomp_base; bool write_old_base_state = true; if (write_old_base_state) { - ng = IntVect{1}; - ncomp = 3; + ng_base = IntVect{1}; + ncomp_base = 3; } else { - ng = base_state[lev].nGrowVect(); - ncomp = base_state[lev].nComp(); + ng_base = base_state[lev].nGrowVect(); + ncomp_base = base_state[lev].nComp(); } - MultiFab base(grids[lev],dmap[lev],ncomp,ng); - MultiFab::Copy(base,base_state[lev],0,0,ncomp,ng); + MultiFab base(grids[lev],dmap[lev],ncomp_base,ng_base); + MultiFab::Copy(base,base_state[lev],0,0,ncomp_base,ng_base); VisMF::Write(base, MultiFabFileFullPrefix(lev, checkpointname, "Level_", "BaseState")); if (solverChoice.use_terrain) { // Note that we also write the ghost cells of z_phys_nd - ng = z_phys_nd[lev]->nGrowVect(); + IntVect ng = z_phys_nd[lev]->nGrowVect(); MultiFab z_height(convert(grids[lev],IntVect(1,1,1)),dmap[lev],1,ng); MultiFab::Copy(z_height,*z_phys_nd[lev],0,0,1,ng); VisMF::Write(z_height, MultiFabFileFullPrefix(lev, checkpointname, "Level_", "Z_Phys_nd")); @@ -160,10 +160,10 @@ ERF::WriteCheckpointFile () const micro->Get_Qmoist_Restart_Vars(lev, solverChoice, qmoist_indices, qmoist_names); int qmoist_nvar = qmoist_indices.size(); for (int var = 0; var < qmoist_nvar; var++) { - ng = qmoist[lev][qmoist_indices[var]]->nGrowVect(); - ncomp = 1; - MultiFab moist_vars(grids[lev],dmap[lev],ncomp,ng); - MultiFab::Copy(moist_vars,*(qmoist[lev][qmoist_indices[var]]),0,0,ncomp,ng); + IntVect ng_moist = qmoist[lev][qmoist_indices[var]]->nGrowVect(); + const int ncomp = 1; + MultiFab moist_vars(grids[lev],dmap[lev],ncomp,ng_moist); + MultiFab::Copy(moist_vars,*(qmoist[lev][qmoist_indices[var]]),0,0,ncomp,ng_moist); VisMF::Write(moist_vars, amrex::MultiFabFileFullPrefix(lev, checkpointname, "Level_", qmoist_names[var])); } @@ -171,9 +171,9 @@ ERF::WriteCheckpointFile () const if(solverChoice.windfarm_type == WindFarmType::Fitch or solverChoice.windfarm_type == WindFarmType::EWP or solverChoice.windfarm_type == WindFarmType::SimpleAD){ - ng = Nturb[lev].nGrowVect(); - MultiFab mf_Nturb(grids[lev],dmap[lev],1,ng); - MultiFab::Copy(mf_Nturb,Nturb[lev],0,0,1,ng); + IntVect ng_turb = Nturb[lev].nGrowVect(); + MultiFab mf_Nturb(grids[lev],dmap[lev],1,ng_turb); + MultiFab::Copy(mf_Nturb,Nturb[lev],0,0,1,ng_turb); VisMF::Write(mf_Nturb, amrex::MultiFabFileFullPrefix(lev, checkpointname, "Level_", "NumTurb")); } #endif @@ -182,7 +182,7 @@ ERF::WriteCheckpointFile () const for (int mvar(0); mvarboxArray(); DistributionMapping dm = lsm_data[lev][mvar]->DistributionMap(); - ng = lsm_data[lev][mvar]->nGrowVect(); + IntVect ng = lsm_data[lev][mvar]->nGrowVect(); int nvar = lsm_data[lev][mvar]->nComp(); MultiFab lsm_vars(ba,dm,nvar,ng); MultiFab::Copy(lsm_vars,*(lsm_data[lev][mvar]),0,0,nvar,ng); @@ -197,7 +197,7 @@ ERF::WriteCheckpointFile () const } BoxArray ba2d(std::move(bl2d)); - ng = mapfac_m[lev]->nGrowVect(); + IntVect ng = mapfac_m[lev]->nGrowVect(); MultiFab mf_m(ba2d,dmap[lev],1,ng); MultiFab::Copy(mf_m,*mapfac_m[lev],0,0,1,ng); VisMF::Write(mf_m, MultiFabFileFullPrefix(lev, checkpointname, "Level_", "MapFactor_m")); @@ -414,18 +414,18 @@ ERF::ReadCheckpointFile () // The original base state only had 3 components and 1 ghost cell -- we read this // here to be consistent with the old style - IntVect ng; int ncomp; + IntVect ng_base; int ncomp_base; bool read_old_base_state = true; if (read_old_base_state) { - ng = IntVect{1}; - ncomp = 3; + ng_base = IntVect{1}; + ncomp_base = 3; } else { - ng = base_state[lev].nGrowVect(); - ncomp = base_state[lev].nComp(); + ng_base = base_state[lev].nGrowVect(); + ncomp_base = base_state[lev].nComp(); } - MultiFab base(grids[lev],dmap[lev],ncomp,ng); + MultiFab base(grids[lev],dmap[lev],ncomp_base,ng_base); VisMF::Read(base, MultiFabFileFullPrefix(lev, restart_chkfile, "Level_", "BaseState")); - MultiFab::Copy(base_state[lev],base,0,0,ncomp,ng); + MultiFab::Copy(base_state[lev],base,0,0,ncomp_base,ng_base); if (read_old_base_state) { for (MFIter mfi(base_state[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi) { @@ -442,7 +442,7 @@ ERF::ReadCheckpointFile () if (solverChoice.use_terrain) { // Note that we also read the ghost cells of z_phys_nd - ng = z_phys_nd[lev]->nGrowVect(); + IntVect ng = z_phys_nd[lev]->nGrowVect(); MultiFab z_height(convert(grids[lev],IntVect(1,1,1)),dmap[lev],1,ng); VisMF::Read(z_height, MultiFabFileFullPrefix(lev, restart_chkfile, "Level_", "Z_Phys_nd")); MultiFab::Copy(*z_phys_nd[lev],z_height,0,0,1,ng); @@ -455,11 +455,11 @@ ERF::ReadCheckpointFile () micro->Get_Qmoist_Restart_Vars(lev, solverChoice, qmoist_indices, qmoist_names); int qmoist_nvar = qmoist_indices.size(); for (int var = 0; var < qmoist_nvar; var++) { - ng = qmoist[lev][qmoist_indices[var]]->nGrowVect(); - ncomp = 1; - MultiFab moist_vars(grids[lev],dmap[lev],ncomp,ng); + IntVect ng_moist = qmoist[lev][qmoist_indices[var]]->nGrowVect(); + const int ncomp_moist = 1; + MultiFab moist_vars(grids[lev],dmap[lev],ncomp_moist,ng_moist); VisMF::Read(moist_vars, amrex::MultiFabFileFullPrefix(lev, restart_chkfile, "Level_", qmoist_names[var])); - MultiFab::Copy(*(qmoist[lev][qmoist_indices[var]]),moist_vars,0,0,ncomp,ng); + MultiFab::Copy(*(qmoist[lev][qmoist_indices[var]]),moist_vars,0,0,ncomp_moist,ng_moist); } #if defined(ERF_USE_WINDFARM) @@ -477,7 +477,7 @@ ERF::ReadCheckpointFile () for (int mvar(0); mvarboxArray(); DistributionMapping dm = lsm_data[lev][mvar]->DistributionMap(); - ng = lsm_data[lev][mvar]->nGrowVect(); + IntVect ng = lsm_data[lev][mvar]->nGrowVect(); int nvar = lsm_data[lev][mvar]->nComp(); MultiFab lsm_vars(ba,dm,nvar,ng); VisMF::Read(lsm_vars, MultiFabFileFullPrefix(lev, restart_chkfile, "Level_", "LsmVars")); @@ -492,7 +492,8 @@ ERF::ReadCheckpointFile () } BoxArray ba2d(std::move(bl2d)); - ng = mapfac_m[lev]->nGrowVect(); + { + IntVect ng = mapfac_m[lev]->nGrowVect(); MultiFab mf_m(ba2d,dmap[lev],1,ng); VisMF::Read(mf_m, MultiFabFileFullPrefix(lev, restart_chkfile, "Level_", "MapFactor_m")); MultiFab::Copy(*mapfac_m[lev],mf_m,0,0,1,ng); @@ -506,6 +507,7 @@ ERF::ReadCheckpointFile () MultiFab mf_v(convert(ba2d,IntVect(0,1,0)),dmap[lev],1,ng); VisMF::Read(mf_v, MultiFabFileFullPrefix(lev, restart_chkfile, "Level_", "MapFactor_v")); MultiFab::Copy(*mapfac_v[lev],mf_v,0,0,1,ng); + } } #ifdef ERF_USE_PARTICLES diff --git a/Source/IO/ERF_NCColumnFile.cpp b/Source/IO/ERF_NCColumnFile.cpp index eba2f2931..68720d353 100644 --- a/Source/IO/ERF_NCColumnFile.cpp +++ b/Source/IO/ERF_NCColumnFile.cpp @@ -121,10 +121,16 @@ ERF::writeToNCColumnFile (const int lev, // Need data in one grow cell for interpolation // Note that vars_new is what's filled here; rU_new/rV_new/rW_new are just used as scratch space - FillPatch(lev, t_new[lev], {&vars_new[lev][Vars::cons], &vars_new[lev][Vars::xvel], - &vars_new[lev][Vars::yvel], &vars_new[lev][Vars::zvel]}, - {&vars_new[lev][Vars::cons], &rU_new[lev], - &rV_new[lev], &rW_new[lev]}); + if (lev == 0) { + FillPatch(lev, t_new[lev], {&vars_new[lev][Vars::cons], &vars_new[lev][Vars::xvel], + &vars_new[lev][Vars::yvel], &vars_new[lev][Vars::zvel]}); + } else { + FillPatch(lev, t_new[lev], {&vars_new[lev][Vars::cons], &vars_new[lev][Vars::xvel], + &vars_new[lev][Vars::yvel], &vars_new[lev][Vars::zvel]}, + {&vars_new[lev][Vars::cons], &rU_new[lev], + &rV_new[lev], &rW_new[lev]}, + base_state[lev], base_state[lev]); + } MultiFab& S_new = vars_new[lev][Vars::cons]; MultiFab& U_new = vars_new[lev][Vars::xvel]; diff --git a/Source/IO/ERF_Plotfile.cpp b/Source/IO/ERF_Plotfile.cpp index ad08b384c..c8b7294e9 100644 --- a/Source/IO/ERF_Plotfile.cpp +++ b/Source/IO/ERF_Plotfile.cpp @@ -185,12 +185,6 @@ ERF::PlotFileVarNames (Vector plot_var_names ) void ERF::WritePlotFile (int which, PlotFileType plotfile_type, Vector plot_var_names) { - if (plotfile_type == PlotFileType::Amrex) { - amrex::Print() << "WRITE AMREX " << std::endl; - } else if (plotfile_type == PlotFileType::Netcdf) { - amrex::Print() << "WRITE NETCDF " << std::endl; - } - const Vector varnames = PlotFileVarNames(plot_var_names); const int ncomp_mf = varnames.size(); @@ -202,12 +196,18 @@ ERF::WritePlotFile (int which, PlotFileType plotfile_type, Vector p // which require ghost cells to be filled. We do not need to call FillPatcher // because we don't need to set interior fine points. // NOTE: the momenta here are only used as scratch space, the momenta themselves are not fillpatched - for (int lev = 0; lev <= finest_level; ++lev) { + + // Level 0 FilLPatch + FillPatch(0, t_new[0], {&vars_new[0][Vars::cons], &vars_new[0][Vars::xvel], + &vars_new[0][Vars::yvel], &vars_new[0][Vars::zvel]}); + + for (int lev = 1; lev <= finest_level; ++lev) { bool fillset = false; FillPatch(lev, t_new[lev], {&vars_new[lev][Vars::cons], &vars_new[lev][Vars::xvel], &vars_new[lev][Vars::yvel], &vars_new[lev][Vars::zvel]}, - {&vars_new[lev][Vars::cons], &rU_new[lev], - &rV_new[lev], &rW_new[lev]}, fillset); + {&vars_new[lev][Vars::cons], + &rU_new[lev], &rV_new[lev], &rW_new[lev]}, + base_state[lev], base_state[lev], fillset); } // Get qmoist pointers if using moisture diff --git a/Source/Initialization/ERF_init1d.cpp b/Source/Initialization/ERF_init1d.cpp index 1f6bc28f5..5bbe00d32 100644 --- a/Source/Initialization/ERF_init1d.cpp +++ b/Source/Initialization/ERF_init1d.cpp @@ -124,8 +124,6 @@ ERF::initHSE (int lev) // but inside the domain have already been filled in the call above to InterpFromCoarseLevel // (*physbcs_base[lev])(base_state[lev],0,base_state[lev].nComp(),base_state[lev].nGrowVect()); - - base_state[lev].FillBoundary(geom[lev].periodicity()); } void @@ -158,8 +156,6 @@ ERF::erf_enforce_hse (int lev, const auto geomdata = geom[lev].data(); const Real dz = geomdata.CellSize(2); - const Box& domain = geom[lev].Domain(); - for ( MFIter mfi(dens, TileNoZ()); mfi.isValid(); ++mfi ) { // Create a flat box with same horizontal extent but only one cell in vertical @@ -167,8 +163,10 @@ ERF::erf_enforce_hse (int lev, int klo = tbz.smallEnd(2); int khi = tbz.bigEnd(2); + // Note we only grow by 1 because that is how big z_cc is. Box b2d = tbz; // Copy constructor - b2d.grow(0,1); b2d.grow(1,1); // Grow by one in the lateral directions + b2d.grow(0,1); + b2d.grow(1,1); b2d.setRange(2,0); // We integrate to the first cell (and below) by using rho in this cell @@ -191,13 +189,11 @@ ERF::erf_enforce_hse (int lev, const Real rdOcp = solverChoice.rdOcp; - IntVect ngvect = base_state[lev].nGrowVect(); - int ng_z = ngvect[2]; - ParallelFor(b2d, [=] AMREX_GPU_DEVICE (int i, int j, int) { // Set value at surface from Newton iteration for rho - if (klo == 0) { + if (klo == 0) + { // Physical height of the terrain at cell center Real hz; if (l_use_terrain) { @@ -210,17 +206,16 @@ ERF::erf_enforce_hse (int lev, pi_arr(i,j,klo) = getExnergivenP(pres_arr(i,j,klo), rdOcp); th_arr(i,j,klo) =getRhoThetagivenP(pres_arr(i,j,klo)) / rho_arr(i,j,klo); + // // Set ghost cell with dz and rho at boundary + // (We will set the rest of the ghost cells in the boundary condition routine) + // pres_arr(i,j,klo-1) = p_0 + hz * rho_arr(i,j,klo) * l_gravity; pi_arr(i,j,klo-1) = getExnergivenP(pres_arr(i,j,klo-1), rdOcp); th_arr(i,j,klo-1) = getRhoThetagivenP(pres_arr(i,j,klo-1)) / rho_arr(i,j,klo-1); - for (int kk = 2; kk <= ng_z; kk++) { - pres_arr(i,j,klo-kk) = pres_arr(i,j,klo-1); - pi_arr(i,j,klo-kk) = pi_arr(i,j,klo-1); - th_arr(i,j,klo-kk) = th_arr(i,j,klo-1); - } } else { + // If level > 0 and klo > 0, we need to use the value of pres_arr(i,j,klo-1) which was // filled from FillPatch-ing it. Real dz_loc; @@ -236,11 +231,8 @@ ERF::erf_enforce_hse (int lev, pi_arr(i,j,klo ) = getExnergivenP(pres_arr(i,j,klo ), rdOcp); th_arr(i,j,klo ) = getRhoThetagivenP(pres_arr(i,j,klo )) / rho_arr(i,j,klo ); - for (int kk = 1; kk <= ng_z; kk++) { - pi_arr(i,j,klo-kk) = getExnergivenP(pres_arr(i,j,klo-kk), rdOcp); - th_arr(i,j,klo-kk) = getRhoThetagivenP(pres_arr(i,j,klo-kk)) / rho_arr(i,j,klo-kk); - } - + pi_arr(i,j,klo-1) = getExnergivenP(pres_arr(i,j,klo-1), rdOcp); + th_arr(i,j,klo-1) = getRhoThetagivenP(pres_arr(i,j,klo-1)) / rho_arr(i,j,klo-1); } Real dens_interp; @@ -262,68 +254,7 @@ ERF::erf_enforce_hse (int lev, } }); - bool is_periodic_in_x = geom[lev].isPeriodic(0); - bool is_periodic_in_y = geom[lev].isPeriodic(1); - - int domlo_x = domain.smallEnd(0); int domhi_x = domain.bigEnd(0); - int domlo_y = domain.smallEnd(1); int domhi_y = domain.bigEnd(1); - - if (!is_periodic_in_x) { - if (pres[mfi].box().smallEnd(0) < domlo_x) { - Box bx = mfi.nodaltilebox(2); - int ng = ngvect[0]; - bx.setSmall(0,domlo_x-1); - bx.setBig(0,domlo_x-ng); - ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) - { - pres_arr(i,j,k) = pres_arr(domlo_x,j,k); - pi_arr(i,j,k) = pi_arr(domlo_x,j,k); - th_arr(i,j,k) = th_arr(domlo_x,j,k); - }); - } - - if (pres[mfi].box().bigEnd(0) > domhi_x) { - Box bx = mfi.nodaltilebox(2); - int ng = ngvect[0]; - bx.setSmall(0,domhi_x+1); - bx.setBig(0,domhi_x+ng); - ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) - { - pres_arr(i,j,k) = pres_arr(domhi_x,j,k); - pi_arr(i,j,k) = pi_arr(domhi_x,j,k); - th_arr(i,j,k) = th_arr(domhi_x,j,k); - }); - } - } - - if (!is_periodic_in_y) { - if (pres[mfi].box().smallEnd(1) < domlo_y) { - Box bx = mfi.nodaltilebox(2); - int ng = ngvect[1]; - bx.setSmall(1,domlo_y-ng); - bx.setBig(1,domlo_y-1); - ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) - { - pres_arr(i,j,k) = pres_arr(i,domlo_y,k); - pi_arr(i,j,k) = pi_arr(i,domlo_y,k); - th_arr(i,j,k) = th_arr(i,domlo_y,k); - }); - } - - if (pres[mfi].box().bigEnd(1) > domhi_y) { - Box bx = mfi.nodaltilebox(2); - int ng = ngvect[1]; - bx.setSmall(1,domhi_y+1); - bx.setBig(1,domhi_y+ng); - ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) - { - pres_arr(i,j,k) = pres_arr(i,domhi_y,k); - pi_arr(i,j,k) = pi_arr(i,domhi_y,k); - th_arr(i,j,k) = th_arr(i,domhi_y,k); - }); - } - } - } + } // mfi dens.FillBoundary(geom[lev].periodicity()); pres.FillBoundary(geom[lev].periodicity()); diff --git a/Source/TimeIntegration/ERF_Advance.cpp b/Source/TimeIntegration/ERF_Advance.cpp index 87fc9ee79..2a8a73497 100644 --- a/Source/TimeIntegration/ERF_Advance.cpp +++ b/Source/TimeIntegration/ERF_Advance.cpp @@ -95,8 +95,14 @@ ERF::Advance (int lev, Real time, Real dt_lev, int iteration, int /*ncycle*/) // // NOTE: the momenta here are not fillpatched (they are only used as scratch space) // - FillPatch(lev, time, {&S_old, &U_old, &V_old, &W_old}, - {&S_old, &rU_old[lev], &rV_old[lev], &rW_old[lev]}); + if (lev == 0) { + FillPatch(lev, time, {&S_old, &U_old, &V_old, &W_old}); + } else { + FillPatch(lev, time, {&S_old, &U_old, &V_old, &W_old}, + {&S_old, &rU_old[lev], &rV_old[lev], &rW_old[lev]}, + base_state[lev], base_state[lev]); + } + // // So we must convert the fillpatched to momenta, including the ghost values // @@ -192,6 +198,20 @@ ERF::Advance (int lev, Real time, Real dt_lev, int iteration, int /*ncycle*/) evolveTracers( lev, dt_lev, vars_new, z_phys_nd ); #endif + // *********************************************************************************************** + // Impose domain boundary conditions here so that in FillPatching the fine data we won't + // need to re-fill these + // *********************************************************************************************** + if (lev < finest_level) { + IntVect ngvect_vels = vars_new[lev][Vars::xvel].nGrowVect(); + (*physbcs_cons[lev])(vars_new[lev][Vars::cons],0,vars_new[lev][Vars::cons].nComp(), + vars_new[lev][Vars::cons].nGrowVect(),time,BCVars::cons_bc,true); + (*physbcs_u[lev])(vars_new[lev][Vars::xvel],0,1,ngvect_vels,time,BCVars::xvel_bc,true); + (*physbcs_v[lev])(vars_new[lev][Vars::yvel],0,1,ngvect_vels,time,BCVars::yvel_bc,true); + (*physbcs_w[lev])(vars_new[lev][Vars::zvel], vars_new[lev][Vars::xvel], vars_new[lev][Vars::yvel], + ngvect_vels,time,BCVars::zvel_bc,true); + } + // ************************************************************************************** // Register old and new coarse data if we are at a level less than the finest level // ************************************************************************************** diff --git a/Source/Utils/ERF_AverageDown.cpp b/Source/Utils/ERF_AverageDown.cpp new file mode 100644 index 000000000..1741b011f --- /dev/null +++ b/Source/Utils/ERF_AverageDown.cpp @@ -0,0 +1,162 @@ +/** + * \file ERF_AverageDown.cpp + */ + +/** + * Main class in ERF code, instantiated from main.cpp +*/ + +#include +#include + +using namespace amrex; + +// Set covered coarse cells to be the average of overlying fine cells for all levels +void +ERF::AverageDown () +{ + AMREX_ALWAYS_ASSERT(solverChoice.coupling_type == CouplingType::TwoWay); + int src_comp = 0; + int num_comp = vars_new[0][Vars::cons].nComp(); + for (int lev = finest_level-1; lev >= 0; --lev) + { + AverageDownTo(lev,src_comp,num_comp); + } +} + +// Set covered coarse cells to be the average of overlying fine cells at level crse_lev +void +ERF::AverageDownTo (int crse_lev, int scomp, int ncomp) // NOLINT +{ + AMREX_ALWAYS_ASSERT(scomp == 0); + AMREX_ALWAYS_ASSERT(ncomp == vars_new[crse_lev][Vars::cons].nComp()); + AMREX_ALWAYS_ASSERT(solverChoice.coupling_type == CouplingType::TwoWay); + + // ****************************************************************************************** + // First do cell-centered quantities + // The quantity that is conserved is not (rho S), but rather (rho S / m^2) where + // m is the map scale factor at cell centers + // Here we pre-divide (rho S) by m^2 before average down + // ****************************************************************************************** + for (int lev = crse_lev; lev <= crse_lev+1; lev++) { + for (MFIter mfi(vars_new[lev][Vars::cons], TilingIfNotGPU()); mfi.isValid(); ++mfi) { + const Box& bx = mfi.tilebox(); + const Array4< Real> cons_arr = vars_new[lev][Vars::cons].array(mfi); + const Array4 mapfac_arr = mapfac_m[lev]->const_array(mfi); + if (solverChoice.use_terrain) { + const Array4 detJ_arr = detJ_cc[lev]->const_array(mfi); + ParallelFor(bx, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept + { + cons_arr(i,j,k,scomp+n) *= detJ_arr(i,j,k) / (mapfac_arr(i,j,0)*mapfac_arr(i,j,0)); + }); + } else { + ParallelFor(bx, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept + { + cons_arr(i,j,k,scomp+n) /= (mapfac_arr(i,j,0)*mapfac_arr(i,j,0)); + }); + } + } // mfi + } // lev + + int fine_lev = crse_lev+1; + + if (interpolation_type == StateInterpType::Perturbational) { + // Make the fine rho and (rho theta) be perturbational + MultiFab::Divide(vars_new[fine_lev][Vars::cons],vars_new[fine_lev][Vars::cons], + Rho_comp,RhoTheta_comp,1,IntVect{0}); + MultiFab::Subtract(vars_new[fine_lev][Vars::cons],base_state[fine_lev], + BaseState::r0_comp,Rho_comp,1,IntVect{0}); + MultiFab::Subtract(vars_new[fine_lev][Vars::cons],base_state[fine_lev], + BaseState::th0_comp,RhoTheta_comp,1,IntVect{0}); + + // Make the crse rho and (rho theta) be perturbational + MultiFab::Divide(vars_new[crse_lev][Vars::cons],vars_new[crse_lev][Vars::cons], + Rho_comp,RhoTheta_comp,1,IntVect{0}); + MultiFab::Subtract(vars_new[crse_lev][Vars::cons],base_state[crse_lev], + BaseState::r0_comp,Rho_comp,1,IntVect{0}); + MultiFab::Subtract(vars_new[crse_lev][Vars::cons],base_state[crse_lev], + BaseState::th0_comp,RhoTheta_comp,1,IntVect{0}); + } + + average_down(vars_new[crse_lev+1][Vars::cons],vars_new[crse_lev ][Vars::cons], + scomp, ncomp, refRatio(crse_lev)); + + if (interpolation_type == StateInterpType::Perturbational) { + // Restore the fine data to what it was + MultiFab::Add(vars_new[fine_lev][Vars::cons],base_state[fine_lev], + BaseState::r0_comp,Rho_comp,1,IntVect{0}); + MultiFab::Add(vars_new[fine_lev][Vars::cons],base_state[fine_lev], + BaseState::th0_comp,RhoTheta_comp,1,IntVect{0}); + MultiFab::Multiply(vars_new[fine_lev][Vars::cons],vars_new[fine_lev][Vars::cons], + Rho_comp,RhoTheta_comp,1,IntVect{0}); + + // Make the crse data be full values not perturbational + MultiFab::Add(vars_new[crse_lev][Vars::cons],base_state[crse_lev], + BaseState::r0_comp,Rho_comp,1,IntVect{0}); + MultiFab::Add(vars_new[crse_lev][Vars::cons],base_state[crse_lev], + BaseState::th0_comp,RhoTheta_comp,1,IntVect{0}); + MultiFab::Multiply(vars_new[crse_lev][Vars::cons],vars_new[crse_lev][Vars::cons], + Rho_comp,RhoTheta_comp,1,IntVect{0}); + } + + vars_new[crse_lev][Vars::cons].FillBoundary(geom[crse_lev].periodicity()); + + // Here we multiply (rho S) by m^2 after average down + for (int lev = crse_lev; lev <= crse_lev+1; lev++) { + for (MFIter mfi(vars_new[lev][Vars::cons], TilingIfNotGPU()); mfi.isValid(); ++mfi) { + const Box& bx = mfi.tilebox(); + const Array4< Real> cons_arr = vars_new[lev][Vars::cons].array(mfi); + const Array4 mapfac_arr = mapfac_m[lev]->const_array(mfi); + if (solverChoice.use_terrain) { + const Array4 detJ_arr = detJ_cc[lev]->const_array(mfi); + ParallelFor(bx, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept + { + cons_arr(i,j,k,scomp+n) *= (mapfac_arr(i,j,0)*mapfac_arr(i,j,0)) / detJ_arr(i,j,k); + }); + } else { + ParallelFor(bx, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept + { + cons_arr(i,j,k,scomp+n) *= (mapfac_arr(i,j,0)*mapfac_arr(i,j,0)); + }); + } + } // mfi + } // lev + + // ****************************************************************************************** + // Now average down momenta. + // Note that vars_new holds velocities not momenta, but we want to do conservative + // averaging so we first convert to momentum, then average down, then convert + // back to velocities -- only on the valid region + // ****************************************************************************************** + for (int lev = crse_lev; lev <= crse_lev+1; lev++) + { + // FillBoundary for density so we can go back and forth between velocity and momentum + vars_new[lev][Vars::cons].FillBoundary(geom[lev].periodicity()); + + VelocityToMomentum(vars_new[lev][Vars::xvel], IntVect(0,0,0), + vars_new[lev][Vars::yvel], IntVect(0,0,0), + vars_new[lev][Vars::zvel], IntVect(0,0,0), + vars_new[lev][Vars::cons], + rU_new[lev], + rV_new[lev], + rW_new[lev], + Geom(lev).Domain(), + domain_bcs_type); + } + + average_down_faces(rU_new[crse_lev+1], rU_new[crse_lev], refRatio(crse_lev), geom[crse_lev]); + average_down_faces(rV_new[crse_lev+1], rV_new[crse_lev], refRatio(crse_lev), geom[crse_lev]); + average_down_faces(rW_new[crse_lev+1], rW_new[crse_lev], refRatio(crse_lev), geom[crse_lev]); + + for (int lev = crse_lev; lev <= crse_lev+1; lev++) { + MomentumToVelocity(vars_new[lev][Vars::xvel], + vars_new[lev][Vars::yvel], + vars_new[lev][Vars::zvel], + vars_new[lev][Vars::cons], + rU_new[lev], + rV_new[lev], + rW_new[lev], + Geom(lev).Domain(), + domain_bcs_type); + } +} diff --git a/Source/Utils/Make.package b/Source/Utils/Make.package index f757c99f9..238eccb8c 100644 --- a/Source/Utils/Make.package +++ b/Source/Utils/Make.package @@ -16,6 +16,7 @@ CEXE_headers += ERF_Sat_methods.H CEXE_headers += ERF_Water_vapor_saturation.H CEXE_headers += ERF_DirectionSelector.H +CEXE_sources += ERF_AverageDown.cpp CEXE_sources += ERF_ChopGrids.cpp CEXE_sources += ERF_MomentumToVelocity.cpp CEXE_sources += ERF_VelocityToMomentum.cpp