Add support for source terms to cells from particles. (#520)

* particle clang-tidy cleanup

* working on coupled particle solver

* first draft of coupled particle solver

* changed the coupling to be prestep

* prototype code

* added test for coupled particles

* working with the api changes

* adding more fields to the test

* updated for vol mass matrix in FVM field

* added restart test

* version bump

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

src/domain/subDomain.hpp

@@ -257,7 +257,7 @@ class SubDomain : public io::Serializable {
      * Returns raw access to the global dm
      * @return
      */
-    inline DM& GetDM() const noexcept { return domain.GetDM(); }
+    [[nodiscard]] inline DM& GetDM() const noexcept { return domain.GetDM(); }
 
     /**
      * Returns the dm describing the aux fields living in this subdomain.  The dm is defined across

src/particles/CMakeLists.txt

@@ -2,11 +2,13 @@ target_sources(ablateLibrary
         PRIVATE
         particleSolver.cpp
         fieldDescription.cpp
+        coupledParticleSolver.cpp
 
         PUBLIC
         field.hpp
         fieldDescription.hpp
         particleSolver.hpp
+        coupledParticleSolver.hpp
         )
 
 add_subdirectory(initializers)

src/particles/accessors/CMakeLists.txt

@@ -8,4 +8,5 @@ target_sources(ablateLibrary
         swarmAccessor.hpp
         rhsAccessor.hpp
         eulerianAccessor.hpp
+        eulerianSourceAccessor.hpp
         )

src/particles/accessors/accessor.hpp

@@ -29,7 +29,7 @@ class Accessor {
     /**
      * Keep a vector of destructors to call
      */
-    std::vector<std::function<void()>> destructors;
+    std::vector<std::function<void()>> destructors{};
 
    protected:
     /**
@@ -40,7 +40,7 @@ class Accessor {
     virtual Data<DataType> CreateData(const std::string& fieldName) = 0;
 
    public:
-    explicit Accessor(bool cachePointData) : cachePointData(cachePointData), destructors(cachePointData ? 5 : 0) {}
+    explicit Accessor(bool cachePointData) : cachePointData(cachePointData) {}
 
     virtual ~Accessor() {
         for (auto& function : destructors) {
@@ -62,11 +62,10 @@ class Accessor {
      */
     inline Data<DataType> GetData(const std::string& fieldName) {
         if (cachePointData) {
-            if (dataCache.count(fieldName)) {
+            if (!dataCache.count(fieldName)) {
                 dataCache[fieldName] = CreateData(fieldName);
             }
             return dataCache.at(fieldName);
-
         } else {
             return CreateData(fieldName);
         }
@@ -84,4 +83,4 @@ class Accessor {
     Accessor(const Accessor&) = delete;
 };
 }  // namespace ablate::particles::accessors
-#endif  // ABLATELIBRARY_SWARMACCESSOR_HPP
+#endif  // ABLATELIBRARY_PARTICLEACCESSOR_HPP

src/particles/accessors/eulerianAccessor.cpp

@@ -6,11 +6,11 @@
 ablate::particles::accessors::EulerianAccessor::EulerianAccessor(bool cachePointData, std::shared_ptr<ablate::domain::SubDomain> subDomain, SwarmAccessor& swarm, PetscReal currentTime)
     : Accessor(cachePointData), subDomain(std::move(subDomain)), currentTime(currentTime), np(swarm.GetNumberParticles()) {
     // Resize and copy over the coordinates
-    auto coodinatesField = swarm.GetData(ablate::particles::ParticleSolver::ParticleCoordinates);
+    auto coordinatesField = swarm.GetData(ablate::particles::ParticleSolver::ParticleCoordinates);
 
     // Size up and copy the coordinates
-    coordinates.resize(np * coodinatesField.numberComponents);
-    coodinatesField.CopyAll(coordinates.data(), np);
+    coordinates.resize(np * coordinatesField.numberComponents);
+    coordinatesField.CopyAll(coordinates.data(), np);
 }
 
 ablate::particles::accessors::ConstPointData ablate::particles::accessors::EulerianAccessor::CreateData(const std::string& fieldName) {

src/particles/accessors/eulerianAccessor.hpp

@@ -44,7 +44,7 @@ class EulerianAccessor : public Accessor<const PetscReal> {
     /**
      * Get the number of dimensions
      */
-    inline PetscInt GetDimensions() const { return subDomain->GetDimensions(); }
+    [[nodiscard]] inline PetscInt GetDimensions() const { return subDomain->GetDimensions(); }
 };
 }  // namespace ablate::particles::accessors
-#endif  // ABLATELIBRARY_SWARMDATA_HPP
+#endif  // ABLATELIBRARY_EULERIANDATA_HPP

src/particles/accessors/eulerianSourceAccessor.hpp

@@ -0,0 +1,62 @@
+#ifndef ABLATELIBRARY_EULERIANSOURCEACCESSOR_HPP
+#define ABLATELIBRARY_EULERIANSOURCEACCESSOR_HPP
+
+#include <petsc.h>
+#include <map>
+#include "accessor.hpp"
+#include "domain/subDomain.hpp"
+#include "particles/field.hpp"
+#include "swarmAccessor.hpp"
+#include "utilities/petscUtilities.hpp"
+
+namespace ablate::particles::accessors {
+/**
+ * Allows pushing source terms to the particle source array based upon variable/component name.
+ * The coupled solver is used to push back to the flow field ts
+ */
+class EulerianSourceAccessor : public Accessor<PetscReal> {
+   public:
+    // A string to hold the coupled source terms name
+    inline static const char CoupledSourceTermPostfix[] = "_CoupledSourceTerm";
+
+   private:
+    //! borrowed reference to
+    const DM& swarmDm;
+
+    //! a map of fields for easy field lookup
+    const std::map<std::string, Field>& fieldsMap;
+
+   public:
+    EulerianSourceAccessor(bool cachePointData, const DM& swarmDm, const std::map<std::string, Field>& fieldsMap) : Accessor(cachePointData), swarmDm(swarmDm), fieldsMap(fieldsMap) {}
+
+    /**
+     * Create point data from the source field in the DM
+     * @param fieldName
+     * @return
+     */
+    PointData CreateData(const std::string& fieldName) override {
+        const auto& field = fieldsMap.at(fieldName + CoupledSourceTermPostfix);
+        if (field.location == domain::FieldLocation::SOL) {
+            throw std::invalid_argument("Eulerian Source Fields should not be SOL fields");
+        } else {
+            // get the field from the dm
+            PetscScalar* values;
+            DMSwarmGetField(swarmDm, field.name.c_str(), nullptr, nullptr, (void**)&values) >> utilities::PetscUtilities::checkError;
+
+            // Register the cleanup
+            RegisterCleanupFunction([=]() {
+                const std::string name = field.name;
+                DMSwarmRestoreField(swarmDm, name.c_str(), nullptr, nullptr, (void**)&values) >> utilities::PetscUtilities::checkError;
+            });
+
+            return {values, field};
+        }
+    }
+
+    /**
+     * prevent copy of this class
+     */
+    EulerianSourceAccessor(const EulerianSourceAccessor&) = delete;
+};
+}  // namespace ablate::particles::accessors
+#endif  // ABLATELIBRARY_EULERIANSOURCEACCESSOR_HPP

src/particles/accessors/pointData.hpp

@@ -27,7 +27,7 @@ struct Data {
     /**
      * empty default constructor
      */
-    Data() {}
+    Data() = default;
 
     /**
      * The default constructor
@@ -125,4 +125,4 @@ using ConstPointData = Data<const PetscReal>;
 using PointData = Data<PetscReal>;
 
 }  // namespace ablate::particles::accessors
-#endif  // ABLATELIBRARY_SWARMDATA_HPP
+#endif  // ABLATELIBRARY_POINTDATA_HPP

src/particles/accessors/rhsAccessor.hpp

@@ -20,7 +20,7 @@ class RhsAccessor : public Accessor<PetscReal> {
     Vec rhsVec;
 
     //! the array for the rhs values
-    PetscScalar* rhsValues;
+    PetscScalar* rhsValues{};
 
    public:
     RhsAccessor(bool cachePointData, const std::map<std::string, Field>& fieldsMap, Vec rhsVec) : Accessor(cachePointData), fieldsMap(fieldsMap), rhsVec(rhsVec) {
@@ -42,7 +42,7 @@ class RhsAccessor : public Accessor<PetscReal> {
         try {
             const auto& field = fieldsMap.at(fieldName);
             if (field.location == domain::FieldLocation::SOL) {
-                return PointData(rhsValues, field);
+                return {rhsValues, field};
             } else {
                 throw std::invalid_argument("The field " + std::string(fieldName) + " is not a solution variable");
             }
@@ -57,4 +57,4 @@ class RhsAccessor : public Accessor<PetscReal> {
     RhsAccessor(const RhsAccessor&) = delete;
 };
 }  // namespace ablate::particles::accessors
-#endif  // ABLATELIBRARY_SWARMDATA_HPP
+#endif  // ABLATELIBRARY_RHSACCESSOR_HPP

src/particles/accessors/swarmAccessor.hpp

@@ -23,7 +23,7 @@ class SwarmAccessor : public Accessor<const PetscReal> {
     Vec solutionVec;
 
     //! the array for the solution values
-    const PetscScalar* solutionValues;
+    const PetscScalar* solutionValues{};
 
    public:
     SwarmAccessor(bool cachePointData, const DM& swarmDm, const std::map<std::string, Field>& fieldsMap, Vec solutionVec)
@@ -35,15 +35,20 @@ class SwarmAccessor : public Accessor<const PetscReal> {
     ~SwarmAccessor() override { VecRestoreArrayRead(solutionVec, &solutionValues) >> utilities::PetscUtilities::checkError; }
 
     /**
-     * Returns the local size of the particlesdestination
+     * Returns the local size of the particles destination
      * @return
      */
-    inline PetscInt GetNumberParticles() const {
+    [[nodiscard]] inline PetscInt GetNumberParticles() const {
         PetscInt size;
         DMSwarmGetLocalSize(swarmDm, &size) >> utilities::PetscUtilities::checkError;
         return size;
     }
 
+    /**
+     * prevent copy of this class
+     */
+    SwarmAccessor(const SwarmAccessor&) = delete;
+
    protected:
     /**
      * Create point data from the solution field or swarmdm
@@ -53,7 +58,7 @@ class SwarmAccessor : public Accessor<const PetscReal> {
     ConstPointData CreateData(const std::string& fieldName) override {
         const auto& field = fieldsMap.at(fieldName);
         if (field.location == domain::FieldLocation::SOL) {
-            return ConstPointData(solutionValues, field);
+            return {solutionValues, field};
         } else {
             // get the field from the dm
             PetscScalar* values;
@@ -65,14 +70,9 @@ class SwarmAccessor : public Accessor<const PetscReal> {
                 DMSwarmRestoreField(swarmDm, name.c_str(), nullptr, nullptr, (void**)&values) >> utilities::PetscUtilities::checkError;
             });
 
-            return ConstPointData(values, field);
+            return {values, field};
         }
     }
-
-    /**
-     * prevent copy of this class
-     */
-    SwarmAccessor(const SwarmAccessor&) = delete;
 };
 }  // namespace ablate::particles::accessors
 #endif  // ABLATELIBRARY_SWARMACCESSOR_HPP