eb_aux_ Class Reference

Stores face-centered EB geometry reconstructed from a cell-centered factory. More...

#include <ERF_EBAux.H>

Collaboration diagram for eb_aux_:

Public Member Functions
	eb_aux_ ()
	Construct an undefined auxiliary EB factory. More...
	~eb_aux_ ()
	Destroy the auxiliary EB factory. More...
void	define (int const &a_level, int const &a_idim, amrex::Geometry const &a_geom, amrex::BoxArray const &a_grids, amrex::DistributionMapping const &a_dmap, amrex::Vector< int > const &a_ngrow, amrex::EBFArrayBoxFactory const *a_factory)
	Define face-centered EB geometry for one staggered direction. More...
void	set_verbose ()
	Enable verbose diagnostic output during factory construction. More...
const amrex::FabArray< amrex::EBCellFlagFab > &	getMultiEBCellFlagFab () const
	Return the reconstructed EB cell flags. More...
const amrex::MultiFab &	getVolFrac () const
	Return the reconstructed volume fractions. More...
const amrex::MultiFab &	getCentroid () const
	Return the reconstructed volume centroids. More...
const amrex::MultiFab &	getBndryArea () const
	Return the reconstructed EB boundary areas. More...
const amrex::MultiFab &	getBndryCent () const
	Return the reconstructed EB boundary centroids. More...
const amrex::MultiFab &	getBndryNormal () const
	Return the reconstructed EB boundary normals. More...
amrex::Array< const amrex::MultiFab *, AMREX_SPACEDIM >	getAreaFrac () const
	Return reconstructed area fractions for all coordinate directions. More...
amrex::Array< const amrex::MultiFab *, AMREX_SPACEDIM >	getFaceCent () const
	Return reconstructed face centroids for all coordinate directions. More...

Private Attributes
int	m_verbose
amrex::FabArray< amrex::EBCellFlagFab > *	m_cellflags = nullptr
amrex::MultiFab *	m_volfrac = nullptr
amrex::MultiFab *	m_volcent = nullptr
amrex::MultiFab *	m_bndryarea = nullptr
amrex::MultiFab *	m_bndrycent = nullptr
amrex::MultiFab *	m_bndrynorm = nullptr
amrex::Array< amrex::MultiFab *, AMREX_SPACEDIM >	m_areafrac {{AMREX_D_DECL(nullptr, nullptr, nullptr)}}
amrex::Array< amrex::MultiFab *, AMREX_SPACEDIM >	m_facecent {{AMREX_D_DECL(nullptr, nullptr, nullptr)}}

Detailed Description

Stores face-centered EB geometry reconstructed from a cell-centered factory.

This class mirrors the subset of the AMReX EB factory interface needed by ERF staggered momentum fields when native face-centered EB factories are disabled.

Constructor & Destructor Documentation

eb_aux_()

eb_aux_::eb_aux_

(

)

Construct an undefined auxiliary EB factory.

  : m_verbose(0)
// ,m_defined(0)
{}

~eb_aux_()

eb_aux_::~eb_aux_

(

)

Destroy the auxiliary EB factory.

{
}

Member Function Documentation

define()

void eb_aux_::define	(	int const &	a_level,
		int const &	a_idim,
		amrex::Geometry const &	a_geom,
		amrex::BoxArray const &	a_grids,
		amrex::DistributionMapping const &	a_dmap,
		amrex::Vector< int > const &	a_ngrow,
		amrex::EBFArrayBoxFactory const *	a_factory
	)

Define face-centered EB geometry for one staggered direction.

Parameters

a_level	AMR level index.
a_idim	Staggered direction for the target velocity grid.
a_geom	Cell-centered level geometry.
a_grids	Cell-centered grid boxes.
a_dmap	Distribution mapping for the grids.
a_ngrow	Ghost-cell counts for flags, volumes, and centroids.
a_factory	Source cell-centered EB factory.

{
  // Box dbox(a_geom.Domain());
 
  // small_volfrac
  Real small_volfrac = Real(1.e-14);
  ParmParse pp("eb2");
  pp.queryAdd("small_volfrac", small_volfrac);
  const Real small_value = Real(1.e-15);
 
  const IntVect vdim(IntVect::TheDimensionVector(a_idim));
 
  const BoxArray& my_grids = amrex::convert(a_grids, vdim);
 
  m_cellflags = new FabArray<EBCellFlagFab>(my_grids, a_dmap, 1, a_ngrow[0], MFInfo(),
                                            DefaultFabFactory<EBCellFlagFab>());
 
  // Set m_cellflags type to singlevalued
  m_cellflags->setVal(EBCellFlag::TheDefaultCell());
  for (MFIter mfi(*m_cellflags, false); mfi.isValid(); ++mfi) {
    auto& fab = (*m_cellflags)[mfi];
    fab.setType(FabType::singlevalued);
  }
 
  m_volfrac = new MultiFab(my_grids, a_dmap, 1, a_ngrow[1], MFInfo(), FArrayBoxFactory());
  m_volcent = new MultiFab(my_grids, a_dmap, AMREX_SPACEDIM, a_ngrow[2], MFInfo(), FArrayBoxFactory());
 
  for (int idim = 0; idim < AMREX_SPACEDIM; ++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());
  }
 
  m_bndryarea = new MultiFab(my_grids, a_dmap, 1, a_ngrow[2], MFInfo(), FArrayBoxFactory());
  m_bndrycent = new MultiFab(my_grids, a_dmap, AMREX_SPACEDIM, a_ngrow[2], MFInfo(), FArrayBoxFactory());
  m_bndrynorm = new MultiFab(my_grids, a_dmap, AMREX_SPACEDIM, a_ngrow[2], MFInfo(), FArrayBoxFactory());
 
  // Initialize with zeros
  m_volfrac->setVal(0.0);
  m_volcent->setVal(0.0);
 
  for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
    m_areafrac[idim]->setVal(0.0);
    m_facecent[idim]->setVal(0.0);
  }
 
  m_bndryarea->setVal(0.0);
  m_bndrycent->setVal(0.0);
  m_bndrynorm->setVal(0.0);
 
  const auto& FlagFab = a_factory->getMultiEBCellFlagFab(); // EBFArrayBoxFactory, EBDataCollection
 
  for (MFIter mfi(*m_cellflags, false); mfi.isValid(); ++mfi) {
 
    const Box& bx = mfi.validbox();
    const Box& bx_grown = mfi.growntilebox();
    const Box tbx = mfi.nodaltilebox(a_idim);
    const Box domain = surroundingNodes(a_geom.Domain(), a_idim);
 
    GpuArray<Real, AMREX_SPACEDIM> dx = a_geom.CellSizeArray();
    bool l_periodic   = a_geom.isPeriodic(a_idim);
 
    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);
 
    if (FlagFab[mfi].getType(bx) == FabType::covered ) {
 
      ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
      {
        aux_flag(i,j,k).setCovered();
        aux_flag(i,j,k).setDisconnected();
        if (i==bx.bigEnd(0)) {
          aux_flag(i+1,j,k).setCovered();
        }
        if (j==bx.bigEnd(1)) {
          aux_flag(i,j+1,k).setCovered();
        }
        if (k==bx.bigEnd(2)) {
          aux_flag(i,j,k+1).setCovered();
        }
      });
 
    } else if (FlagFab[mfi].getType(bx) == FabType::regular ) {
 
      ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
      {
        aux_flag(i,j,k).setRegular();
        aux_flag(i,j,k).setDisconnected();
        aux_vfrac(i,j,k) = one;
        aux_afrac_x(i,j,k) = one;
        aux_afrac_y(i,j,k) = one;
        aux_afrac_z(i,j,k) = one;
        if (i==bx.bigEnd(0)) {
          aux_flag(i+1,j,k).setRegular();
          aux_vfrac(i+1,j,k) = one;
          aux_afrac_x(i+1,j,k) = one;
        }
        if (j==bx.bigEnd(1)) {
          aux_flag(i,j+1,k).setRegular();
          aux_vfrac(i,j+1,k) = one;
          aux_afrac_y(i,j+1,k) = one;
        }
        if (k==bx.bigEnd(2)) {
          aux_flag(i,j,k+1).setRegular();
          aux_vfrac(i,j,k+1) = one;
          aux_afrac_z(i,j,k+1) = one;
        }
      });
 
    } else if (FlagFab[mfi].getType(bx) == FabType::singlevalued ) {
 
      // Initialization
 
      // CC cell quantities
      Array4<EBCellFlag const> const& flag = FlagFab.const_array(mfi);
      Array4<Real const> const& afrac = (a_factory->getAreaFrac()[a_idim])->const_array(mfi);
      Array4<Real const> const& bnorm = a_factory->getBndryNormal()[mfi].const_array();
      Array4<Real const> const& bcent = a_factory->getBndryCent()[mfi].const_array();
 
      // aux quantities
      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);
 
      // Extended domain in the direction of periodicity
      Box dom_grown = domain;
      for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
        if (a_geom.isPeriodic(idim)) {
          dom_grown.grow(idim, a_ngrow[0]);
        }
      }
 
      const IntVect dom_grown_lo = dom_grown.smallEnd();
      const IntVect dom_grown_hi = dom_grown.bigEnd();
 
      BoxList diffList = boxDiff(bx_grown, bx);
      for (const Box& b : diffList) {
        ParallelFor(b, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
        {
          if ( i < dom_grown_lo[0] || i > dom_grown_hi[0] ||
               j < dom_grown_lo[1] || j > dom_grown_hi[1] ||
               k < dom_grown_lo[2] || k > dom_grown_hi[2] ) {
            aux_flag(i,j,k).setCovered();
            aux_flag(i,j,k).setDisconnected();
          }
        });
      }
 
#ifndef AMREX_USE_GPU
      int const verbose=m_verbose;
#endif
 
      ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
      {
        // defaults to covered and disconnected.
        aux_flag(i,j,k).setCovered();
        aux_flag(i,j,k).setDisconnected();
 
        if (i==bx.bigEnd(0)) {
          aux_flag(i+1,j,k).setCovered();
        }
        if (j==bx.bigEnd(1)) {
          aux_flag(i,j+1,k).setCovered();
        }
        if (k==bx.bigEnd(2)) {
          aux_flag(i,j,k+1).setCovered();
        }
 
        // Index for low and hi cells
        IntVect iv_hi(i,j,k);
        IntVect iv_lo(iv_hi - vdim);
 
        bool lo_isCovered = flag(iv_lo).isCovered();
        bool hi_isCovered = flag(iv_hi).isCovered();
        bool lo_isRegular = flag(iv_lo).isRegular();
        bool hi_isRegular = flag(iv_hi).isRegular();
        bool lo_isSingleValued = flag(iv_lo).isSingleValued();
        bool hi_isSingleValued = flag(iv_hi).isSingleValued();
 
        const bool at_lo_boundary = (!l_periodic && iv_hi[a_idim]==domain.smallEnd(a_idim));
        const bool at_hi_boundary = (!l_periodic && iv_hi[a_idim]==domain.bigEnd(a_idim));
 
        // Treatment of lower boundary
 
        if (at_lo_boundary) {
          if (hi_isCovered) {
            lo_isCovered = true;
            lo_isRegular = false;
            lo_isSingleValued = false;
          } else if (hi_isRegular) {
            lo_isCovered = false;
            lo_isRegular = true;
            lo_isSingleValued = false;
          } else if (hi_isSingleValued) {
            if (almostEqual(afrac(i,j,k),zero)) {
              lo_isCovered = true;
              lo_isRegular = false;
              lo_isSingleValued = false;
            } else if (almostEqual(afrac(i,j,k),one)) {
              lo_isCovered = false;
              lo_isRegular = true;
              lo_isSingleValued = false;
            } else {
              lo_isCovered = false;
              lo_isRegular = false;
              lo_isSingleValued = true;
              iv_lo = iv_hi; // At the lower boundary, low cell takes the values of the high cell.
            }
          }
        }
 
        // Treatment of upper boundary
 
        if (at_hi_boundary) {
          if (lo_isCovered) { // Covered
            hi_isCovered = true;
            hi_isRegular = false;
            hi_isSingleValued = false;
          } else if (lo_isRegular) { // Regular
            hi_isCovered = false;
            hi_isRegular = true;
            hi_isSingleValued = false;
          } else if (lo_isSingleValued) { // SingleValued
            if (almostEqual(afrac(i,j,k),zero)) { //Covered
              hi_isCovered = true;
              hi_isRegular = false;
              hi_isSingleValued = false;
            } else if (almostEqual(afrac(i,j,k),one)) { //Regular
              hi_isCovered = false;
              hi_isRegular = true;
              hi_isSingleValued = false;
            } else { // SingleValued
              hi_isCovered = false;
              hi_isRegular = false;
              hi_isSingleValued = true;
              iv_hi = iv_lo; // At the upper boundary, hi cell takes the values of the low cell.
            }
          }
        }
 
        if ( lo_isCovered && hi_isCovered) {
 
          // defaults to covered and disconnected.
 
        } else if ( lo_isRegular && hi_isRegular) {
 
          aux_flag(i,j,k).setRegular();
          aux_flag(i,j,k).setConnected();
 
          aux_vfrac(i,j,k) = one;
 
          aux_afrac_x(i,j,k) = one;
          aux_afrac_y(i,j,k) = one;
          aux_afrac_z(i,j,k) = one;
 
          if (i==bx.bigEnd(0)) {
            aux_afrac_x(i+1,j,k) = one;
          }
          if (j==bx.bigEnd(1)) {
            aux_afrac_y(i,j+1,k) = one;
          }
          if (k==bx.bigEnd(2)) {
            aux_afrac_z(i,j,k+1) = one;
          }
 
        } else {
 
#ifndef AMREX_USE_GPU
          if (verbose) { Print() << "\ncell: " << amrex::IntVect(i,j,k) << "\n"; }
#endif
          Array<Real,AMREX_SPACEDIM> lo_arr = {-myhalf,-myhalf,-myhalf};
          Array<Real,AMREX_SPACEDIM> hi_arr = { myhalf, myhalf, myhalf};
 
          //-----------------------
          // Low EB cut cell
          //-----------------------
 
          // Map bcent and bnorm to the isoparametric space for anisotropic grids.
          // (This step is needed because bcent in AMReX is isotropically normalized.)
 
          RealVect lo_point (bcent(iv_lo,0), bcent(iv_lo,1), bcent(iv_lo,2));
          RealVect lo_normal(bnorm(iv_lo,0), bnorm(iv_lo,1), bnorm(iv_lo,2));
 
          if (at_lo_boundary) { // At lower boundary
            lo_point[a_idim] += one; // Move the boundary centroid upward in the a_idim direction.
          }
 
          if (lo_isSingleValued ) {
            Real bnorm_x = bnorm(iv_lo,0) * dx[0];
            Real bnorm_y = bnorm(iv_lo,1) * dx[1];
            Real bnorm_z = bnorm(iv_lo,2) * dx[2];
 
            Real norm = std::sqrt( bnorm_x*bnorm_x + bnorm_y*bnorm_y + bnorm_z*bnorm_z);
 
            RealVect bnorm_isoparam ( bnorm_x / norm, bnorm_y / norm, bnorm_z / norm);
 
            lo_normal = bnorm_isoparam;
          }
 
          // High side of low cell
          lo_arr[a_idim] = zero;
          hi_arr[a_idim] = myhalf;
          RealBox lo_rbx(lo_arr.data(), hi_arr.data());
 
          eb_cut_cell_ lo_eb_cc(flag(iv_lo), lo_rbx, lo_point, lo_normal);
 
          // cell iv_lo covered (regular) implies lo_eb_cc is covered (regular)
          // The inverse is not always true.
          AMREX_ASSERT( !lo_isCovered || lo_eb_cc.isCovered() );
          AMREX_ASSERT( !lo_isRegular || lo_eb_cc.isRegular() );
 
          //-----------------------
          // High EB cut cell
          //-----------------------
 
          RealVect hi_point (bcent(iv_hi,0), bcent(iv_hi,1), bcent(iv_hi,2));
          RealVect hi_normal(bnorm(iv_hi,0), bnorm(iv_hi,1), bnorm(iv_hi,2));
 
          if (at_hi_boundary) {
            hi_point[a_idim] += -one; // Move the boundary centroid downward in the a_idim direction.
          }
 
          if (hi_isSingleValued ) {
            Real bnorm_x = bnorm(iv_hi,0) * dx[0];
            Real bnorm_y = bnorm(iv_hi,1) * dx[1];
            Real bnorm_z = bnorm(iv_hi,2) * dx[2];
 
            Real norm = std::sqrt( bnorm_x*bnorm_x + bnorm_y*bnorm_y + bnorm_z*bnorm_z);
 
            RealVect bnorm_isoparam ( bnorm_x / norm, bnorm_y / norm, bnorm_z / norm);
 
            hi_normal = bnorm_isoparam;
          }
 
          // Low side of high cell
          lo_arr[a_idim] = -myhalf;
          hi_arr[a_idim] =  zero;
          RealBox hi_rbx(lo_arr.data(), hi_arr.data());
 
          eb_cut_cell_ hi_eb_cc(flag(iv_hi), hi_rbx, hi_point, hi_normal);
 
          // cell iv_hi covered (regular) implies hi_eb_cc is covered (regular)
          // The inverse is not always true.
          AMREX_ASSERT( !hi_isCovered || hi_eb_cc.isCovered() );
          AMREX_ASSERT( !hi_isRegular || hi_eb_cc.isRegular() );
 
#if 0
#if defined(AMREX_DEBUG) || defined(AMREX_TESTING) || 1
 
          { /***************************** SANITY CHECK ***********************\
            * Perform some basic sanity checks to verify that what we computed *
            * for cell (i,j,k) compares to what we know to be true.           *
            \******************************************************************/
 
            // Compute the cut-cell for the high side of the high cell. This is
            // only needed for sanity checks.
 
            eb_cut_cell_ hi_hi_eb_cc(flag(iv_hi), lo_rbx, hi_point, hi_normal);
 
            // cell iv_hi covered (regular) implies hi_hi_eb_cc is covered (regular)
            // The inverse is not always true.
#ifndef AMREX_USE_GPU
            if ( !(!hi_isRegular || hi_hi_eb_cc.isRegular()) ||
                 !(!hi_isCovered || hi_hi_eb_cc.isCovered()) ) {
              Print() << "flag(iv_hi) and hi_hi_eb_cc flags do not agree\n"
                      << "\n  isRegular() " << hi_isRegular << "  " << hi_hi_eb_cc.isRegular()
                      << "\n  isCovered() " << hi_isCovered << "  " << hi_hi_eb_cc.isCovered()
                      << "\n";
            }
#endif
            // If cell iv_hi is regular or covered, then hi_hi_eb_cc must also
            // be regular or covered. The inverse is not true.
            AMREX_ALWAYS_ASSERT( !hi_isRegular || hi_hi_eb_cc.isRegular() );
            AMREX_ALWAYS_ASSERT( !hi_isCovered || hi_hi_eb_cc.isCovered() );
 
            // The area and volume fractions that are computed for the scalar grid
            // are slightly different than those we compute from the geometric
            // reconstruction using the EB point and normal. However, we expect
            // that the area fractions computed here will give back the same
            // normal we used to compute them.
            if ( hi_isSingleValued ) {
 
              Real const adx = (a_idim == 0)
                             ? (hi_eb_cc.areaLo(0) - hi_hi_eb_cc.areaHi(0)) * dx[1] * dx[2]
                             : (hi_eb_cc.areaLo(0) + hi_hi_eb_cc.areaLo(0)) * dx[1] * dx[2]
                             - (hi_eb_cc.areaHi(0) + hi_hi_eb_cc.areaHi(0)) * dx[1] * dx[2];
 
              Real const ady = (a_idim == 1)
                             ? (hi_eb_cc.areaLo(1) - hi_hi_eb_cc.areaHi(1)) * dx[0] * dx[2]
                             : (hi_eb_cc.areaLo(1) + hi_hi_eb_cc.areaLo(1)) * dx[0] * dx[2]
                             - (hi_eb_cc.areaHi(1) + hi_hi_eb_cc.areaHi(1)) * dx[0] * dx[2];
 
              Real const adz = (a_idim == 2)
                             ? (hi_eb_cc.areaLo(2) - hi_hi_eb_cc.areaHi(2)) * dx[0] * dx[1]
                             : (hi_eb_cc.areaLo(2) + hi_hi_eb_cc.areaLo(2)) * dx[0] * dx[1]
                             - (hi_eb_cc.areaHi(2) + hi_hi_eb_cc.areaHi(2)) * dx[0] * dx[1];
 
              Real const apnorm = std::sqrt(adx*adx + ady*ady + adz*adz);
 
              // EB normal
              Real const apnorminv = one / apnorm;
              RealVect const normal(adx*apnorminv, ady*apnorminv, adz*apnorminv);
              Real const dot_normals = normal.dotProduct(hi_normal);
 
#ifndef AMREX_USE_GPU
              if ( !amrex::almostEqual(dot_normals, one) ) {
                Print() << "\nFail: check-1 dot_normals " << dot_normals
                        << '\n';
 
                hi_eb_cc.debug();
                hi_hi_eb_cc.debug();
 
              } else if (verbose) {
                Print() << "Pass: dot_normals = one\n";
 
              }
#endif
              AMREX_ALWAYS_ASSERT( amrex::almostEqual(dot_normals, one) );
            }
 
            // The a_idim area of hi_eb_cc.areaHi() should equal hi_hi_eb_cc.areaLo()
            {
#ifndef AMREX_USE_GPU
            Real const abs_err = std::abs( hi_eb_cc.areaHi(a_idim) - hi_hi_eb_cc.areaLo(a_idim) );
            Real machine_tol = Real(10.0)*std::numeric_limits<amrex::Real>::epsilon();
            if ( abs_err >= machine_tol ) {
                Print() << "\nFail: check-2 area abs_err: " << abs_err
                        << "\n  hi_eb_cc.areaHi " << hi_eb_cc.areaHi(a_idim)
                        << "\n  hi_hi_eb_cc.areaLo " << hi_hi_eb_cc.areaLo(a_idim)
                        << '\n';
            } else if (verbose) {
                Print() << "Pass: hi_eb_cc.areaHi = hi_hi_eb_cc.areaLo"
                        << "  abs_err: " << abs_err << "\n";
            }
            AMREX_ALWAYS_ASSERT( abs_err < machine_tol );
#endif
            }
 
            // The low-side area of hi_eb_cc should equal a_idim afrac.
            { Real const abs_err = amrex::max(std::abs(lo_eb_cc.areaHi(a_idim) - afrac(iv_hi)),
                                              std::abs(hi_eb_cc.areaLo(a_idim) - afrac(iv_hi)));
              Real compare_tol = Real(5.0e-6);
#ifndef AMREX_USE_GPU
              if ( abs_err >= compare_tol ) {
                //hi_eb_cc.debug();
                Print() << "\nFail: check-3 area abs_err " << abs_err
                        << "\n  hi_eb_cc.areaLo(" << a_idim << ") = " << hi_eb_cc.areaLo(a_idim)
                        << "\n  lo_eb_cc.areaHi(" << a_idim << ") = " << lo_eb_cc.areaHi(a_idim)
                        << "\n  afrac" << iv_hi << " =  " << afrac(iv_hi)
                        << '\n';
              } else if (verbose) {
                Print() << "Pass: hi_eb_cc.areaLo = afrac = " << afrac(iv_hi)
                        << "  abs_err: " << abs_err << "\n";
              }
#endif
              AMREX_ALWAYS_ASSERT( abs_err < compare_tol );
            }
 
            // The combined volumes of hi_eb_cc.areaHi() and hi_hi_eb_cc should
            // equal vfrac(iv_hi).
            { Real const vol = hi_eb_cc.volume() + hi_hi_eb_cc.volume();
              Real const abs_err = amrex::Math::abs(vfrac(iv_hi) - vol);
              Real compare_tol = Real(5.0e-6);
#ifndef AMREX_USE_GPU
              if ( abs_err >= compare_tol ) {
                hi_eb_cc.debug();
                hi_hi_eb_cc.debug();
                amrex::Print() << "\nFail: check-4 volume abs_err: " << abs_err
                               << "\n  point:  " << hi_point
                               << "\n  normal: " << hi_normal
                               << "\n     hi_eb_cc.volume() " <<    hi_eb_cc.volume()
                               << "\n  hi_hi_eb_cc.volume() " << hi_hi_eb_cc.volume()
                               << "\n  vfrac:   " << vfrac(iv_hi)
                               << '\n';
              } else if (verbose) {
                Print() << "Pass: hi_eb_cc + hi_hi_eb_cc = vfrac = " << vfrac(iv_hi)
                        << "  abs_err: " << abs_err << "\n";
              }
#endif
              AMREX_ALWAYS_ASSERT( abs_err < compare_tol );
            }
          } //
#endif
#endif // 0
 
          //-----------------------
          // Fill out aux_ arrays
          //-----------------------
 
          if (lo_eb_cc.isCovered() && hi_eb_cc.isCovered()) {
 
            // defaults to covered and disconnected.
 
          } else if (lo_eb_cc.isRegular() && hi_eb_cc.isRegular()) {
 
            aux_flag(i,j,k).setRegular();
            aux_flag(i,j,k).setConnected();
 
            aux_vfrac(i,j,k) = one;
 
            aux_afrac_x(i,j,k) = one;
            aux_afrac_y(i,j,k) = one;
            aux_afrac_z(i,j,k) = one;
 
            aux_fcent_x(i,j,k,0) = zero; aux_fcent_x(i,j,k,1) = zero;
            aux_fcent_y(i,j,k,0) = zero; aux_fcent_y(i,j,k,1) = zero;
            aux_fcent_z(i,j,k,0) = zero; aux_fcent_z(i,j,k,1) = zero;
 
            if (i==bx.bigEnd(0)) {
              aux_afrac_x(i+1,j,k) = one;
              aux_fcent_x(i+1,j,k,0) = zero; aux_fcent_x(i+1,j,k,1) = zero;
            }
            if (j==bx.bigEnd(1)) {
              aux_afrac_y(i,j+1,k) = one;
              aux_fcent_y(i,j+1,k,0) = zero; aux_fcent_y(i,j+1,k,1) = zero;
            }
            if (k==bx.bigEnd(2)) {
              aux_afrac_z(i,j,k+1) = one;
              aux_fcent_z(i,j,k+1,0) = zero; aux_fcent_z(i,j,k+1,1) = zero;
            }
 
          } else if ( (lo_eb_cc.isRegular() && hi_eb_cc.isCovered())
                   || (lo_eb_cc.isCovered() && hi_eb_cc.isRegular()) ) {
 
            // This is a problematic situation.
#ifndef AMREX_USE_GPU
            Print()<< "eb_aux_ / Check: Regular and Covered cut cells are facing each other." << std::endl;
#endif
 
          } else {
 
            // zero Cell Flag
 
            aux_flag(i,j,k).setSingleValued();
 
            // one Volume Fraction
 
            Real lo_vol {lo_eb_cc.volume()}; AMREX_ASSERT(lo_vol >= zero && lo_vol <= myhalf);
            Real hi_vol {hi_eb_cc.volume()}; AMREX_ASSERT(hi_vol >= zero && hi_vol <= myhalf);
 
            aux_vfrac(i,j,k) = lo_vol + hi_vol;
 
            // two Volume Centroid
 
            /* centVol() returns the coordinates based on m_rbx.
              The coordinates in the a_idim direction are in [zero,myhalf] for the low cell and in [-myhalf,zero] for the hi cell.
              Therefore, they need to be mapped to the eb_aux space, by shifting:
              x' = x - myhalf (low cell), x + myhalf (hi cell) if a_idim = 0
              y' = y - myhalf (low cell), y + myhalf (hi cell) if a_idim = 1
              z' = z - myhalf (low cell), z + myhalf (hi cell) if a_idim = 2
            */
 
            RealVect lo_vcent {lo_eb_cc.centVol()};
            RealVect hi_vcent {hi_eb_cc.centVol()};
 
            lo_vcent[a_idim] = lo_vcent[a_idim] - myhalf;
            hi_vcent[a_idim] = hi_vcent[a_idim] + myhalf;
 
            aux_vcent(i,j,k,0) = ( lo_vol * lo_vcent[0] + hi_vol * hi_vcent[0] ) / aux_vfrac(i,j,k);
            aux_vcent(i,j,k,1) = ( lo_vol * lo_vcent[1] + hi_vol * hi_vcent[1] ) / aux_vfrac(i,j,k);
            aux_vcent(i,j,k,2) = ( lo_vol * lo_vcent[2] + hi_vol * hi_vcent[2] ) / aux_vfrac(i,j,k);
 
            // three Area Fraction
 
            Real lo_areaLo_x {lo_eb_cc.areaLo(0)};
            Real lo_areaLo_y {lo_eb_cc.areaLo(1)};
            Real lo_areaLo_z {lo_eb_cc.areaLo(2)};
 
            Real hi_areaLo_x {hi_eb_cc.areaLo(0)};
            Real hi_areaLo_y {hi_eb_cc.areaLo(1)};
            Real hi_areaLo_z {hi_eb_cc.areaLo(2)};
 
            aux_afrac_x(i,j,k) = (a_idim == 0) ? lo_areaLo_x : lo_areaLo_x + hi_areaLo_x;
            aux_afrac_y(i,j,k) = (a_idim == 1) ? lo_areaLo_y : lo_areaLo_y + hi_areaLo_y;
            aux_afrac_z(i,j,k) = (a_idim == 2) ? lo_areaLo_z : lo_areaLo_z + hi_areaLo_z;
 
            if (i==bx.bigEnd(0)) {
              Real lo_areaHi_x {lo_eb_cc.areaHi(0)};
              Real hi_areaHi_x {hi_eb_cc.areaHi(0)};
              aux_afrac_x(i+1,j,k) = (a_idim == 0) ? hi_areaHi_x : lo_areaHi_x + hi_areaHi_x;
            }
            if (j==bx.bigEnd(1)) {
              Real lo_areaHi_y {lo_eb_cc.areaHi(1)};
              Real hi_areaHi_y {hi_eb_cc.areaHi(1)};
              aux_afrac_y(i,j+1,k) = (a_idim == 1) ? hi_areaHi_y : lo_areaHi_y + hi_areaHi_y;
            }
            if (k==bx.bigEnd(2)) {
              Real lo_areaHi_z {lo_eb_cc.areaHi(2)};
              Real hi_areaHi_z {hi_eb_cc.areaHi(2)};
              aux_afrac_z(i,j,k+1) = (a_idim == 2) ? hi_areaHi_z : lo_areaHi_z + hi_areaHi_z;
            }
 
            // Real(4.) Face Centroid
 
            /* fcentLo returns the coordinates based on m_rbx.
              The coordinates in the a_idim direction are in [zero,myhalf] for the low cell and in [-myhalf,zero] for the hi cell.
              Therefore, they need to be mapped to the eb_aux space, by shifting:
              x' = x - myhalf (low cell), x + myhalf (hi cell) if a_idim = 0
              y' = y - myhalf (low cell), y + myhalf (hi cell) if a_idim = 1
              z' = z - myhalf (low cell), z + myhalf (hi cell) if a_idim = 2
            */
 
            RealVect lo_centLo_x {lo_eb_cc.centLo(0)};
            RealVect lo_centLo_y {lo_eb_cc.centLo(1)};
            RealVect lo_centLo_z {lo_eb_cc.centLo(2)};
 
            RealVect hi_centLo_x {hi_eb_cc.centLo(0)};
            RealVect hi_centLo_y {hi_eb_cc.centLo(1)};
            RealVect hi_centLo_z {hi_eb_cc.centLo(2)};
 
            if (a_idim == 0) {
              aux_fcent_x(i,j,k,0) = lo_centLo_x[1];      // y
              aux_fcent_x(i,j,k,1) = lo_centLo_x[2];      // z
              aux_fcent_y(i,j,k,0) = (aux_afrac_y(i,j,k) > zero)   // x (mapped)
                                    ? ( lo_areaLo_y * (lo_centLo_y[0] - myhalf)
                                      + hi_areaLo_y * (hi_centLo_y[0] + myhalf) ) / aux_afrac_y(i,j,k)
                                    : zero;
              aux_fcent_y(i,j,k,1) = (aux_afrac_y(i,j,k) > zero)   // z
                                    ? ( lo_areaLo_y * lo_centLo_y[2]
                                      + hi_areaLo_y * hi_centLo_y[2] ) / aux_afrac_y(i,j,k)
                                    : zero;
              aux_fcent_z(i,j,k,0) = (aux_afrac_z(i,j,k) > zero)   // x (mapped)
                                    ? ( lo_areaLo_z * (lo_centLo_z[0] - myhalf)
                                      + hi_areaLo_z * (hi_centLo_z[0] + myhalf) ) / aux_afrac_z(i,j,k)
                                    : zero;
              aux_fcent_z(i,j,k,1) = (aux_afrac_z(i,j,k) > zero)   // y
                                    ? ( lo_areaLo_z * lo_centLo_z[1]
                                      + hi_areaLo_z * hi_centLo_z[1] ) / aux_afrac_z(i,j,k)
                                    : zero;
            } else if (a_idim == 1) {
              aux_fcent_x(i,j,k,0) = (aux_afrac_x(i,j,k) > zero)   // y (mapped)
                                    ? ( lo_areaLo_x * (lo_centLo_x[1] - myhalf)
                                      + hi_areaLo_x * (hi_centLo_x[1] + myhalf) ) / aux_afrac_x(i,j,k)
                                    : zero;
              aux_fcent_x(i,j,k,1) = (aux_afrac_x(i,j,k) > zero)   // z
                                    ? ( lo_areaLo_x * lo_centLo_x[2]
                                      + hi_areaLo_x * hi_centLo_x[2] ) / aux_afrac_x(i,j,k)
                                    : zero;
              aux_fcent_y(i,j,k,0) = lo_centLo_y[0];      // x
              aux_fcent_y(i,j,k,1) = lo_centLo_y[2];      // z
              aux_fcent_z(i,j,k,0) = (aux_afrac_z(i,j,k) > zero)   // x
                                    ? ( lo_areaLo_z * lo_centLo_z[0]
                                      + hi_areaLo_z * hi_centLo_z[0] ) / aux_afrac_z(i,j,k)
                                    : zero;
              aux_fcent_z(i,j,k,1) = (aux_afrac_z(i,j,k) > zero)   // y (mapped)
                                    ? ( lo_areaLo_z * (lo_centLo_z[1] - myhalf)
                                      + hi_areaLo_z * (hi_centLo_z[1] + myhalf) ) / aux_afrac_z(i,j,k)
                                    : zero;
            } else if (a_idim == 2) {
              aux_fcent_x(i,j,k,0) = (aux_afrac_x(i,j,k) > zero)   // y
                                    ? ( lo_areaLo_x * lo_centLo_x[1]
                                      + hi_areaLo_x * hi_centLo_x[1] ) / aux_afrac_x(i,j,k)
                                    : zero;
              aux_fcent_x(i,j,k,1) = (aux_afrac_x(i,j,k) > zero)   // z (mapped)
                                    ? ( lo_areaLo_x * (lo_centLo_x[2] - myhalf)
                                      + hi_areaLo_x * (hi_centLo_x[2] + myhalf) ) / aux_afrac_x(i,j,k)
                                    : zero;
              aux_fcent_y(i,j,k,0) = (aux_afrac_y(i,j,k) > zero)   // x
                                    ? ( lo_areaLo_y * lo_centLo_y[0]
                                      + hi_areaLo_y * hi_centLo_y[0] ) / aux_afrac_y(i,j,k)
                                    : zero;
              aux_fcent_y(i,j,k,1) = (aux_afrac_y(i,j,k) > zero)   // z (mapped)
                                    ? ( lo_areaLo_y * (lo_centLo_y[2] - myhalf)
                                      + hi_areaLo_y * (hi_centLo_y[2] + myhalf) ) / aux_afrac_y(i,j,k)
                                    : zero;
              aux_fcent_z(i,j,k,0) = lo_centLo_z[0];      // x
              aux_fcent_z(i,j,k,1) = lo_centLo_z[1];      // y
            }
 
            if (i==bx.bigEnd(0)) {
              Real lo_areaHi_x {lo_eb_cc.areaHi(0)};
              Real hi_areaHi_x {hi_eb_cc.areaHi(0)};
              RealVect lo_centHi_x {lo_eb_cc.centHi(0)};
              RealVect hi_centHi_x {hi_eb_cc.centHi(0)};
              if (a_idim == 0) {
                aux_fcent_x(i+1,j,k,0) = hi_centHi_x[1];      // y
                aux_fcent_x(i+1,j,k,1) = hi_centHi_x[2];      // z
              } else if (a_idim == 1) {
                aux_fcent_x(i+1,j,k,0) = (aux_afrac_x(i+1,j,k) > zero)   // y (mapped)
                                      ? ( lo_areaHi_x * (lo_centHi_x[1] - myhalf)
                                        + hi_areaHi_x * (hi_centHi_x[1] + myhalf) ) / aux_afrac_x(i+1,j,k)
                                      : zero;
                aux_fcent_x(i+1,j,k,1) = (aux_afrac_x(i+1,j,k) > zero)   // z
                                      ? ( lo_areaHi_x * lo_centHi_x[2]
                                        + hi_areaHi_x * hi_centHi_x[2] ) / aux_afrac_x(i+1,j,k)
                                      : zero;
              } else if (a_idim == 2) {
                aux_fcent_x(i+1,j,k,0) = (aux_afrac_x(i+1,j,k) > zero)   // y
                                      ? ( lo_areaHi_x * lo_centHi_x[1]
                                        + hi_areaHi_x * hi_centHi_x[1] ) / aux_afrac_x(i+1,j,k)
                                      : zero;
                aux_fcent_x(i+1,j,k,1) = (aux_afrac_x(i+1,j,k) > zero)   // z (mapped)
                                      ? ( lo_areaHi_x * (lo_centHi_x[2] - myhalf)
                                        + hi_areaHi_x * (hi_centHi_x[2] + myhalf) ) / aux_afrac_x(i+1,j,k)
                                      : zero;
              }
            }
            if (j==bx.bigEnd(1)) {
              Real lo_areaHi_y {lo_eb_cc.areaHi(1)};
              Real hi_areaHi_y {hi_eb_cc.areaHi(1)};
              RealVect lo_centHi_y {lo_eb_cc.centHi(1)};
              RealVect hi_centHi_y {hi_eb_cc.centHi(1)};
              if (a_idim == 0) {
                aux_fcent_y(i,j+1,k,0) = (aux_afrac_y(i,j+1,k) > zero)   // x (mapped)
                                      ? ( lo_areaHi_y * (lo_centHi_y[0] - myhalf)
                                        + hi_areaHi_y * (hi_centHi_y[0] + myhalf) ) / aux_afrac_y(i,j+1,k)
                                      : zero;
                aux_fcent_y(i,j+1,k,1) = (aux_afrac_y(i,j+1,k) > zero)   // z
                                      ? ( lo_areaHi_y * lo_centHi_y[2]
                                        + hi_areaHi_y * hi_centHi_y[2] ) / aux_afrac_y(i,j+1,k)
                                      : zero;
              } else if (a_idim == 1) {
                aux_fcent_y(i,j+1,k,0) = lo_centHi_y[0];      // x
                aux_fcent_y(i,j+1,k,1) = lo_centHi_y[2];      // z
              } else if (a_idim == 2) {
                aux_fcent_y(i,j+1,k,0) = (aux_afrac_y(i,j+1,k) > zero)   // x
                                      ? ( lo_areaHi_y * lo_centHi_y[0]
                                        + hi_areaHi_y * hi_centHi_y[0] ) / aux_afrac_y(i,j+1,k)
                                      : zero;
                aux_fcent_y(i,j+1,k,1) = (aux_afrac_y(i,j+1,k) > zero)   // z (mapped)
                                      ? ( lo_areaHi_y * (lo_centHi_y[2] - myhalf)
                                        + hi_areaHi_y * (hi_centHi_y[2] + myhalf) ) / aux_afrac_y(i,j+1,k)
                                      : zero;
              }
            }
            if (k==bx.bigEnd(2)) {
              Real lo_areaHi_z {lo_eb_cc.areaHi(2)};
              Real hi_areaHi_z {hi_eb_cc.areaHi(2)};
              RealVect lo_centHi_z {lo_eb_cc.centHi(2)};
              RealVect hi_centHi_z {hi_eb_cc.centHi(2)};
              if (a_idim == 0) {
                aux_fcent_z(i,j,k+1,0) = (aux_afrac_z(i,j,k+1) > zero)   // x (mapped)
                                      ? ( lo_areaHi_z * (lo_centHi_z[0] - myhalf)
                                        + hi_areaHi_z * (hi_centHi_z[0] + myhalf) ) / aux_afrac_z(i,j,k+1)
                                      : zero;
                aux_fcent_z(i,j,k+1,1) = (aux_afrac_z(i,j,k+1) > zero)   // y
                                      ? ( lo_areaHi_z * lo_centHi_z[1]
                                        + hi_areaHi_z * hi_centHi_z[1] ) / aux_afrac_z(i,j,k+1)
                                      : zero;
              } else if (a_idim == 1) {
                aux_fcent_z(i,j,k+1,0) = (aux_afrac_z(i,j,k+1) > zero)   // x
                                      ? ( lo_areaHi_z * lo_centHi_z[0]
                                        + hi_areaHi_z * hi_centHi_z[0] ) / aux_afrac_z(i,j,k+1)
                                      : zero;
                aux_fcent_z(i,j,k+1,1) = (aux_afrac_z(i,j,k+1) > zero)   // y (mapped)
                                      ? ( lo_areaHi_z * (lo_centHi_z[1] - myhalf)
                                        + hi_areaHi_z * (hi_centHi_z[1] + myhalf) ) / aux_afrac_z(i,j,k+1)
                                      : zero;
              } else if (a_idim == 2) {
                aux_fcent_z(i,j,k+1,0) = lo_centHi_z[0];      // x
                aux_fcent_z(i,j,k+1,1) = lo_centHi_z[1];      // y
              }
            }
 
            // Real(5.) Boundary Area
 
            Real lo_areaBoun {lo_eb_cc.areaBoun()};
            Real hi_areaBoun {hi_eb_cc.areaBoun()};
 
            aux_barea(i,j,k) = lo_areaBoun + hi_areaBoun;
 
            // Real(6.) Boundary Centroid
 
            RealVect lo_centBoun {lo_eb_cc.centBoun()};
            RealVect hi_centBoun {hi_eb_cc.centBoun()};
 
            if (a_idim == 0) {
              aux_bcent(i,j,k,0) = ( lo_areaBoun * (lo_centBoun[0]-myhalf) + hi_areaBoun * (hi_centBoun[0]+myhalf) ) / aux_barea(i,j,k);  // x (mapped)
              aux_bcent(i,j,k,1) = ( lo_areaBoun * lo_centBoun[1] + hi_areaBoun * hi_centBoun[1] ) / aux_barea(i,j,k);              // y
              aux_bcent(i,j,k,2) = ( lo_areaBoun * lo_centBoun[2] + hi_areaBoun * hi_centBoun[2] ) / aux_barea(i,j,k);              // z
            } else if (a_idim == 1) {
              aux_bcent(i,j,k,0) = ( lo_areaBoun * lo_centBoun[0] + hi_areaBoun * hi_centBoun[0] ) / aux_barea(i,j,k);              // x
              aux_bcent(i,j,k,1) = ( lo_areaBoun * (lo_centBoun[1]-myhalf) + hi_areaBoun * (hi_centBoun[1]+myhalf) ) / aux_barea(i,j,k);  // y (mapped)
              aux_bcent(i,j,k,2) = ( lo_areaBoun * lo_centBoun[2] + hi_areaBoun * hi_centBoun[2] ) / aux_barea(i,j,k);              // z
            } else if (a_idim == 2) {
              aux_bcent(i,j,k,0) = ( lo_areaBoun * lo_centBoun[0] + hi_areaBoun * hi_centBoun[0] ) / aux_barea(i,j,k);              // x
              aux_bcent(i,j,k,1) = ( lo_areaBoun * lo_centBoun[1] + hi_areaBoun * hi_centBoun[1] ) / aux_barea(i,j,k);              // y
              aux_bcent(i,j,k,2) = ( lo_areaBoun * (lo_centBoun[2]-myhalf) + hi_areaBoun * (hi_centBoun[2]+myhalf) ) / aux_barea(i,j,k);  // z (mapped)
            }
 
            // Real(7.) Boundary Normal
 
            RealVect eb_normal = ( lo_areaBoun * lo_normal + hi_areaBoun * hi_normal )/ aux_barea(i,j,k);
 
            aux_bnorm(i,j,k,0) = eb_normal[0];
            aux_bnorm(i,j,k,1) = eb_normal[1];
            aux_bnorm(i,j,k,2) = eb_normal[2];
 
          }
 
        } // flag(iv_lo) and flag(iv_hi)
 
      });
 
      ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
      {
        if (aux_vfrac(i,j,k) < small_volfrac)
        {
          aux_vfrac(i,j,k)   = zero;
        }
      });
 
    } // if (FlagFab[mfi].getType(bx) == FabType::singlevalued )
 
  } // MFIter
 
  // We FillBoundary volfrac here so that we can use tests on volfrac in ghost cells below
  m_volfrac->FillBoundary(a_geom.periodicity());
 
  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);
    Array4<Real>       const& aux_afrac_idim = m_areafrac[a_idim]->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 ) {
 
      // Corrections for small cells
      Box my_xbx(bx); my_xbx.growHi(0,1);
      int xbx_lo = my_xbx.smallEnd(0);
      int xbx_hi = my_xbx.bigEnd(0);
      ParallelFor(my_xbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
      {
        if ((i == xbx_lo && aux_vfrac(i,j,k) < small_volfrac) ||
            (i == xbx_hi && aux_vfrac(i-1,j,k) < small_volfrac) ||
            (i > xbx_lo && i < xbx_hi &&
            (aux_vfrac(i,j,k) < small_volfrac || aux_vfrac(i-1,j,k) < small_volfrac))) {
            aux_afrac_x(i,j,k) = zero;
        }
      });
 
      Box my_ybx(bx); my_ybx.growHi(1,1);
      int ybx_lo = my_ybx.smallEnd(1);
      int ybx_hi = my_ybx.bigEnd(1);
      ParallelFor(my_ybx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
      {
        if ((j == ybx_lo && aux_vfrac(i,j,k) < small_volfrac) ||
            (j == ybx_hi && aux_vfrac(i,j-1,k) < small_volfrac) ||
            (j > ybx_lo && j < ybx_hi &&
            (aux_vfrac(i,j,k) < small_volfrac || aux_vfrac(i,j-1,k) < small_volfrac))) {
            aux_afrac_y(i,j,k) = zero;
        }
      });
 
      Box my_zbx(bx); my_zbx.growHi(2,1);
      int zbx_lo = my_zbx.smallEnd(2);
      int zbx_hi = my_zbx.bigEnd(2);
      ParallelFor(my_zbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
      {
        if ((k == zbx_lo && aux_vfrac(i,j,k) < small_volfrac) ||
            (k == zbx_hi && aux_vfrac(i,j,k-1) < small_volfrac) ||
            (k > zbx_lo && k < zbx_hi &&
            (aux_vfrac(i,j,k) < small_volfrac || aux_vfrac(i,j,k-1) < small_volfrac))) {
            aux_afrac_z(i,j,k) = zero;
        }
      });
 
      ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
      {
        if (aux_vfrac(i,j,k) < small_volfrac)
        {
          aux_vcent(i,j,k,0) = zero;
          aux_vcent(i,j,k,1) = zero;
          aux_vcent(i,j,k,2) = zero;
 
          aux_fcent_x(i  ,j  ,k  ,0) = zero;
          aux_fcent_x(i  ,j  ,k  ,1) = zero;
          aux_fcent_x(i+1,j  ,k  ,0) = zero;
          aux_fcent_x(i+1,j  ,k  ,1) = zero;
 
          aux_fcent_y(i  ,j  ,k  ,0) = zero;
          aux_fcent_y(i  ,j  ,k  ,1) = zero;
          aux_fcent_y(i  ,j+1,k  ,0) = zero;
          aux_fcent_y(i  ,j+1,k  ,1) = zero;
 
          aux_fcent_z(i  ,j  ,k  ,0) = zero;
          aux_fcent_z(i  ,j  ,k  ,1) = zero;
          aux_fcent_z(i  ,j  ,k+1,0) = zero;
          aux_fcent_z(i  ,j  ,k+1,1) = zero;
 
          aux_barea(i,j,k) = zero;
 
          aux_bcent(i,j,k,0) = zero;
          aux_bcent(i,j,k,1) = zero;
          aux_bcent(i,j,k,2) = zero;
 
          aux_bnorm(i,j,k,0) = zero;
          aux_bnorm(i,j,k,1) = zero;
          aux_bnorm(i,j,k,2) = zero;
 
          aux_flag(i,j,k).setCovered();
        }
 
        if (std::abs(aux_vcent(i,j,k,0)) < small_value) aux_vcent(i,j,k,0) = zero;
        if (std::abs(aux_vcent(i,j,k,1)) < small_value) aux_vcent(i,j,k,1) = zero;
        if (std::abs(aux_vcent(i,j,k,2)) < small_value) aux_vcent(i,j,k,2) = zero;
        if (std::abs(aux_bcent(i,j,k,0)) < small_value) aux_bcent(i,j,k,0) = zero;
        if (std::abs(aux_bcent(i,j,k,1)) < small_value) aux_bcent(i,j,k,1) = zero;
        if (std::abs(aux_bcent(i,j,k,2)) < small_value) aux_bcent(i,j,k,2) = zero;
      });
 
      // Area fraction MultiFab has one more slice at bigEnd(idim),
      // and this slice is not filled by fillBoundary(), for higher levels.
      // (Lower level might be filled by fillBoundary().)
      // Fill the ghost region for the last slice at bigEnd(idim)
      // by the value of the nearest point. And let fillBoundary() overwrite it.
 
      Box upper_slab = makeSlab(bx_grown, a_idim, bx.bigEnd(a_idim)+1);
      Box bx_grown_1 = bx; bx_grown_1.grow(a_idim,1);
      BoxList slab_diffList = boxDiff(upper_slab, bx_grown_1);
 
      for (const Box& b : slab_diffList) {
        ParallelFor(b, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
        {
          IntVect iv(AMREX_D_DECL(i,j,k));
          IntVect iv_nearest = iv;
          for (int d=0; d<AMREX_SPACEDIM; ++d) {
              iv_nearest[d] = Clamp(iv[d], bx_grown_1.smallEnd(d), bx_grown_1.bigEnd(d));
          }
          aux_afrac_idim(iv) = aux_afrac_idim(iv_nearest);
        });
      }
 
    } // if (FlagFab[mfi].getType(bx) == FabType::singlevalued )
 
  } // MFIter
 
  // Fill Boundary
 
  // The FB call for volfrac is done above
  // m_volfrac->FillBoundary(a_geom.periodicity());
 
  m_volcent->FillBoundary(a_geom.periodicity());
  for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
    m_areafrac[idim]->FillBoundary(a_geom.periodicity());
    m_facecent[idim]->FillBoundary(a_geom.periodicity());
  }
  m_bndryarea->FillBoundary(a_geom.periodicity());
  m_bndrycent->FillBoundary(a_geom.periodicity());
  m_bndrynorm->FillBoundary(a_geom.periodicity());
 
  // Set Connectivities
  for (MFIter mfi(*m_cellflags, false); mfi.isValid(); ++mfi) {
 
    const Box& bx = mfi.validbox();
    const Box domain = surroundingNodes(a_geom.Domain(), a_idim);
 
    if (FlagFab[mfi].getType(bx) == FabType::singlevalued ) {
 
      Array4<EBCellFlag> const& aux_flag  = m_cellflags->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);
 
      ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
      {
        EB2::build_cellflag_from_ap (i, j, k, aux_flag, aux_afrac_x, aux_afrac_y, aux_afrac_z);
      });
 
      // Set disconnected non-periodicfaces
 
      bool l_periodic_x = a_geom.isPeriodic(0);
      bool l_periodic_y = a_geom.isPeriodic(1);
      bool l_periodic_z = a_geom.isPeriodic(2);
 
      if (!l_periodic_x) {
        const Box dom_grown = grow(grow(domain,1,1),2,1);
        const Box bx_grown  = grow(grow(    bx,1,1),2,1);
        const Box bx_face_x_lo = bx_grown & makeSlab(dom_grown,0,domain.smallEnd(0));
        const Box bx_face_x_hi = bx_grown & makeSlab(dom_grown,0,domain.bigEnd(0));
 
        ParallelFor(bx_face_x_lo, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
        {
          for(int kk(-1); kk<=1; kk++) {
          for(int jj(-1); jj<=1; jj++) {
            aux_flag(i,j,k).setDisconnected(-1,jj,kk);
          }}
        });
        ParallelFor(bx_face_x_hi, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
        {
          for(int kk(-1); kk<=1; kk++) {
          for(int jj(-1); jj<=1; jj++) {
            aux_flag(i,j,k).setDisconnected( 1,jj,kk);
          }}
        });
      }
 
      if (!l_periodic_y) {
        const Box dom_grown = grow(grow(domain,0,1),2,1);
        const Box bx_grown  = grow(grow(    bx,0,1),2,1);
        const Box bx_face_y_lo = bx_grown & makeSlab(dom_grown,1,domain.smallEnd(1));
        const Box bx_face_y_hi = bx_grown & makeSlab(dom_grown,1,domain.bigEnd(1));
 
        ParallelFor(bx_face_y_lo, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
        {
          for(int kk(-1); kk<=1; kk++) {
          for(int ii(-1); ii<=1; ii++) {
            aux_flag(i,j,k).setDisconnected(ii,-1,kk);
          }}
        });
        ParallelFor(bx_face_y_hi, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
        {
          for(int kk(-1); kk<=1; kk++) {
          for(int ii(-1); ii<=1; ii++) {
            aux_flag(i,j,k).setDisconnected(ii, 1,kk);
          }}
        });
      }
 
      if (!l_periodic_z) {
        const Box dom_grown = grow(grow(domain,0,1),1,1);
        const Box bx_grown  = grow(grow(    bx,0,1),1,1);
        const Box bx_face_z_lo = bx_grown & makeSlab(dom_grown,2,domain.smallEnd(2));
        const Box bx_face_z_hi = bx_grown & makeSlab(dom_grown,2,domain.bigEnd(2));
 
        ParallelFor(bx_face_z_lo, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
        {
          for(int jj(-1); jj<=1; jj++) {
          for(int ii(-1); ii<=1; ii++) {
            aux_flag(i,j,k).setDisconnected(ii,jj,-1);
          }}
        });
        ParallelFor(bx_face_z_hi, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
        {
          for(int jj(-1); jj<=1; jj++) {
          for(int ii(-1); ii<=1; ii++) {
            aux_flag(i,j,k).setDisconnected(ii,jj, 1);
          }}
        });
      }
 
    } // FabType::singlevalued
 
  } // MFIter
 
  // Set disconnected zero-volume-fraction cells
  // (equivalent to eb_::set_connection_flags for CC grids)
 
  for (MFIter mfi(*m_cellflags, false); mfi.isValid(); ++mfi) {
 
    const Box& bx = mfi.validbox();
    const Box gbx = amrex::grow(bx, m_cellflags->nGrow()-1); // Leave one cell layer
 
    Array4<EBCellFlag> const& aux_flag  = m_cellflags->array(mfi);
    Array4<Real>       const& aux_vfrac = m_volfrac->array(mfi);
 
    ParallelFor(gbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
    {
      for(int kk(-1); kk<=1; kk++) {
      for(int jj(-1); jj<=1; jj++) {
      for(int ii(-1); ii<=1; ii++)
      {
        if (aux_vfrac(i+ii,j+jj,k+kk) == zero) {
            aux_flag(i,j,k).setDisconnected(ii,jj,kk);
        }
      }}}
    });
 
    ParallelFor(gbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
    {
        if (aux_vfrac(i,j,k)==zero) {
            aux_flag(i,j,k).setCovered();
        }
    });
 
  } // MFIter
 
  // Fill Boundary
 
  m_cellflags->FillBoundary(a_geom.periodicity());
 
#if 0
  // We leave these here for debugging if necessary.
  // If you uncomment these, make sure to uncomment AMReX_VisMF include above
  if (a_idim == 0) {
      amrex::VisMF::Write(*m_volfrac,"UVOL");
      amrex::VisMF::Write(*m_areafrac[0],"UAREAX");
      amrex::VisMF::Write(*m_areafrac[1],"UAREAY");
      amrex::VisMF::Write(*m_areafrac[2],"UAREAZ");
  } else if (a_idim == 1) {
      amrex::VisMF::Write(*m_volfrac,"VVOL");
      amrex::VisMF::Write(*m_areafrac[0],"VAREAX");
      amrex::VisMF::Write(*m_areafrac[1],"VAREAY");
      amrex::VisMF::Write(*m_areafrac[2],"VAREAZ");
  } else {
      amrex::VisMF::Write(*m_volfrac,"WVOL");
      amrex::VisMF::Write(*m_areafrac[0],"WAREAX");
      amrex::VisMF::Write(*m_areafrac[1],"WAREAY");
      amrex::VisMF::Write(*m_areafrac[2],"WAREAZ");
  }
#endif
}

Referenced by eb_::make_all_factories().

Here is the call graph for this function:

Here is the caller graph for this function:

getAreaFrac()

Array< const MultiFab *, AMREX_SPACEDIM > eb_aux_::getAreaFrac

(

)

const

Return reconstructed area fractions for all coordinate directions.

{
    AMREX_ASSERT(m_areafrac[0] != nullptr);
    return {AMREX_D_DECL(m_areafrac[0], m_areafrac[1], m_areafrac[2])};
}

getBndryArea()

const MultiFab & eb_aux_::getBndryArea

(

)

const

Return the reconstructed EB boundary areas.

{
    AMREX_ASSERT(m_bndryarea != nullptr);
    return *m_bndryarea;
}

getBndryCent()

const MultiFab & eb_aux_::getBndryCent

(

)

const

Return the reconstructed EB boundary centroids.

{
    AMREX_ASSERT(m_bndrycent != nullptr);
    return *m_bndrycent;
}

getBndryNormal()

const MultiFab & eb_aux_::getBndryNormal

(

)

const

Return the reconstructed EB boundary normals.

{
    AMREX_ASSERT(m_bndrynorm != nullptr);
    return *m_bndrynorm;
}

getCentroid()

const MultiFab & eb_aux_::getCentroid

(

)

const

Return the reconstructed volume centroids.

{
    AMREX_ASSERT(m_volcent != nullptr);
    return *m_volcent;
}

getFaceCent()

Array< const MultiFab *, AMREX_SPACEDIM > eb_aux_::getFaceCent

(

)

const

Return reconstructed face centroids for all coordinate directions.

{
    AMREX_ASSERT(m_facecent[0] != nullptr);
    return {AMREX_D_DECL(m_facecent[0], m_facecent[1], m_facecent[2])};
}

getMultiEBCellFlagFab()

const FabArray< EBCellFlagFab > & eb_aux_::getMultiEBCellFlagFab

(

)

const

Return the reconstructed EB cell flags.

{
    AMREX_ASSERT(m_cellflags != nullptr);
    return *m_cellflags;
}

Referenced by AdvectionSrcForMom_EB(), compute_gradp_xy(), compute_gradp_z(), and DiffusionSrcForMom_EB().

Here is the caller graph for this function:

getVolFrac()

const MultiFab & eb_aux_::getVolFrac

(

)

const

Return the reconstructed volume fractions.

{
    AMREX_ASSERT(m_volfrac != nullptr);
    return *m_volfrac;
}

set_verbose()

void eb_aux_::set_verbose ( )

inline

Enable verbose diagnostic output during factory construction.

{ m_verbose = 1; }

Member Data Documentation

m_areafrac

amrex::Array<amrex::MultiFab*,AMREX_SPACEDIM> eb_aux_::m_areafrac {{AMREX_D_DECL(nullptr, nullptr, nullptr)}}

private

Referenced by define(), and getAreaFrac().

m_bndryarea

amrex::MultiFab* eb_aux_::m_bndryarea = nullptr

private

Referenced by define(), and getBndryArea().

m_bndrycent

amrex::MultiFab* eb_aux_::m_bndrycent = nullptr

private

Referenced by define(), and getBndryCent().

m_bndrynorm

amrex::MultiFab* eb_aux_::m_bndrynorm = nullptr

private

Referenced by define(), and getBndryNormal().

m_cellflags

amrex::FabArray<amrex::EBCellFlagFab>* eb_aux_::m_cellflags = nullptr

private

Referenced by define(), and getMultiEBCellFlagFab().

m_facecent

amrex::Array<amrex::MultiFab*,AMREX_SPACEDIM> eb_aux_::m_facecent {{AMREX_D_DECL(nullptr, nullptr, nullptr)}}

private

Referenced by define(), and getFaceCent().

m_verbose

int eb_aux_::m_verbose

private

Referenced by define(), and set_verbose().

m_volcent

amrex::MultiFab* eb_aux_::m_volcent = nullptr

private

Referenced by define(), and getCentroid().

m_volfrac

amrex::MultiFab* eb_aux_::m_volfrac = nullptr

private

Referenced by define(), and getVolFrac().

---

The documentation for this class was generated from the following files:

• Source/EB/ERF_EBAux.H
• Source/EB/ERF_EBAux.cpp