Mcgurn/accessors (#521)

* added protype accessors to the domain

* adding missing header

* version bump

* added missing header to cmake

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.26)
+project(ablateLibrary VERSION 0.12.27)
 
 # Load the Required 3rd Party Libaries
 pkg_check_modules(PETSc REQUIRED IMPORTED_TARGET GLOBAL PETSc)

src/boundarySolver/physics/sublimation.cpp

@@ -1,5 +1,6 @@
 #include "sublimation.hpp"
 #include <utility>
+#include "domain/constAccessor.hpp"
 #include "domain/dynamicRange.hpp"
 #include "finiteVolume/compressibleFlowFields.hpp"
 #include "io/interval/fixedInterval.hpp"
@@ -377,70 +378,50 @@ PetscErrorCode ablate::boundarySolver::physics::Sublimation::SublimationPreRHS(B
 
 void ablate::boundarySolver::physics::Sublimation::UpdateSpecies(TS ts, ablate::solver::Solver &solver) {
     // Get the density and densityYi field info
-    const auto &eulerFieldInfo = solver.GetSubDomain().GetField(finiteVolume::CompressibleFlowFields::EULER_FIELD);
-    const auto &densityYiFieldInfo = solver.GetSubDomain().GetField(finiteVolume::CompressibleFlowFields::DENSITY_YI_FIELD);
-    const auto &yiFieldInfo = solver.GetSubDomain().GetField(finiteVolume::CompressibleFlowFields::YI_FIELD);
-
-    // Get the solution vec and dm
-    auto dm = solver.GetSubDomain().GetDM();
-    auto solVec = solver.GetSubDomain().GetSolutionVector();
-    auto auxDm = solver.GetSubDomain().GetAuxDM();
-    auto auxVec = solver.GetSubDomain().GetAuxVector();
+    auto eulerFieldAccessor = solver.GetSubDomain().GetConstSolutionAccessor(finiteVolume::CompressibleFlowFields::EULER_FIELD);
+    auto densityYiAccessor = solver.GetSubDomain().GetSolutionAccessor(finiteVolume::CompressibleFlowFields::DENSITY_YI_FIELD);
+    auto yiFieldAccessor = solver.GetSubDomain().GetAuxAccessor(finiteVolume::CompressibleFlowFields::YI_FIELD);
 
     // get the time from the ts
     PetscReal time;
     TSGetTime(ts, &time) >> utilities::PetscUtilities::checkError;
 
-    // Get the array vector
-    PetscScalar *solutionArray;
-    VecGetArray(solVec, &solutionArray) >> utilities::PetscUtilities::checkError;
-    PetscScalar *auxArray;
-    VecGetArray(auxVec, &auxArray) >> utilities::PetscUtilities::checkError;
-
     // March over each cell in this domain
     ablate::domain::Range cellRange;
     solver.GetCellRange(cellRange);
     auto dim = solver.GetSubDomain().GetDimensions();
 
     // Get the cell geom vec
     Vec cellGeomVec;
-    const PetscScalar *cellGeomArray;
-    DM cellGeomDm;
-    DMPlexGetDataFVM(dm, nullptr, &cellGeomVec, nullptr, nullptr) >> utilities::PetscUtilities::checkError;
-    VecGetDM(cellGeomVec, &cellGeomDm) >> utilities::PetscUtilities::checkError;
-    VecGetArrayRead(cellGeomVec, &cellGeomArray) >> utilities::PetscUtilities::checkError;
+    DMPlexGetDataFVM(solver.GetSubDomain().GetDM(), nullptr, &cellGeomVec, nullptr, nullptr) >> utilities::PetscUtilities::checkError;
+
+    // create an accessor the cell data
+    domain::ConstAccessor<PetscFVCellGeom> cellGeomAccessor(cellGeomVec);
 
     for (PetscInt c = cellRange.start; c < cellRange.end; ++c) {
         PetscInt cell = cellRange.points ? cellRange.points[c] : c;
 
         // Get the euler and density field
-        const PetscScalar *euler = nullptr;
-        DMPlexPointGlobalFieldRead(dm, cell, eulerFieldInfo.id, solutionArray, &euler) >> utilities::PetscUtilities::checkError;
-        PetscScalar *densityYi;
-        DMPlexPointGlobalFieldRef(dm, cell, densityYiFieldInfo.id, solutionArray, &densityYi) >> utilities::PetscUtilities::checkError;
-        PetscScalar *yi;
-        DMPlexPointLocalFieldRead(auxDm, cell, yiFieldInfo.id, auxArray, &yi) >> utilities::PetscUtilities::checkError;
-        PetscFVCellGeom *cellGeom;
-        DMPlexPointLocalRead(cellGeomDm, cell, cellGeomArray, &cellGeom) >> utilities::PetscUtilities::checkError;
+        const PetscScalar *euler = eulerFieldAccessor[cell];
+        PetscScalar *densityYi = densityYiAccessor[cell];
+        PetscScalar *yi = yiFieldAccessor[cell];
+        const PetscFVCellGeom *cellGeom = cellGeomAccessor[cell];
 
         // compute the mass fractions on the boundary
-        massFractionsFunction(dim, time, cellGeom->centroid, yiFieldInfo.numberComponents, yi, massFractionsContext);
+        massFractionsFunction(dim, time, cellGeom->centroid, yiFieldAccessor.GetField().numberComponents, yi, massFractionsContext);
 
         // Only update if in the global vector
         if (euler) {
             // Get density
             const PetscScalar density = euler[finiteVolume::CompressibleFlowFields::RHO];
 
-            for (PetscInt sp = 0; sp < densityYiFieldInfo.numberComponents; sp++) {
+            for (PetscInt sp = 0; sp < densityYiAccessor.GetField().numberComponents; sp++) {
                 densityYi[sp] = yi[sp] * density;
             }
         }
     }
 
     // cleanup
-    VecRestoreArrayRead(cellGeomVec, &cellGeomArray) >> utilities::PetscUtilities::checkError;
-    VecRestoreArray(auxVec, &auxArray) >> utilities::PetscUtilities::checkError;
-    VecRestoreArray(solVec, &solutionArray) >> utilities::PetscUtilities::checkError;
     solver.RestoreRange(cellRange);
 }
 

src/domain/CMakeLists.txt

@@ -41,6 +41,10 @@ target_sources(ablateLibrary
         hdf5Initializer.hpp
         initializerList.hpp
         meshGenerator.hpp
+        fieldAccessor.hpp
+        constFieldAccessor.hpp
+        accessor.hpp
+        constAccessor.hpp
 )
 
 add_subdirectory(modifiers)

src/domain/accessor.hpp

@@ -0,0 +1,75 @@
+#ifndef ABLATELIBRARY_ACCESSOR_HPP
+#define ABLATELIBRARY_ACCESSOR_HPP
+
+#include <petsc.h>
+#include "field.hpp"
+
+namespace ablate::domain {
+/**
+ * Class responsible for computing point data locations for dm plex integration
+ */
+template <class DataType, bool isLocal = true>
+class Accessor {
+   private:
+    /**
+     * The petsc vector used to hold the data
+     */
+    Vec dataVector;
+
+    /**
+     * The extracted data from the vector
+     */
+    PetscScalar* dataArray;
+
+    /**
+     * Store the DM for the vector to find memory locations
+     */
+    DM dataDM;
+
+   public:
+    Accessor(Vec dataVectorIn, DM dm = nullptr) : dataVector(dataVectorIn), dataDM(dm) {
+        // Get read/write access to the vector
+        VecGetArray(dataVector, &dataArray) >> ablate::utilities::PetscUtilities::checkError;
+
+        // Get the dm from the vector if not specified
+        if (!dataDM) {
+            VecGetDM(dataVector, &dataDM) >> ablate::utilities::PetscUtilities::checkError;
+        }
+    }
+
+    ~Accessor() {
+        // Put back the vector
+        VecRestoreArray(dataVector, &dataArray) >> ablate::utilities::PetscUtilities::checkError;
+    };
+
+    //! Inline function to compute the memory address at this point
+    template <class IndexType>
+    inline DataType* operator[](IndexType point) {
+        DataType* field;
+        if constexpr (isLocal) {
+            DMPlexPointLocalRef(dataDM, point, dataArray, &field) >> ablate::utilities::PetscUtilities::checkError;
+        } else {
+            DMPlexPointGlobalRef(dataDM, point, dataArray, &field) >> ablate::utilities::PetscUtilities::checkError;
+        }
+        return field;
+    }
+
+    //! Inline function to compute the memory address at this point using PETSC style return error
+    template <class IndexType>
+    inline PetscErrorCode operator()(IndexType point, DataType** field) {
+        PetscFunctionBeginHot;
+        if constexpr (isLocal) {
+            PetscCall(DMPlexPointLocalRef(dataDM, point, dataArray, field));
+        } else {
+            PetscCall(DMPlexPointGlobalRef(dataDM, point, dataArray, field));
+        }
+        PetscFunctionReturn(PETSC_SUCCESS);
+    }
+
+    /**
+     * prevent copy of this class
+     */
+    Accessor(const Accessor&) = delete;
+};
+}  // namespace ablate::domain
+#endif  // ABLATELIBRARY_FIELDACCESSOR_HPP

src/domain/constAccessor.hpp

@@ -0,0 +1,75 @@
+#ifndef ABLATELIBRARY_CONSTACCESSOR_HPP
+#define ABLATELIBRARY_CONSTACCESSOR_HPP
+
+#include <petsc.h>
+#include "field.hpp"
+
+namespace ablate::domain {
+/**
+ * Class responsible for computing point data locations for dm plex integration
+ */
+template <class DataType, bool isLocal = true>
+class ConstAccessor {
+   private:
+    /**
+     * The petsc vector used to hold the data
+     */
+    Vec dataVector;
+
+    /**
+     * The extracted data from the vector
+     */
+    const PetscScalar* dataArray;
+
+    /**
+     * Store the DM for the vector to find memory locations
+     */
+    DM dataDM;
+
+   public:
+    explicit ConstAccessor(Vec dataVectorIn, DM dm = nullptr) : dataVector(dataVectorIn), dataDM(dm) {
+        // Get read/write access to the vector
+        VecGetArrayRead(dataVector, &dataArray) >> ablate::utilities::PetscUtilities::checkError;
+
+        // Get the dm from the vector if not specified
+        if (!dataDM) {
+            VecGetDM(dataVector, &dataDM) >> ablate::utilities::PetscUtilities::checkError;
+        }
+    }
+
+    ~ConstAccessor() {
+        // Put back the vector
+        VecRestoreArrayRead(dataVector, &dataArray) >> ablate::utilities::PetscUtilities::checkError;
+    };
+
+    //! Inline function to compute the memory address at this point
+    template <class IndexType>
+    inline const DataType* operator[](IndexType point) const {
+        const DataType* field;
+        if constexpr (isLocal) {
+            DMPlexPointLocalRead(dataDM, point, dataArray, &field) >> ablate::utilities::PetscUtilities::checkError;
+        } else {
+            DMPlexPointGlobalRead(dataDM, point, dataArray, &field) >> ablate::utilities::PetscUtilities::checkError;
+        }
+        return field;
+    }
+
+    //! Inline function to compute the memory address at this point using PETSC style return error
+    template <class IndexType>
+    inline PetscErrorCode operator()(IndexType point, const DataType** field) const {
+        PetscFunctionBeginHot;
+        if constexpr (isLocal) {
+            PetscCall(DMPlexPointLocalRead(dataDM, point, dataArray, field));
+        } else {
+            PetscCall(DMPlexPointGlobalRead(dataDM, point, dataArray, field));
+        }
+        PetscFunctionReturn(PETSC_SUCCESS);
+    }
+
+    /**
+     * prevent copy of this class
+     */
+    ConstAccessor(const ConstAccessor&) = delete;
+};
+}  // namespace ablate::domain
+#endif  // ABLATELIBRARY_FIELDACCESSOR_HPP

src/domain/constFieldAccessor.hpp

@@ -0,0 +1,83 @@
+#ifndef ABLATELIBRARY_CONSTFIELDACCESSOR_HPP
+#define ABLATELIBRARY_CONSTFIELDACCESSOR_HPP
+
+#include <petsc.h>
+#include "field.hpp"
+
+namespace ablate::domain {
+/**
+ * Class responsible for computing const point data locations for dm plex integration
+ */
+template <class DataType, bool isLocal = true>
+class ConstFieldAccessor {
+   private:
+    /**
+     * The petsc vector used to hold the data
+     */
+    Vec dataVector;
+
+    /**
+     * The field information needed to extract the data
+     */
+    const Field& dataField;
+
+    /**
+     * The extracted data from the vector
+     */
+    const PetscScalar* dataArray{};
+
+    /**
+     * Store the DM for the vector to find memory locations
+     */
+    DM dataDM;
+
+   public:
+    ConstFieldAccessor(Vec dataVectorIn, const Field& dataFieldIn, DM dm = nullptr) : dataVector(dataVectorIn), dataField(dataFieldIn), dataDM(dm) {
+        // Get read/write access to the vector
+        VecGetArrayRead(dataVector, &dataArray) >> ablate::utilities::PetscUtilities::checkError;
+
+        // Get the dm from the vector if not specified
+        if (!dataDM) {
+            VecGetDM(dataVector, &dataDM) >> ablate::utilities::PetscUtilities::checkError;
+        }
+    }
+
+    ~ConstFieldAccessor() {
+        // Put back the vector
+        VecRestoreArrayRead(dataVector, &dataArray) >> ablate::utilities::PetscUtilities::checkError;
+    };
+
+    //! provide easy access to the field information
+    [[nodiscard]] inline const Field& GetField() const { return dataField; }
+
+    //! Inline function to compute the memory address at this point
+    template <class IndexType>
+    inline const DataType* operator[](IndexType point) const {
+        const DataType* field;
+        if constexpr (isLocal) {
+            DMPlexPointLocalFieldRead(dataDM, point, dataField.id, dataArray, &field) >> ablate::utilities::PetscUtilities::checkError;
+        } else {
+            DMPlexPointGlobalFieldRead(dataDM, point, dataField.id, dataArray, &field) >> ablate::utilities::PetscUtilities::checkError;
+        }
+        return field;
+    }
+
+    //! Inline function to compute the memory address at this point using PETSC style return error
+    template <class IndexType>
+    inline PetscErrorCode operator()(IndexType point, const DataType** field) const {
+        PetscFunctionBeginHot;
+        if constexpr (isLocal) {
+            PetscCall(DMPlexPointLocalFieldRead(dataDM, point, dataField.id, dataArray, field));
+        } else {
+            PetscCall(DMPlexPointGlobalFieldRead(dataDM, point, dataField.id, dataArray, field));
+        }
+        PetscFunctionReturn(PETSC_SUCCESS);
+    }
+
+    /**
+     * prevent copy of this class
+     */
+    ConstFieldAccessor(const ConstFieldAccessor&) = delete;
+};
+}  // namespace ablate::domain
+#endif  // ABLATELIBRARY_CONSTFIELDACCESSOR_HPP

src/domain/field.cpp

@@ -53,4 +53,4 @@ ablate::domain::Field ablate::domain::Field::CreateSubField(PetscInt newSubId, P
 ablate::domain::Field ablate::domain::Field::Rename(std::string newName) const {
     return ablate::domain::Field{
         .name = std::move(newName), .numberComponents = numberComponents, .components = components, .id = id, .subId = subId, .offset = offset, .location = location, .type = type, .tags = tags};
-}
+}