Generalize Interpolator classes #176

TaigoFr · 2021-05-27T16:55:13Z

This PR 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 already to use class. That is the purpose of this PR.

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). One can simply take a look at the InterpolatorTest to see how to use them.

Use case examples:

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

KAClough · 2021-06-02T11:35:57Z

Again @mirenradia is probably the best person to review this, but I'm ok with the principle of it.

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)

mirenradia

I know why you wanted to do this but at the moment it feels a little hacky. However, I think it is possible to salvage (see my comments below).

mirenradia · 2021-08-01T17:36:59Z

Source/AMRInterpolator/AMRInterpolator.hpp

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


I think I'm ok with moving this here but we should delete the InterpolationLayout.hpp file too.

mirenradia · 2021-08-01T17:43:41Z

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


Might as well fix the typo while we're at it

Suggested change

// Abstrace base class to identify points in a grid structure

// Abstract base class to identify points in a grid structure

mirenradia · 2021-08-01T17:45:15Z

Source/AMRInterpolator/InterpSource.hpp

 {
  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;


Can you add virtual and override specifiers to the corresponding override in GRAMRLevel?

mirenradia · 2021-08-01T17:51:21Z

Source/AMRInterpolator/InterpSource.hpp

-    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;


It looks like this function has to be defined even if the specialised class represents an InterpSource with grid spacings that are different in each direction which is not ideal. Perhaps only the one below should be pure virtual and this can be removed?

mirenradia · 2021-08-01T18:00:04Z

Source/AMRInterpolator/Lagrange.impl.hpp

+    // TF: assume nearest point is at 'std::round(eval_index[dim])'
+    // this used to be an argument, but I think this assumption should always
+    // hold
+    int candidate = std::round(eval_index[dim]);


Won't this be problematic when the eval_index[dim] exactly coincides with the edge of the high boundary of a box? Unless I'm wrong, I'd prefer to keep the nearest argument. It also makes the implementation of the nearest neighbour interpolation algorithm more straightforward.

mirenradia · 2021-08-01T22:02:55Z

Source/AMRInterpolator/SimpleArrayBox.hpp

+    // IntVect will be CH_SPACEDIM dimension, but ignore the lasts components if
+    // N_DIMS < CH_SPACEDIM
+    // the 'comp' argument doesn't matter, always assume 'm_f' is what matters
+    Real get(const IntVect &a_iv, int a_comp) const


What if I want to do N_DIMS > CH_SPACEDIM? I'm not convinced by the use of IntVects for this interface.

mirenradia · 2021-08-01T22:03:55Z

Source/AMRInterpolator/SimpleArrayBox.hpp

+    SimpleArrayBox(std::array<int, N_DIMS> a_points_per_dir,
+                   std::vector<double> a_f,
+                   std::array<bool, N_DIMS> a_is_periodic = {false})
+    {
+        set(a_points_per_dir, a_f, a_is_periodic);
+    }
+
+    void set(std::array<int, N_DIMS> a_points_per_dir, std::vector<double> a_f,
+             std::array<bool, N_DIMS> a_is_periodic = {false})


Maybe pass at least the std::vector<double>s by reference?

mirenradia · 2021-08-01T22:35:14Z

Source/AMRInterpolator/SimpleInterpSource.hpp

+    const LevelData<FArrayBox> &
+    getLevelData(const VariableType var_type = VariableType::evolution) const
+    {
+        CH_assert("This should not be needed.");
+        return *m_fake;
+    }


I am not a fan of this. It's quite hacky. Also I'm not sure whether a string evaluates to true or false! Since the InterpAlgo classes only need the get_dxs() and contains() members of InterpSource, maybe it should be refactored into a base (needed by the InterpAlgo classes) and child (needed by the AMRInterpolator class)? For example

class LevelDataInterpSource : public InterpSource

Alternatively, I'm not even sure of the point of AMRInterpolator having a generic InterpSource anyway so we could just make it explicitly use GRAMRLevel and remove the getLevelData() function from InterpSource?

mirenradia · 2021-08-01T22:44:44Z

Tests/InterpolatorTest/InterpolatorTest.cpp

+    mainSetup(argc, argv);
+
+    int status = runInterpolatorTest(argc, argv);
+
+    if (status == 0)
+        std::cout << "Interpolator test passed." << endl;
+    else
+        std::cout << "Interpolator test failed with return code " << status
+                  << endl;
+
+    mainFinalize();


Are mainSetup() and mainFinalize() necessary here? If not, then we can also get rid of the InterpolatorTest.inputs, SimulationParameters.hpp and UserVariables.hpp files too.

TaigoFr requested review from KAClough and mirenradia May 27, 2021 16:55

mirenradia added the enhancement Modification of existing feature/general improvement label May 27, 2021

mirenradia assigned TaigoFr May 27, 2021

TaigoFr force-pushed the enhancement/interpolator_usability branch from c8294a8 to ff9401a Compare June 4, 2021 16:33

TaigoFr force-pushed the enhancement/interpolator_usability branch from ff9401a to 85c9231 Compare July 21, 2021 10:57

mirenradia requested changes Aug 1, 2021

View reviewed changes

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Generalize Interpolator classes #176

Generalize Interpolator classes #176

TaigoFr commented May 27, 2021 •

edited

Loading

KAClough commented Jun 2, 2021

mirenradia left a comment

mirenradia Aug 1, 2021

mirenradia Aug 1, 2021

mirenradia Aug 1, 2021

mirenradia Aug 1, 2021

mirenradia Aug 1, 2021 •

edited

Loading

mirenradia Aug 1, 2021

mirenradia Aug 1, 2021

mirenradia Aug 1, 2021

mirenradia Aug 1, 2021

	// Abstrace base class to identify points in a grid structure
	// Abstract base class to identify points in a grid structure

Generalize Interpolator classes #176

Are you sure you want to change the base?

Generalize Interpolator classes #176

Conversation

TaigoFr commented May 27, 2021 • edited Loading

KAClough commented Jun 2, 2021

mirenradia left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

mirenradia Aug 1, 2021 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

TaigoFr commented May 27, 2021 •

edited

Loading

mirenradia Aug 1, 2021 •

edited

Loading