fix bcs with terrain and clean up time integrator

CMake/BuildERFExe.cmake

@@ -203,10 +203,10 @@ function(build_erf_lib erf_lib_name)
        ${SRC_DIR}/Microphysics/Kessler/ERF_Init_Kessler.cpp
        ${SRC_DIR}/Microphysics/Kessler/ERF_Kessler.cpp
        ${SRC_DIR}/Microphysics/Kessler/ERF_Update_Kessler.cpp
-	   ${SRC_DIR}/WindFarmParametrization/Fitch/ERF_AdvanceFitch.cpp
-	   ${SRC_DIR}/WindFarmParametrization/EWP/ERF_AdvanceEWP.cpp
-	   ${SRC_DIR}/WindFarmParametrization/SimpleActuatorDisk/ERF_AdvanceSimpleAD.cpp
-	   ${SRC_DIR}/WindFarmParametrization/GeneralActuatorDisk/ERF_AdvanceGeneralAD.cpp
+       ${SRC_DIR}/WindFarmParametrization/Fitch/ERF_AdvanceFitch.cpp
+       ${SRC_DIR}/WindFarmParametrization/EWP/ERF_AdvanceEWP.cpp
+       ${SRC_DIR}/WindFarmParametrization/SimpleActuatorDisk/ERF_AdvanceSimpleAD.cpp
+       ${SRC_DIR}/WindFarmParametrization/GeneralActuatorDisk/ERF_AdvanceGeneralAD.cpp
        ${SRC_DIR}/LandSurfaceModel/SLM/ERF_SLM.cpp
        ${SRC_DIR}/LandSurfaceModel/MM5/ERF_MM5.cpp
   )

Source/BoundaryConditions/ERF_BoundaryConditions_basestate.cpp

@@ -259,7 +259,7 @@ void ERFPhysBCFunct_base::impose_vertical_basestate_bcs (const Array4<Real>& des
         {
             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::r0_comp) * hz * CONST_GRAV;
             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);
         }

Source/BoundaryConditions/ERF_BoundaryConditions_cons.cpp

@@ -261,7 +261,7 @@ void ERFPhysBCFunct_cons::impose_lateral_cons_bcs (const Array4<Real>& dest_arr,
 void ERFPhysBCFunct_cons::impose_vertical_cons_bcs (const Array4<Real>& dest_arr, const Box& bx, const Box& domain,
                                                     const Array4<Real const>& z_phys_nd,
                                                     const GpuArray<Real,AMREX_SPACEDIM> dxInv,
-                                                    int icomp, int ncomp)
+                                                    int icomp, int ncomp, bool do_terrain_adjustment)
 {
     BL_PROFILE_VAR("impose_vertical_cons_bcs()",impose_vertical_cons_bcs);
     const auto& dom_lo = lbound(domain);
@@ -416,7 +416,7 @@ void ERFPhysBCFunct_cons::impose_vertical_cons_bcs (const Array4<Real>& dest_arr
         );
     }
 
-    if (m_z_phys_nd) {
+    if (do_terrain_adjustment &&  m_z_phys_nd) {
         const auto&  bx_lo = lbound(bx);
         const auto&  bx_hi = ubound(bx);
         const BCRec* bc_ptr_h = bcrs.data();

Source/BoundaryConditions/ERF_FillIntermediatePatch.cpp

@@ -61,11 +61,13 @@ 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::Print() << "LEVEL " << lev << " CONS ONLY   " << cons_only <<
+    //                   " 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);
+    if (!cons_only) {
+        AMREX_ALWAYS_ASSERT(mfs_mom.size() == IntVars::NumTypes);
+        AMREX_ALWAYS_ASSERT(mfs_vel.size() == Vars::NumTypes);
+    }
 
     // Enforce no penetration for thin immersed body
     if (!cons_only) {
@@ -122,7 +124,9 @@ ERF::FillIntermediatePatch (int lev, Real time,
         Vector<Real> ftime    = {time,time};
 
         // 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);
+        bool do_fb = true; bool do_terrain_adjustment = false;
+        (*physbcs_cons[lev])(*mfs_vel[Vars::cons],icomp_cons,ncomp_cons,IntVect{ng_cons},time,BCVars::cons_bc,
+                             do_fb, do_terrain_adjustment);
 
         if ( (icomp_cons+ncomp_cons > 1) && (interpolation_type == StateInterpType::Perturbational) )
         {

Source/BoundaryConditions/ERF_PhysBCFunct.H

@@ -53,7 +53,7 @@ public:
     */
     void operator() (amrex::MultiFab& mf, int icomp, int ncomp,
                      amrex::IntVect const& nghost, const amrex::Real time, int bccomp_cons,
-                     bool do_fb);
+                     bool do_fb = true, bool do_terrain_adjustment = true);
 
     void impose_lateral_cons_bcs (const amrex::Array4<amrex::Real>& dest_arr,
                                   const amrex::Box& bx, const amrex::Box& domain,
@@ -62,7 +62,7 @@ public:
                                    const amrex::Box& bx, const amrex::Box& domain,
                                    const amrex::Array4<amrex::Real const>& z_nd,
                                    const amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> dxInv,
-                                   int icomp, int ncomp);
+                                   int icomp, int ncomp, bool do_terrain_adjustment = true);
 
 private:
     int                  m_lev;

Source/BoundaryConditions/ERF_PhysBCFunct.cpp

@@ -15,7 +15,7 @@ using namespace amrex;
 
 void ERFPhysBCFunct_cons::operator() (MultiFab& mf, int icomp, int ncomp,
                                       IntVect const& nghost, const Real /*time*/, int /*bccomp*/,
-                                      bool do_fb)
+                                      bool do_fb, bool do_terrain_adjustment)
 {
     BL_PROFILE("ERFPhysBCFunct_cons::()");
 
@@ -32,21 +32,6 @@ void ERFPhysBCFunct_cons::operator() (MultiFab& mf, int icomp, int ncomp,
         }
     }
 
-    MultiFab z_nd_mf_loc;
-    if (m_z_phys_nd) {
-        m_z_phys_nd->FillBoundary(m_geom.periodicity());
-        BoxList bl_z_phys = convert(mf.boxArray(),IntVect(1,1,1)).boxList();
-        for (auto& b : bl_z_phys) {
-            b.setSmall(2,0);
-            b.setBig(2,1);
-        }
-        BoxArray ba_z(std::move(bl_z_phys));
-
-        z_nd_mf_loc.define(ba_z,mf.DistributionMap(),1,IntVect(nghost[0],nghost[1],0));
-        z_nd_mf_loc.ParallelCopy(*m_z_phys_nd,0,0,1,m_z_phys_nd->nGrowVect(),
-                                                    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.
@@ -77,7 +62,7 @@ void ERFPhysBCFunct_cons::operator() (MultiFab& mf, int icomp, int ncomp,
 
             if (m_z_phys_nd)
             {
-                z_nd_arr = z_nd_mf_loc.const_array(mfi);
+                z_nd_arr = m_z_phys_nd->const_array(mfi);
             }
 
             if (!gdomain.contains(cbx2))
@@ -91,7 +76,7 @@ void ERFPhysBCFunct_cons::operator() (MultiFab& mf, int icomp, int ncomp,
                 }
 
                 // We send the full FAB box with ghost cells
-                impose_vertical_cons_bcs(cons_arr,cbx2,domain,z_nd_arr,dxInv,icomp,ncomp);
+                impose_vertical_cons_bcs(cons_arr,cbx2,domain,z_nd_arr,dxInv,icomp,ncomp,do_terrain_adjustment);
             }
 
         } // MFIter
@@ -116,20 +101,6 @@ void ERFPhysBCFunct_u::operator() (MultiFab& mf, int /*icomp*/, int /*ncomp*/,
         }
     }
 
-    MultiFab z_nd_mf_loc;
-    if (m_z_phys_nd) {
-        m_z_phys_nd->FillBoundary(m_geom.periodicity());
-        BoxList bl_z_phys = convert(mf.boxArray(),IntVect(1,1,1)).boxList();
-        for (auto& b : bl_z_phys) {
-            b.setSmall(2,0);
-            b.setBig(2,1);
-        }
-        BoxArray ba_z(std::move(bl_z_phys));
-        z_nd_mf_loc.define(ba_z,mf.DistributionMap(),1,IntVect(nghost[0]+1,nghost[1],0));
-        z_nd_mf_loc.ParallelCopy(*m_z_phys_nd,0,0,1,m_z_phys_nd->nGrowVect(),
-                                                    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.
@@ -163,7 +134,7 @@ void ERFPhysBCFunct_u::operator() (MultiFab& mf, int /*icomp*/, int /*ncomp*/,
 
             if (m_z_phys_nd)
             {
-                z_nd_arr = z_nd_mf_loc.const_array(mfi);
+                z_nd_arr = m_z_phys_nd->const_array(mfi);
             }
 
             if (!gdomainx.contains(xbx2))
@@ -202,20 +173,6 @@ void ERFPhysBCFunct_v::operator() (MultiFab& mf, int /*icomp*/, int /*ncomp*/,
         }
     }
 
-    MultiFab z_nd_mf_loc;
-    if (m_z_phys_nd) {
-        m_z_phys_nd->FillBoundary(m_geom.periodicity());
-        BoxList bl_z_phys = convert(mf.boxArray(),IntVect(1,1,1)).boxList();
-        for (auto& b : bl_z_phys) {
-            b.setSmall(2,0);
-            b.setBig(2,1);
-        }
-        BoxArray ba_z(std::move(bl_z_phys));
-        z_nd_mf_loc.define(ba_z,mf.DistributionMap(),1,IntVect(nghost[0],nghost[1]+1,0));
-        z_nd_mf_loc.ParallelCopy(*m_z_phys_nd,0,0,1,m_z_phys_nd->nGrowVect(),
-                                                    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.
@@ -249,7 +206,7 @@ void ERFPhysBCFunct_v::operator() (MultiFab& mf, int /*icomp*/, int /*ncomp*/,
 
             if (m_z_phys_nd)
             {
-                z_nd_arr = z_nd_mf_loc.const_array(mfi);
+                z_nd_arr = m_z_phys_nd->const_array(mfi);
             }
 
             if (!gdomainy.contains(ybx2))
@@ -293,21 +250,6 @@ void ERFPhysBCFunct_w::operator() (MultiFab& mf, MultiFab& xvel, MultiFab& yvel,
     //
     if (gdomainz.smallEnd(2) == 0) gdomainz.setSmall(2,1);
 
-    Box ndomain  = convert(domain,IntVect(1,1,1));
-
-    MultiFab z_nd_mf_loc;
-    if (m_z_phys_nd) {
-        BoxList bl_z_phys = convert(mf.boxArray(),IntVect(1,1,1)).boxList();
-        for (auto& b : bl_z_phys) {
-            b &= ndomain;
-        }
-        BoxArray ba_z(std::move(bl_z_phys));
-        z_nd_mf_loc.define(ba_z,mf.DistributionMap(),1,IntVect(nghost[0],nghost[1],0));
-        z_nd_mf_loc.ParallelCopy(*m_z_phys_nd,0,0,1,m_z_phys_nd->nGrowVect(),
-                                                    z_nd_mf_loc.nGrowVect());
-    }
-    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.
@@ -339,7 +281,7 @@ void ERFPhysBCFunct_w::operator() (MultiFab& mf, MultiFab& xvel, MultiFab& yvel,
 
             if (m_z_phys_nd)
             {
-                z_nd_arr = z_nd_mf_loc.const_array(mfi);
+                z_nd_arr = m_z_phys_nd->const_array(mfi);
             }
 
             //

Source/ERF_make_new_level.cpp

@@ -47,7 +47,9 @@ 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;
+    if (verbose) {
+        amrex::Print() <<" BA FROM SCRATCH AT LEVEL " << lev << " " << ba << std::endl;
+    }
 
     if (lev == 0) init_bcs();
 
@@ -179,7 +181,9 @@ ERF::MakeNewLevelFromCoarse (int lev, Real time, const BoxArray& ba,
 {
     AMREX_ALWAYS_ASSERT(lev > 0);
 
-    // amrex::Print() <<" NEW BA FROM COARSE AT LEVEL " << lev << " " << ba << std::endl;
+    if (verbose) {
+        amrex::Print() <<" NEW BA FROM COARSE AT LEVEL " << lev << " " << ba << std::endl;
+    }
 
     //********************************************************************************************
     // This allocates all kinds of things, including but not limited to: solution arrays,
@@ -283,15 +287,19 @@ 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;
+    if (verbose) {
+        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);
 
     BoxArray            ba_old(vars_new[lev][Vars::cons].boxArray());
     DistributionMapping dm_old(vars_new[lev][Vars::cons].DistributionMap());
 
-    // amrex::Print() <<"               OLD BA AT LEVEL " << lev << " " << ba_old << std::endl;
+    if (verbose) {
+        amrex::Print() <<"               OLD BA AT LEVEL " << lev << " " << ba_old << std::endl;
+    }
 
     int     ncomp_cons  = vars_new[lev][Vars::cons].nComp();
     IntVect ngrow_state = vars_new[lev][Vars::cons].nGrowVect();

Source/IO/ERF_Plotfile.cpp

@@ -1410,10 +1410,15 @@ ERF::WritePlotFile (int which, PlotFileType plotfile_type, Vector<std::string> p
 
             int lev0 = 0;
             int desired_ratio = std::max(std::max(ref_ratio[lev0][0],ref_ratio[lev0][1]),ref_ratio[lev0][2]);
-            bool any_ratio_one = ( ( (ref_ratio[lev0][0] == 0) || (ref_ratio[lev0][1] == 0) ) ||
-                                     (ref_ratio[lev0][2] == 0) );
+            bool any_ratio_one = ( ( (ref_ratio[lev0][0] == 1) || (ref_ratio[lev0][1] == 1) ) ||
+                                     (ref_ratio[lev0][2] == 1) );
+            for (int lev = 1; lev < finest_level; lev++) {
+                any_ratio_one = any_ratio_one ||
+                                     ( ( (ref_ratio[lev][0] == 1) || (ref_ratio[lev][1] == 1) ) ||
+                                         (ref_ratio[lev][2] == 1) );
+            }
 
-            if (any_ratio_one == 1 && m_expand_plotvars_to_unif_rr)
+            if (any_ratio_one && m_expand_plotvars_to_unif_rr)
             {
                 Vector<IntVect>   r2(finest_level);
                 Vector<Geometry>  g2(finest_level+1);
@@ -1435,9 +1440,9 @@ ERF::WritePlotFile (int which, PlotFileType plotfile_type, Vector<std::string> p
 
                 for (int lev = 1; lev <= finest_level; ++lev) {
                     if (lev > 1) {
-                        r2[lev-1][0] *= desired_ratio / ref_ratio[lev-1][0];
-                        r2[lev-1][1] *= desired_ratio / ref_ratio[lev-1][1];
-                        r2[lev-1][2] *= desired_ratio / ref_ratio[lev-1][2];
+                        r2[lev-1][0] = r2[lev-2][0] * desired_ratio / ref_ratio[lev-1][0];
+                        r2[lev-1][1] = r2[lev-2][1] * desired_ratio / ref_ratio[lev-1][1];
+                        r2[lev-1][2] = r2[lev-2][2] * desired_ratio / ref_ratio[lev-1][2];
                     }
 
                     mf2[lev].define(refine(grids[lev],r2[lev-1]), dmap[lev], ncomp_mf, 0);