Salac/rbf interp cell (#522)

* Allowing for RBF interpolation if the cell is already known. The DMPlexGetContainingCell is REALLY slow (even with the hash). Need to see if there is an alternative that's more efficient. * Version bump * Code cleanup * Format update * Updated tests to account for the new factorization * Forgot to check formatting * Another unit test tolerance update * This runs fine on the laptop. * Now we're on to ARM64 * Another arm64 update * More updates * Fixed bug in rbf->evalder * More ARM64 test updates * Format update
UBCHREST · Mar 13, 2024 · 452903b · 452903b
1 parent e94dfe2
commit 452903b
Show file tree

Hide file tree

Showing 4 changed files with 116 additions and 75 deletions.
diff --git a/CMakeLists.txt b/CMakeLists.txt
@@ -4,7 +4,7 @@ cmake_minimum_required(VERSION 3.18.4)
 include(config/petscCompilers.cmake)
 
 # Set the project details
-project(ablateLibrary VERSION 0.12.27)
+project(ablateLibrary VERSION 0.12.28)
 
 # Load the Required 3rd Party Libaries
 pkg_check_modules(PETSc REQUIRED IMPORTED_TARGET GLOBAL PETSc)

diff --git a/src/domain/RBF/rbf.cpp b/src/domain/RBF/rbf.cpp
@@ -1,5 +1,6 @@
 #include "rbf.hpp"
 #include <petsc/private/dmpleximpl.h>
+#include "utilities/petscSupport.hpp"
 
 using namespace ablate::domain::rbf;
 
@@ -139,7 +140,7 @@ void RBF::Matrix(const PetscInt c) {
     const DM dm = RBF::subDomain->GetSubDM();
 
     // Get the list of neighbor cells
-    DMPlexGetNeighbors(dm, c, -1, -1.0, RBF::minNumberCells, RBF::useCells, RBF::returnNeighborVertices, &nCells, &list) >> utilities::PetscUtilities::checkError;
+    DMPlexGetNeighbors(dm, c, -1, -1.0, RBF::minNumberCells, RBF::useCells, (PetscBool)(RBF::returnNeighborVertices), &nCells, &list) >> utilities::PetscUtilities::checkError;
 
     if (numPoly >= nCells) {
         throw std::invalid_argument("Number of surrounding cells, " + std::to_string(nCells) + ", can not support a requested polynomial order of " + std::to_string(p) + " which requires " +
@@ -228,19 +229,34 @@ void RBF::Matrix(const PetscInt c) {
     MatViewFromOptions(A, NULL, "-ablate::domain::rbf::RBF::A_view") >> utilities::PetscUtilities::checkError;
 
     // Factor the matrix
-    MatLUFactor(A, NULL, NULL, NULL) >> utilities::PetscUtilities::checkError;
+    Mat F;
+
+    MatFactorInfo info;
+    MatFactorInfoInitialize(&info) >> utilities::PetscUtilities::checkError;
+
+    MatGetFactor(A, MATSOLVERPETSC, MAT_FACTOR_QR, &F) >> utilities::PetscUtilities::checkError;
+    MatQRFactorSymbolic(F, A, NULL, &info) >> utilities::PetscUtilities::checkError;
+    MatQRFactorNumeric(F, A, &info) >> utilities::PetscUtilities::checkError;
+
+    //    MatGetFactor(A, MATSOLVERPETSC, MAT_FACTOR_LU, &F) >> utilities::PetscUtilities::checkError;
+    //    info.shifttype = (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE;
+    //    info.shiftamount = 0.1;
+    //    MatLUFactorSymbolic(F, A, NULL, NULL, &info) >> utilities::PetscUtilities::checkError;
+    //    MatLUFactorNumeric(F, A, &info) >> utilities::PetscUtilities::checkError;
+
+    MatDestroy(&A) >> utilities::PetscUtilities::checkError;
 
     PetscFree(xp) >> utilities::PetscUtilities::checkError;
 
     // Assign output
-    RBF::RBFMatrix[c] = A;
+    RBF::RBFMatrix[c] = F;
     RBF::stencilXLocs[c] = x;
     RBF::nStencil[c] = nCells;
     PetscMalloc1(nCells, &(RBF::stencilList[c])) >> utilities::PetscUtilities::checkError;
     PetscArraycpy(RBF::stencilList[c], list, nCells) >> utilities::PetscUtilities::checkError;
 
     // Return the work arrays
-    DMPlexRestoreNeighbors(dm, c, -1, -1.0, RBF::minNumberCells, RBF::useCells, RBF::returnNeighborVertices, &nCells, &list) >> utilities::PetscUtilities::checkError;
+    DMPlexRestoreNeighbors(dm, c, -1, -1.0, RBF::minNumberCells, RBF::useCells, (PetscBool)(RBF::returnNeighborVertices), &nCells, &list) >> utilities::PetscUtilities::checkError;
 }
 
 /************ Begin Derivative Code **********************/
@@ -264,9 +280,9 @@ void RBF::SetDerivatives(PetscInt numDer, PetscInt dx[], PetscInt dy[], PetscInt
 }
 
 /**
- * Set derivatives, defaulting to using vertices
+ * Set derivatives, defaulting to using common vertices
  */
-void RBF::SetDerivatives(PetscInt numDer, PetscInt dx[], PetscInt dy[], PetscInt dz[]) { RBF::SetDerivatives(numDer, dx, dy, dz, PETSC_FALSE); }
+void RBF::SetDerivatives(PetscInt numDer, PetscInt dx[], PetscInt dy[], PetscInt dz[]) { RBF::SetDerivatives(numDer, dx, dy, dz, PETSC_TRUE); }
 
 // Compute the RBF weights at the cell center of p using a cell-list
 // c - The center cell in cellRange ordering
@@ -389,17 +405,15 @@ void RBF::SetupDerivativeStencils() {
 PetscReal RBF::EvalDer(const ablate::domain::Field *field, PetscInt c, PetscInt dx, PetscInt dy, PetscInt dz) {
     RBF::CheckField(field);
 
-    return RBF::EvalDer(field, RBF::subDomain->GetVec(*field), c, dx, dy, dz);
+    return RBF::EvalDer(RBF::subDomain->GetFieldDM(*field), RBF::subDomain->GetVec(*field), field->id, c, dx, dy, dz);
 }
 
-PetscReal RBF::EvalDer(const ablate::domain::Field *field, Vec vec, PetscInt c, PetscInt dx, PetscInt dy, PetscInt dz) {
+PetscReal RBF::EvalDer(DM dm, Vec vec, const PetscInt fid, PetscInt c, PetscInt dx, PetscInt dy, PetscInt dz) {
     PetscReal *wt = nullptr;
     PetscScalar val = 0.0, *f;
     const PetscScalar *array;
     PetscInt nCells = -1, *lst = nullptr;
     PetscInt derID = -1, numDer = RBF::nDer;
-    DM dm = RBF::subDomain->GetFieldDM(*field);
-    const PetscInt fid = field->id;
 
     PetscBool hasKey;
     PetscInt derKey = RBF::derivativeKey(dx, dy, dz);
@@ -408,7 +422,7 @@ PetscReal RBF::EvalDer(const ablate::domain::Field *field, Vec vec, PetscInt c,
     PetscHMapIGet(RBF::hash, derKey, &derID);
 
     // If the stencil hasn't been setup yet do so
-    if (RBF::nStencil[c] < 1) {
+    if (RBF::stencilWeights[c] == nullptr) {
         RBF::SetupDerivativeStencils(c);
     }
 
@@ -435,15 +449,21 @@ PetscReal RBF::EvalDer(const ablate::domain::Field *field, Vec vec, PetscInt c,
 /************ End Derivative Code **********************/
 
 /************ Begin Interpolation Code **********************/
+PetscReal RBF::Interpolate(const ablate::domain::Field *field, Vec f, PetscReal xEval[3]) {
+    DM dm = RBF::subDomain->GetFieldDM(*field);
 
-PetscReal RBF::Interpolate(const ablate::domain::Field *field, PetscReal xEval[3]) {
-    RBF::CheckField(field);
+    PetscInt c;
+    DMPlexGetContainingCell(dm, xEval, &c) >> utilities::PetscUtilities::checkError;
+    if (c < 0) {
+        throw std::runtime_error("ablate::domain::RBF::Interpolate could not determine the location of (" + std::to_string(xEval[0]) + ", " + std::to_string(xEval[1]) + ", " +
+                                 std::to_string(xEval[2]) + ").");
+    }
 
-    return RBF::Interpolate(field, RBF::subDomain->GetVec(*field), xEval);
+    return RBF::Interpolate(field, f, c, xEval);
 }
 
-PetscReal RBF::Interpolate(const ablate::domain::Field *field, Vec f, PetscReal xEval[3]) {
-    PetscInt i, c, nCells, *lst;
+PetscReal RBF::Interpolate(const ablate::domain::Field *field, Vec f, const PetscInt c, PetscReal xEval[3]) {
+    PetscInt i, nCells, *lst;
     PetscScalar *vals, *v;
     const PetscScalar *fvals;
     PetscReal *x, x0[3];
@@ -452,12 +472,7 @@ PetscReal RBF::Interpolate(const ablate::domain::Field *field, Vec f, PetscReal
     DM dm = RBF::subDomain->GetFieldDM(*field);
     const PetscInt fid = field->id;
 
-    DMPlexGetContainingCell(dm, xEval, &c) >> utilities::PetscUtilities::checkError;
-
-    if (c < 0) {
-        throw std::runtime_error("ablate::domain::RBF::Interpolate could not determine the location of (" + std::to_string(xEval[0]) + ", " + std::to_string(xEval[1]) + ", " +
-                                 std::to_string(xEval[2]) + ").");
-    }
+    RBF::CheckField(field);
 
     if (RBF::RBFMatrix[c] == nullptr) {
         RBF::Matrix(c);
@@ -478,8 +493,6 @@ PetscReal RBF::Interpolate(const ablate::domain::Field *field, Vec f, PetscReal
     VecGetArray(rhs, &vals) >> utilities::PetscUtilities::checkError;
 
     for (i = 0; i < nCells; ++i) {
-        // DMPlexPointLocalFieldRead isn't behaving like I would expect. If I don't make f a pointer then it just returns zero.
-        //    Additionally, it looks like it allows for the editing of the value.
         if (fid >= 0) {
             DMPlexPointLocalFieldRead(dm, lst[i], fid, fvals, &v) >> utilities::PetscUtilities::checkError;
         } else {
@@ -559,24 +572,36 @@ PetscReal RBF::Interpolate(const ablate::domain::Field *field, Vec f, PetscReal
 /************ End Interpolation Code **********************/
 
 /************ Constructor, Setup, and Initialization Code **********************/
-RBF::RBF(int polyOrder, bool hasDerivatives, bool hasInterpolation) : polyOrder(polyOrder), hasDerivatives(hasDerivatives), hasInterpolation(hasInterpolation) {}
+RBF::RBF(int polyOrder, bool hasDerivatives, bool hasInterpolation, bool returnNeighborVertices)
+    : polyOrder(polyOrder), returnNeighborVertices(returnNeighborVertices), hasDerivatives(hasDerivatives), hasInterpolation(hasInterpolation) {}
 
 RBF::~RBF() {
+    RBF::FreeStencilData();
+
+    if (dxyz) {
+        PetscFree(dxyz);
+    }
+    if (hash) {
+        PetscHMapIDestroy(&hash);
+    }
+}
+
+void RBF::FreeStencilData() {
     if ((RBF::cEnd - RBF::cStart) > 0) {
         for (PetscInt c = RBF::cStart; c < RBF::cEnd; ++c) {
             PetscFree(RBF::stencilList[c]);
             if (RBF::RBFMatrix[c]) MatDestroy(&(RBF::RBFMatrix[c]));
             PetscFree(RBF::stencilWeights[c]);
             PetscFree(RBF::stencilXLocs[c]);
         }
+        RBF::cellList += cStart;
+        RBF::nStencil += cStart;
+        RBF::stencilList += cStart;
+        RBF::RBFMatrix += cStart;
+        RBF::stencilXLocs += cStart;
+        RBF::stencilWeights += cStart;
         PetscFree6(RBF::cellList, RBF::nStencil, RBF::stencilList, RBF::RBFMatrix, RBF::stencilXLocs, RBF::stencilWeights) >> utilities::PetscUtilities::checkError;
     }
-    if (dxyz) {
-        PetscFree(dxyz);
-    }
-    if (hash) {
-        PetscHMapIDestroy(&hash);
-    }
 }
 
 void RBF::CheckField(const ablate::domain::Field *field) {  // Checks whether the field is SOL or AUX
@@ -669,24 +694,21 @@ void RBF::Setup(std::shared_ptr<ablate::domain::SubDomain> subDomainIn) {
             dz[numDer++] = 2;
         }
 
-        SetDerivatives(numDer, dx, dy, dz);
+        SetDerivatives(numDer, dx, dy, dz, PETSC_TRUE);
     }
 }
 
-void RBF::Initialize(ablate::domain::Range cellRange) {
+void RBF::Initialize() {
+    ablate::domain::Range range;
+
+    // Grab the range of cells from the subDomain
+    RBF::subDomain->GetCellRange(nullptr, range);
+
     // If this is called due to a grid change then release the old memory. In this case cEnd - cStart will be greater than zero.
-    if ((RBF::cEnd - RBF::cStart) > 0) {
-        for (PetscInt c = RBF::cStart; c < RBF::cEnd; ++c) {
-            PetscFree(RBF::stencilList[c]);
-            if (RBF::RBFMatrix[c]) MatDestroy(&(RBF::RBFMatrix[c]));
-            PetscFree(RBF::stencilWeights[c]);
-            PetscFree(RBF::stencilXLocs[c]);
-        }
-        PetscFree6(RBF::cellList, RBF::nStencil, RBF::stencilList, RBF::RBFMatrix, RBF::stencilXLocs, RBF::stencilWeights) >> utilities::PetscUtilities::checkError;
-    }
+    RBF::FreeStencilData();
 
-    RBF::cStart = cellRange.start;
-    RBF::cEnd = cellRange.end;
+    RBF::cStart = range.start;
+    RBF::cEnd = range.end;
 
     // Both interpolation and derivatives need the list of points
     PetscInt nCells = RBF::cEnd - RBF::cStart;
@@ -702,7 +724,7 @@ void RBF::Initialize(ablate::domain::Range cellRange) {
     RBF::stencilWeights -= cStart;
 
     for (PetscInt c = cStart; c < cEnd; ++c) {
-        RBF::cellList[c] = cellRange.GetPoint(c);
+        RBF::cellList[c] = range.GetPoint(c);
 
         RBF::nStencil[c] = -1;
         RBF::stencilList[c] = nullptr;
@@ -711,4 +733,6 @@ void RBF::Initialize(ablate::domain::Range cellRange) {
         RBF::stencilXLocs[c] = nullptr;
         RBF::stencilWeights[c] = nullptr;
     }
+
+    subDomain->RestoreRange(range);
 }
diff --git a/src/domain/RBF/rbf.hpp b/src/domain/RBF/rbf.hpp
@@ -4,7 +4,6 @@
 #include <petsc/private/hashmapi.h>
 #include "domain/range.hpp"  // For domain::Range
 #include "domain/subDomain.hpp"
-#include "utilities/petscSupport.hpp"
 
 #define __RBF_DEFAULT_POLYORDER 3
 
@@ -17,10 +16,10 @@ class RBF {
     // Radial Basis Function type and parameters
     const int polyOrder = 4;
 
-    PetscInt nPoly = -1;                             // The number of polynomial components to include
-    PetscInt minNumberCells = -1;                    // Minimum number of cells-vertices needed to compute the RBF
-    PetscBool useCells = PETSC_FALSE;                // Use vertices or edges/faces when computing neighbor cells/vertices
-    PetscBool returnNeighborVertices = PETSC_FALSE;  // If it is true, it returns neighbor vertices, else it returns neighbor cells
+    PetscInt nPoly = -1;                // The number of polynomial components to include
+    PetscInt minNumberCells = -1;       // Minimum number of cells-vertices needed to compute the RBF
+    PetscBool useCells = PETSC_FALSE;   // Use vertices or edges/faces when computing neighbor cells/vertices
+    const bool returnNeighborVertices;  // If it is true formulates the RBF based on vertices surrounding a cell, otherwise will use cells surrounding a cell
 
     // Information from the subDomain cell range
     PetscInt cStart = 0, cEnd = 0;  // The cell range
@@ -50,18 +49,20 @@ class RBF {
 
     void CheckField(const ablate::domain::Field *field);  // Checks whether the field is SOL or AUX
 
+    void FreeStencilData();
+
    protected:
     PetscReal DistanceSquared(PetscInt dim, PetscReal x[], PetscReal y[]);
     PetscReal DistanceSquared(PetscInt dim, PetscReal x[]);
     void Loc3D(PetscInt dim, PetscReal xIn[], PetscReal x[3]);
 
    public:
-    explicit RBF(int polyOrder = 4, bool hasDerivatives = true, bool hasInterpolation = true);
+    explicit RBF(int polyOrder = 4, bool hasDerivatives = true, bool hasInterpolation = true, bool returnNeighborVertices = false);
 
     virtual ~RBF();
 
     /** SubDomain Register and Setup **/
-    void Initialize(ablate::domain::Range cellRange);
+    void Initialize();
     void Setup(std::shared_ptr<ablate::domain::SubDomain> subDomain);
 
     // Derivative stuff
@@ -95,14 +96,25 @@ class RBF {
 
     /**
      * Return the derivative of a field at a given location
-     * @param field - The field to take the derivative of
-     * @param f - The local vector containing the data
+     * @param dm - The mesh
+     * @param vec - Vector containing the data
+     * @param fid - Field id.
      * @param c - The location in ablate::domain::Range
      * @param dx, dy, dz - The derivative
      */
-    PetscReal EvalDer(const ablate::domain::Field *field, Vec f, PetscInt c, PetscInt dx, PetscInt dy, PetscInt dz);  // Evaluate a derivative
+    PetscReal EvalDer(DM dm, Vec vec, const PetscInt fid, PetscInt c, PetscInt dx, PetscInt dy, PetscInt dz);  // Evaluate a derivative
 
     // Interpolation stuff
+
+    /**
+     * Return the interpolation of a field at a given location
+     * @param field - The field to interpolate
+     * @param f - The local vector containing the data
+     * @param c - Cell containing the evaluation point
+     * @param xEval - The location where to perform the interpolation
+     */
+    PetscReal Interpolate(const ablate::domain::Field *field, Vec f, const PetscInt c, PetscReal xEval[3]);
+
     /**
      * Return the interpolation of a field at a given location
      * @param field - The field to interpolate
@@ -116,7 +128,7 @@ class RBF {
      * @param field - The field to interpolate
      * @param xEval - The location where to perform the interpolation
      */
-    PetscReal Interpolate(const ablate::domain::Field *field, PetscReal xEval[3]);
+    [[nodiscard]] inline PetscReal Interpolate(const ablate::domain::Field *field, PetscReal xEval[3]) { return RBF::Interpolate(field, RBF::subDomain->GetVec(*field), xEval); }
 
     // These will be overwritten in the derived classes
     /**
@@ -139,6 +151,11 @@ class RBF {
      * The RBF kernel type
      */
     virtual std::string_view type() const = 0;
+
+    /**
+     *
+     */
+    inline PetscInt GetDimensions() { return RBF::subDomain->GetDimensions(); }
 };
 
 }  // namespace ablate::domain::rbf