Skip to content

Commit

Permalink
EB Data outside domain
Browse files Browse the repository at this point in the history 
We do not guarantee there is valid EB information outside the domain, unless
the user specifies a large enough number of ghost cells when building
EB. However, this makes things difficult for application codes with AMR,
because the number of ghost cells needed depends on the number of AMR
levels. In this commit, we try to fill some basic EB information
(EBCellFlag, volume fraction and area fraction) for all cells outside the
domain.
  • Loading branch information
@WeiqunZhang
WeiqunZhang committed Sep 17, 2023
1 parent 6eb91be commit 03afc90
Show file tree
Hide file tree
Showing 3 changed files with 152 additions and 1 deletion.
2 changes: 2 additions & 0 deletions Src/EB/AMReX_EB2_Level.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -62,6 +62,8 @@ public:
const Geometry& Geom () const noexcept { return m_geom; }
IndexSpace const* getEBIndexSpace () const noexcept { return m_parent; }

IntVect const& nGrowVect () const noexcept { return m_ngrow; }

void write_to_chkpt_file (const std::string& fname, bool extend_domain_face, int max_grid_size) const;

bool hasEBInfo () const noexcept { return m_has_eb_info; }
Expand Down
3 changes: 3 additions & 0 deletions Src/EB/AMReX_EBDataCollection.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -42,6 +42,9 @@ public:
[[nodiscard]] Array<const MultiCutFab*, AMREX_SPACEDIM> getEdgeCent () const;
[[nodiscard]] const iMultiFab* getCutCellMask () const;

// public for cuda
void extendDataOutsideDomain (IntVect const& level_ng);

private:

Vector<int> m_ngrow;
Expand Down
148 changes: 147 additions & 1 deletion Src/EB/AMReX_EBDataCollection.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,7 @@
#include <AMReX_MultiCutFab.H>

#include <AMReX_EB2_Level.H>
#include <algorithm>
#include <utility>

namespace amrex {
Expand All @@ -13,7 +14,7 @@ EBDataCollection::EBDataCollection (const EB2::Level& a_level,
const Geometry& a_geom,
const BoxArray& a_ba_in,
const DistributionMapping& a_dm,
Vector<int> a_ngrow, EBSupport a_support)
Vector<int> a_ngrow, EBSupport a_support)
: m_ngrow(std::move(a_ngrow)),
m_support(a_support),
m_geom(a_geom)
Expand All @@ -33,6 +34,8 @@ EBDataCollection::EBDataCollection (const EB2::Level& a_level,

if (m_support >= EBSupport::volume)
{
AMREX_ALWAYS_ASSERT(m_ngrow[1] <= m_ngrow[0]);

m_volfrac = new MultiFab(a_ba, a_dm, 1, m_ngrow[1], MFInfo(), FArrayBoxFactory());
a_level.fillVolFrac(*m_volfrac, m_geom);

Expand All @@ -42,6 +45,8 @@ EBDataCollection::EBDataCollection (const EB2::Level& a_level,

if (m_support == EBSupport::full)
{
AMREX_ALWAYS_ASSERT(m_ngrow[2] <= m_ngrow[0]);

const int ng = m_ngrow[2];

m_bndrycent = new MultiCutFab(a_ba, a_dm, AMREX_SPACEDIM, ng, *m_cellflags);
Expand Down Expand Up @@ -72,6 +77,147 @@ EBDataCollection::EBDataCollection (const EB2::Level& a_level,
DefaultFabFactory<IArrayBox>());
a_level.fillCutCellMask(*m_cutcellmask, m_geom);
}

extendDataOutsideDomain(a_level.nGrowVect());
}

void EBDataCollection::extendDataOutsideDomain (IntVect const& level_ng)
{
if (m_cellflags == nullptr) { return; }

int const ngrow = m_ngrow[0];
Box const& data_domain = amrex::grow(m_geom.Domain(), ngrow);

Box level_domain = m_geom.Domain();
for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
if (m_geom.isPeriodic(idim)) {
level_domain.grow(idim, ngrow);
} else {
level_domain.grow(idim, level_ng[idim]);
}
}

if (level_domain.contains(data_domain)) { return; }

Box const& level_nodal_domain = amrex::surroundingNodes(level_domain);
Array<Box,AMREX_SPACEDIM> lev_ap_domain
{AMREX_D_DECL(amrex::surroundingNodes(level_domain,0),
amrex::surroundingNodes(level_domain,1),
amrex::surroundingNodes(level_domain,2))};

#ifdef AMREX_USE_OMP
#pragma omp parallel if (Gpu::notInLaunchRegion())
#endif
for (MFIter mfi(*m_cellflags); mfi.isValid(); ++mfi) {
Box const& bx = mfi.fabbox();
if (! level_domain.contains(bx)) {
Box const& nbx = amrex::surroundingNodes(bx);
auto const& ls_a = m_levelset->array(mfi);
auto const& flag_a = m_cellflags->array(mfi);
Array4<Real> vfrc_a;
if (m_volfrac) {
vfrc_a = m_volfrac->array(mfi);
}
amrex::ParallelFor(nbx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
{
if (! level_nodal_domain.contains(i,j,k)) {
int ii = std::clamp(i, level_nodal_domain.smallEnd(0),
level_nodal_domain.bigEnd(0));
int jj = std::clamp(j, level_nodal_domain.smallEnd(1),
level_nodal_domain.bigEnd(1));
#if (AMREX_SPACEDIM > 2)
int kk = std::clamp(k, level_nodal_domain.smallEnd(2),
level_nodal_domain.bigEnd(2));
#else
int kk = 0;
#endif
ls_a(i,j,k) = ls_a(ii,jj,kk);
}
});
amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
{
if (! level_domain.contains(i,j,k)) {
EBCellFlag flag;
int ii = std::clamp(i, level_domain.smallEnd(0),
level_domain.bigEnd(0));
int jj = std::clamp(j, level_domain.smallEnd(1),
level_domain.bigEnd(1));
#if (AMREX_SPACEDIM > 2)
int kk = std::clamp(k, level_domain.smallEnd(2),
level_domain.bigEnd(2));
#else
int kk = 0;
#endif
if (flag_a(ii,jj,kk).isCovered()) {
flag.setCovered();
} else if (flag_a(ii,jj,kk).isRegular()) {
flag.setRegular();
} else {
int ncov = 0;
#if (AMREX_SPACEDIM > 2)
for (int ko = k; ko <= k+1; ++ko) {
#else
for (int ko = k; ko <= k ; ++ko) {
#endif
for (int jo = j; jo <= j+1; ++jo) {
for (int io = i; io <= i+1; ++io) {
if (ls_a(io,jo,ko) >= Real(0.0)) { ++ncov; }
}}}
constexpr int nnodes = (AMREX_SPACEDIM == 2) ? 4 : 8;
if (ncov == nnodes) {
flag.setCovered();
} else if (ncov == 0) {
flag.setRegular();
} else {
flag.setSingleValued();
}
}
if (flag.isCovered()) {
flag_a(i,j,k).setCovered();
if (vfrc_a && vfrc_a.contains(i,j,k)) {
vfrc_a(i,j,k) = Real(0.0);
}
} else if (flag.isRegular()) {
flag_a(i,j,k).setRegular();
if (vfrc_a && vfrc_a.contains(i,j,k)) {
vfrc_a(i,j,k) = Real(1.0);
}
} else {
flag_a(i,j,k).setSingleValued();
if (vfrc_a && vfrc_a.contains(i,j,k)) {
vfrc_a(i,j,k) = vfrc_a(ii,jj,kk);
}
}
}
});

for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
if (m_areafrac[idim]) {
auto const& ap = m_areafrac[idim]->array(mfi);
Box apbx = Box(ap);
if (apbx.smallEnd(idim) == nbx.smallEnd(idim)) {
apbx.growLo(idim,-1);
}
if (apbx.bigEnd(idim) == nbx.bigEnd(idim)) {
apbx.growHi(idim,-1);
}
auto lev_apidim_domain = lev_ap_domain[idim];
Dim3 off = IntVect::TheDimensionVector(idim).dim3();
amrex::ParallelFor(apbx,
[=] AMREX_GPU_DEVICE (int i, int j, int k)
{
if (! lev_apidim_domain.contains(i,j,k)) {
if (flag_a(i-off.x,j-off.y,k-off.z).isCovered() ||
flag_a(i,j,k).isCovered())
{
ap(i,j,k) = Real(0.0);
}
}
});
}
}
}
}
}

EBDataCollection::~EBDataCollection ()
Expand Down

0 comments on commit 03afc90

Please sign in to comment.