From c333708d59f01ab664363ba426e85ad24c1fb23d Mon Sep 17 00:00:00 2001 From: "Aaron M. Lattanzi" <103702284+AMLattanzi@users.noreply.github.com> Date: Mon, 21 Oct 2024 12:35:31 -0700 Subject: [PATCH] Bounded sampling (#4195) 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](https://github.com/erf-model/ERF/pull/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. --- Src/Base/AMReX_MultiFabUtil.H | 43 ++++++++++++++++++++++----------- Src/Base/AMReX_MultiFabUtil.cpp | 19 +++++++++------ 2 files changed, 41 insertions(+), 21 deletions(-) diff --git a/Src/Base/AMReX_MultiFabUtil.H b/Src/Base/AMReX_MultiFabUtil.H index 228070a13c..0740069c68 100644 --- a/Src/Base/AMReX_MultiFabUtil.H +++ b/Src/Base/AMReX_MultiFabUtil.H @@ -169,7 +169,8 @@ namespace amrex std::unique_ptr 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 @@ -188,7 +189,7 @@ namespace amrex * specified by a direction and a cell. */ template ::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 @@ -996,8 +997,10 @@ Vector get_cell_data (MF const& mf, IntVect const& cell } template ::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(ba.size()); @@ -1005,17 +1008,29 @@ MF get_line_data (MF const& mf, int dir, IntVect const& cell) BoxList bl(ba.ixType()); Vector procmap; Vector 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); + } } } diff --git a/Src/Base/AMReX_MultiFabUtil.cpp b/Src/Base/AMReX_MultiFabUtil.cpp index 86a1e29054..721919f509 100644 --- a/Src/Base/AMReX_MultiFabUtil.cpp +++ b/Src/Base/AMReX_MultiFabUtil.cpp @@ -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(std::floor((real_box.lo(0) - geom.ProbLo(0))/geom.CellSize(0)));, @@ -24,12 +24,11 @@ namespace { } - std::unique_ptr allocateSlice(int dir, const MultiFab& cell_centered_data, - int ncomp, const Geometry& geom, Real dir_coord, - Vector& slice_to_full_ba_map) { + std::unique_ptr allocateSlice (int dir, const MultiFab& cell_centered_data, + int ncomp, const Geometry& geom, Real dir_coord, + Vector& 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); @@ -550,7 +549,7 @@ namespace amrex return amrex::cast > > (imf); } - std::unique_ptr get_slice_data(int dir, Real coord, const MultiFab& cc, const Geometry& geom, int start_comp, int ncomp, bool interpolate) { + std::unique_ptr 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"); @@ -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 slice_to_full_ba_map; - std::unique_ptr slice = allocateSlice(dir, cc, ncomp, geom, coord, slice_to_full_ba_map); + std::unique_ptr slice = allocateSlice(dir, cc, ncomp, geom, coord, slice_to_full_ba_map, real_slice); if (!slice) { return nullptr;