FillPatcher class (AMReX-Codes#2972)

This adds a class FillPatcher for filling fine level data. It's not as
general as the various FillPatch functions (e.g., FillPatchTwoLevels).
However, it can reduce the amount of communication data. Suppose we use
RK2 with subcycling and the refinement ratio is 2. For each step on
level 0, there are two steps on level 1. With RK2, each fine step needs
to call FillPatch twice. So the total number of FillPatch calls is 4 in
the two fine steps. Using the free function, one ParallelCopy per
FillPatch call is needed for copying coarse data for spatial
interpolation. With the FillPatcher class, two ParallelCopy calls will
be done to copy old and new coarse data. Then these data will be used in
the four FillPatcher::fill calls. This new approach saves two
ParallelCopy calls per coarse step for a two levels run. It could save
more if the time stepping requires more substeps or the refinement ratio
is higher. Note that many of our AMReX codes use a time stepping
algorithm that needs only one FillPatch call per step. For those codes,
this new approach will not save any communication for a refinement ratio
of 2. However, it will save communication when the refinement ratio is
4.

Src/Amr/AMReX_AmrLevel.H

@@ -15,6 +15,7 @@
 #include <AMReX_StateDescriptor.H>
 #include <AMReX_StateData.H>
 #include <AMReX_VisMF.H>
+#include <AMReX_FillPatcher.H>
 #ifdef AMREX_USE_EB
 #include <AMReX_EBSupport.H>
 #endif
@@ -243,12 +244,14 @@ public:
     Long countCells () const noexcept;
 
     //! Get the area not to tag.
-    const BoxArray& getAreaNotToTag() noexcept;
-    const Box& getAreaToTag() noexcept;
+    const BoxArray& getAreaNotToTag () noexcept;
+    const Box& getAreaToTag () noexcept;
     //! Construct the area not to tag.
-    void constructAreaNotToTag();
+    void constructAreaNotToTag ();
     //! Set the area not to tag.
-    void setAreaNotToTag(BoxArray& ba) noexcept;
+    void setAreaNotToTag (BoxArray& ba) noexcept;
+
+    void resetFillPatcher ();
 
     /**
     * \brief Error estimation for regridding. This is a pure virtual
@@ -365,6 +368,20 @@ public:
     virtual void particle_redistribute (int /*lbase*/ = 0, bool /*a_init*/ = false) {;}
 #endif
 
+    /**
+     * \brief Fill with FillPatcher on level > 0 and AmrLevel::FillPatch on level 0.
+     *
+     * \param mf     destination MultiFab
+     * \param dcomp  starting component for the destination
+     * \param ncomp  number of component to fill
+     * \param nghost number of ghost cells to fill
+     * \param time   time
+     * \param state_index StateData index
+     * \param scomp  starting component in the StateData
+     */
+    void FillPatcherFill (amrex::MultiFab& mf, int dcomp, int ncomp, int nghost,
+                          amrex::Real time, int state_index, int scomp);
+
     static void FillPatch (AmrLevel& amrlevel,
                            MultiFab& leveldata,
                            int       boxGrow,
@@ -425,7 +442,7 @@ protected:
     IntVect               fine_ratio;   // Refinement ratio to finer level.
     static DeriveList     derive_lst;   // List of derived quantities.
     static DescriptorList desc_lst;     // List of state variables.
-    Vector<StateData>      state;        // Array of state data.
+    Vector<StateData>     state;        // Array of state data.
 
     BoxArray              m_AreaNotToTag; //Area which shouldn't be tagged on this level.
     Box                   m_AreaToTag;    //Area which is allowed to be tagged on this level.
@@ -436,6 +453,8 @@ protected:
 
     std::unique_ptr<FabFactory<FArrayBox> > m_factory;
 
+    Vector<std::unique_ptr<FillPatcher<MultiFab>>> m_fillpatcher;
+
 private:
 
     mutable BoxArray      edge_grids[AMREX_SPACEDIM];  // face-centered grids

Src/Amr/AMReX_AmrLevel.cpp

@@ -102,6 +102,7 @@ AmrLevel::AmrLevel (Amr&            papa,
     }
 
     state.resize(desc_lst.size());
+    m_fillpatcher.resize(desc_lst.size());
 
 #ifdef AMREX_USE_EB
     if (EB2::TopIndexSpaceIfPresent()) {
@@ -451,6 +452,8 @@ AmrLevel::restart (Amr&          papa,
         }
     }
 
+    m_fillpatcher.resize(ndesc);
+
     if (parent->useFixedCoarseGrids()) constructAreaNotToTag();
 
     post_step_regrid = 0;
@@ -2096,6 +2099,63 @@ void AmrLevel::constructAreaNotToTag ()
     }
 }
 
+void
+AmrLevel::resetFillPatcher ()
+{
+    for (auto& fp : m_fillpatcher) {
+        fp.reset();
+    }
+}
+
+void
+AmrLevel::FillPatcherFill (MultiFab& mf, int dcomp, int ncomp, int nghost,
+                           Real time, int state_index, int scomp)
+{
+    if (level == 0) {
+        FillPatch(*this, mf, nghost, time, state_index, scomp, ncomp, dcomp);
+    } else {
+        AmrLevel& fine_level = *this;
+        AmrLevel& crse_level = parent->getLevel(level-1);
+        const Geometry& geom_fine = fine_level.geom;
+        const Geometry& geom_crse = crse_level.geom;
+
+        Vector<MultiFab*> smf_crse;
+        Vector<Real> stime_crse;
+        StateData& statedata_crse = crse_level.state[state_index];
+        statedata_crse.getData(smf_crse,stime_crse,time);
+        StateDataPhysBCFunct physbcf_crse(statedata_crse,scomp,geom_crse);
+
+        Vector<MultiFab*> smf_fine;
+        Vector<Real> stime_fine;
+        StateData& statedata_fine = fine_level.state[state_index];
+        statedata_fine.getData(smf_fine,stime_fine,time);
+        StateDataPhysBCFunct physbcf_fine(statedata_fine,scomp,geom_fine);
+
+        const StateDescriptor& desc = AmrLevel::desc_lst[state_index];
+
+        if (level > 1 &&!amrex::ProperlyNested(fine_level.crse_ratio,
+                                               parent->blockingFactor(fine_level.level),
+                                               nghost, mf.ixType(),
+                                               desc.interp(scomp))) {
+            amrex::Abort("FillPatcherFill: Grids are not properly nested.  Must increase blocking factor.");
+        }
+
+        auto& fillpatcher = m_fillpatcher[state_index];
+        if (fillpatcher == nullptr) {
+            fillpatcher = std::make_unique<FillPatcher<MultiFab>>
+                (parent->boxArray(level), parent->DistributionMap(level), geom_fine,
+                 parent->boxArray(level-1), parent->DistributionMap(level-1), geom_crse,
+                 IntVect(nghost), desc.nComp(), desc.interp(scomp));
+        }
+
+        fillpatcher->fill(mf, IntVect(nghost), time,
+                          smf_crse, stime_crse, smf_fine, stime_fine,
+                          scomp, dcomp, ncomp,
+                          physbcf_crse, scomp, physbcf_fine, scomp,
+                          desc.getBCs(), scomp);
+    }
+}
+
 void
 AmrLevel::FillPatch (AmrLevel& amrlevel,
                      MultiFab& leveldata,

Src/AmrCore/AMReX_FillPatchUtil.H

@@ -28,12 +28,17 @@
 namespace amrex
 {
 
-    template <typename FAB>
+    template <typename MFFAB>
     struct NullInterpHook
     {
-        void operator() (FAB& /*fab*/, const Box& /*bx*/, int /*icomp*/, int /*ncomp*/) const {}
+        template <class F=MFFAB, std::enable_if_t<IsBaseFab<F>::value,int> = 0>
+        void operator() (MFFAB& /*fab*/, const Box& /*bx*/, int /*icomp*/, int /*ncomp*/) const {}
 
-        void operator() (Array<FAB*, AMREX_SPACEDIM> /*fab*/, const Box& /*bx*/, int /*icomp*/, int /*ncomp*/) const {}
+        template <class F=MFFAB, std::enable_if_t<IsBaseFab<F>::value,int> = 0>
+        void operator() (Array<MFFAB*, AMREX_SPACEDIM> /*fab*/, const Box& /*bx*/, int /*icomp*/, int /*ncomp*/) const {}
+
+        template <class F=MFFAB, std::enable_if_t<IsFabArray<F>::value,int> = 0>
+        void operator() (MFFAB& /*mf*/, int /*icomp*/, int /*ncomp*/) const {}
     };
 
     template <typename Interp>