more fixes for multilevel (#2077)

* more fixes for multilevel * we need to impose bc's on base state when moving terrain as well * fix for moving geometry which I don't quite understand
erf-model · Jan 22, 2025 · 51783b2 · 51783b2
1 parent fbd2925
commit 51783b2
Show file tree

Hide file tree

Showing 3 changed files with 22 additions and 7 deletions.
diff --git a/Source/BoundaryConditions/ERF_FillCoarsePatch.cpp b/Source/BoundaryConditions/ERF_FillCoarsePatch.cpp
@@ -77,6 +77,11 @@ ERF::FillCoarsePatch (int lev, Real time)
                           geom[lev-1], geom[lev],
                           refRatio(lev-1), mapper_c, domain_bcs_type, BCVars::cons_bc);
 
+    // ***************************************************************************
+    // Physical bc's for cell centered variables at domain boundary
+    // ***************************************************************************
+    (*physbcs_cons[lev])(vars_new[lev][Vars::cons],0,ncomp_cons,ngvect_cons,time,BCVars::cons_bc,true);
+
     //
     //************************************************************************************************
     // Interpolate x-momentum from coarse to fine level
@@ -127,7 +132,6 @@ ERF::FillCoarsePatch (int lev, Real time)
     // ***************************************************************************
     IntVect ngvect_vels = vars_new[lev][Vars::xvel].nGrowVect();
 
-    (*physbcs_cons[lev])(vars_new[lev][Vars::cons],0,ncomp_cons,ngvect_cons,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],

diff --git a/Source/BoundaryConditions/ERF_PhysBCFunct.cpp b/Source/BoundaryConditions/ERF_PhysBCFunct.cpp
@@ -314,11 +314,6 @@ 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
@@ -358,7 +353,9 @@ void ERFPhysBCFunct_base::operator() (MultiFab& mf, int /*icomp*/, int ncomp, In
                 const Array4<Real> base_arr = mf.array(mfi);
 
                 impose_lateral_basestate_bcs(base_arr,cbx1,domain,ncomp,nghost);
-                impose_vertical_basestate_bcs(base_arr,cbx2,domain,ncomp,nghost);
+                if (!m_moving_terrain) { // TODO: I don't know why the CI test fails if this is called
+                   impose_vertical_basestate_bcs(base_arr,cbx2,domain,ncomp,nghost);
+                }
             }
 
         } // MFIter

diff --git a/Source/Initialization/ERF_Init1D.cpp b/Source/Initialization/ERF_Init1D.cpp
@@ -76,6 +76,11 @@ ERF::initHSE (int lev)
 
         erf_enforce_hse(lev, r_hse, p_hse, pi_hse, th_hse, z_phys_cc[lev]);
 
+        //
+        // Impose physical bc's on the base state
+        //
+        (*physbcs_base[lev])(base_state[lev],0,base_state[lev].nComp(),base_state[lev].nGrowVect());
+
     } else {
 
         BoxArray ba_new(domain);
@@ -114,6 +119,11 @@ ERF::initHSE (int lev)
 
         erf_enforce_hse(lev, new_r_hse, new_p_hse, new_pi_hse, new_th_hse, new_z_phys_cc);
 
+        //
+        // Impose physical bc's on the base state
+        //
+        (*physbcs_base[lev])(new_base_state,0,new_base_state.nComp(),new_base_state.nGrowVect());
+
         // Now copy back into the original arrays
         base_state[lev].ParallelCopy(new_base_state,0,0,base_state[lev].nComp(),
                                      base_state[lev].nGrowVect(),base_state[lev].nGrowVect());
@@ -169,6 +179,10 @@ ERF::erf_enforce_hse (int lev,
         b2d.grow(1,1);
         b2d.setRange(2,0);
 
+        // Intersect this box with the domain
+        Box zdomain = convert(geom[lev].Domain(),tbz.ixType());
+        b2d &= zdomain;
+
         // We integrate to the first cell (and below) by using rho in this cell
         // If gravity == 0 this is constant pressure
         // If gravity != 0, hence this is a wall, this gives gp0 = dens[0] * gravity