Use cell-centered grid for EB area fraction and face centroid (#2702)

* Apply FillBoundary for momenta in slow_rhs_post.

* check on domain boundaries not box boundaries

* replace bx lo/hi by domain lo/hi

* revert incorrect change in last commit

* fix test on small cells

* Added wrappers for non-const EB factory members.

* Use cell-centered grids for eb_aux_ area fraction and face centroids.

* Corrected rayleigh damping thickness in inputs_FittedMesh.

* Remove FB from slow_rhs_post, since we do this in erf_slow_rhs_pre.

---------

Co-authored-by: Ann Almgren <[email protected]>

Author: skang67

Exec/DryRegTests/WitchOfAgnesi/inputs_FittedMesh

@@ -64,8 +64,8 @@ erf.abl_driver_type = "PressureGradient"
 erf.abl_pressure_grad = -0.02 0. 0.
 
 erf.rayleigh_damp_W = true
-erf.rayleigh_zdamp = 5000.0
-erf.rayleigh_dampcoef = 0.2
+erf.rayleigh_zdamp = 50.0
+erf.rayleigh_dampcoef = 0.25
 
 #erf.init_type = "input_sounding"
 #erf.init_sounding_ideal = true

Source/EB/ERF_EB.H

@@ -99,8 +99,44 @@ class eb_ {
     inline amrex::FabArray<amrex::EBCellFlagFab>&
     getNonConstEBCellFlags(const amrex::EBFArrayBoxFactory& ebfact)
     {
-        const amrex::FabArray<amrex::EBCellFlagFab>& flags_const = ebfact.getMultiEBCellFlagFab();
-        return const_cast<amrex::FabArray<amrex::EBCellFlagFab>&>(flags_const);
+      const amrex::FabArray<amrex::EBCellFlagFab>& flags_const = ebfact.getMultiEBCellFlagFab();
+      return const_cast<amrex::FabArray<amrex::EBCellFlagFab>&>(flags_const);
+    }
+
+    inline amrex::MultiFab&
+    getNonConstVolFrac(const amrex::EBFArrayBoxFactory& ebfact)
+    {
+      const amrex::MultiFab& vfrac_const = ebfact.getVolFrac();
+      return const_cast<amrex::MultiFab&>(vfrac_const);
+    }
+
+    inline amrex::MultiCutFab&
+    getNonConstCentroid(const amrex::EBFArrayBoxFactory& ebfact)
+    {
+      const amrex::MultiCutFab& vcent_const = ebfact.getCentroid();
+      return const_cast<amrex::MultiCutFab&>(vcent_const);
+    }
+
+    inline amrex::Array<amrex::MultiCutFab*, AMREX_SPACEDIM>
+    getNonConstAreaFrac(const amrex::EBFArrayBoxFactory& ebfact)
+    {
+      auto afrac_const = ebfact.getAreaFrac();
+      amrex::Array<amrex::MultiCutFab*, AMREX_SPACEDIM> afrac;
+      for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) {
+          afrac[dir] = const_cast<amrex::MultiCutFab*>(afrac_const[dir]);
+      }
+      return afrac;
+    }
+
+    inline amrex::Array<amrex::MultiCutFab*, AMREX_SPACEDIM>
+    getNonConstFaceCent(const amrex::EBFArrayBoxFactory& ebfact)
+    {
+      auto fcent_const = ebfact.getFaceCent();
+      amrex::Array<amrex::MultiCutFab*, AMREX_SPACEDIM> fcent;
+      for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) {
+          fcent[dir] = const_cast<amrex::MultiCutFab*>(fcent_const[dir]);
+      }
+      return fcent;
     }
 
 };

Source/EB/ERF_EBAux.cpp

@@ -53,14 +53,8 @@ define( [[maybe_unused]] int const& a_level,
   m_volcent = new MultiFab(my_grids, a_dmap, AMREX_SPACEDIM, a_ngrow[2], MFInfo(), FArrayBoxFactory());
 
   for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
-      const BoxArray& faceba = amrex::convert(a_grids, IntVect::TheDimensionVector(idim));
-      if (idim == a_idim) {
-          m_areafrac[idim] = new MultiFab(a_grids, a_dmap,                1, a_ngrow[1]+1, MFInfo(), FArrayBoxFactory());
-          m_facecent[idim] = new MultiFab(a_grids, a_dmap, AMREX_SPACEDIM-1, a_ngrow[2], MFInfo(), FArrayBoxFactory());
-      } else {
-          m_areafrac[idim] = new MultiFab(faceba, a_dmap, 1, a_ngrow[1], MFInfo(), FArrayBoxFactory());
-          m_facecent[idim] = new MultiFab(faceba, a_dmap, AMREX_SPACEDIM-1, a_ngrow[2], MFInfo(), FArrayBoxFactory());
-      }
+    m_areafrac[idim] = new MultiFab(a_grids, a_dmap,                1, a_ngrow[1]+1, MFInfo(), FArrayBoxFactory());
+    m_facecent[idim] = new MultiFab(a_grids, a_dmap, AMREX_SPACEDIM-1, a_ngrow[2], MFInfo(), FArrayBoxFactory());
   }
 
   m_bndryarea = new MultiFab(my_grids, a_dmap, 1, a_ngrow[2], MFInfo(), FArrayBoxFactory());
@@ -848,23 +842,22 @@ define( [[maybe_unused]] int const& a_level,
 
   for (MFIter mfi(*m_cellflags, false); mfi.isValid(); ++mfi) {
 
-      const Box& bx = mfi.validbox();
-      const Box& bx_grown = mfi.growntilebox();
-
-      Array4<EBCellFlag> const& aux_flag  = m_cellflags->array(mfi);
-      Array4<Real>       const& aux_vfrac = m_volfrac->array(mfi);
-      Array4<Real>       const& aux_afrac_x = m_areafrac[0]->array(mfi);
-      Array4<Real>       const& aux_afrac_y = m_areafrac[1]->array(mfi);
-      Array4<Real>       const& aux_afrac_z = m_areafrac[2]->array(mfi);
+    const Box& bx = mfi.validbox();
+    const Box& bx_grown = mfi.growntilebox();
 
-      Array4<Real>       const& aux_vcent = m_volcent->array(mfi);
-      Array4<Real>       const& aux_fcent_x = m_facecent[0]->array(mfi);
-      Array4<Real>       const& aux_fcent_y = m_facecent[1]->array(mfi);
-      Array4<Real>       const& aux_fcent_z = m_facecent[2]->array(mfi);
-      Array4<Real>       const& aux_barea = m_bndryarea->array(mfi);
-      Array4<Real>       const& aux_bcent = m_bndrycent->array(mfi);
-      Array4<Real>       const& aux_bnorm = m_bndrynorm->array(mfi);
+    Array4<EBCellFlag> const& aux_flag  = m_cellflags->array(mfi);
+    Array4<Real>       const& aux_vfrac = m_volfrac->array(mfi);
+    Array4<Real>       const& aux_afrac_x = m_areafrac[0]->array(mfi);
+    Array4<Real>       const& aux_afrac_y = m_areafrac[1]->array(mfi);
+    Array4<Real>       const& aux_afrac_z = m_areafrac[2]->array(mfi);
 
+    Array4<Real>       const& aux_vcent = m_volcent->array(mfi);
+    Array4<Real>       const& aux_fcent_x = m_facecent[0]->array(mfi);
+    Array4<Real>       const& aux_fcent_y = m_facecent[1]->array(mfi);
+    Array4<Real>       const& aux_fcent_z = m_facecent[2]->array(mfi);
+    Array4<Real>       const& aux_barea = m_bndryarea->array(mfi);
+    Array4<Real>       const& aux_bcent = m_bndrycent->array(mfi);
+    Array4<Real>       const& aux_bnorm = m_bndrynorm->array(mfi);
 
     if (FlagFab[mfi].getType(bx) == FabType::singlevalued ) {