Skip to content

Commit

Permalink
Bounded sampling (AMReX-Codes#4195)
Browse files Browse the repository at this point in the history
The `get_line_data` and `get_slice_data` routines do not currently allow
for a natural way to bound the returned MultiFab. In the former case,
the line length is dictated by the bottom and top of the boxarray at
that level. In the latter case, the plane is dictated by the problem
domain at that level.

This PR adds an optional `Box` argument to the `get_line_data` call and
an optional `RealBox` argument to the `get_slice_data` call. These
optional arguments bound the returned MuliFab. This feature has
application in the ERF code in [PR
1897](erf-model/ERF#1897) where it is of great
value to sample velocities in planes/lines just upstream of a wind
turbine. Since the wind turbine is small in relation to the problem
domain, the bounding feature ensures the returned data is actually what
the user wants.
  • Loading branch information
AMLattanzi authored Oct 21, 2024
1 parent b00c828 commit c333708
Show file tree
Hide file tree
Showing 2 changed files with 41 additions and 21 deletions.
43 changes: 29 additions & 14 deletions Src/Base/AMReX_MultiFabUtil.H
Original file line number Diff line number Diff line change
Expand Up @@ -169,7 +169,8 @@ namespace amrex
std::unique_ptr<MultiFab> get_slice_data(int dir, Real coord,
const MultiFab& cc,
const Geometry& geom, int start_comp, int ncomp,
bool interpolate=false);
bool interpolate=false,
RealBox const& bnd_rbx = RealBox());

/**
* \brief Get data in a cell of MultiFab/FabArray
Expand All @@ -188,7 +189,7 @@ namespace amrex
* specified by a direction and a cell.
*/
template <typename MF, std::enable_if_t<IsFabArray<MF>::value,int> FOO = 0>
MF get_line_data (MF const& mf, int dir, IntVect const& cell);
MF get_line_data (MF const& mf, int dir, IntVect const& cell, Box const& bnd_bx = Box());

//! Return an iMultiFab that has the same BoxArray and DistributionMapping
//! as the coarse MultiFab cmf. Cells covered by the coarsened fine grids
Expand Down Expand Up @@ -996,26 +997,40 @@ Vector<typename MF::value_type> get_cell_data (MF const& mf, IntVect const& cell
}

template <typename MF, std::enable_if_t<IsFabArray<MF>::value,int> FOO>
MF get_line_data (MF const& mf, int dir, IntVect const& cell)
MF get_line_data (MF const& mf, int dir, IntVect const& cell, Box const& bnd_bx)
{
bool do_bnd = (!bnd_bx.isEmpty());

BoxArray const& ba = mf.boxArray();
DistributionMapping const& dm = mf.DistributionMap();
const auto nboxes = static_cast<int>(ba.size());

BoxList bl(ba.ixType());
Vector<int> procmap;
Vector<int> index_map;
for (int i = 0; i < nboxes; ++i) {
Box const& b = ba[i];
IntVect lo = cell;
lo[dir] = b.smallEnd(dir);
if (b.contains(lo)) {
IntVect hi = lo;
hi[dir] = b.bigEnd(dir);
Box b1d(lo,hi,b.ixType());
bl.push_back(b1d);
procmap.push_back(dm[i]);
index_map.push_back(i);
if (!do_bnd) {
for (int i = 0; i < nboxes; ++i) {
Box const& b = ba[i];
IntVect lo = cell;
lo[dir] = b.smallEnd(dir);
if (b.contains(lo)) {
IntVect hi = lo;
hi[dir] = b.bigEnd(dir);
Box b1d(lo,hi,b.ixType());
bl.push_back(b1d);
procmap.push_back(dm[i]);
index_map.push_back(i);
}
}
} else {
for (int i = 0; i < nboxes; ++i) {
Box const& b = ba[i];
Box const& b1d = bnd_bx & b;
if (b1d.ok()) {
bl.push_back(b1d);
procmap.push_back(dm[i]);
index_map.push_back(i);
}
}
}

Expand Down
19 changes: 12 additions & 7 deletions Src/Base/AMReX_MultiFabUtil.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,7 @@ namespace {
using namespace amrex;

Box
getIndexBox(const RealBox& real_box, const Geometry& geom) {
getIndexBox (const RealBox& real_box, const Geometry& geom) {
IntVect slice_lo, slice_hi;

AMREX_D_TERM(slice_lo[0]=static_cast<int>(std::floor((real_box.lo(0) - geom.ProbLo(0))/geom.CellSize(0)));,
Expand All @@ -24,12 +24,11 @@ namespace {
}


std::unique_ptr<MultiFab> allocateSlice(int dir, const MultiFab& cell_centered_data,
int ncomp, const Geometry& geom, Real dir_coord,
Vector<int>& slice_to_full_ba_map) {
std::unique_ptr<MultiFab> allocateSlice (int dir, const MultiFab& cell_centered_data,
int ncomp, const Geometry& geom, Real dir_coord,
Vector<int>& slice_to_full_ba_map, RealBox real_slice) {

// Get our slice and convert to index space
RealBox real_slice = geom.ProbDomain();
real_slice.setLo(dir, dir_coord);
real_slice.setHi(dir, dir_coord);
Box slice_box = getIndexBox(real_slice, geom);
Expand Down Expand Up @@ -550,7 +549,7 @@ namespace amrex
return amrex::cast<FabArray<BaseFab<Long> > > (imf);
}

std::unique_ptr<MultiFab> get_slice_data(int dir, Real coord, const MultiFab& cc, const Geometry& geom, int start_comp, int ncomp, bool interpolate) {
std::unique_ptr<MultiFab> get_slice_data(int dir, Real coord, const MultiFab& cc, const Geometry& geom, int start_comp, int ncomp, bool interpolate, RealBox const& bnd_rbx) {

BL_PROFILE("amrex::get_slice_data");

Expand All @@ -559,9 +558,15 @@ namespace amrex
}

const auto geomdata = geom.data();
RealBox real_slice;
if (bnd_rbx.ok()) {
real_slice = bnd_rbx;
} else {
real_slice = geom.ProbDomain();
}

Vector<int> slice_to_full_ba_map;
std::unique_ptr<MultiFab> slice = allocateSlice(dir, cc, ncomp, geom, coord, slice_to_full_ba_map);
std::unique_ptr<MultiFab> slice = allocateSlice(dir, cc, ncomp, geom, coord, slice_to_full_ba_map, real_slice);

if (!slice) {
return nullptr;
Expand Down

0 comments on commit c333708

Please sign in to comment.