Generalize Interpolator classes

This is an attempt to allow the interpolation algorithms (namely Lagrange) to
be used outside of the AMRInterpolator to anyone that needs it and to the
AHFinder. The algorithm we have is prepared for any dimension, but we only use
it for CH_SPACEDIM. If we want to interpolate something else, this should be a
ready to use class. That is the purpose of this commit. The algorithm is
actually "ready" to even do non-uniform grids ('dx' not constant), but changing
that required the structure to change too much.
For these alternative interpolations, the classes 'SimpleArrayBox' and
'SimpleInterpSource' were added, mimicking 'FArrayBox' and 'GRAMRLevel'. They
represent a D-dimensional grid that can be provided based on a flattened 1D
array, potentially with periodic BCs and constant 'dx' (per direction).

Use cases:
1) interpolating the AHFinder AH centers from the 'stats' output when frequency
of solving the AHFinder changes at a restart (1D interpolation)
2) interpolating the AH surface when 'num_points_u' or 'num_points_v' change at
a restart (2D interpolation)
3) interpolating spherically symmetric initial data into the 3D grid (1D
interpolation)

Source/AMRInterpolator/AMRInterpolator.hpp

@@ -12,15 +12,12 @@
 
 // Chombo includes
 #include "AMRLevel.H"
+#include "LoHiSide.H"
 
 // Our includes
 #include "BoundaryConditions.hpp"
 #include "GRAMR.hpp"
-#include "InterpSource.hpp"
-#include "InterpolationAlgorithm.hpp"
-#include "InterpolationLayout.hpp"
 #include "InterpolationQuery.hpp"
-
 #include "MPIContext.hpp"
 #include "UserVariables.hpp"
 
@@ -31,6 +28,19 @@
 
 template <typename InterpAlgo> class AMRInterpolator
 {
+    struct InterpolationLayout
+    {
+        std::vector<int> rank;
+        std::vector<int> level_idx;
+        std::vector<int> box_idx;
+
+        InterpolationLayout(int num_points)
+            : rank(num_points, -1), level_idx(num_points, -1),
+              box_idx(num_points, -1)
+        {
+        }
+    };
+
   public:
     // constructor for backward compatibility
     // (adds an artificial BC with only periodic BC)
@@ -56,15 +66,17 @@ template <typename InterpAlgo> class AMRInterpolator
     void limit_num_levels(unsigned int num_levels);
     void interp(InterpolationQuery &query);
     const AMR &getAMR() const;
-    const std::array<double, CH_SPACEDIM> &get_coarsest_dx();
-    const std::array<double, CH_SPACEDIM> &get_coarsest_origin();
+    const std::array<double, CH_SPACEDIM> &get_coarsest_dx() const;
+    const std::array<double, CH_SPACEDIM> &get_coarsest_origin() const;
+    bool get_boundary_reflective(Side::LoHiSide a_side, int a_dir) const;
+    bool get_boundary_periodic(int a_dir) const;
 
   private:
     void computeLevelLayouts();
     InterpolationLayout findBoxes(InterpolationQuery &query);
 
     void prepareMPI(InterpolationQuery &query,
-                    const InterpolationLayout layout);
+                    const InterpolationLayout &layout);
     void exchangeMPIQuery();
     void calculateAnswers(InterpolationQuery &query);
     void exchangeMPIAnswer();

Source/AMRInterpolator/AMRInterpolator.impl.hpp

@@ -3,6 +3,10 @@
  * Please refer to LICENSE in GRChombo's root directory.
  */
 
+#if !defined(AMRINTERPOLATOR_HPP_)
+#error "This file should only be included through AMRInterpolator.hpp"
+#endif
+
 #ifndef AMRINTERPOLATOR_IMPL_HPP_
 #define AMRINTERPOLATOR_IMPL_HPP_
 
@@ -41,9 +45,7 @@ void AMRInterpolator<InterpAlgo>::refresh(const bool a_fill_ghosts)
 {
     CH_TIME("AMRInterpolator::refresh");
 
-    const Vector<AMRLevel *> &levels =
-        const_cast<GRAMR &>(m_gr_amr).getAMRLevels();
-    m_num_levels = levels.size();
+    m_num_levels = const_cast<GRAMR &>(m_gr_amr).getAMRLevels().size();
 
     m_mem_level.clear();
     m_mem_box.clear();
@@ -73,17 +75,32 @@ const AMR &AMRInterpolator<InterpAlgo>::getAMR() const
 
 template <typename InterpAlgo>
 const std::array<double, CH_SPACEDIM> &
-AMRInterpolator<InterpAlgo>::get_coarsest_dx()
+AMRInterpolator<InterpAlgo>::get_coarsest_dx() const
 {
     return m_coarsest_dx;
 }
 template <typename InterpAlgo>
 const std::array<double, CH_SPACEDIM> &
-AMRInterpolator<InterpAlgo>::get_coarsest_origin()
+AMRInterpolator<InterpAlgo>::get_coarsest_origin() const
 {
     return m_coarsest_origin;
 }
 
+template <typename InterpAlgo>
+bool AMRInterpolator<InterpAlgo>::get_boundary_reflective(Side::LoHiSide a_side,
+                                                          int a_dir) const
+{
+    if (a_side == Side::Lo)
+        return m_lo_boundary_reflective[a_dir];
+    else
+        return m_hi_boundary_reflective[a_dir];
+}
+template <typename InterpAlgo>
+bool AMRInterpolator<InterpAlgo>::get_boundary_periodic(int a_dir) const
+{
+    return m_bc_params.is_periodic[a_dir];
+}
+
 template <typename InterpAlgo>
 void AMRInterpolator<InterpAlgo>::limit_num_levels(unsigned int num_levels)
 {
@@ -239,10 +256,14 @@ void AMRInterpolator<InterpAlgo>::computeLevelLayouts()
         if (m_verbosity >= 2)
         {
             _pout << "    Level " << level_idx << "\t"
-                  << "dx=(" << m_dx[level_idx][0] << "," << m_dx[level_idx][1]
-                  << "," << m_dx[level_idx][2] << ")\t"
-                  << "grid_origin=(" << m_origin[level_idx][0] << ","
-                  << m_origin[level_idx][1] << "," << m_origin[level_idx][2]
+                  << "dx=("
+                  << D_TERM(m_dx[level_idx][0], << "," << m_dx[level_idx][1],
+                                                << "," << m_dx[level_idx][2])
+                  << ")\t"
+                  << "grid_origin=("
+                  << D_TERM(m_origin[level_idx][0],
+                            << "," << m_origin[level_idx][1],
+                            << "," << m_origin[level_idx][2])
                   << ")" << endl;
         }
 
@@ -257,7 +278,7 @@ void AMRInterpolator<InterpAlgo>::computeLevelLayouts()
 }
 
 template <typename InterpAlgo>
-InterpolationLayout
+typename AMRInterpolator<InterpAlgo>::InterpolationLayout
 AMRInterpolator<InterpAlgo>::findBoxes(InterpolationQuery &query)
 {
     CH_TIME("AMRInterpolator::findBoxes");
@@ -293,7 +314,7 @@ AMRInterpolator<InterpAlgo>::findBoxes(InterpolationQuery &query)
     //    const AMRLevel& level = *levels[level_idx];
 
     //    const LevelData<FArrayBox>& level_data = dynamic_cast<const
-    //    InterpSource&>(level).getLevelData(); const DisjointBoxLayout&
+    //    InterpSource<>&>(level).getLevelData(); const DisjointBoxLayout&
     //    box_layout = level_data.disjointBoxLayout(); const Box& domain_box =
     //    level.problemDomain().domainBox();
 
@@ -381,7 +402,7 @@ AMRInterpolator<InterpAlgo>::findBoxes(InterpolationQuery &query)
         const AMRLevel &level = *levels[level_idx];
 
         const LevelData<FArrayBox> &level_data =
-            dynamic_cast<const InterpSource &>(level).getLevelData(
+            dynamic_cast<const InterpSource<> &>(level).getLevelData(
                 VariableType::evolution);
         const DisjointBoxLayout &box_layout = level_data.disjointBoxLayout();
         const Box &domain_box = level.problemDomain().domainBox();
@@ -495,7 +516,7 @@ AMRInterpolator<InterpAlgo>::findBoxes(InterpolationQuery &query)
 
 template <typename InterpAlgo>
 void AMRInterpolator<InterpAlgo>::prepareMPI(InterpolationQuery &query,
-                                             const InterpolationLayout layout)
+                                             const InterpolationLayout &layout)
 {
     CH_TIME("AMRInterpolator::prepareMPI");
 
@@ -631,15 +652,15 @@ void AMRInterpolator<InterpAlgo>::calculateAnswers(InterpolationQuery &query)
     const int num_answers = m_mpi.totalAnswerCount();
 
     std::array<double, CH_SPACEDIM> grid_coord;
-    IntVect nearest;
 
     for (int answer_idx = 0; answer_idx < num_answers; ++answer_idx)
     {
         const int box_idx = m_answer_box[answer_idx];
         const int level_idx = m_answer_level[answer_idx];
 
         const AMRLevel &level = *levels[level_idx];
-        const InterpSource &source = dynamic_cast<const InterpSource &>(level);
+        const InterpSource<> &source =
+            dynamic_cast<const InterpSource<> &>(level);
         const LevelData<FArrayBox> *const evolution_level_data_ptr =
             &source.getLevelData(VariableType::evolution);
         const LevelData<FArrayBox> *diagnostics_level_data_ptr;
@@ -657,10 +678,6 @@ void AMRInterpolator<InterpAlgo>::calculateAnswers(InterpolationQuery &query)
                 &diagnostics_level_data_ptr->disjointBoxLayout();
         }
 
-        const Box &domain_box = level.problemDomain().domainBox();
-        const IntVect &small_end = domain_box.smallEnd();
-        const IntVect &big_end = domain_box.bigEnd();
-
         // Convert the LayoutIndex to DataIndex
         const DataIndex evolution_data_idx(
             evolution_box_layout_ptr->layoutIterator()[box_idx]);
@@ -696,22 +713,6 @@ void AMRInterpolator<InterpAlgo>::calculateAnswers(InterpolationQuery &query)
                   << (box.bigEnd()[i] + 0.5) << "]";
                 MayDay::Abort(s.str().c_str());
             }
-
-            // point lies beyond the "small end" of the whole domain, but still
-            // within the boundary cell
-            if (grid_coord[i] < small_end[i] &&
-                grid_coord[i] >= small_end[i] - 0.5)
-                nearest[i] = small_end[i];
-
-            // point lies beyond the "big end" of the whole domain, but still
-            // within the boundary cell
-            else if (grid_coord[i] > big_end[i] &&
-                     grid_coord[i] <= big_end[i] + 0.5)
-                nearest[i] = big_end[i];
-
-            // otherwise we round to nearest grid point
-            else
-                nearest[i] = (int)ceil(grid_coord[i] - 0.5);
         }
 
         if (m_verbosity >= 2)
@@ -737,7 +738,7 @@ void AMRInterpolator<InterpAlgo>::calculateAnswers(InterpolationQuery &query)
 
             typedef std::vector<typename InterpolationQuery::out_t> comps_t;
             comps_t &comps = deriv_it->second;
-            algo.setup(deriv, m_dx[level_idx], grid_coord, nearest);
+            algo.setup(deriv, grid_coord);
 
             for (typename comps_t::iterator it = comps.begin();
                  it != comps.end(); ++it)

Source/AMRInterpolator/InterpSource.hpp

@@ -16,14 +16,21 @@
 // Chombo namespace
 #include "UsingNamespace.H"
 
-// Abstrace base class to get the FABs out of an AMRLevel
-class InterpSource
+// Abstrace base class to identify points in a grid structure
+template <int N_DIMS = CH_SPACEDIM> class InterpSource
 {
   public:
     virtual const LevelData<FArrayBox> &getLevelData(
         const VariableType var_type = VariableType::evolution) const = 0;
-    virtual bool
-    contains(const std::array<double, CH_SPACEDIM> &point) const = 0;
+    virtual bool contains(const std::array<double, N_DIMS> &point) const = 0;
+    virtual double get_dx(int dir) const = 0;
+    virtual std::array<double, N_DIMS> get_dxs() const
+    {
+        std::array<double, N_DIMS> dxs;
+        for (int i = 0; i < N_DIMS; ++i)
+            dxs[i] = get_dx(i);
+        return dxs;
+    }
 };
 
 #endif /* INTERPSOURCE_H_ */

Source/AMRInterpolator/Lagrange.hpp

@@ -9,9 +9,9 @@
 #include "InterpSource.hpp"
 #include <utility>
 
-template <int Order> class Lagrange
+template <int Order, int N_DIMS = CH_SPACEDIM> class Lagrange
 {
-    const InterpSource &m_source;
+    const InterpSource<N_DIMS> &m_source;
     bool m_verbosity;
 
     struct Stencil;
@@ -41,10 +41,10 @@ template <int Order> class Lagrange
     // Helper function to generate tensor product weights
     // Argument 'dim' is used for recursion over dimensions.
     pair<std::vector<IntVect>, std::vector<double>>
-    generateStencil(const std::array<int, CH_SPACEDIM> &deriv,
-                    const std::array<double, CH_SPACEDIM> &dx,
-                    const std::array<double, CH_SPACEDIM> &evalCoord,
-                    const IntVect &nearest, int dim = CH_SPACEDIM - 1);
+    generateStencil(const std::array<int, N_DIMS> &deriv,
+                    const std::array<double, N_DIMS> &dx,
+                    const std::array<double, N_DIMS> &eval_index,
+                    int dim = N_DIMS - 1);
 
     std::vector<IntVect> m_interp_points;
     std::vector<double> m_interp_weights;
@@ -56,13 +56,16 @@ template <int Order> class Lagrange
     multiset<double> m_interp_pos;
 
   public:
-    Lagrange(const InterpSource &source, bool verbosity = false);
+    Lagrange(const InterpSource<N_DIMS> &source, bool verbosity = false);
 
-    void setup(const std::array<int, CH_SPACEDIM> &deriv,
-               const std::array<double, CH_SPACEDIM> &dx,
-               const std::array<double, CH_SPACEDIM> &evalCoord,
-               const IntVect &nearest);
-    double interpData(const FArrayBox &fab, int comp);
+    // eval_index is in 'index' coordinates, not physical coordinates
+    void setup(const std::array<int, N_DIMS> &deriv,
+               const std::array<double, N_DIMS> &eval_index);
+
+    // any class with a method:
+    // Real get(const IntVect &a_iv, int a_comp) const
+    template <class GeneralArrayBox>
+    double interpData(const GeneralArrayBox &fab, int comp = 0);
 
     const static string TAG;
 };