FillPatcher class

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,10 @@ public:
     virtual void particle_redistribute (int /*lbase*/ = 0, bool /*a_init*/ = false) {;}
 #endif
 
+    // Fill with FillPatcher on level > 0 and AmrLevel::FillPatch on level 0.
+    void FillPatcherFill (amrex::MultiFab& mf, int nghost, amrex::Real time,
+                          int state_index);
+
     static void FillPatch (AmrLevel& amrlevel,
                            MultiFab& leveldata,
                            int       boxGrow,
@@ -425,7 +432,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 +443,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,59 @@ void AmrLevel::constructAreaNotToTag ()
     }
 }
 
+void
+AmrLevel::resetFillPatcher ()
+{
+    for (auto& fp : m_fillpatcher) {
+        fp.reset();
+    }
+}
+
+void
+AmrLevel::FillPatcherFill (MultiFab& mf, int nghost, Real time, int state_index)
+{
+    if (level == 0) {
+        FillPatch(*this, mf, nghost, time, state_index, 0, mf.nComp());
+    } 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,0,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,0,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(0))) {
+            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.interp(0));
+        }
+
+        fillpatcher->fill(mf, time, smf_crse, stime_crse, smf_fine, stime_fine,
+                          physbcf_crse, physbcf_fine, desc.getBCs());
+    }
+}
+
 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>