faceInterpolant.hpp

#ifndef ABLATELIBRARY_FACEINTERPOLANT_HPP
#define ABLATELIBRARY_FACEINTERPOLANT_HPP
 
#include <memory>
#include "domain/range.hpp"
#include "domain/subDomain.hpp"
#include "stencils/stencil.hpp"
 
namespace ablate::finiteVolume {
 
class FaceInterpolant {
   private:
    std::shared_ptr<ablate::domain::SubDomain> subDomain;
 
    PetscInt solTotalSize = 0;
 
    PetscInt auxTotalSize = 0;
 
    PetscInt globalFaceStart = -1;
 
    DM faceSolutionDm = nullptr;
 
    DM faceSolutionGradDm = nullptr;
 
    DM faceAuxDm = nullptr;
 
    DM faceAuxGradDm = nullptr;
 
    static void CreateFaceDm(PetscInt totalDim, DM dm, DM& newDm);
 
    std::vector<stencil::Stencil> stencils;
 
    template <class I, class T>
    static inline void AddToArray(I size, const T* input, T* sum, T factor) {
        for (I d = 0; d < size; d++) {
            sum[d] += factor * input[d];
        }
    }
 
   public:
    FaceInterpolant(const std::shared_ptr<ablate::domain::SubDomain>& subDomain, const std::shared_ptr<domain::Region> solverRegion, Vec faceGeomVec, Vec cellGeomVec);
    ~FaceInterpolant();
 
    using ContinuousFluxFunction = PetscErrorCode (*)(PetscInt dim, const PetscFVFaceGeom* fg, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar field[], const PetscScalar grad[],
                                                      const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar aux[], const PetscScalar gradAux[], PetscScalar flux[], void* ctx);
 
    struct ContinuousFluxFunctionDescription {
        ContinuousFluxFunction function;
        void* context;
 
        PetscInt field;
        std::vector<PetscInt> inputFields;
        std::vector<PetscInt> auxFields;
    };
 
    void ComputeRHS(PetscReal time, Vec locXVec, Vec locAuxVec, Vec locFVec, const std::shared_ptr<domain::Region>& solverRegion,
                    std::vector<FaceInterpolant::ContinuousFluxFunctionDescription>& rhsFunctions, const ablate::domain::Range& faceRange, Vec cellGeomVec, Vec faceGeomVec);
 
    void GetInterpolatedFaceVectors(Vec solutionVec, Vec auxVec, Vec& faceSolutionVec, Vec& faceAuxVec, Vec& faceSolutionGradVec, Vec& faceAuxGradVec);
 
    void RestoreInterpolatedFaceVectors(Vec solutionVec, Vec auxVec, Vec& faceSolutionVec, Vec& faceAuxVec, Vec& faceSolutionGradVec, Vec& faceAuxGradVec);
};
 
}  // namespace ablate::finiteVolume
#endif  // ABLATELIBRARY_FACEINTERPOLANT_HPP