particleSolver.hpp

#ifndef ABLATELIBRARY_PARTICLESOLVER_HPP
#define ABLATELIBRARY_PARTICLESOLVER_HPP
 
#include "field.hpp"
#include "fieldDescription.hpp"
#include "initializers/initializer.hpp"
#include "processes/process.hpp"
#include "solver/solver.hpp"
 
namespace ablate::particles {
 
class ParticleSolver : public solver::Solver, public io::Serializable {
   public:
    inline static const char ParticleVelocity[] = "ParticleVelocity";
    inline static const char ParticleDiameter[] = "ParticleDiameter";
    inline static const char ParticleDensity[] = "ParticleDensity";
    inline static const char PackedSolution[] = "PackedSolution";
    inline static const char ParticleInitialLocation[] = "InitialLocation";
 
    inline static const char ParticleCoordinates[] = "coordinates";
 
   protected:
    DM swarmDm = nullptr;
 
    TS particleTs = nullptr;
 
    PetscReal timeInitial = 0.0;
 
    PetscReal timeFinal = 0.0;
 
    PetscInt ndims = 0;
 
    bool dmChanged = false;
 
    std::vector<FieldDescription> fieldsDescriptions;
 
    std::vector<Field> fields;
 
    std::map<std::string, Field> fieldsMap;
 
    std::vector<std::shared_ptr<processes::Process>> processes;
 
    std::shared_ptr<particles::initializers::Initializer> initializer = nullptr;
 
    const std::vector<std::shared_ptr<mathFunctions::FieldFunction>> fieldInitialization;
 
    const std::vector<std::shared_ptr<mathFunctions::FieldFunction>> exactSolutions;
 
   public:
    ParticleSolver(std::string solverId, std::shared_ptr<domain::Region>, std::shared_ptr<parameters::Parameters> options, std::vector<FieldDescription> fields,
                   std::vector<std::shared_ptr<processes::Process>> processes, std::shared_ptr<initializers::Initializer> initializer,
                   std::vector<std::shared_ptr<mathFunctions::FieldFunction>> fieldInitialization, std::vector<std::shared_ptr<mathFunctions::FieldFunction>> exactSolutions = {});
 
    ParticleSolver(std::string solverId, std::shared_ptr<domain::Region>, std::shared_ptr<parameters::Parameters> options, const std::vector<std::shared_ptr<FieldDescription>>& fields,
                   std::vector<std::shared_ptr<processes::Process>> processes, std::shared_ptr<initializers::Initializer> initializer,
                   std::vector<std::shared_ptr<mathFunctions::FieldFunction>> fieldInitialization, std::vector<std::shared_ptr<mathFunctions::FieldFunction>> exactSolutions = {});
 
    ~ParticleSolver() override;
 
    void Setup() override;
 
    /*** Set up mesh dependent initialization, this may be called multiple times if the mesh changes **/
    void Initialize() override;
 
    inline DM GetParticleDM() { return swarmDm; }
 
    inline TS GetParticleTS() { return particleTs; }
 
    static PetscErrorCode ComputeParticleExactSolution(TS particleTS, Vec);
 
    [[nodiscard]] const std::string& GetId() const override { return GetSolverId(); }
 
    PetscErrorCode Save(PetscViewer viewer, PetscInt steps, PetscReal time) override;
 
    PetscErrorCode Restore(PetscViewer viewer, PetscInt steps, PetscReal time) override;
 
   protected:
    virtual void MacroStepParticles(TS macroTS, bool swarmMigrate);
 
    void RegisterParticleField(const FieldDescription& fieldDescriptor);
 
    void StoreInitialParticleLocations();
 
    static PetscErrorCode ComputeParticleError(TS particleTS, Vec u, Vec e);
 
    void ProjectFunction(const std::shared_ptr<mathFunctions::FieldFunction>& fieldFunction, PetscReal time = 0.0);
 
    void SwarmMigrate();
 
    void CoordinatesToSolutionVector();
 
    void CoordinatesFromSolutionVector();
 
   protected:
    template <class T>
    void GetField(const Field& field, T** values) {
        if (field.location == domain::FieldLocation::SOL) {
            // Get the solution vector
            DMSwarmGetField(swarmDm, PackedSolution, nullptr, nullptr, (void**)values) >> utilities::PetscUtilities::checkError;
        } else {
            // get the raw field
            DMSwarmGetField(swarmDm, field.name.c_str(), nullptr, nullptr, (void**)values) >> utilities::PetscUtilities::checkError;
        }
    }
 
    template <class T>
    void RestoreField(const Field& field, T** values) {
        if (field.location == domain::FieldLocation::SOL) {
            // Get the solution vector
            DMSwarmRestoreField(swarmDm, PackedSolution, nullptr, nullptr, (void**)values) >> utilities::PetscUtilities::checkError;
        } else {
            // get the raw field
            DMSwarmRestoreField(swarmDm, field.name.c_str(), nullptr, nullptr, (void**)values) >> utilities::PetscUtilities::checkError;
        }
    }
 
    template <class T>
    const Field& GetField(const std::string& fieldName, T** values) {
        const auto& field = GetField(fieldName);
        GetField(field, values);
        return field;
    }
 
    template <class T>
    void RestoreField(const std::string& fieldName, T** values) {
        const auto& field = GetField(fieldName);
        RestoreField(field, values);
    }
 
    [[nodiscard]] const Field& GetField(const std::string& fieldName) const { return fieldsMap.at(fieldName); }
 
    static PetscErrorCode ComputeParticleRHS(TS ts, PetscReal t, Vec X, Vec F, void* ctx);
};
 
}  // namespace ablate::particles
#endif  // ABLATELIBRARY_PARTICLESOLVER_HPP