Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

EB outside domain #1201

Merged
merged 2 commits into from
Jul 27, 2020
Merged

EB outside domain #1201

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
bc8b72a
EB outside domain
WeiqunZhang Jul 26, 2020
22d87aa
fix fillFab for EB_levelset
WeiqunZhang Jul 27, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
7 changes: 7 additions & 0 deletions Src/Base/AMReX_Algorithm.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -55,6 +55,13 @@ namespace amrex
t1 = std::move(t2);
t2 = std::move(temp);
}

template <typename T>
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
constexpr const T& Clamp (const T& v, const T& lo, const T& hi)
{
return (v < lo) ? lo : (hi < v) ? hi : v;
}
}

#endif
11 changes: 8 additions & 3 deletions Src/EB/AMReX_EB2.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -64,7 +64,8 @@ public:

IndexSpaceImp (const G& gshop, const Geometry& geom,
int required_coarsening_level, int max_coarsening_level,
int ngrow, bool build_coarse_level_by_coarsening);
int ngrow, bool build_coarse_level_by_coarsening,
bool extend_domain_face);

IndexSpaceImp (IndexSpaceImp<G> const&) = delete;
IndexSpaceImp (IndexSpaceImp<G> &&) = delete;
Expand Down Expand Up @@ -92,17 +93,21 @@ private:

#include <AMReX_EB2_IndexSpaceI.H>

bool ExtendDomainFace ();

template <typename G>
void
Build (const G& gshop, const Geometry& geom,
int required_coarsening_level, int max_coarsening_level,
int ngrow = 4, bool build_coarse_level_by_coarsening = true)
int ngrow = 4, bool build_coarse_level_by_coarsening = true,
bool extend_domain_face = ExtendDomainFace())
{
BL_PROFILE("EB2::Initialize()");
IndexSpace::push(new IndexSpaceImp<G>(gshop, geom,
required_coarsening_level,
max_coarsening_level,
ngrow, build_coarse_level_by_coarsening));
ngrow, build_coarse_level_by_coarsening,
extend_domain_face));
}

void Build (const Geometry& geom,
Expand Down
7 changes: 7 additions & 0 deletions Src/EB/AMReX_EB2.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -18,11 +18,13 @@ namespace amrex { namespace EB2 {
Vector<std::unique_ptr<IndexSpace> > IndexSpace::m_instance;

int max_grid_size = 64;
bool extend_domain_face = false;

void Initialize ()
{
ParmParse pp("eb2");
pp.query("max_grid_size", max_grid_size);
pp.query("extend_domain_face", extend_domain_face);

amrex::ExecOnFinalize(Finalize);
}
Expand All @@ -32,6 +34,11 @@ void Finalize ()
IndexSpace::clear();
}

bool ExtendDomainFace ()
{
return extend_domain_face;
}

void
IndexSpace::push (IndexSpace* ispace)
{
Expand Down
56 changes: 34 additions & 22 deletions Src/EB/AMReX_EB2_GeometryShop.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -299,43 +299,52 @@ public:
static constexpr bool isGPUable () noexcept { return false; }

template <class U=F, typename std::enable_if<IsGPUable<U>::value>::type* FOO = nullptr >
void fillFab (BaseFab<Real>& levelset, const Geometry& geom, RunOn run_on) const noexcept
void fillFab (BaseFab<Real>& levelset, const Geometry& geom, RunOn run_on,
Box const& bounding_box) const noexcept
{
const auto problo = geom.ProbLoArray();
const auto dx = geom.CellSizeArray();
const Box& bx = levelset.box();
const auto& a = levelset.array();
const auto blo = amrex::lbound(bounding_box);
const auto bhi = amrex::ubound(bounding_box);
auto f = m_f;
AMREX_HOST_DEVICE_FOR_3D_FLAG(run_on, bx, i, j, k,
{
a(i,j,k) = f(AMREX_D_DECL(problo[0]+i*dx[0],
problo[1]+j*dx[1],
problo[2]+k*dx[2]));
a(i,j,k) = f(AMREX_D_DECL(problo[0]+amrex::Clamp(i,blo.x,bhi.x)*dx[0],
problo[1]+amrex::Clamp(j,blo.y,bhi.y)*dx[1],
problo[2]+amrex::Clamp(k,blo.z,bhi.z)*dx[2]));
});
}

template <class U=F, typename std::enable_if<!IsGPUable<U>::value>::type* BAR = nullptr >
void fillFab (BaseFab<Real>& levelset, const Geometry& geom, RunOn) const noexcept
void fillFab (BaseFab<Real>& levelset, const Geometry& geom, RunOn,
Box const& bounding_box) const noexcept
{
const auto problo = geom.ProbLoArray();
const auto dx = geom.CellSizeArray();
const Box& bx = levelset.box();
const auto& a = levelset.array();
const auto blo = amrex::lbound(bounding_box);
const auto bhi = amrex::ubound(bounding_box);
amrex::LoopOnCpu(bx, [&] (int i, int j, int k) noexcept
{
a(i,j,k) = m_f(RealArray{AMREX_D_DECL(problo[0]+i*dx[0],
problo[1]+j*dx[1],
problo[2]+k*dx[2])});
a(i,j,k) = m_f(RealArray{AMREX_D_DECL(problo[0]+amrex::Clamp(i,blo.x,bhi.x)*dx[0],
problo[1]+amrex::Clamp(j,blo.y,bhi.y)*dx[1],
problo[2]+amrex::Clamp(k,blo.z,bhi.z)*dx[2])});
});
}

template <class U=F, typename std::enable_if<IsGPUable<U>::value>::type* FOO = nullptr >
void getIntercept (Array<BaseFab<Real>,AMREX_SPACEDIM>& inter_fab,
Array<BaseFab<Type_t>,AMREX_SPACEDIM> const& type_fab,
Geometry const& geom, RunOn run_on) const noexcept
Geometry const& geom, RunOn run_on,
Box const& bounding_box) const noexcept
{
auto const& dx = geom.CellSizeArray();
auto const& problo = geom.ProbLoArray();
const auto blo = amrex::lbound(bounding_box);
const auto bhi = amrex::ubound(bounding_box);
auto f = m_f;
for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
Array4<Real> const& inter = inter_fab[idim].array();
Expand All @@ -349,12 +358,12 @@ public:
IntVect ivhi(AMREX_D_DECL(i,j,k));
ivhi[dir] += 1;
inter(i,j,k) = BrentRootFinder
({AMREX_D_DECL(problo[0]+ivlo[0]*dx[0],
problo[1]+ivlo[1]*dx[1],
problo[2]+ivlo[2]*dx[2])},
{AMREX_D_DECL(problo[0]+ivhi[0]*dx[0],
problo[1]+ivhi[1]*dx[1],
problo[2]+ivhi[2]*dx[2])},
({AMREX_D_DECL(problo[0]+amrex::Clamp(ivlo[0],blo.x,bhi.x)*dx[0],
problo[1]+amrex::Clamp(ivlo[1],blo.y,bhi.y)*dx[1],
problo[2]+amrex::Clamp(ivlo[2],blo.z,bhi.z)*dx[2])},
{AMREX_D_DECL(problo[0]+amrex::Clamp(ivhi[0],blo.x,bhi.x)*dx[0],
problo[1]+amrex::Clamp(ivhi[1],blo.y,bhi.y)*dx[1],
problo[2]+amrex::Clamp(ivhi[2],blo.z,bhi.z)*dx[2])},
dir, f);
} else {
inter(i,j,k) = std::numeric_limits<Real>::quiet_NaN();
Expand All @@ -366,10 +375,13 @@ public:
template <class U=F, typename std::enable_if<!IsGPUable<U>::value>::type* BAR = nullptr >
void getIntercept (Array<BaseFab<Real>,AMREX_SPACEDIM>& inter_fab,
Array<BaseFab<Type_t>,AMREX_SPACEDIM> const& type_fab,
Geometry const& geom, RunOn) const noexcept
Geometry const& geom, RunOn,
Box const& bounding_box) const noexcept
{
auto const& dx = geom.CellSizeArray();
auto const& problo = geom.ProbLoArray();
const auto blo = amrex::lbound(bounding_box);
const auto bhi = amrex::ubound(bounding_box);
for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
Array4<Real> const& inter = inter_fab[idim].array();
Array4<Type_t const> const& type = type_fab[idim].array();
Expand All @@ -382,12 +394,12 @@ public:
IntVect ivhi(AMREX_D_DECL(i,j,k));
ivhi[dir] += 1;
inter(i,j,k) = BrentRootFinder
({AMREX_D_DECL(problo[0]+ivlo[0]*dx[0],
problo[1]+ivlo[1]*dx[1],
problo[2]+ivlo[2]*dx[2])},
{AMREX_D_DECL(problo[0]+ivhi[0]*dx[0],
problo[1]+ivhi[1]*dx[1],
problo[2]+ivhi[2]*dx[2])},
({AMREX_D_DECL(problo[0]+amrex::Clamp(ivlo[0],blo.x,bhi.x)*dx[0],
problo[1]+amrex::Clamp(ivlo[1],blo.y,bhi.y)*dx[1],
problo[2]+amrex::Clamp(ivlo[2],blo.z,bhi.z)*dx[2])},
{AMREX_D_DECL(problo[0]+amrex::Clamp(ivhi[0],blo.x,bhi.x)*dx[0],
problo[1]+amrex::Clamp(ivhi[1],blo.y,bhi.y)*dx[1],
problo[2]+amrex::Clamp(ivhi[2],blo.z,bhi.z)*dx[2])},
dir, m_f);
} else {
inter(i,j,k) = std::numeric_limits<Real>::quiet_NaN();
Expand Down
7 changes: 4 additions & 3 deletions Src/EB/AMReX_EB2_IndexSpaceI.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,8 @@ template <typename G>
IndexSpaceImp<G>::IndexSpaceImp (const G& gshop, const Geometry& geom,
int required_coarsening_level,
int max_coarsening_level,
int ngrow, bool build_coarse_level_by_coarsening)
int ngrow, bool build_coarse_level_by_coarsening,
bool extend_domain_face)
{
// build finest level (i.e., level 0) first
AMREX_ALWAYS_ASSERT(required_coarsening_level >= 0 && required_coarsening_level <= 30);
Expand All @@ -19,7 +20,7 @@ IndexSpaceImp<G>::IndexSpaceImp (const G& gshop, const Geometry& geom,
m_domain.push_back(geom.Domain());
m_ngrow.push_back(ngrow_finest);
m_gslevel.reserve(max_coarsening_level+1);
m_gslevel.emplace_back(this, gshop, geom, EB2::max_grid_size, ngrow_finest);
m_gslevel.emplace_back(this, gshop, geom, EB2::max_grid_size, ngrow_finest, extend_domain_face);

for (int ilev = 1; ilev <= max_coarsening_level; ++ilev)
{
Expand All @@ -43,7 +44,7 @@ IndexSpaceImp<G>::IndexSpaceImp (const G& gshop, const Geometry& geom,
if (build_coarse_level_by_coarsening) {
amrex::Abort("Failed to build required coarse EB level "+std::to_string(ilev));
} else {
m_gslevel.emplace_back(this, gshop, cgeom, EB2::max_grid_size, ng);
m_gslevel.emplace_back(this, gshop, cgeom, EB2::max_grid_size, ng, extend_domain_face);
}
} else {
break;
Expand Down
20 changes: 14 additions & 6 deletions Src/EB/AMReX_EB2_Level.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -94,14 +94,14 @@ class GShopLevel
: public Level
{
public:
GShopLevel (IndexSpace const* is, G const& gshop, const Geometry& geom, int max_grid_size, int ngrow);
GShopLevel (IndexSpace const* is, G const& gshop, const Geometry& geom, int max_grid_size, int ngrow, bool extend_domain_face);
GShopLevel (IndexSpace const* is, int ilev, int max_grid_size, int ngrow,
const Geometry& geom, GShopLevel<G>& fineLevel);
};

template <typename G>
GShopLevel<G>::GShopLevel (IndexSpace const* is, G const& gshop, const Geometry& geom,
int max_grid_size, int ngrow)
int max_grid_size, int ngrow, bool extend_domain_face)
: Level(is, geom)
{
if (std::is_same<typename G::FunctionType, AllRegularIF>::value) {
Expand Down Expand Up @@ -130,6 +130,8 @@ GShopLevel<G>::GShopLevel (IndexSpace const* is, G const& gshop, const Geometry&
}
}
domain_grown.grow(m_ngrow);
Box bounding_box = (extend_domain_face) ? geom.Domain() : domain_grown;
bounding_box.surroundingNodes();

m_grids.define(domain_grown);
m_grids.maxSize(max_grid_size);
Expand All @@ -142,7 +144,7 @@ GShopLevel<G>::GShopLevel (IndexSpace const* is, G const& gshop, const Geometry&
{
const Box& vbx = mfi.validbox();
const Box& gbx = amrex::surroundingNodes(amrex::grow(vbx,1));
int box_type = gshop.getBoxType(gbx, geom, RunOn::Gpu);
int box_type = gshop.getBoxType(gbx & bounding_box, geom, RunOn::Gpu);
if (box_type == gshop.allcovered) {
covered_boxes.push_back(vbx);
} else if (box_type == gshop.mixedcells) {
Expand Down Expand Up @@ -194,6 +196,12 @@ GShopLevel<G>::GShopLevel (IndexSpace const* is, G const& gshop, const Geometry&
const auto dx = geom.CellSizeArray();
const auto problo = geom.ProbLoArray();

for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
if (!extend_domain_face || geom.isPeriodic(idim)) {
bounding_box.grow(idim,GFab::ng);
}
}

RunOn gshop_run_on = (Gpu::inLaunchRegion() && gshop.isGPUable())
? RunOn::Gpu : RunOn::Cpu;

Expand All @@ -213,7 +221,7 @@ GShopLevel<G>::GShopLevel (IndexSpace const* is, G const& gshop, const Geometry&
const Box& vbx = gfab.validbox();

auto& levelset = gfab.getLevelSet();
gshop.fillFab(levelset, geom, gshop_run_on);
gshop.fillFab(levelset, geom, gshop_run_on, bounding_box);

if (hybrid) levelset.prefetchToDevice();

Expand Down Expand Up @@ -264,7 +272,7 @@ GShopLevel<G>::GShopLevel (IndexSpace const* is, G const& gshop, const Geometry&
}
}

gshop.getIntercept(intercept, edgetype, geom, gshop_run_on);
gshop.getIntercept(intercept, edgetype, geom, gshop_run_on, bounding_box);

if (hybrid) {
for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
Expand Down Expand Up @@ -313,7 +321,7 @@ GShopLevel<G>::GShopLevel (IndexSpace const* is, G const& gshop, const Geometry&
}
}

gshop.getIntercept(intercept, facetype, geom, gshop_run_on);
gshop.getIntercept(intercept, facetype, geom, gshop_run_on, bounding_box);

if (hybrid) {
for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
Expand Down
14 changes: 12 additions & 2 deletions Src/EB/AMReX_EB_levelset.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -392,11 +392,21 @@ class GShopLSFactory {
std::unique_ptr<MultiFab> mf_impfunc = std::unique_ptr<MultiFab>(new MultiFab);
mf_impfunc->define(m_ls_ba, m_ls_dm, 1, m_ls_pad);

Box bounding_box = m_geom.Domain();
bounding_box.surroundingNodes();
bool extend_domain_face = EB2::ExtendDomainFace();
for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
if (!extend_domain_face || m_geom.isPeriodic(idim)) {
bounding_box.grow(m_ls_pad);
}
}

#ifdef _OPENMP
#pragma omp parallel if (Gpu::notInLaunchRegion())
#endif
for(MFIter mfi(* mf_impfunc); mfi.isValid(); ++ mfi)
m_gshop.fillFab((* mf_impfunc)[mfi], m_geom, RunOn::Gpu);
for(MFIter mfi(* mf_impfunc); mfi.isValid(); ++ mfi) {
m_gshop.fillFab((* mf_impfunc)[mfi], m_geom, RunOn::Gpu, bounding_box);
}

mf_impfunc->FillBoundary(m_geom.periodicity());

Expand Down