Add MultiFab::sum(region) and sum_unique(region)

Src/Base/AMReX_MultiFab.H

@@ -238,17 +238,28 @@ public:
     * \brief Returns the sum of component "comp" over the MultiFab -- no ghost cells are included.
     */
     [[nodiscard]] Real sum (int comp = 0, bool local = false) const;
+    /**
+    * \brief Returns the sum of component "comp" in the given "region". -- no ghost cells are included.
+    */
+    [[nodiscard]] Real sum (Box const& region, int comp = 0, bool local = false) const;
 
     using FabArray<FArrayBox>::sum;
 
     /**
     * \brief Same as sum with \p local =false, but for non-cell-centered data, this
-    * skips non-unique points that are owned by multiple boxes.
+    * only adds non-unique points that are owned by multiple boxes once.
     */
     [[nodiscard]] Real sum_unique (int comp = 0,
                      bool local = false,
                      const Periodicity& period = Periodicity::NonPeriodic()) const;
     /**
+    * \brief Returns the unique sum of component "comp" in the given
+    * region. Non-unique points owned by multiple boxes in the MultiFab are
+    * only added once. No ghost cells are included. This function does not take
+    * periodicity into account in the determination of uniqueness of points.
+    */
+    [[nodiscard]] Real sum_unique (Box const& region, int comp = 0, bool local = false) const;
+    /**
     * \brief Adds the scalar value \p val to the value of each cell in the
     * specified subregion of the MultiFab.
     *

Src/Base/AMReX_MultiFab.cpp

@@ -1209,14 +1209,56 @@ MultiFab::sum (int comp, bool local) const
     return FabArray<FArrayBox>::sum(comp, IntVect(0), local);
 }
 
+Real
+MultiFab::sum (Box const& region, int comp, bool local) const
+{
+    BL_PROFILE("MultiFab::sum(region)");
+
+    auto sm = 0.0_rt;
+
+#ifdef AMREX_USE_GPU
+    if (Gpu::inLaunchRegion()) {
+        auto const& ma = this->const_arrays();
+        sm = ParReduce(TypeList<ReduceOpSum>{}, TypeList<Real>{}, *this, IntVect(0),
+        [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) noexcept -> GpuTuple<Real>
+        {
+            return (region.contains(i,j,k)) ? ma[box_no](i,j,k,comp) : 0.0_rt;
+        });
+    } else
+#endif
+    {
+#ifdef AMREX_USE_OMP
+#pragma omp parallel if (!system::regtest_reduction) reduction(+:sm)
+#endif
+        for (MFIter mfi(*this,true); mfi.isValid(); ++mfi) {
+            Box const& bx = mfi.tilebox() & region;
+            if (bx.ok()) {
+                auto const& a = this->const_array(mfi);
+                auto tmp = 0.0_rt;
+                AMREX_LOOP_3D(bx, i, j, k,
+                {
+                    tmp += a(i,j,k,comp);
+                });
+                sm += tmp; // Do it this way so that it does not break regression tests.
+            }
+        }
+    }
+
+    if (!local) {
+        ParallelAllReduce::Sum(sm, ParallelContext::CommunicatorSub());
+    }
+
+    return sm;
+}
+
 Real
 MultiFab::sum_unique (int comp,
                       bool local,
                       const Periodicity& period) const
 {
     BL_PROFILE("MultiFab::sum_unique()");
 
-    // no duplicatly distributed points if cell centered
+    // no duplicately distributed points if cell centered
     if (ixType().cellCentered()) {
         return this->sum(comp, local);
     }
@@ -1262,6 +1304,60 @@ MultiFab::sum_unique (int comp,
     return sm;
 }
 
+Real
+MultiFab::sum_unique (Box const& region, int comp, bool local) const
+{
+    BL_PROFILE("MultiFab::sum_unique(region)");
+
+    // no duplicately distributed points if cell centered
+    if (ixType().cellCentered()) {
+        return this->sum(region, comp, local);
+    }
+
+    // Owner is the grid with the lowest grid number containing the data
+    std::unique_ptr<iMultiFab> owner_mask = OwnerMask();
+
+    auto sm = 0.0_rt;
+#ifdef AMREX_USE_GPU
+    if (Gpu::inLaunchRegion()) {
+        auto const& ma = this->const_arrays();
+        auto const& msk = owner_mask->const_arrays();
+        sm = ParReduce(TypeList<ReduceOpSum>{}, TypeList<Real>{}, *this, IntVect(0),
+        [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) noexcept
+                       -> GpuTuple<Real>
+        {
+            return (region.contains(i,j,k) && msk[box_no](i,j,k))
+                ? ma[box_no](i,j,k,comp) : 0.0_rt;
+        });
+    } else
+#endif
+    {
+#ifdef AMREX_USE_OMP
+#pragma omp parallel if (!system::regtest_reduction) reduction(+:sm)
+#endif
+        for (MFIter mfi(*this,true); mfi.isValid(); ++mfi)
+        {
+            Box const& bx = mfi.tilebox() & region;
+            if (bx.ok()) {
+                Array4<Real const> const& a = this->const_array(mfi);
+                Array4<int const> const& msk = owner_mask->const_array(mfi);
+                auto tmp = 0.0_rt;
+                AMREX_LOOP_3D(bx, i, j, k,
+                {
+                    if (msk(i,j,k)) { tmp += a(i,j,k,comp); }
+                });
+                sm += tmp; // Do it this way so that it does not break regression tests.
+            }
+        }
+    }
+
+    if (!local) {
+        ParallelAllReduce::Sum(sm, ParallelContext::CommunicatorSub());
+    }
+
+    return sm;
+}
+
 void
 MultiFab::minus (const MultiFab& mf, int strt_comp, int num_comp, int nghost)
 {