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](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.

Author: AMLattanzi

Src/Base/AMReX_MultiFabUtil.H

@@ -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
@@ -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
@@ -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);
+            }
         }
     }
 

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<int>(std::floor((real_box.lo(0) - geom.ProbLo(0))/geom.CellSize(0)));,
@@ -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);
@@ -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");
 
@@ -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;