pnsolver.h 4.84 KB
Newer Older
1
2
3
#ifndef PNSOLVER_H
#define PNSOLVER_H

4
#include "solverbase.h"
5
6
7
8

class PNSolver : public Solver
{
  public:
9
10
    /*! @brief PNSolver constructor
     *  @param settings stores all needed information
11
12
13
     */
    PNSolver( Config* settings );

14
15
16
    virtual void Solve() override;              /*! @brief Solve functions runs main time loop. (Run Solver) */
    void Save() const override;                 /*! @brief Save Output solution to VTK file */
    void Save( int currEnergy ) const override; /*! @brief Save Output solution at given energy (pseudo time) to VTK file */
steffen.schotthoefer's avatar
steffen.schotthoefer committed
17

18
  protected:
19
20
    unsigned _nTotalEntries; /*! @brief: total number of equations in the system */
    unsigned _LMaxDegree;    /*! @brief: maximal degree of the spherical harmonics basis*/
21

22
    VectorVector _sigmaA; /*!  @brief: Absorption coefficient for all energies*/
23
24
    // System Matrix for x, y and z flux
    //    ==> not needed after computation of A+ and A- ==> maybe safe only temporarly and remove as member?
25
26
27
    SymMatrix _Ax; /*! @brief:  Flux Jacbioan in x direction */
    SymMatrix _Ay; /*! @brief:  Flux Jacbioan in x direction */
    SymMatrix _Az; /*! @brief:  Flux Jacbioan in x direction */
28
29

    // Upwinding Matrices
30
31
32
33
34
35
36
37
38
39
40
41
42
43
    Matrix _AxPlus;  /*! @brief:  Flux Jacbioan in x direction, positive part */
    Matrix _AxMinus; /*! @brief:  Flux Jacbioan in x direction, negative part */
    Matrix _AxAbs;   /*! @brief:  Flux Jacbioan in x direction, absolute part */
    Matrix _AyPlus;  /*! @brief:  Flux Jacbioan in y direction, positive part */
    Matrix _AyMinus; /*! @brief:  Flux Jacbioan in y direction, negative part */
    Matrix _AyAbs;   /*! @brief:  Flux Jacbioan in y direction, absolute part */
    Matrix _AzPlus;  /*! @brief:  Flux Jacbioan in z direction, positive part */
    Matrix _AzMinus; /*! @brief:  Flux Jacbioan in z direction, negative part */
    Matrix _AzAbs;   /*! @brief:  Flux Jacbioan in z direction, absolute part */

    // double _combinedSpectralRadius; /*! @brief:  Combined spectral radius of sum of flux jacobians*/

    Vector _scatterMatDiag; /*! @brief: diagonal of the scattering matrix (its a diagonal matrix by construction) */

44
    // Output related members
45
46
47
    std::vector<std::vector<std::vector<double>>> _outputFields; /*! @brief: Solver Output: dimensions (GroupID,FieldID,CellID). !Protoype output for
                                                                    multiple output fields. Will replace _solverOutput */
    std::vector<std::vector<std::string>> _outputFieldNames;     /*! @brief: Names of the outputFields: dimensions (GroupID,FieldID) */
48

49
50
51
52
53
54
55
56
57
58
59
60
    // ---- Member functions ----

    // IO
    /*! @brief Initializes the output groups and fields of this solver and names the fields */
    void PrepareOutputFields();

    /*! @brief Function that prepares VTK export and csv export of the current solver iteration
        @returns: Mass of current iteration
    */
    double WriteOutputFields();

    // Solver
61
62
63
64
    /*! @brief: parameter functions for setting up system matrix
     *  @param: degree l, it must hold: 0 <= l <=_nq
     *  @param : order k, it must hold: -l <=k <= l
     */
65
66
67
68
69
70
71
72
73
74
75
76
    double AParam( int l, int k ) const;
    double BParam( int l, int k ) const;
    double CParam( int l, int k ) const;
    double DParam( int l, int k ) const;
    double EParam( int l, int k ) const;
    double FParam( int l, int k ) const;

    double CTilde( int l, int k ) const;
    double DTilde( int l, int k ) const;
    double ETilde( int l, int k ) const;
    double FTilde( int l, int k ) const;

77
78
79
    /*! @brief: mathematical + index functions. Helper functions for setting up system matrix.
     *  @param: k: arbitrary integer
     */
80
81
82
    int Sgn( int k ) const;
    int kPlus( int k ) const;
    int kMinus( int k ) const;
steffen.schotthoefer's avatar
steffen.schotthoefer committed
83
84
85
86
87
88

    /*! @brief : computes the global index of the moment corresponding to basis function (l,k)
     *  @param : degree l, it must hold: 0 <= l <=_nq
     *  @param : order k, it must hold: -l <=k <= l
     *  @returns : global index
     */
89
    int GlobalIndex( int l, int k ) const;
90

steffen.schotthoefer's avatar
steffen.schotthoefer committed
91
92
93
94
95
    /*! @brief : Checks, if index invariant for global index holds
     *  @returns : True, if invariant holds, else false
     */
    bool CheckIndex( int l, int k ) const;

96
    /*! @brief: function for computing and setting up system matrices */
97
    void ComputeSystemMatrices();
98
    /*! @brief:  function for computing and setting up flux matrices for upwinding */
99
    void ComputeFluxComponents();
100
    /*! @brief:  fucntion for computing and setting up diagonal matrix for scattering kernel */
101
    void ComputeScatterMatrix();
102
    /*! @brief:  Computes Legedre polinomial of oder l at point x */
103
    double LegendrePoly( double x, int l );
104
105
106
107

    /*! @brief: Sets Entries of FluxMatrices to zero, if they are below double precision,
     *          to prevent floating point inaccuracies later in the solver
     */
108
    void CleanFluxMatrices();
109
};
110
111

#endif    // PNSOLVER_H