solverbase.h 7.02 KB
Newer Older
Jonas Kusch's avatar
Jonas Kusch committed
1
2
3
4
#ifndef SOLVER_H
#define SOLVER_H

// include Matrix, Vector definitions
5
6
#include "common/globalconstants.h"
#include "common/typedef.h"
7

8
9
10
// externals
#include "spdlog/spdlog.h"

steffen.schotthoefer's avatar
steffen.schotthoefer committed
11
// Forward Declarations
12
13
14
15
16
class NumericalFlux;
class Mesh;
class Config;
class ProblemBase;
class QuadratureBase;
17
class Reconstructor;
steffen.schotthoefer's avatar
steffen.schotthoefer committed
18

Jonas Kusch's avatar
Jonas Kusch committed
19
20
class Solver
{
21
  protected:
22
23
24
25
    Mesh* _mesh;           /*! @brief mesh object for writing out information */
    NumericalFlux* _g;     /*! @brief class for numerical flux */
    Config* _settings;     /*! @brief config class for global information */
    ProblemBase* _problem; /*! @brief problem class for initial conditions */
steffen.schotthoefer's avatar
steffen.schotthoefer committed
26

27
28
    // --------- Often used variables of member classes for faster access ----

jannick.wolters's avatar
jannick.wolters committed
29
30
31
32
33
34
    unsigned _nEnergies;             /*! @brief number of energy/time steps, number of nodal energy values for CSD */
    double _dE;                      /*! @brief energy/time step size */
    Vector _energies;                // energy groups used in the simulation [keV]
    std::vector<double> _density;    // patient density, dim(_density) = _nCells
    Vector _s;                       // stopping power, dim(_s) = _nTimeSteps
    std::vector<VectorVector> _Q;    /*!  @brief  external source term */
35
36
37
38
39

    VectorVector _sigmaS; /*!  @brief scattering cross section for all energies */
    VectorVector _sigmaT; /*!  @brief total cross section for all energies */

    // quadrature related numbers
40
41
42
43
44
    QuadratureBase* _quadrature; /*! @brief quadrature to create members below */
    unsigned _nq;                /*! @brief number of quadrature points */

    // VectorVector _quadPoints;    /*!  @brief quadrature points, dim(_quadPoints) = (_nSystem,spatialDim) */
    // Vector _weights;             /*!  @brief quadrature weights, dim(_weights) = (_NCells) */
Jonas Kusch's avatar
Jonas Kusch committed
45

46
47
48
49
50
51
52
53
54
    // Mesh related members
    unsigned _nCells;                          /*! @brief number of spatial cells */
    std::vector<BOUNDARY_TYPE> _boundaryCells; /*! boundary type for all cells, dim(_boundary) = (_NCells) */
    std::vector<double> _areas;                /*! @brief surface area of all spatial cells, dim(_areas) = _NCells */
    /*! @brief edge normals multiplied by edge length, dim(_normals) = (_NCells,nEdgesPerCell,spatialDim) */
    std::vector<std::vector<Vector>> _normals;
    /*! @brief edge neighbor cell ids, dim(_neighbors) = (_NCells,nEdgesPerCell) */
    std::vector<std::vector<unsigned>> _neighbors;

Tianbai Xiao's avatar
Tianbai Xiao committed
55
    // slope related params
56
57
58
59
60
61
    Reconstructor* _reconstructor; /*! @brief reconstructor object for high-order scheme */
    unsigned _reconsOrder; /*! @brief reconstruction order (current: 1 & 2) */
    VectorVector _psiDx; /*! @brief slope of solutions in X direction */
    VectorVector _psiDy; /*! @brief slope of solutions in Y direction */
    VectorVector _cellMidPoints; /*! @brief middle point locations of elements */
    std::vector<std::vector<Vector>> _interfaceMidPoints; /*! @brief middle point locations of edges */
62

63
64
    // Solution related members
    VectorVector _sol;                 /*! @brief solution of the PDE, e.g. angular flux or moments */
65
    std::vector<double> _solverOutput; /*! @brief LEGACY: Outputfield for solver ==> Will be replaced by _outputFields in the near future */
66

67
    // Output related members
68
    std::vector<std::vector<std::vector<double>>> _outputFields; /*! @brief: Solver Output: dimensions (GroupID,FieldID,CellID).*/
69
    std::vector<std::vector<std::string>> _outputFieldNames;     /*! @brief: Names of the outputFields: dimensions (GroupID,FieldID) */
70
    // we will have to add a further dimension for quadPoints and weights once we start with multilevel SN
71

72
73
74
75
    // Output related members
    std::vector<double> _screenOutputFields;          /*! @brief: Solver Output: dimensions (FieldID). */
    std::vector<std::string> _screenOutputFieldNames; /*! @brief: Names of the outputFields: dimensions (FieldID) */

76
77
78
79
    // Output related members
    std::vector<double> _historyOutputFields;          /*! @brief: Solver Output: dimensions (FieldID). */
    std::vector<std::string> _historyOutputFieldNames; /*! @brief: Names of the outputFields: dimensions (FieldID) */

80
81
82
83
84
    // Internal Members
    VectorVector _solNew; /*! @brief: VectorVector to store the new flux and later the new solution per iteration */    // REPLACES psiNEW
    Vector _fluxNew; /*! @brief: Vector to store the new Flux */
    Vector _flux;    /*! @brief: Vector to store the old Flux */

85
86
    // ---- Member functions ----

87
    // Solver
88
89
    /*! @brief Performs preprocessing for the current solver iteration */
    virtual void IterPreprocessing() = 0;
90
91
    /*! @brief Performs postprocessing for the current solver iteration */
    virtual void IterPostprocessing() = 0;
92
    /*! @brief Constructs  the flux update for the current iteration and stores it in psiNew*/
93
    virtual void FluxUpdate() = 0;
94
    /*! @brief Computes the finite Volume update step for the current iteration */
95
    virtual void FVMUpdate( unsigned idx_energy ) = 0;
96

97
    // Helper
98
    /*! @brief ComputeTimeStep calculates the maximal stable time step */
99
    double ComputeTimeStep( double cfl ) const;
100
101
102
    /*! @brief: Computes the flux of the solution to check conservation properties */
    virtual void ComputeRadFlux() = 0;

103
104
105
106
107
    // IO
    /*! @brief Initializes the output groups and fields of this solver and names the fields */
    virtual void PrepareVolumeOutput() = 0;
    /*! @brief Function that prepares VTK export and csv export of the current solver iteration  */
    virtual void WriteVolumeOutput( unsigned iteration ) = 0;
Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
108
109
    /*! @brief Save Output solution at given energy (pseudo time) to VTK file. Write frequency is given by
               option VOLUME_OUTPUT_FREQUENCY. Always prints last iteration without iteration affix.*/
110
    void PrintVolumeOutput( int currEnergy ) const;
Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
111
    /*! @brief: Initialized the output fields and their Names for the screenoutput */
112
    void PrepareScreenOutput();
113
114
    /*! @brief Function that writes screen and history output fields */
    void WriteScalarOutput( unsigned iteration );
Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
115
116
    /*! @brief Prints ScreenOutputFields to Screen and to logger. Write frequency is given by
               option SCREEN_OUTPUT_FREQUENCY. Always prints last iteration. */
117
    void PrintScreenOutput( unsigned iteration );
Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
118
119
    /*! @brief: Initialized the historyOutputFields and their Names for history output. Write frequency is given by
               option HISTORY_OUTPUT_FREQUENCY. Always prints last iteration. */
120
    void PrepareHistoryOutput();
121
122
    /*! @brief Prints HistoryOutputFields to logger */
    void PrintHistoryOutput( unsigned iteration );
123

124
  public:
125
126
    /*! @brief Solver constructor
     *  @param settings stores all needed information */
127
    Solver( Config* settings );
Jonas Kusch's avatar
Jonas Kusch committed
128

129
130
    ~Solver();

131
132
133
    /*! @brief Create constructor
     *  @param settings stores all needed information
     *  @return pointer to Solver */
134
    static Solver* Create( Config* settings );
Jonas Kusch's avatar
Jonas Kusch committed
135

136
    /*! @brief Solve functions runs main time loop */
137
    virtual void Solve();
138

139
    /*! @brief Save Output solution to VTK file */
Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
140
    void PrintVolumeOutput() const;    // Only for debugging purposes.
Jonas Kusch's avatar
Jonas Kusch committed
141
};
142
#endif    // SOLVER_H