boundarySolver.hpp

#ifndef ABLATELIBRARY_BOUNDARYSOLVER_HPP
#define ABLATELIBRARY_BOUNDARYSOLVER_HPP
 
#include <memory>
#include "solver/cellSolver.hpp"
#include "solver/timeStepper.hpp"
 
namespace ablate::boundarySolver {
 
// forward declare the boundaryProcess
class BoundaryProcess;
 
class BoundarySolver : public solver::CellSolver, public solver::RHSFunction, private utilities::Loggable<BoundarySolver>, public io::Serializable {
   public:
    typedef struct {
        PetscInt faceId;       /* local id for this face.  For merged faces this is the first face merged **/
        PetscReal normal[3];   /* normals (pointing into the boundary from the other region) */
        PetscReal areas[3];    /* Area-scaled normals */
        PetscReal centroid[3]; /* Location of centroid (quadrature point) */
    } BoundaryFVFaceGeom;
 
    using BoundarySourceFunction = PetscErrorCode (*)(PetscInt dim, const BoundaryFVFaceGeom* fg, const PetscFVCellGeom* boundaryCell, const PetscInt uOff[], const PetscScalar* boundaryValues,
                                                      const PetscScalar* stencilValues[], const PetscInt aOff[], const PetscScalar* auxValues, const PetscScalar* stencilAuxValues[],
                                                      PetscInt stencilSize, const PetscInt stencil[], const PetscScalar stencilWeights[], const PetscInt sOff[], PetscScalar source[], void* ctx);
 
    using BoundaryUpdateFunction = PetscErrorCode (*)(PetscInt dim, const BoundaryFVFaceGeom* fg, const PetscFVCellGeom* boundaryCell, const PetscInt uOff[], PetscScalar* boundaryValues,
                                                      const PetscScalar* stencilValues, const PetscInt aOff[], PetscScalar* auxValues, const PetscScalar* stencilAuxValues, void* ctx);
 
    using BoundaryPreRHSFunctionDefinition = PetscErrorCode (*)(BoundarySolver&, TS ts, PetscReal time, bool initialStage, Vec locX, void* ctx);
 
    enum class BoundarySourceType {
        Point,       
        Distributed, 
        Flux,        
        Face         
    };
 
    static void ComputeGradient(PetscInt dim, PetscScalar boundaryValue, PetscInt stencilSize, const PetscScalar* stencilValues, const PetscScalar* stencilWeights, PetscScalar* grad);
 
    static void ComputeGradientAlongNormal(PetscInt dim, const BoundaryFVFaceGeom* fg, PetscScalar boundaryValue, PetscInt stencilSize, const PetscScalar* stencilValues,
                                           const PetscScalar* stencilWeights, PetscScalar& dPhiDNorm);
 
    using BoundaryPreRHSPointFunction = PetscErrorCode (*)(PetscReal time, PetscReal dt, PetscInt dim, const BoundaryFVFaceGeom* fg, const PetscFVCellGeom* boundaryCell, const PetscInt uOff[],
                                                           PetscScalar* boundaryValues, const PetscScalar* stencilValues[], const PetscInt aOff[], PetscScalar* auxValues,
                                                           const PetscScalar* stencilAuxValues[], PetscInt stencilSize, const PetscInt stencil[], const PetscScalar stencilWeights[], void* ctx);
 
    struct BoundaryPreRHSPointFunctionDefinition {
        BoundaryPreRHSPointFunction function;
        void* context;
 
        std::vector<PetscInt> inputFieldsOffset;
        std::vector<PetscInt> auxFieldsOffset;
    };
 
    PetscErrorCode ComputeBoundaryPreRHSPointFunction(PetscReal time, PetscReal dt, Vec locXVec, const BoundaryPreRHSPointFunctionDefinition& boundaryPreRhsPointFunction);
 
   private:
    struct GradientStencil {
        PetscInt cellId;
        BoundaryFVFaceGeom geometry;
        std::vector<PetscInt> stencil;
        PetscInt stencilSize;
        std::vector<PetscScalar> gradientWeights;
        std::vector<PetscScalar> distributionWeights;
        std::vector<PetscScalar> volumes;
    };
 
    struct BoundarySourceFunctionDescription {
        BoundarySourceFunction function;
        void* context;
        BoundarySourceType type;
 
        std::vector<PetscInt> sourceFieldsOffset;
        std::vector<PetscInt> inputFieldsOffset;
        std::vector<PetscInt> auxFieldsOffset;
    };
 
    struct BoundaryUpdateFunctionDescription {
        BoundaryUpdateFunction function;
        void* context;
 
        std::vector<PetscInt> inputFields;
        std::vector<PetscInt> auxFields;
    };
 
    // Hold the region used to define the boundary faces
    const std::shared_ptr<domain::Region> fieldBoundary;
 
    // hold the update functions for flux and point sources
    std::vector<BoundarySourceFunctionDescription> boundarySourceFunctions;
 
    // boundary output functions that can be used for
    std::vector<BoundarySourceFunctionDescription> boundaryOutputFunctions;
 
    // keep track of the output field components
    std::vector<std::string> outputComponents;
 
    // hold the update functions for flux and point sources
    std::vector<BoundaryUpdateFunctionDescription> boundaryUpdateFunctions;
 
    // Hold a list of boundaryProcesses that contribute to this solver
    std::vector<std::shared_ptr<BoundaryProcess>> boundaryProcesses;
 
    // allow the use of any arbitrary pre rhs functions
    std::vector<std::pair<BoundaryPreRHSFunctionDefinition, void*>> preRhsFunctions;
 
   protected:
    // Hold a list of GradientStencils, this order corresponds to the face order
    std::vector<GradientStencil> gradientStencils;
 
   private:
    // keep track of maximumStencilSize
    PetscInt maximumStencilSize = 0;
 
    // The PetscFV (usually the least squares method) is used to compute the gradient weights
    PetscFV gradientCalculator = nullptr;
 
    // Determine if multiple faces should be merged for a single cell
    const bool mergeFaces;
 
    void CreateGradientStencil(PetscInt cellId, const BoundaryFVFaceGeom& geometry, const std::vector<PetscInt>& stencil, DM cellDM, const PetscScalar* cellGeomArray);
 
    void UpdateVariablesPreStep(TS ts, ablate::solver::Solver&);
 
   public:
    BoundarySolver(std::string solverId, std::shared_ptr<domain::Region> region, std::shared_ptr<domain::Region> fieldBoundary, std::vector<std::shared_ptr<BoundaryProcess>> boundaryProcesses,
                   std::shared_ptr<parameters::Parameters> options, bool mergeFaces = false);
    ~BoundarySolver() override;
 
    void Setup() override;
    void Initialize() override;
 
    void RegisterFunction(BoundarySourceFunction function, void* context, const std::vector<std::string>& sourceFields, const std::vector<std::string>& inputFields,
                          const std::vector<std::string>& auxFields, BoundarySourceType type = BoundarySourceType::Point);
 
    void RegisterFunction(BoundaryUpdateFunction function, void* context, const std::vector<std::string>& inputFields, const std::vector<std::string>& auxFields);
 
    void RegisterPreRHSFunction(BoundaryPreRHSFunctionDefinition function, void* context);
 
    PetscErrorCode ComputeRHSFunction(PetscReal time, Vec locXVec, Vec locFVec) override;
 
    PetscErrorCode ComputeRHSFunction(PetscReal time, Vec locXVec, Vec locFVec, const std::vector<BoundarySourceFunctionDescription>& boundarySourceFunctions);
 
    void InsertFieldFunctions(const std::vector<std::shared_ptr<mathFunctions::FieldFunction>>& initialization, PetscReal time = 0.0);
 
    [[nodiscard]] std::vector<GradientStencil> GetBoundaryGeometry(PetscInt cell) const;
 
    const std::vector<GradientStencil>& GetBoundaryGeometry() const { return gradientStencils; }
 
    inline const std::vector<std::string>& GetOutputComponents() { return outputComponents; }
 
    inline const std::vector<BoundarySourceFunctionDescription>& GetOutputFunctions() { return boundaryOutputFunctions; }
 
    PetscErrorCode PreRHSFunction(TS ts, PetscReal time, bool initialStage, Vec locX) override;
 
    [[nodiscard]] SerializerType Serialize() const override;
 
    [[nodiscard]] const std::string& GetId() const override { return GetSolverId(); }
 
    PetscErrorCode Save(PetscViewer viewer, PetscInt sequenceNumber, PetscReal time) override;
 
    PetscErrorCode Restore(PetscViewer viewer, PetscInt sequenceNumber, PetscReal time) override;
};
 
std::istream& operator>>(std::istream& is, BoundarySolver::BoundarySourceType& v);
 
}  // namespace ablate::boundarySolver
#endif  // ABLATELIBRARY_BOUNDARYSOLVER_HPP