RadiationDiffusionFDDiscretization.h

Go to the documentation of this file.

#ifndef RAD_DIF_FD_DISCRETIZATION
#define RAD_DIF_FD_DISCRETIZATION
 
#include "AMP/IO/AsciiWriter.h"
#include "AMP/IO/PIO.h"
#include "AMP/discretization/MultiDOF_Manager.h"
#include "AMP/discretization/boxMeshDOFManager.h"
#include "AMP/geometry/shapes/Box.h"
#include "AMP/matrices/MatrixBuilder.h"
#include "AMP/mesh/structured/BoxMesh.h"
#include "AMP/mesh/structured/structuredMeshElement.h"
#include "AMP/operators/LinearOperator.h"
#include "AMP/operators/Operator.h"
#include "AMP/operators/OperatorParameters.h"
#include "AMP/utils/AMPManager.h"
#include "AMP/vectors/MultiVector.h"
#include "AMP/vectors/Vector.h"
#include "AMP/vectors/VectorBuilder.h"
 
#include <optional>
#include <variant>
 
 
namespace AMP::Operator {
 
// Foward declaration of classes declared in this file
class RadDifCoefficients;      // Static class defining PDE coefficients
class FDMeshOps;               // Static class for creating mesh-related data and FD coefficients
class FDMeshGlobalIndexingOps; // Class for index conversions on global DOFs
class FDBoundaryUtils;         // Static class for boundary-related utilities
//
class RadDifOp;               // The main, nonlinear operator
class RadDifOpPJac;           // Its linearization
class RadDifOpPJacParameters; // Operator parameters for linearization
struct RadDifOpPJacData;      // Data structure for storing the linearization data
 
class RadDifCoefficients
{
 
public:
    // Hack. This should really be templated against. Realistically the operator RadDifOp, and its
    // linearization, should be templated against a coefficient class, and a non-templated base
    // class added. The realistic use case is for the nonlinear coefficients, but the linear
    // coefficients can be useful for testing, debugging, etc.
    static constexpr bool IsNonlinear = true;
 
    RadDifCoefficients() = delete;
 
    static double diffusionE( double k11, double T, double zatom );
 
    static double diffusionT( double k21, double T );
 
    static void reaction( double k12,
                          double k22,
                          double T,
                          double zatom,
                          double &REE,
                          double &RET,
                          double &RTE,
                          double &RTT );
};
 
 
class FDMeshOps
{
 
private:
    static constexpr size_t WEST   = 0;
    static constexpr size_t ORIGIN = 1;
    static constexpr size_t EAST   = 2;
 
public:
    // Hack. This should really be templated against. That said, for applications flux limiting is
    // usually required.
    static constexpr bool IsFluxLimited = true;
 
    FDMeshOps() = delete;
 
    static void
    createMeshData( std::shared_ptr<AMP::Mesh::Mesh> mesh,
                    std::shared_ptr<AMP::Mesh::BoxMesh> &d_BoxMesh,
                    size_t &d_dim,
                    AMP::Mesh::GeomType &d_CellCenteredGeom,
                    std::shared_ptr<AMP::Discretization::DOFManager> &d_scalarDOFMan,
                    std::shared_ptr<AMP::Discretization::multiDOFManager> &d_multiDOFMan,
                    std::shared_ptr<AMP::Mesh::BoxMesh::Box> &d_globalBox,
                    std::shared_ptr<AMP::Mesh::BoxMesh::Box> &d_localBox,
                    std::shared_ptr<AMP::ArraySize> &d_localArraySize,
                    std::vector<double> &d_h,
                    std::vector<double> &d_rh2 );
 
    template<bool computeE, bool computeT>
    static void FaceDiffusionCoefficients( std::array<double, 3> &ELoc3,
                                           std::array<double, 3> &TLoc3,
                                           double k11,
                                           double k21,
                                           double zatom,
                                           double h,
                                           double *Dr_WO,
                                           double *Dr_OE,
                                           double *DT_WO,
                                           double *DT_OE );
 
    // Helper function
private:
    static void
    createDOFManagers( const AMP::Mesh::GeomType &Geom,
                       std::shared_ptr<AMP::Mesh::BoxMesh> &mesh,
                       std::shared_ptr<AMP::Discretization::DOFManager> &scalarDOFMan,
                       std::shared_ptr<AMP::Discretization::multiDOFManager> &multiDOFMan );
};
 
 
class FDMeshGlobalIndexingOps
{
 
    // Data
private:
    static constexpr size_t WEST   = 0;
    static constexpr size_t ORIGIN = 1;
    static constexpr size_t EAST   = 2;
 
    std::shared_ptr<AMP::Mesh::BoxMesh> d_BoxMesh;
    AMP::Mesh::GeomType d_geom;
    std::shared_ptr<AMP::Discretization::DOFManager> d_scalarDOFMan;
    std::shared_ptr<AMP::Discretization::multiDOFManager> d_multiDOFMan;
 
 
public:
    FDMeshGlobalIndexingOps( std::shared_ptr<AMP::Mesh::BoxMesh> BoxMesh,
                             AMP::Mesh::GeomType &geom,
                             std::shared_ptr<AMP::Discretization::DOFManager> scalarDOFMan,
                             std::shared_ptr<AMP::Discretization::multiDOFManager> multiDOFMan );
 
    size_t gridIndsToScalarDOF( const std::array<size_t, 3> &ijk ) const;
 
    std::array<size_t, 3> scalarDOFToGridInds( size_t dof ) const;
 
    AMP::Mesh::structuredMeshElement
    gridIndsToMeshElement( const std::array<size_t, 3> &ijk ) const;
};
 
 
class FDBoundaryUtils
{
 
private:
    static constexpr std::string_view a_keys[] = { "a1", "a2", "a3", "a4", "a5", "a6" };
    static constexpr std::string_view b_keys[] = { "b1", "b2", "b3", "b4", "b5", "b6" };
    static constexpr std::string_view r_keys[] = { "r1", "r2", "r3", "r4", "r5", "r6" };
    static constexpr std::string_view n_keys[] = { "n1", "n2", "n3", "n4", "n5", "n6" };
 
    //
public:
    FDBoundaryUtils() = delete;
 
    enum class BoundarySide { WEST, EAST };
 
    static size_t getBoundaryIDFromDim( size_t dim, BoundarySide side );
 
    static size_t getDimFromBoundaryID( size_t boundaryID );
 
    static void getBCConstantsFromDB( const AMP::Database &db,
                                      size_t boundaryID,
                                      double &ak,
                                      double &bk,
                                      double &rk,
                                      double &nk );
 
    static void ghostValuesSolve( double a,
                                  double b,
                                  const std::function<double( double T )> &cHandle,
                                  double r,
                                  double n,
                                  double h,
                                  double Eint,
                                  double Tint,
                                  double &Eg,
                                  double &Tg );
 
    static double ghostValueSolveT( double n, double h, double Tint );
 
    static double ghostValueSolveE( double a, double b, double r, double c, double h, double Eint );
};
 
 
class RadDifOp : public AMP::Operator::Operator
{
 
    // Data
public:
    std::shared_ptr<AMP::Database> d_db = nullptr;
 
 
    // Methods
public:
    RadDifOp( std::shared_ptr<const AMP::Operator::OperatorParameters> params );
    virtual ~RadDifOp(){};
 
    std::string type() const override { return "RadDifOp"; };
 
    void apply( std::shared_ptr<const AMP::LinearAlgebra::Vector> ET,
                std::shared_ptr<AMP::LinearAlgebra::Vector> LET ) override;
 
    bool isValidVector( std::shared_ptr<const AMP::LinearAlgebra::Vector> ET ) override;
 
    std::shared_ptr<AMP::LinearAlgebra::Vector> createInputVector() const override;
 
    const std::vector<double> &getMeshSize() const;
 
    AMP::Mesh::GeomType getGeomType() const;
 
    std::shared_ptr<const AMP::Discretization::DOFManager> getScalarDOFManager() const;
 
    void setBoundaryFunctionE(
        const std::function<double( size_t boundaryID, const AMP::Mesh::Point &boundaryPoint )>
            &fn_ );
 
    void setBoundaryFunctionT(
        const std::function<double( size_t boundaryID, const AMP::Mesh::Point &boundaryPoint )>
            &fn_ );
 
    // Data
private:
    const double d_k11;
    const double d_k12;
    const double d_k21;
    const double d_k22;
 
    std::array<double, 6> d_ak;
    std::array<double, 6> d_bk;
    std::array<double, 6> d_rk;
    std::array<double, 6> d_nk;
 
private:
    std::shared_ptr<AMP::Discretization::multiDOFManager> d_multiDOFMan;
    std::shared_ptr<AMP::Discretization::DOFManager> d_scalarDOFMan;
    std::shared_ptr<AMP::Mesh::BoxMesh> d_BoxMesh;
    std::shared_ptr<AMP::Mesh::BoxMesh::Box> d_globalBox = nullptr;
    std::shared_ptr<AMP::Mesh::BoxMesh::Box> d_localBox = nullptr;
    std::shared_ptr<AMP::ArraySize> d_localArraySize = nullptr;
    AMP::Mesh::GeomType CellCenteredGeom;
    size_t d_dim = (size_t) -1;
    std::vector<double> d_h;
    std::vector<double> d_rh2;
 
    static constexpr size_t WEST   = 0;
    static constexpr size_t ORIGIN = 1;
    static constexpr size_t EAST   = 2;
 
    std::shared_ptr<FDMeshGlobalIndexingOps> d_meshIndexingOps = nullptr;
 
 
private:
    void applyInterior( std::shared_ptr<const AMP::LinearAlgebra::Vector> E_vec,
                        std::shared_ptr<const AMP::LinearAlgebra::Vector> T_vec,
                        std::shared_ptr<AMP::LinearAlgebra::Vector> LE_vec,
                        std::shared_ptr<AMP::LinearAlgebra::Vector> LT_vec );
 
    void applyBoundary( std::shared_ptr<const AMP::LinearAlgebra::Vector> E_vec,
                        std::shared_ptr<const AMP::LinearAlgebra::Vector> T_vec,
                        std::shared_ptr<AMP::LinearAlgebra::Vector> LE_vec,
                        std::shared_ptr<AMP::LinearAlgebra::Vector> LT_vec );
 
 
    void ghostValuesSolveWrapper( size_t boundaryID,
                                  const AMP::Mesh::Point &boundaryPoint,
                                  double Eint,
                                  double Tint,
                                  double &Eg,
                                  double &Tg );
 
    void getNNDataBoundary( std::shared_ptr<const AMP::LinearAlgebra::Vector> E_vec,
                            std::shared_ptr<const AMP::LinearAlgebra::Vector> T_vec,
                            std::array<size_t, 3> &ijk,
                            size_t dim,
                            std::array<double, 3> &ELoc3,
                            std::array<double, 3> &TLoc3 );
 
 
    std::function<double( size_t boundaryID, const AMP::Mesh::Point &boundaryPoint )>
        d_robinFunctionE;
 
    std::function<double( size_t boundaryID, const AMP::Mesh::Point &boundaryPoint )>
        d_pseudoNeumannFunctionT;
 
    //
protected:
    std::shared_ptr<AMP::Operator::OperatorParameters>
    getJacobianParameters( AMP::LinearAlgebra::Vector::const_shared_ptr u_in ) override;
};
// End of RadDifOp
 
 
class RadDifOpPJac : public AMP::Operator::LinearOperator
{
 
    //
public:
    std::shared_ptr<AMP::Database> d_db = nullptr;
    std::shared_ptr<RadDifOpPJacData> d_data = nullptr;
    AMP::LinearAlgebra::Vector::shared_ptr d_frozenVec = nullptr;
 
    RadDifOpPJac( std::shared_ptr<const AMP::Operator::OperatorParameters> params_ );
 
    virtual ~RadDifOpPJac(){};
 
    std::shared_ptr<AMP::LinearAlgebra::Vector> createInputVector() const override;
 
    static std::unique_ptr<AMP::Operator::Operator>
    create( std::shared_ptr<AMP::Operator::OperatorParameters> params )
    {
        return std::make_unique<RadDifOpPJac>( params );
    };
 
    void reset( std::shared_ptr<const AMP::Operator::OperatorParameters> params ) override;
 
    std::string type() const override { return "RadDifOpPJac"; };
 
    void apply( std::shared_ptr<const AMP::LinearAlgebra::Vector> ET,
                std::shared_ptr<AMP::LinearAlgebra::Vector> LET ) override;
 
    void applyWithOverwrittenDataIsValid() { d_applyWithOverwrittenDataIsValid = true; };
 
 
    // Data
private:
    const double d_k11;
    const double d_k12;
    const double d_k21;
    const double d_k22;
 
    std::array<double, 6> d_ak;
    std::array<double, 6> d_bk;
    std::array<double, 6> d_rk;
    std::array<double, 6> d_nk;
 
    bool d_applyWithOverwrittenDataIsValid = false;
 
 
private:
    std::shared_ptr<AMP::Discretization::multiDOFManager> d_multiDOFMan;
    std::shared_ptr<AMP::Discretization::DOFManager> d_scalarDOFMan;
    std::shared_ptr<AMP::Mesh::BoxMesh> d_BoxMesh;
    std::shared_ptr<AMP::Mesh::BoxMesh::Box> d_globalBox = nullptr;
    std::shared_ptr<AMP::Mesh::BoxMesh::Box> d_localBox = nullptr;
    std::shared_ptr<AMP::ArraySize> d_localArraySize = nullptr;
    AMP::Mesh::GeomType CellCenteredGeom;
    size_t d_dim = (size_t) -1;
    std::vector<double> d_h;
    std::vector<double> d_rh2;
 
    std::array<size_t, 5> d_ijk;
    std::array<double, 3> d_ELoc3;
    std::array<double, 3> d_TLoc3;
    std::array<size_t, 3> d_dofsLoc3;
 
    static constexpr size_t WEST   = 0;
    static constexpr size_t ORIGIN = 1;
    static constexpr size_t EAST   = 2;
 
    std::shared_ptr<FDMeshGlobalIndexingOps> d_meshIndexingOps = nullptr;
    //
private:
    void applyFromData( std::shared_ptr<const AMP::LinearAlgebra::Vector> ET_,
                        std::shared_ptr<AMP::LinearAlgebra::Vector> LET_ );
 
    void setData();
 
    void setDataReaction( std::shared_ptr<const AMP::LinearAlgebra::Vector> T_vec );
 
    template<size_t Component>
    void fillDiffusionMatrixWithData( std::shared_ptr<AMP::LinearAlgebra::Matrix> matrix );
 
    template<size_t Component>
    void getCSRDataDiffusionMatrix( std::shared_ptr<const AMP::LinearAlgebra::Vector> E_vec,
                                    std::shared_ptr<const AMP::LinearAlgebra::Vector> T_vec,
                                    const double *E_rawData,
                                    const double *T_rawData,
                                    size_t row,
                                    std::vector<size_t> &cols,
                                    std::vector<double> &data );
 
    template<size_t Component>
    void getCSRDataDiffusionMatrixInterior( const double *E_rawData,
                                            const double *T_rawData,
                                            size_t rowLocal,
                                            std::array<size_t, 5> &ijkLocal,
                                            std::vector<size_t> &colsLocal,
                                            std::vector<double> &data );
 
    template<size_t Component>
    void getCSRDataDiffusionMatrixBoundary( std::shared_ptr<const AMP::LinearAlgebra::Vector> E_vec,
                                            std::shared_ptr<const AMP::LinearAlgebra::Vector> T_vec,
                                            size_t row,
                                            std::vector<size_t> &cols,
                                            std::vector<double> &data );
 
 
    double PicardCorrectionCoefficient( size_t component, size_t boundaryID, double ck ) const;
 
    void ghostValuesSolveWrapper( size_t boundaryID,
                                  const AMP::Mesh::Point &boundaryPoint,
                                  double Eint,
                                  double Tint,
                                  double &Eg,
                                  double &Tg );
 
    void getNNDataBoundary(
        std::shared_ptr<const AMP::LinearAlgebra::Vector> E_vec,
        std::shared_ptr<const AMP::LinearAlgebra::Vector> T_vec,
        std::array<size_t, 3> &ijk, // is modified locally, but returned in same state
        size_t dim,
        std::array<double, 3> &ELoc3,
        std::array<double, 3> &TLoc3,
        std::array<size_t, 3> &dofsLoc3,
        std::optional<FDBoundaryUtils::BoundarySide> &boundaryIntersection );
 
    std::function<double( size_t boundaryID, const AMP::Mesh::Point &boundaryPoint )>
        d_robinFunctionE;
 
    std::function<double( size_t boundaryID, const AMP::Mesh::Point &boundaryPoint )>
        d_pseudoNeumannFunctionT;
};
// End of RadDifOpPJac
 
 
struct RadDifOpPJacData {
    friend class RadDifOpPJac;
 
public:
    std::tuple<std::shared_ptr<AMP::LinearAlgebra::Matrix>,
               std::shared_ptr<AMP::LinearAlgebra::Matrix>,
               std::shared_ptr<AMP::LinearAlgebra::Vector>,
               std::shared_ptr<AMP::LinearAlgebra::Vector>,
               std::shared_ptr<AMP::LinearAlgebra::Vector>,
               std::shared_ptr<AMP::LinearAlgebra::Vector>>
    get();
 
private:
    bool d_dataMaybeOverwritten = false;
 
    std::shared_ptr<AMP::LinearAlgebra::Matrix> d_E  = nullptr;
    std::shared_ptr<AMP::LinearAlgebra::Matrix> d_T  = nullptr;
    std::shared_ptr<AMP::LinearAlgebra::Vector> r_EE = nullptr;
    std::shared_ptr<AMP::LinearAlgebra::Vector> r_ET = nullptr;
    std::shared_ptr<AMP::LinearAlgebra::Vector> r_TE = nullptr;
    std::shared_ptr<AMP::LinearAlgebra::Vector> r_TT = nullptr;
};
 
 
class RadDifOpPJacParameters : public AMP::Operator::OperatorParameters
{
public:
    // Constructor
    explicit RadDifOpPJacParameters( std::shared_ptr<AMP::Database> db )
        : OperatorParameters( db ){};
    virtual ~RadDifOpPJacParameters(){};
 
    AMP::LinearAlgebra::Vector::shared_ptr d_frozenSolution = nullptr;
 
    std::function<double( size_t boundaryID, const AMP::Mesh::Point &boundaryPoint )>
        d_robinFunctionE;
    std::function<double( size_t boundaryID, const AMP::Mesh::Point &boundaryPoint )>
        d_pseudoNeumannFunctionT;
};
 
} // namespace AMP::Operator
 
 
#endif