Skip to content

Commit

Permalink
New EB optimization parameter: eb2.num_coarsen_opt
Browse files Browse the repository at this point in the history 
At the beginning of EB generation, we chop the entire finest domain into
boxes and find out the type of the boxes.  We then collect the completely
covered boxes and cut boxes into two BoxArrays.  This process can be costly
because of the number of calls to the implicit functions.  In this commit,
we have introduced a new ParmParse parameter, eb2.num_coarsen_opt with a
default value of zero.  If for instance it is set to 3, we start the box
type categorization at a resolution that is coarsened by a factor of 2^3.
For the provisional cut boxes, we refine them by a factor of 2, Then we chop
them into small boxes and categorize the new boxes.  This process is
performed recursively until we are at the original finest resolution.

The users should be aware that, if eb2.num_coaren_opt is too big, this could
produce in erroneous results because evaluating the implicit function on
coarse boxes could miss fine structures in the EB.

Thank Robert Marskar for sharing this algorithm.
  • Loading branch information
@WeiqunZhang
WeiqunZhang committed Jul 7, 2022
1 parent 557aae8 commit 057cf52
Show file tree
Hide file tree
Showing 10 changed files with 119 additions and 67 deletions.
13 changes: 9 additions & 4 deletions Src/EB/AMReX_EB2.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -66,7 +66,7 @@ public:
IndexSpaceImp (const G& gshop, const Geometry& geom,
int required_coarsening_level, int max_coarsening_level,
int ngrow, bool build_coarse_level_by_coarsening,
bool extend_domain_face);
bool extend_domain_face, int num_coarsen_opt);

IndexSpaceImp (IndexSpaceImp<G> const&) = delete;
IndexSpaceImp (IndexSpaceImp<G> &&) = delete;
Expand Down Expand Up @@ -95,27 +95,32 @@ private:
#include <AMReX_EB2_IndexSpaceI.H>

bool ExtendDomainFace ();
int NumCoarsenOpt ();

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,
bool extend_domain_face = ExtendDomainFace())
bool extend_domain_face = ExtendDomainFace(),
int num_coarsen_opt = NumCoarsenOpt())
{
BL_PROFILE("EB2::Initialize()");
IndexSpace::push(new IndexSpaceImp<G>(gshop, geom,
required_coarsening_level,
max_coarsening_level,
ngrow, build_coarse_level_by_coarsening,
extend_domain_face));
extend_domain_face,
num_coarsen_opt));
}

void Build (const Geometry& geom,
int required_coarsening_level,
int max_coarsening_level,
int ngrow = 4,
bool build_coarse_level_by_coarsening = true);
bool build_coarse_level_by_coarsening = true,
bool extend_domain_face = ExtendDomainFace(),
int num_coarsen_opt = NumCoarsenOpt());

int maxCoarseningLevel (const Geometry& geom);
int maxCoarseningLevel (IndexSpace const* ebis, const Geometry& geom);
Expand Down
34 changes: 25 additions & 9 deletions Src/EB/AMReX_EB2.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,12 +21,14 @@ AMREX_EXPORT Vector<std::unique_ptr<IndexSpace> > IndexSpace::m_instance;

AMREX_EXPORT int max_grid_size = 64;
AMREX_EXPORT bool extend_domain_face = true;
AMREX_EXPORT int num_coarsen_opt = 0;

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

amrex::ExecOnFinalize(Finalize);
}
Expand All @@ -41,6 +43,11 @@ bool ExtendDomainFace ()
return extend_domain_face;
}

int NumCoarsenOpt ()
{
return num_coarsen_opt;
}

void
IndexSpace::push (IndexSpace* ispace)
{
Expand Down Expand Up @@ -74,7 +81,8 @@ const IndexSpace* TopIndexSpaceIfPresent() noexcept {

void
Build (const Geometry& geom, int required_coarsening_level,
int max_coarsening_level, int ngrow, bool build_coarse_level_by_coarsening)
int max_coarsening_level, int ngrow, bool build_coarse_level_by_coarsening,
bool a_extend_domain_face, int a_num_coarsen_opt)
{
ParmParse pp("eb2");
std::string geom_type;
Expand All @@ -85,7 +93,8 @@ Build (const Geometry& geom, int required_coarsening_level,
EB2::AllRegularIF rif;
EB2::GeometryShop<EB2::AllRegularIF> gshop(rif);
EB2::Build(gshop, geom, required_coarsening_level,
max_coarsening_level, ngrow, build_coarse_level_by_coarsening);
max_coarsening_level, ngrow, build_coarse_level_by_coarsening,
a_extend_domain_face, a_num_coarsen_opt);
}
else if (geom_type == "box")
{
Expand All @@ -102,7 +111,8 @@ Build (const Geometry& geom, int required_coarsening_level,

EB2::GeometryShop<EB2::BoxIF> gshop(bf);
EB2::Build(gshop, geom, required_coarsening_level,
max_coarsening_level, ngrow, build_coarse_level_by_coarsening);
max_coarsening_level, ngrow, build_coarse_level_by_coarsening,
a_extend_domain_face, a_num_coarsen_opt);
}
else if (geom_type == "cylinder")
{
Expand All @@ -127,7 +137,8 @@ Build (const Geometry& geom, int required_coarsening_level,

EB2::GeometryShop<EB2::CylinderIF> gshop(cf);
EB2::Build(gshop, geom, required_coarsening_level,
max_coarsening_level, ngrow, build_coarse_level_by_coarsening);
max_coarsening_level, ngrow, build_coarse_level_by_coarsening,
a_extend_domain_face, a_num_coarsen_opt);
}
else if (geom_type == "plane")
{
Expand All @@ -141,7 +152,8 @@ Build (const Geometry& geom, int required_coarsening_level,

EB2::GeometryShop<EB2::PlaneIF> gshop(pf);
EB2::Build(gshop, geom, required_coarsening_level,
max_coarsening_level, ngrow, build_coarse_level_by_coarsening);
max_coarsening_level, ngrow, build_coarse_level_by_coarsening,
a_extend_domain_face, a_num_coarsen_opt);
}
else if (geom_type == "sphere")
{
Expand All @@ -158,7 +170,8 @@ Build (const Geometry& geom, int required_coarsening_level,

EB2::GeometryShop<EB2::SphereIF> gshop(sf);
EB2::Build(gshop, geom, required_coarsening_level,
max_coarsening_level, ngrow, build_coarse_level_by_coarsening);
max_coarsening_level, ngrow, build_coarse_level_by_coarsening,
a_extend_domain_face, a_num_coarsen_opt);
}
else if (geom_type == "torus")
{
Expand All @@ -177,7 +190,8 @@ Build (const Geometry& geom, int required_coarsening_level,

EB2::GeometryShop<EB2::TorusIF> gshop(sf);
EB2::Build(gshop, geom, required_coarsening_level,
max_coarsening_level, ngrow, build_coarse_level_by_coarsening);
max_coarsening_level, ngrow, build_coarse_level_by_coarsening,
a_extend_domain_face, a_num_coarsen_opt);
}
else if (geom_type == "parser")
{
Expand All @@ -188,7 +202,8 @@ Build (const Geometry& geom, int required_coarsening_level,
EB2::ParserIF pif(parser.compile<3>());
EB2::GeometryShop<EB2::ParserIF,Parser> gshop(pif,parser);
EB2::Build(gshop, geom, required_coarsening_level,
max_coarsening_level, ngrow, build_coarse_level_by_coarsening);
max_coarsening_level, ngrow, build_coarse_level_by_coarsening,
a_extend_domain_face, a_num_coarsen_opt);
}
else if (geom_type == "stl")
{
Expand All @@ -206,7 +221,8 @@ Build (const Geometry& geom, int required_coarsening_level,
geom, required_coarsening_level,
max_coarsening_level, ngrow,
build_coarse_level_by_coarsening,
extend_domain_face));
a_extend_domain_face,
a_num_coarsen_opt));
}
else
{
Expand Down
8 changes: 5 additions & 3 deletions Src/EB/AMReX_EB2_IndexSpaceI.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,7 @@ 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,
bool extend_domain_face)
bool extend_domain_face, int num_coarsen_opt)
{
// build finest level (i.e., level 0) first
AMREX_ALWAYS_ASSERT(required_coarsening_level >= 0 && required_coarsening_level <= 30);
Expand All @@ -20,7 +20,8 @@ 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, extend_domain_face);
m_gslevel.emplace_back(this, gshop, geom, EB2::max_grid_size, ngrow_finest, extend_domain_face,
num_coarsen_opt);

for (int ilev = 1; ilev <= max_coarsening_level; ++ilev)
{
Expand All @@ -44,7 +45,8 @@ 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, extend_domain_face);
m_gslevel.emplace_back(this, gshop, cgeom, EB2::max_grid_size, ng, extend_domain_face,
num_coarsen_opt-ilev);
}
} else {
break;
Expand Down
2 changes: 1 addition & 1 deletion Src/EB/AMReX_EB2_IndexSpace_STL.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,7 +19,7 @@ public:
const Geometry& geom, int required_coarsening_level,
int max_coarsening_level, int ngrow,
bool build_coarse_level_by_coarsening,
bool extend_domain_face);
bool extend_domain_face, int num_coarsen_opt);

IndexSpaceSTL (IndexSpaceSTL const&) = delete;
IndexSpaceSTL (IndexSpaceSTL &&) = delete;
Expand Down
6 changes: 3 additions & 3 deletions Src/EB/AMReX_EB2_IndexSpace_STL.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@ IndexSpaceSTL::IndexSpaceSTL (const std::string& stl_file, Real stl_scale,
const Geometry& geom, int required_coarsening_level,
int max_coarsening_level, int ngrow,
bool build_coarse_level_by_coarsening,
bool extend_domain_face)
bool extend_domain_face, int num_coarsen_opt)
{
Gpu::LaunchSafeGuard lsg(true); // Always use GPU

Expand All @@ -29,7 +29,7 @@ IndexSpaceSTL::IndexSpaceSTL (const std::string& stl_file, Real stl_scale,
m_ngrow.push_back(ngrow_finest);
m_stllevel.reserve(max_coarsening_level+1);
m_stllevel.emplace_back(this, stl_tools, geom, EB2::max_grid_size, ngrow_finest,
extend_domain_face);
extend_domain_face, num_coarsen_opt);

for (int ilev = 1; ilev <= max_coarsening_level; ++ilev)
{
Expand All @@ -54,7 +54,7 @@ IndexSpaceSTL::IndexSpaceSTL (const std::string& stl_file, Real stl_scale,
amrex::Abort("Failed to build required coarse EB level "+std::to_string(ilev));
} else {
m_stllevel.emplace_back(this, stl_tools, cgeom, EB2::max_grid_size, ng,
extend_domain_face);
extend_domain_face, num_coarsen_opt-ilev);
}
} else {
break;
Expand Down
115 changes: 71 additions & 44 deletions Src/EB/AMReX_EB2_Level.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -98,12 +98,13 @@ class GShopLevel
: public Level
{
public:
GShopLevel (IndexSpace const* is, G const& gshop, const Geometry& geom, int max_grid_size, int ngrow, bool extend_domain_face);
GShopLevel (IndexSpace const* is, G const& gshop, const Geometry& geom, int max_grid_size,
int ngrow, bool extend_domain_face, int num_crse_opt);
GShopLevel (IndexSpace const* is, int ilev, int max_grid_size, int ngrow,
const Geometry& geom, GShopLevel<G>& fineLevel);
GShopLevel (IndexSpace const* is, const Geometry& geom);
void define_fine (G const& gshop, const Geometry& geom,
int max_grid_size, int ngrow, bool extend_domain_face);
int max_grid_size, int ngrow, bool extend_domain_face, int num_crse_opt);
};

template <typename G>
Expand All @@ -113,7 +114,7 @@ GShopLevel<G>::GShopLevel (IndexSpace const* is, const Geometry& geom)

template <typename G>
GShopLevel<G>::GShopLevel (IndexSpace const* is, G const& gshop, const Geometry& geom,
int max_grid_size, int ngrow, bool extend_domain_face)
int max_grid_size, int ngrow, bool extend_domain_face, int num_crse_opt)
: Level(is, geom)
{
if (std::is_same<typename G::FunctionType, AllRegularIF>::value) {
Expand All @@ -122,13 +123,13 @@ GShopLevel<G>::GShopLevel (IndexSpace const* is, G const& gshop, const Geometry&
return;
}

define_fine(gshop, geom, max_grid_size, ngrow, extend_domain_face);
define_fine(gshop, geom, max_grid_size, ngrow, extend_domain_face, num_crse_opt);
}

template <typename G>
void
GShopLevel<G>::define_fine (G const& gshop, const Geometry& geom,
int max_grid_size, int ngrow, bool extend_domain_face)
int max_grid_size, int ngrow, bool extend_domain_face, int num_crse_opt)
{
if (amrex::Verbose() > 0 && extend_domain_face == false) {
amrex::Print() << "AMReX WARNING: extend_domain_face=false is not recommended!\n";
Expand Down Expand Up @@ -166,65 +167,91 @@ GShopLevel<G>::define_fine (G const& gshop, const Geometry& geom,
Box bounding_box = (extend_domain_face) ? domain : domain_grown;
bounding_box.surroundingNodes();

BoxList bl(domain);
bl.maxSize(max_grid_size);
if (m_ngrow != 0) {
const IntVect& domlo = domain.smallEnd();
const IntVect& domhi = domain.bigEnd();
for (auto& b : bl) {
for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
if (m_ngrow[idim] != 0) {
if (b.smallEnd(idim) == domlo[idim]) {
b.growLo(idim,m_ngrow[idim]);
}
if (b.bigEnd(idim) == domhi[idim]) {
b.growHi(idim,m_ngrow[idim]);
}
}
BoxList cut_boxes;
BoxList covered_boxes;

const int nprocs = ParallelDescriptor::NProcs();
const int iproc = ParallelDescriptor::MyProc();

num_crse_opt = std::max(0,std::min(8,num_crse_opt));
for (int clev = num_crse_opt; clev >= 0; --clev) {
IntVect crse_ratio(1 << clev);
if (domain.coarsenable(crse_ratio)) {
Box const& crse_bounding_box = amrex::coarsen(bounding_box, crse_ratio);
Geometry const& crse_geom = amrex::coarsen(geom, crse_ratio);
BoxList test_boxes;
if (cut_boxes.isEmpty()) {
covered_boxes.clear();
test_boxes = BoxList(crse_geom.Domain());
test_boxes.maxSize(max_grid_size);
} else {
test_boxes.swap(cut_boxes);
test_boxes.coarsen(crse_ratio);
test_boxes.maxSize(max_grid_size);
}
}
}

m_grids.define(std::move(bl));
m_dmap.define(m_grids);

Vector<Box> cut_boxes;
Vector<Box> covered_boxes;
const Long nboxes = test_boxes.size();
const auto& boxes = test_boxes.data();
for (Long i = iproc; i < nboxes; i += nprocs) {
const Box& vbx = boxes[i];
const Box& gbx = amrex::surroundingNodes(amrex::grow(vbx,1));
auto box_type = gshop.getBoxType(gbx&crse_bounding_box,crse_geom,RunOn::Gpu);
if (box_type == gshop.allcovered) {
covered_boxes.push_back(amrex::refine(vbx, crse_ratio));
} else if (box_type == gshop.mixedcells) {
cut_boxes.push_back(amrex::refine(vbx, crse_ratio));
}
}

for (MFIter mfi(m_grids, m_dmap); mfi.isValid(); ++mfi)
{
const Box& vbx = mfi.validbox();
const Box& gbx = amrex::surroundingNodes(amrex::grow(vbx,1));
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) {
cut_boxes.push_back(vbx);
amrex::AllGatherBoxes(cut_boxes.data());
amrex::AllGatherBoxes(covered_boxes.data());
}
}

amrex::AllGatherBoxes(cut_boxes);
amrex::AllGatherBoxes(covered_boxes);
if (m_ngrow != 0) {
auto grow_at_domain_boundary = [&] (BoxList& bl)
{
const IntVect& domlo = domain.smallEnd();
const IntVect& domhi = domain.bigEnd();
for (auto& b : bl) {
for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
if (m_ngrow[idim] != 0) {
if (b.smallEnd(idim) == domlo[idim]) {
b.growLo(idim,m_ngrow[idim]);
}
if (b.bigEnd(idim) == domhi[idim]) {
b.growHi(idim,m_ngrow[idim]);
}
}
}
}
};
grow_at_domain_boundary(covered_boxes);
grow_at_domain_boundary(cut_boxes);
}

if ( cut_boxes.empty() &&
!covered_boxes.empty())
if ( cut_boxes.isEmpty() &&
!covered_boxes.isEmpty())
{
amrex::Abort("AMReX_EB2_Level.H: Domain is completely covered");
}

if (!covered_boxes.empty()) {
m_covered_grids = BoxArray(BoxList(std::move(covered_boxes)));
if (!covered_boxes.isEmpty()) {
if (num_crse_opt > 2) { // don't want the box too big
covered_boxes.maxSize(max_grid_size*4);
}
m_covered_grids = BoxArray(std::move(covered_boxes));
}

if (cut_boxes.empty()) {
if (cut_boxes.isEmpty()) {
m_grids = BoxArray();
m_dmap = DistributionMapping();
m_allregular = true;
m_ok = true;
return;
}

m_grids = BoxArray(BoxList(std::move(cut_boxes)));
m_grids = BoxArray(std::move(cut_boxes));
m_dmap = DistributionMapping(m_grids);

m_mgf.define(m_grids, m_dmap);
Expand Down
Loading

0 comments on commit 057cf52

Please sign in to comment.