Correct z_phys and z_ref in MOST for EB.

Source/BoundaryConditions/ERF_MOSTAverage.H

@@ -64,6 +64,9 @@ public:
     // Populate a 2D iMF k_indx (w/o terrain)
     void set_k_indices_N (const int& lev);
 
+    // Populate a 2D iMF k_indx (for EB)
+    void set_k_indices_EB (const int& lev);
+
     // Populate a 2D iMF k_indx (w/ terrain)
     void set_k_indices_T (const int& lev);
 

Source/BoundaryConditions/ERF_MOSTAverage.cpp

@@ -97,6 +97,8 @@ MOSTAverage::make_MOSTAverage_at_level (const int& lev,
     bool use_terrain_fitted_coords = ( (m_terrain_type == TerrainType::StaticFittedMesh) ||
                                        (m_terrain_type == TerrainType::MovingFittedMesh) );
 
+    bool use_eb = (m_terrain_type == TerrainType::EB);
+
     { // Nodal in x
         auto& mf = *vars_old[Vars::xvel];
         // Create a 2D ba, dm, & ghost cells
@@ -208,6 +210,8 @@ MOSTAverage::make_MOSTAverage_at_level (const int& lev,
         set_norm_indices_T(lev);
     } else if (use_terrain_fitted_coords) {                    // Terrain
         set_k_indices_T(lev);
+    } else if (use_eb) {                                       // EB
+        set_k_indices_EB(lev);
     } else {                                                   // No Terrain
         set_k_indices_N(lev);
     }
@@ -449,6 +453,60 @@ MOSTAverage::set_k_indices_N (const int& lev)
     }
 }
 
+/**
+ * Function to set K indices for EB.
+ *
+ */
+void
+MOSTAverage::set_k_indices_EB (const int& lev)
+{
+    ParmParse pp(m_pp_prefix);
+    Real zref_tmp = zref_default;
+    auto read_z = pp.query("most.zref",zref_tmp);
+    auto read_k = pp.queryarr("most.k_arr_in",m_k_in);
+
+    ParmParse pp_eb2("eb2");
+    std::string geometry;
+    pp_eb2.query("geometry", geometry);
+    AMREX_ALWAYS_ASSERT_WITH_MESSAGE(geometry == "plane", "Only plane geometry supported for EB MOST");
+
+    RealArray plane_point{zero, zero, zero};
+    pp_eb2.query("plane_point", plane_point);
+    Real z_eb = plane_point[2];
+
+    // Default behavior is to use the first cell center
+    if (!read_z && !read_k) {
+        Real m_dz  = m_geom[0].CellSize(2);
+        zref_tmp = myhalf * m_dz;
+        m_zref[lev]->setVal( zref_tmp );
+        Print() << "Reference height for MOST set to " << zref_tmp << std::endl;
+        read_z = true;
+    }
+
+    // Specify z_ref & compute k_indx (z_ref takes precedence)
+    if (read_z) {
+        Real m_zlo = m_geom[lev].ProbLo(2);
+        Real m_zhi = m_geom[lev].ProbHi(2);
+        Real m_dz  = m_geom[lev].CellSize(2);
+
+        amrex::ignore_unused(m_zhi);
+
+        AMREX_ALWAYS_ASSERT_WITH_MESSAGE(zref_tmp >= m_zlo + myhalf * m_dz,
+                                         "Query point must be past first z-cell!");
+
+        AMREX_ALWAYS_ASSERT_WITH_MESSAGE(zref_tmp <= m_zhi - myhalf * m_dz,
+                                         "Query point must be below the last z-cell!");
+
+        int lk = static_cast<int>(floor((zref_tmp + z_eb - m_zlo) / m_dz - myhalf));
+        int lk_phys = static_cast<int>(floor((zref_tmp - m_zlo) / m_dz - myhalf));
+
+        m_zref[lev]->setVal( (lk_phys + 0.5) * m_dz + m_zlo );
+
+        AMREX_ALWAYS_ASSERT(lk >= m_radius);
+
+        m_k_indx[lev]->setVal(lk);
+    }
+}
 
 /**
  * Function to set K indices with terrain (w/o terrain normals or interpolation).

Source/ERF_MakeNewArrays.cpp

@@ -53,8 +53,21 @@ ERF::init_stuff (int lev, const BoxArray& ba, const DistributionMapping& dm,
     int ngrow = ComputeGhostCells(solverChoice) + 2;
     tmp_zphys_nd = std::make_unique<MultiFab>(ba_nd,dm,1,IntVect(ngrow,ngrow,ngrow));
 
+    // Offset z-coordinate for interpolate_1d when plane EB is used.
+    Real z_offset = zero;
+    if (solverChoice.terrain_type == TerrainType::EB) {
+        ParmParse pp_eb2("eb2");
+        std::string geometry;
+        pp_eb2.query("geometry", geometry);
+        if (geometry == "plane") {
+            RealArray plane_point{zero, zero, zero};
+            pp_eb2.query("plane_point", plane_point);
+            z_offset = plane_point[2];
+        }
+    }
+
     z_phys_cc[lev] = std::make_unique<MultiFab>(ba,dm,1,2);
-    init_default_zphys(lev, geom[lev], *tmp_zphys_nd, *z_phys_cc[lev]);
+    init_default_zphys(lev, geom[lev], *tmp_zphys_nd, *z_phys_cc[lev], z_offset);
 
     if (solverChoice.terrain_type == TerrainType::MovingFittedMesh)
     {
@@ -656,6 +669,15 @@ ERF::update_diffusive_arrays (int lev, const BoxArray& ba, const DistributionMap
 void
 ERF::init_zphys (int lev, Real elapsed_time)
 {
+    // For EB, z_phys_nd was already initialized with the correct z_offset by init_default_zphys.
+    // The terrain-fitting (BTF) done below is irrelevant for a flat EB mesh and would clobber
+    // the offset, so return early here.
+    if (solverChoice.terrain_type == TerrainType::EB) {
+        Real dzmin = get_dzmin_terrain(*z_phys_nd[lev]);
+        micro->Set_dzmin(lev, dzmin);
+        return;
+    }
+
     if (solverChoice.init_type != InitType::WRFInput && solverChoice.init_type != InitType::Metgrid)
     {
         if (lev > 0) {

Source/Utils/ERF_TerrainMetrics.H

@@ -12,7 +12,8 @@
  */
 void
 init_default_zphys (int lev, const amrex::Geometry& geom,
-                    amrex::MultiFab& z_phys_nd, amrex::MultiFab& z_phys_cc);
+                    amrex::MultiFab& z_phys_nd, amrex::MultiFab& z_phys_cc,
+                    amrex::Real z_offset = zero);
 
 /**
  * Utility routines for constructing terrain metric terms

Source/Utils/ERF_TerrainMetrics.cpp

@@ -12,7 +12,8 @@ using namespace amrex;
  * This will be over-written if we use z_levels, or grid stretching, or terrain-fitted grids
  */
 void
-init_default_zphys (int /*lev*/, const Geometry& geom, MultiFab& z_phys_nd, MultiFab& z_phys_cc)
+init_default_zphys (int /*lev*/, const Geometry& geom, MultiFab& z_phys_nd, MultiFab& z_phys_cc,
+                    Real z_offset)
 {
     const auto& dx = geom.CellSize();
     Real dz = dx[2];
@@ -23,7 +24,7 @@ init_default_zphys (int /*lev*/, const Geometry& geom, MultiFab& z_phys_nd, Mult
         const Array4< Real> z_nd_arr = z_phys_nd.array(mfi);
         ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
         {
-            z_nd_arr(i,j,k) = k * dz;
+            z_nd_arr(i,j,k) = k * dz - z_offset;
         });
     }
 
@@ -33,7 +34,7 @@ init_default_zphys (int /*lev*/, const Geometry& geom, MultiFab& z_phys_nd, Mult
         const Array4< Real> z_cc_arr = z_phys_cc.array(mfi);
         ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
         {
-            z_cc_arr(i,j,k) = (k + myhalf) * dz;
+            z_cc_arr(i,j,k) = (k + myhalf) * dz - z_offset;
         });
     }
 }