From 88d2efd917e457227a5a614f3015336500148ea2 Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Mon, 14 Oct 2024 18:49:31 -0700 Subject: [PATCH] This version runs both compressible and anelastic with factor 2 refinement of the DC problem to 900s --- .../ERF_BoundaryConditions_cons.cpp | 20 ++--- Source/BoundaryConditions/ERF_FillPatch.cpp | 78 +++---------------- Source/BoundaryConditions/ERF_PhysBCFunct.H | 2 +- Source/BoundaryConditions/ERF_PhysBCFunct.cpp | 73 +++++++++++++---- Source/DataStructs/ERF_DataStruct.H | 42 +++++----- Source/ERF.H | 3 - Source/ERF.cpp | 15 +--- Source/ERF_make_new_level.cpp | 7 +- Source/IO/ERF_Plotfile.cpp | 9 ++- .../Microphysics/ERF_EulerianMicrophysics.H | 3 - Source/TimeIntegration/ERF_Advance.cpp | 18 +++-- Source/TimeIntegration/ERF_TI_utils.H | 19 ++++- Source/TimeIntegration/ERF_slow_rhs_pre.cpp | 5 +- 13 files changed, 141 insertions(+), 153 deletions(-) diff --git a/Source/BoundaryConditions/ERF_BoundaryConditions_cons.cpp b/Source/BoundaryConditions/ERF_BoundaryConditions_cons.cpp index b46c5f466..738d69fa5 100644 --- a/Source/BoundaryConditions/ERF_BoundaryConditions_cons.cpp +++ b/Source/BoundaryConditions/ERF_BoundaryConditions_cons.cpp @@ -15,7 +15,7 @@ using namespace amrex; */ void ERFPhysBCFunct_cons::impose_lateral_cons_bcs (const Array4& dest_arr, const Box& bx, const Box& domain, - int icomp, int ncomp, int ngz) + int icomp, int ncomp, IntVect ng) { BL_PROFILE_VAR("impose_lateral_cons_bcs()",impose_lateral_cons_bcs); const auto& dom_lo = lbound(domain); @@ -70,6 +70,7 @@ void ERFPhysBCFunct_cons::impose_lateral_cons_bcs (const Array4& dest_arr, { Box bx_xlo(bx); bx_xlo.setBig (0,dom_lo.x-1); Box bx_xhi(bx); bx_xhi.setSmall(0,dom_hi.x+1); + ParallelFor( bx_xlo, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) { @@ -108,6 +109,7 @@ void ERFPhysBCFunct_cons::impose_lateral_cons_bcs (const Array4& dest_arr, { Box bx_ylo(bx); bx_ylo.setBig (1,dom_lo.y-1); Box bx_yhi(bx); bx_yhi.setSmall(1,dom_hi.y+1); + ParallelFor( bx_ylo, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) { @@ -147,10 +149,10 @@ void ERFPhysBCFunct_cons::impose_lateral_cons_bcs (const Array4& dest_arr, // 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,ngz); - if (bx_xlo.bigEnd(2) != domain.bigEnd(2)) bx_xlo.growHi(2,ngz); - if (bx_xhi.smallEnd(2) != domain.smallEnd(2)) bx_xhi.growLo(2,ngz); - if (bx_xhi.bigEnd(2) != domain.bigEnd(2)) bx_xhi.growHi(2,ngz); + if (bx_xlo.smallEnd(2) != domain.smallEnd(2)) bx_xlo.growLo(2,ng[2]); + if (bx_xlo.bigEnd(2) != domain.bigEnd(2)) bx_xlo.growHi(2,ng[2]); + if (bx_xhi.smallEnd(2) != domain.smallEnd(2)) bx_xhi.growLo(2,ng[2]); + if (bx_xhi.bigEnd(2) != domain.bigEnd(2)) bx_xhi.growHi(2,ng[2]); ParallelFor( bx_xlo, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) { @@ -196,10 +198,10 @@ void ERFPhysBCFunct_cons::impose_lateral_cons_bcs (const Array4& dest_arr, // Populate ghost cells on lo-y and hi-y domain boundaries 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,ngz); - if (bx_ylo.bigEnd(2) != domain.bigEnd(2)) bx_ylo.growHi(2,ngz); - if (bx_yhi.smallEnd(2) != domain.smallEnd(2)) bx_yhi.growLo(2,ngz); - if (bx_yhi.bigEnd(2) != domain.bigEnd(2)) bx_yhi.growHi(2,ngz); + if (bx_ylo.smallEnd(2) != domain.smallEnd(2)) bx_ylo.growLo(2,ng[2]); + if (bx_ylo.bigEnd(2) != domain.bigEnd(2)) bx_ylo.growHi(2,ng[2]); + if (bx_yhi.smallEnd(2) != domain.smallEnd(2)) bx_yhi.growLo(2,ng[2]); + if (bx_yhi.bigEnd(2) != domain.bigEnd(2)) bx_yhi.growHi(2,ng[2]); ParallelFor( bx_ylo, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) { diff --git a/Source/BoundaryConditions/ERF_FillPatch.cpp b/Source/BoundaryConditions/ERF_FillPatch.cpp index 6836eba2e..352ee01a5 100644 --- a/Source/BoundaryConditions/ERF_FillPatch.cpp +++ b/Source/BoundaryConditions/ERF_FillPatch.cpp @@ -84,24 +84,18 @@ ERF::FillPatch (int lev, Real time, FillPatchSingleLevel(*mfs_vel[Vars::cons], ngvect_cons, time, fmf, IntVect(0,0,0), ftime, 0, 0, ncomp, geom[lev]); - (*physbcs_cons[lev])(*mfs_vel[Vars::cons],0,ncomp,ngvect_cons,time,BCVars::cons_bc); - 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]); - (*physbcs_u[lev])(*mfs_vel[Vars::xvel],0,1,ngvect_vels,time,BCVars::xvel_bc); 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]); - (*physbcs_v[lev])(*mfs_vel[Vars::yvel],0,1,ngvect_vels,time,BCVars::xvel_bc); 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]); - (*physbcs_w[lev])(*mfs_vel[Vars::zvel],*mfs_vel[Vars::xvel],*mfs_vel[Vars::yvel], - ngvect_vels,time,BCVars::zvel_bc); } // !cons_only } else { @@ -115,13 +109,10 @@ ERF::FillPatch (int lev, Real time, 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 (*physbcs_cons[lev-1])(vars_old[lev-1][Vars::cons],0,mf_c.nComp(),ngvect_cons,time,BCVars::cons_bc); (*physbcs_cons[lev-1])(vars_new[lev-1][Vars::cons],0,mf_c.nComp(),ngvect_cons,time,BCVars::cons_bc); - // Make sure internal ghost cells are filled as well - vars_old[lev-1][Vars::cons].FillBoundary(geom[lev-1].periodicity()); - vars_new[lev-1][Vars::cons].FillBoundary(geom[lev-1].periodicity()); - // Call FillPatchTwoLevels which ASSUMES that all ghost cells have already been filled FillPatchTwoLevels(mf_c, ngvect_cons, IntVect(0,0,0), time, cmf, ctime, fmf, ftime, @@ -129,9 +120,6 @@ ERF::FillPatch (int lev, Real time, refRatio(lev-1), mapper, domain_bcs_type, BCVars::cons_bc); - // Impose physical bc's on fine data - (*physbcs_cons[lev])(mf_c,0,mf_c.nComp(),ngvect_cons,time,BCVars::cons_bc); - if (!cons_only) { mapper = &face_cons_linear_interp; @@ -143,13 +131,10 @@ ERF::FillPatch (int lev, Real time, // ********************************************************************** // Impose physical bc's on coarse data (note time and 0 are not used) + // Note that we call FillBoundary inside the physbcs call (*physbcs_u[lev-1])(vars_old[lev-1][Vars::xvel],0,1,ngvect_vels,time,BCVars::xvel_bc); (*physbcs_u[lev-1])(vars_new[lev-1][Vars::xvel],0,1,ngvect_vels,time,BCVars::xvel_bc); - // Make sure internal ghost cells are filled as well - vars_old[lev-1][Vars::xvel].FillBoundary(geom[lev-1].periodicity()); - vars_new[lev-1][Vars::xvel].FillBoundary(geom[lev-1].periodicity()); - fmf = {&vars_old[lev ][Vars::xvel], &vars_new[lev ][Vars::xvel]}; cmf = {&vars_old[lev-1][Vars::xvel], &vars_new[lev-1][Vars::xvel]}; @@ -160,19 +145,13 @@ ERF::FillPatch (int lev, Real time, refRatio(lev-1), mapper, domain_bcs_type, BCVars::xvel_bc); - // Impose physical bc's on fine data - (*physbcs_u[lev])(vars_new[lev][Vars::xvel],0,mf_u.nComp(),ngvect_vels,time,BCVars::xvel_bc); - // ********************************************************************** // Impose physical bc's on coarse data (note time and 0 are not used) + // Note that we call FillBoundary inside the physbcs call (*physbcs_v[lev-1])(vars_old[lev-1][Vars::yvel],0,1,ngvect_vels,time,BCVars::yvel_bc); (*physbcs_v[lev-1])(vars_new[lev-1][Vars::yvel],0,1,ngvect_vels,time,BCVars::yvel_bc); - // Make sure internal ghost cells are filled as well - vars_old[lev-1][Vars::yvel].FillBoundary(geom[lev-1].periodicity()); - vars_new[lev-1][Vars::yvel].FillBoundary(geom[lev-1].periodicity()); - fmf = {&vars_old[lev ][Vars::yvel], &vars_new[lev ][Vars::yvel]}; cmf = {&vars_old[lev-1][Vars::yvel], &vars_new[lev-1][Vars::yvel]}; @@ -183,12 +162,10 @@ ERF::FillPatch (int lev, Real time, refRatio(lev-1), mapper, domain_bcs_type, BCVars::yvel_bc); - // Impose physical bc's on fine data - (*physbcs_v[lev])(vars_new[lev][Vars::yvel],0,1,ngvect_vels,time,BCVars::yvel_bc); - // ********************************************************************** // Impose physical bc's on coarse data (note time and 0 are not used) + // Note that we call FillBoundary inside the physbcs call (*physbcs_w[lev-1])(vars_old[lev-1][Vars::zvel], vars_old[lev-1][Vars::xvel], vars_old[lev-1][Vars::yvel], @@ -198,10 +175,6 @@ ERF::FillPatch (int lev, Real time, vars_new[lev-1][Vars::yvel], ngvect_vels,time,BCVars::zvel_bc); - // Make sure internal ghost cells are filled as well - vars_old[lev-1][Vars::zvel].FillBoundary(geom[lev-1].periodicity()); - vars_new[lev-1][Vars::zvel].FillBoundary(geom[lev-1].periodicity()); - fmf = {&vars_old[lev ][Vars::zvel], &vars_new[lev ][Vars::zvel]}; cmf = {&vars_old[lev-1][Vars::zvel], &vars_new[lev-1][Vars::zvel]}; @@ -211,12 +184,6 @@ ERF::FillPatch (int lev, Real time, 0, 0, 1, geom[lev-1], geom[lev], refRatio(lev-1), mapper, domain_bcs_type, BCVars::zvel_bc); - - - // Impose physical bc's on fine data -- note the u and v have been filled above - (*physbcs_w[lev])(*mfs_vel[Vars::zvel],*mfs_vel[Vars::xvel],*mfs_vel[Vars::yvel], - ngvect_vels,time,BCVars::zvel_bc); - } // !cons_only } // lev > 0 @@ -236,6 +203,7 @@ ERF::FillPatch (int lev, Real time, if (m_r2d) fill_from_bndryregs(mfs_vel,time); // We call these even if init_type == 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); if (!cons_only) { (*physbcs_u[lev])(*mfs_vel[Vars::xvel],0,1,ngvect_vels,time,BCVars::xvel_bc); @@ -245,37 +213,6 @@ ERF::FillPatch (int lev, Real time, } } -/* - * Fill ghost cells of qmoist - * - * @param[in] lev level of refinement at which to fill the data - * @param[in] time time at which the data should be filled - * @param[out] mf MultiFab to be filled (qmoist[lev]) - */ -void -ERF::FillPatchMoistVars (int lev, MultiFab& mf) -{ - BL_PROFILE_VAR("ERF::FillPatchMoistVars()",ERF_FillPatchMoistVars); - // *************************************************************************** - // Physical bc's at domain boundary - // *************************************************************************** - int icomp_cons = 0; - int ncomp_cons = 1; // We only fill qv, the first component - - // Note that we are filling qv, stored in qmoist[lev], with the input data (if there is any), stored - // in RhoQ1_comp. - - if (!use_real_bcs) { - Real time = Real(0.0); - IntVect ngvect_cons = mf.nGrowVect(); - int bccomp_cons = BCVars::RhoQ1_bc_comp; - - (*physbcs_cons[lev])(mf,icomp_cons,ncomp_cons,ngvect_cons,time,bccomp_cons); - } - - mf.FillBoundary(geom[lev].periodicity()); -} - /* * Fill valid and ghost data * This version fills mfs in valid regions with the values in "mfs" when it is passed in; @@ -549,7 +486,10 @@ ERF::FillIntermediatePatch (int lev, Real time, domain_bcs_type); } - mfs_mom[Vars::cons]->FillBoundary(geom[lev].periodicity()); + mfs_mom[IntVars::cons]->FillBoundary(geom[lev].periodicity()); + mfs_mom[IntVars::xmom]->FillBoundary(geom[lev].periodicity()); + mfs_mom[IntVars::ymom]->FillBoundary(geom[lev].periodicity()); + mfs_mom[IntVars::zmom]->FillBoundary(geom[lev].periodicity()); } /* diff --git a/Source/BoundaryConditions/ERF_PhysBCFunct.H b/Source/BoundaryConditions/ERF_PhysBCFunct.H index ea5ae9b00..ff09ff93a 100644 --- a/Source/BoundaryConditions/ERF_PhysBCFunct.H +++ b/Source/BoundaryConditions/ERF_PhysBCFunct.H @@ -56,7 +56,7 @@ public: void impose_lateral_cons_bcs (const amrex::Array4& dest_arr, const amrex::Box& bx, const amrex::Box& domain, - int icomp, int ncomp, int ngz); + int icomp, int ncomp, amrex::IntVect ng); void impose_vertical_cons_bcs (const amrex::Array4& dest_arr, const amrex::Box& bx, const amrex::Box& domain, const amrex::Array4& z_nd, diff --git a/Source/BoundaryConditions/ERF_PhysBCFunct.cpp b/Source/BoundaryConditions/ERF_PhysBCFunct.cpp index 06becebfc..a5c6048ca 100644 --- a/Source/BoundaryConditions/ERF_PhysBCFunct.cpp +++ b/Source/BoundaryConditions/ERF_PhysBCFunct.cpp @@ -46,6 +46,12 @@ void ERFPhysBCFunct_cons::operator() (MultiFab& mf, int icomp, int ncomp, z_nd_mf_loc.nGrowVect()); } + // + // 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 @@ -77,9 +83,11 @@ void ERFPhysBCFunct_cons::operator() (MultiFab& mf, int icomp, int ncomp, if (!m_use_real_bcs) { - impose_lateral_cons_bcs(cons_arr,cbx1,domain,icomp,ncomp,nghost[2]); + // We send a box with ghost cells in the lateral directions only + impose_lateral_cons_bcs(cons_arr,cbx1,domain,icomp,ncomp,nghost); } + // We send the full FAB box with ghost cells impose_vertical_cons_bcs(cons_arr,cbx2,domain,z_nd_arr,dxInv,icomp,ncomp); } @@ -118,6 +126,12 @@ void ERFPhysBCFunct_u::operator() (MultiFab& mf, int /*icomp*/, int /*ncomp*/, z_nd_mf_loc.nGrowVect()); } + // + // 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 @@ -133,7 +147,10 @@ void ERFPhysBCFunct_u::operator() (MultiFab& mf, int /*icomp*/, int /*ncomp*/, // // These are the boxes we use to test on relative to the domain // - Box xbx1 = surroundingNodes(bx,0); xbx1.grow(IntVect(nghost[0],nghost[1],0)); + Box xbx1 = surroundingNodes(bx,0); xbx1.grow(nghost); + if(xbx1.smallEnd(2) < domain.smallEnd(2)) xbx1.setSmall(2,domain.smallEnd(2)); + if(xbx1.bigEnd(2) > domain.bigEnd(2)) xbx1.setBig(2,domain.bigEnd(2)); + Box xbx2 = surroundingNodes(bx,0); xbx2.grow(nghost); Array4 z_nd_arr; @@ -157,7 +174,6 @@ void ERFPhysBCFunct_u::operator() (MultiFab& mf, int /*icomp*/, int /*ncomp*/, impose_vertical_xvel_bcs(velx_arr,xbx2,domain,z_nd_arr,dxInv,bccomp,time); } - } // MFIter } // OpenMP } // operator() @@ -193,6 +209,12 @@ void ERFPhysBCFunct_v::operator() (MultiFab& mf, int /*icomp*/, int /*ncomp*/, z_nd_mf_loc.nGrowVect()); } + // + // 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 @@ -208,7 +230,10 @@ void ERFPhysBCFunct_v::operator() (MultiFab& mf, int /*icomp*/, int /*ncomp*/, // // These are the boxes we use to test on relative to the domain // - Box ybx1 = surroundingNodes(bx,1); ybx1.grow(IntVect(nghost[0],nghost[1],0)); + Box ybx1 = surroundingNodes(bx,1); ybx1.grow(nghost); + if (ybx1.smallEnd(2) < domain.smallEnd(2)) ybx1.setSmall(2,domain.smallEnd(2)); + if (ybx1.bigEnd(2) > domain.bigEnd(2)) ybx1.setBig(2,domain.bigEnd(2)); + Box ybx2 = surroundingNodes(bx,1); ybx2.grow(nghost); Array4 z_nd_arr; @@ -272,6 +297,12 @@ void ERFPhysBCFunct_w::operator() (MultiFab& mf, MultiFab& xvel, MultiFab& yvel, } z_nd_mf_loc.FillBoundary(m_geom.periodicity()); + // + // 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 @@ -287,8 +318,17 @@ void ERFPhysBCFunct_w::operator() (MultiFab& mf, MultiFab& xvel, MultiFab& yvel, // // These are the boxes we use to test on relative to the domain // - Box zbx = surroundingNodes(bx,2); zbx.grow(0,nghost[0]); - zbx.grow(1,nghost[1]); + Box zbx1 = surroundingNodes(bx,2); zbx1.grow(nghost); + if (zbx1.smallEnd(2) < domain.smallEnd(2)) zbx1.setSmall(2,domain.smallEnd(2)); + if (zbx1.bigEnd(2) > domain.bigEnd(2)) zbx1.setBig(2,domain.bigEnd(2)+1); + + Box zbx2 = surroundingNodes(bx,2); zbx2.grow(nghost); + + // + // These are the boxes we use to test on relative to the domain + // + // Box zbx = surroundingNodes(bx,2); zbx.grow(0,nghost[0]); + // zbx.grow(1,nghost[1]); Array4 z_nd_arr; if (m_z_phys_nd) @@ -296,20 +336,21 @@ void ERFPhysBCFunct_w::operator() (MultiFab& mf, MultiFab& xvel, MultiFab& yvel, z_nd_arr = z_nd_mf_loc.const_array(mfi); } - Array4 const& velx_arr = xvel.const_array(mfi); - Array4 const& vely_arr = yvel.const_array(mfi); - Array4< Real> const& velz_arr = mf.array(mfi); - - if (!m_use_real_bcs) + if (!gdomainz.contains(zbx2)) { - if (!gdomainz.contains(zbx)) + Array4 const& velx_arr = xvel.const_array(mfi); + Array4 const& vely_arr = yvel.const_array(mfi); + Array4< Real> const& velz_arr = mf.array(mfi); + + if (!m_use_real_bcs) { - impose_lateral_zvel_bcs(velz_arr,velx_arr,vely_arr,zbx,domain,z_nd_arr,dxInv,bccomp_w); + if (!gdomainz.contains(zbx1)) + { + impose_lateral_zvel_bcs(velz_arr,velx_arr,vely_arr,zbx1,domain,z_nd_arr,dxInv,bccomp_w); + } } - } // m_use_real_bcs - if (!gdomainz.contains(zbx)) { - impose_vertical_zvel_bcs(velz_arr,velx_arr,vely_arr,zbx,domain,z_nd_arr,dxInv, + impose_vertical_zvel_bcs(velz_arr,velx_arr,vely_arr,zbx2,domain,z_nd_arr,dxInv, bccomp_u, bccomp_v, bccomp_w, m_terrain_type); } } // MFIter diff --git a/Source/DataStructs/ERF_DataStruct.H b/Source/DataStructs/ERF_DataStruct.H index 9e7cc1470..5a94ec286 100644 --- a/Source/DataStructs/ERF_DataStruct.H +++ b/Source/DataStructs/ERF_DataStruct.H @@ -430,35 +430,35 @@ struct SolverChoice { amrex::Print() << "use_coriolis : " << use_coriolis << std::endl; amrex::Print() << "use_gravity : " << use_gravity << std::endl; - if (coupling_type == CouplingType::TwoWay) { - amrex::Print() << "Using two-way coupling " << std::endl; - } else if (coupling_type == CouplingType::OneWay) { - amrex::Print() << "Using one-way coupling " << std::endl; - } - + amrex::Print() << "Terrain Type: " << std::endl; if (terrain_type == TerrainType::Static) { - amrex::Print() << "Using static terrain " << std::endl; + amrex::Print() << " Static" << std::endl; } else if (terrain_type == TerrainType::Moving) { - amrex::Print() << "Using moving terrain " << std::endl; + amrex::Print() << " Moving" << std::endl; } else { - amrex::Print() << "No terrain " << std::endl; + amrex::Print() << " None" << std::endl; } + amrex::Print() << "ABL Driver Type: " << std::endl; if (abl_driver_type == ABLDriverType::None) { - amrex::Print() << "ABL Driver Type: " << "None" << std::endl; - amrex::Print() << "No ABL driver selected " << std::endl; + amrex::Print() << " None" << std::endl; } else if (abl_driver_type == ABLDriverType::PressureGradient) { - amrex::Print() << "ABL Driver Type: " << "PressureGradient" << std::endl; - amrex::Print() << "Driving abl_pressure_grad: ("; - for (int i = 0; i < AMREX_SPACEDIM; ++i) - amrex::Print() << abl_pressure_grad[i] << " "; - amrex::Print() << ")" << std::endl; + amrex::Print() << " Pressure Gradient " + << amrex::RealVect(abl_pressure_grad[0],abl_pressure_grad[1],abl_pressure_grad[2]) + << std::endl; } else if (abl_driver_type == ABLDriverType::GeostrophicWind) { - amrex::Print() << "ABL Driver Type: " << "GeostrophicWind" << std::endl; - amrex::Print() << "Driving abl_geo_forcing: ("; - for (int i = 0; i < AMREX_SPACEDIM; ++i) - amrex::Print() << abl_geo_forcing[i] << " "; - amrex::Print() << ")" << std::endl; + amrex::Print() << " Geostrophic Wind " + << amrex::RealVect(abl_geo_forcing[0],abl_geo_forcing[1],abl_geo_forcing[2]) + << std::endl; + } + + if (max_level > 0) { + amrex::Print() << "Coupling Type: " << std::endl; + if (coupling_type == CouplingType::TwoWay) { + amrex::Print() << " Two-way" << std::endl; + } else if (coupling_type == CouplingType::OneWay) { + amrex::Print() << " One-way" << std::endl; + } } amrex::Print() << "Buoyancy_type : " << buoyancy_type << std::endl; diff --git a/Source/ERF.H b/Source/ERF.H index 0a8230110..e9ed7ad52 100644 --- a/Source/ERF.H +++ b/Source/ERF.H @@ -492,9 +492,6 @@ private: const amrex::Vector& mfs_mom, bool fillset=true, bool cons_only=false); - // Compute a new MultiFab by copying from valid region and filling ghost cells - - void FillPatchMoistVars (int lev, amrex::MultiFab& mf); - // Compute new multifabs by copying data from valid region and filling ghost cells. // Unlike FillPatch, FillIntermediatePatch will use the supplied multifabs instead of fine level data. // This is to support filling boundary cells at an intermediate time between old/new times diff --git a/Source/ERF.cpp b/Source/ERF.cpp index db10aefb4..427582722 100644 --- a/Source/ERF.cpp +++ b/Source/ERF.cpp @@ -719,17 +719,6 @@ ERF::InitData_post () for (int lev(0); lev <= max_level; ++lev) { make_physbcs(lev); } - - // TODO: Check if this is needed. I don't think it is since we now - // advect all the scalars... - - // Need to fill ghost cells here since we will use this qmoist in advance - if (solverChoice.moisture_type != MoistureType::None) - { - for (int lev = 0; lev <= finest_level; lev++) { - if (qmoist[lev].size() > 0) FillPatchMoistVars(lev, *(qmoist[lev][0])); // qv component - } - } } #ifdef ERF_USE_PARTICLES @@ -1697,7 +1686,6 @@ ERF::MakeHorizontalAverages () auto fab_arr = mf.array(mfi); auto const hse_arr = base_state[lev].const_array(mfi); auto const cons_arr = vars_new[lev][Vars::cons].const_array(mfi); - auto const qv_arr = qmoist[lev][0]->const_array(mfi); int ncomp = vars_new[lev][Vars::cons].nComp(); ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) { @@ -1705,7 +1693,8 @@ ERF::MakeHorizontalAverages () if (is_anelastic) { fab_arr(i,j,k,2) = hse_arr(i,j,k,1); } else { - fab_arr(i, j, k, 2) = getPgivenRTh(cons_arr(i, j, k, RhoTheta_comp), qv_arr(i,j,k)); + Real qv = cons_arr(i, j, k, RhoQ1_comp) / dens; + fab_arr(i, j, k, 2) = getPgivenRTh(cons_arr(i, j, k, RhoTheta_comp), qv); } fab_arr(i, j, k, 3) = (ncomp > RhoQ1_comp ? cons_arr(i, j, k, RhoQ1_comp) / dens : 0.0); fab_arr(i, j, k, 4) = (ncomp > RhoQ2_comp ? cons_arr(i, j, k, RhoQ2_comp) / dens : 0.0); diff --git a/Source/ERF_make_new_level.cpp b/Source/ERF_make_new_level.cpp index 52d59ba27..7fbcae895 100644 --- a/Source/ERF_make_new_level.cpp +++ b/Source/ERF_make_new_level.cpp @@ -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,9 +332,10 @@ 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]}, diff --git a/Source/IO/ERF_Plotfile.cpp b/Source/IO/ERF_Plotfile.cpp index 5bdcf5831..e26ee2ec3 100644 --- a/Source/IO/ERF_Plotfile.cpp +++ b/Source/IO/ERF_Plotfile.cpp @@ -199,6 +199,7 @@ ERF::WritePlotFile (int which, Vector plot_var_names) // We Fillpatch here because some of the derived quantities require derivatives // 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) { bool fillset = false; FillPatch(lev, t_new[lev], {&vars_new[lev][Vars::cons], &vars_new[lev][Vars::xvel], @@ -1511,10 +1512,10 @@ ERF::WriteMultiLevelPlotfileWithTerrain (const std::string& plotfilename, int nl const Vector& mf, const Vector& mf_nd, const Vector& varnames, - const Vector& geom, + const Vector& my_geom, Real time, const Vector& level_steps, - const Vector& ref_ratio, + const Vector& rr, const std::string &versionName, const std::string &levelPrefix, const std::string &mfPrefix, @@ -1523,7 +1524,7 @@ ERF::WriteMultiLevelPlotfileWithTerrain (const std::string& plotfilename, int nl BL_PROFILE("WriteMultiLevelPlotfileWithTerrain()"); AMREX_ALWAYS_ASSERT(nlevels <= mf.size()); - AMREX_ALWAYS_ASSERT(nlevels <= ref_ratio.size()+1); + AMREX_ALWAYS_ASSERT(nlevels <= rr.size()+1); AMREX_ALWAYS_ASSERT(nlevels <= level_steps.size()); AMREX_ALWAYS_ASSERT(mf[0]->nComp() == varnames.size()); @@ -1553,7 +1554,7 @@ ERF::WriteMultiLevelPlotfileWithTerrain (const std::string& plotfilename, int nl std::ofstream::binary); if( ! HeaderFile.good()) FileOpenFailed(HeaderFileName); WriteGenericPlotfileHeaderWithTerrain(HeaderFile, nlevels, boxArrays, varnames, - geom, time, level_steps, ref_ratio, versionName, + my_geom, time, level_steps, rr, versionName, levelPrefix, mfPrefix); }; diff --git a/Source/Microphysics/ERF_EulerianMicrophysics.H b/Source/Microphysics/ERF_EulerianMicrophysics.H index 9892a30fe..cb91854b4 100644 --- a/Source/Microphysics/ERF_EulerianMicrophysics.H +++ b/Source/Microphysics/ERF_EulerianMicrophysics.H @@ -31,14 +31,11 @@ public: a_model_type == MoistureType::SAM_NoIce || a_model_type == MoistureType::SAM_NoPrecip_NoIce) { SetModel(); - amrex::Print() << "SAM moisture model!\n"; } else if (a_model_type == MoistureType::Kessler || a_model_type == MoistureType::Kessler_NoRain) { SetModel(); - amrex::Print() << "Kessler moisture model!\n"; } else if (a_model_type == MoistureType::None) { SetModel(); - amrex::Print() << "No moisture model!\n"; } else { amrex::Abort("EulerianMicrophysics: Dont know this moisture_type!") ; } diff --git a/Source/TimeIntegration/ERF_Advance.cpp b/Source/TimeIntegration/ERF_Advance.cpp index c5c4ca369..db7cbf1a0 100644 --- a/Source/TimeIntegration/ERF_Advance.cpp +++ b/Source/TimeIntegration/ERF_Advance.cpp @@ -92,14 +92,20 @@ ERF::Advance (int lev, Real time, Real dt_lev, int iteration, int /*ncycle*/) V_new.setVal(1.e34,V_new.nGrowVect()); W_new.setVal(1.e34,W_new.nGrowVect()); + // + // 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 (solverChoice.moisture_type != MoistureType::None) { - // TODO: This is only qv - if (qmoist[lev].size() > 0) FillPatchMoistVars(lev, *(qmoist[lev][0])); - } - + // + // So we must convert the fillpatched to momenta, including the ghost values + // + VelocityToMomentum(U_old, rU_old[lev].nGrowVect(), + V_old, rV_old[lev].nGrowVect(), + W_old, rW_old[lev].nGrowVect(), + S_old, rU_old[lev], rV_old[lev], rW_old[lev], + Geom(lev).Domain(), + domain_bcs_type); #if defined(ERF_USE_WINDFARM) if (solverChoice.windfarm_type != WindFarmType::None) { diff --git a/Source/TimeIntegration/ERF_TI_utils.H b/Source/TimeIntegration/ERF_TI_utils.H index 21c330bc9..7870508a4 100644 --- a/Source/TimeIntegration/ERF_TI_utils.H +++ b/Source/TimeIntegration/ERF_TI_utils.H @@ -15,19 +15,32 @@ const Box& gbx = mfi.growntilebox(ng); const Array4& cons_arr = cons_state.array(mfi); const Array4< Real>& prim_arr = S_prim.array(mfi); - const Array4< Real>& pi_stage_arr = pi_stage.array(mfi); - const Real rdOcp = solverChoice.rdOcp; + // + // We may need > one ghost cells of prim in order to compute higher order advective terms + // amrex::ParallelFor(gbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { Real rho = cons_arr(i,j,k,Rho_comp); Real rho_theta = cons_arr(i,j,k,RhoTheta_comp); prim_arr(i,j,k,PrimTheta_comp) = rho_theta / rho; - pi_stage_arr(i,j,k) = getExnergivenRTh(rho_theta, rdOcp); for (int n = 1; n < ncomp_prim; ++n) { prim_arr(i,j,k,PrimTheta_comp + n) = cons_arr(i,j,k,RhoTheta_comp + n) / rho; } }); + + // + // We only use one ghost cell of pi_stage so we only fill one here + // + const Box& gbx1 = mfi.growntilebox(1); + + const Array4< Real>& pi_stage_arr = pi_stage.array(mfi); + const Real rdOcp = solverChoice.rdOcp; + + amrex::ParallelFor(gbx1, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + pi_stage_arr(i,j,k) = getExnergivenRTh(cons_arr(i,j,k,RhoTheta_comp), rdOcp); + }); } // mfi }; diff --git a/Source/TimeIntegration/ERF_slow_rhs_pre.cpp b/Source/TimeIntegration/ERF_slow_rhs_pre.cpp index 7cb16958b..7ec8b392c 100644 --- a/Source/TimeIntegration/ERF_slow_rhs_pre.cpp +++ b/Source/TimeIntegration/ERF_slow_rhs_pre.cpp @@ -373,8 +373,9 @@ void erf_slow_rhs_pre (int level, int finest_level, // Now create Omega with momentum (not velocity) with z_t subtracted if moving terrain if (l_use_terrain) { - Box gbxo_lo = gbxo; gbxo_lo.setBig(2,0); - if (gbxo_lo.smallEnd(2) <= 0) { + Box gbxo_lo = gbxo; gbxo_lo.setBig(2,domain.smallEnd(2)); + int lo_z_face = domain.smallEnd(2); + if (gbxo_lo.smallEnd(2) <= lo_z_face) { ParallelFor(gbxo_lo, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { omega_arr(i,j,k) = 0.; });