Commit f4423675 authored by Steffen Schotthöfer's avatar Steffen Schotthöfer
Browse files

add documentation for mn and pn solver

parent da2db94b
......@@ -125,7 +125,7 @@ class Config
bool _dataGeneratorMode;
unsigned long _tainingSetSize; /*!< @brief Size of training data set for data generator */
unsigned long _maxValFirstMoment; /*!< @brief Size of training data set for data generator */
double _RealizableSetEpsilonU0; /*! @brief Distance to 0 of the sampled moments to the boundary of the realizable set */
double _RealizableSetEpsilonU0; /*!< @brief Distance to 0 of the sampled moments to the boundary of the realizable set */
double _RealizableSetEpsilonU1; /*!< @brief: norm(u_1)/u_0 !< _RealizableSetEpsilonU1 */
bool _normalizedSampling; /*!< @brief: Flag for sampling of normalized moments, i.e. u_0 =1 */
......@@ -330,10 +330,12 @@ class Config
// ---- Setters for option structure
// This section is dangerous
// Quadrature Structure
void SetNQuadPoints( unsigned nq ) { _nQuadPoints = nq; } /*! @brief Never change the nq! This is only for the test framework. */
void SetQuadName( QUAD_NAME quadName ) { _quadName = quadName; } /*! @brief Never change the quadName! This is only for the test framework. */
void SetQuadOrder( unsigned quadOrder ) { _quadOrder = quadOrder; } /*! @brief Never change the quadOrder! This is only for the test framework. */
void SetSNAllGaussPts( bool useall ) { _allGaussPts = useall; } /*! @brief Never change the this! This is only for the test framework. */
void SetNQuadPoints( unsigned nq ) { _nQuadPoints = nq; } /*!< @brief Never change the nq! This is only for the test framework. */
void SetQuadName( QUAD_NAME quadName ) { _quadName = quadName; } /*!< @brief Never change the quadName! This is only for the test framework. */
void SetQuadOrder( unsigned quadOrder ) {
_quadOrder = quadOrder;
} /*!< @brief Never change the quadOrder! This is only for the test framework. */
void SetSNAllGaussPts( bool useall ) { _allGaussPts = useall; } /*!< @brief Never change the this! This is only for the test framework. */
// Mesh Structure
void SetNCells( unsigned nCells ) { _nCells = nCells; }
};
......
......@@ -87,11 +87,11 @@ class ProblemBase
*/
virtual VectorVector SetupIC() = 0;
/*!< @brief: Exact analytical solution for the Line Source Test Case at
@param: x: x coordinate of exact solution
y: y coordinate of exact solution
t: time of the exact solution
scatteringXS: scattering cross section of the exact solution
/*! @brief: Exact analytical solution for the Line Source Test Case at
@param: x: x coordinate of exact solution
y: y coordinate of exact solution
t: time of the exact solution
scatteringXS: scattering cross section of the exact solution
@return: exact solution at x,y,t,scatteringXS
*/ // Default is set to 0. ~> if no analytical solution is available.
double virtual GetAnalyticalSolution( double /*x*/, double /*y*/, double /*t*/, double /*scatteringXS*/ ) { return 0.0; }
......
......@@ -12,7 +12,7 @@ class MNSolver : public SolverBase
public:
/**
* @brief MNSolver constructor
* @param settings stores all needed information
* @param settings : Config class that stores all needed information
*/
MNSolver( Config* settings );
......@@ -21,15 +21,15 @@ class MNSolver : public SolverBase
private:
// --- Private member variables ---
unsigned _nTotalEntries; /*! @brief: Total number of equations in the system */
unsigned short _LMaxDegree; /*! @brief: Max Order of Spherical Harmonics */
unsigned _nSystem; /*! @brief Total number of equations in the system */
unsigned short _polyDegreeBasis; /*! @brief Max polynomial degree of the basis */
// Moment basis
SphericalBase* _basis; /*! @brief: Class to compute and store current spherical harmonics basis */
VectorVector _moments; /*! @brief: Moment Vector pre-computed at each quadrature point: dim= _nq x _nTotalEntries */
SphericalBase* _basis; /*! @brief Class to compute and store current spherical harmonics basis */
VectorVector _moments; /*! @brief Moment Vector pre-computed at each quadrature point: dim= _nq x _nTotalEntries */
// Scattering
Vector _scatterMatDiag; /*! @brief: Diagonal of the scattering matrix (its a diagonal matrix by construction) */
Vector _scatterMatDiag; /*! @brief Diagonal of the scattering matrix (its a diagonal matrix by construction) */
// Quadrature related members
VectorVector _quadPoints; /*! @brief quadrature points, dim(_quadPoints) = (_nq,spatialDim) */
......@@ -37,41 +37,41 @@ class MNSolver : public SolverBase
VectorVector _quadPointsSphere; /*! @brief (my,phi), dim(_quadPoints) = (_nq,2) */
// Entropy Optimization related members
EntropyBase* _entropy; /*! @brief: Class to handle entropy functionals */
VectorVector _alpha; /*! @brief: Lagrange Multipliers for Minimal Entropy problem for each gridCell
Layout: _nCells x _nTotalEntries*/
OptimizerBase* _optimizer; /*! @brief: Class to solve minimal entropy problem */
EntropyBase* _entropy; /*! @brief Class to handle entropy functionals */
VectorVector _alpha; /*! @brief Lagrange Multipliers for Minimal Entropy problem for each gridCell
Layout: _nCells x _nTotalEntries*/
OptimizerBase* _optimizer; /*! @brief Class to solve minimal entropy problem */
VectorVector _solDx; /*! @brief: temporary storage of x-derivatives of alpha */
VectorVector _solDy; /*! @brief: temporary storage of y-derivatives of alpha */
VectorVector _solDx; /*! @brief temporary storage of x-derivatives of alpha */
VectorVector _solDy; /*! @brief temporary storage of y-derivatives of alpha */
// ---- Private Member functions ---
// IO
void PrepareVolumeOutput() override;
void WriteVolumeOutput( unsigned idx_pseudoTime ) override;
void WriteVolumeOutput( unsigned idx_iter ) override;
// Solver
void FVMUpdate( unsigned idx_energy ) override;
void FVMUpdate( unsigned idx_iter ) override;
void FluxUpdate() override;
void IterPreprocessing( unsigned idx_pseudotime ) override;
void IterPostprocessing( unsigned idx_pseudotime );
/*! @brief : Construct flux by computing the Moment of the sum of FVM discretization at the interface of cell
* @param : idx_cell = current cell id
* @returns : sum over all neighbors of flux for all moments at interface of idx_cell, idx_neighbor */
void IterPreprocessing( unsigned /*idx_iter*/ ) override;
void IterPostprocessing( unsigned /*idx_iter*/ ) override;
/*! @brief Construct flux by computing the Moment of the sum of FVM discretization at the interface of cell
* @param idx_cell current cell id
* @returns sum over all neighbors of flux for all moments at interface of idx_cell, idx_neighbor */
Vector ConstructFlux( unsigned idx_cell );
/*! @brief Corrects the solution _sol[idx_cell] to be realizable w.r.t. the reconstructed entropy (eta'(alpha*m))
@param idx_cell = cell where the correction happens*/
@param idx_cell cell where the correction happens*/
void ComputeRealizableSolution( unsigned idx_cell );
// Initialization of the solver
/*! @brief : Pre-Compute Moments at all quadrature points. */
/*! @brief Pre-Compute Moments at all quadrature points. */
void ComputeMoments();
/*! @brief: fucntion for computing and setting up EV matrix for scattering kernel */
/*! @brief Function for computing and setting up EV matrix for scattering kernel */
void ComputeScatterMatrix();
// Helper
/*! @brief: Computes the radiative flux from the solution vector of the moment system */
/*! @brief Computes the radiative flux from the solution vector of the moment system */
void ComputeRadFlux();
};
#endif // MNSOLVER_H
......@@ -15,93 +15,132 @@ class PNSolver : public SolverBase
virtual ~PNSolver() {}
private:
unsigned _nTotalEntries; /*! @brief: total number of equations in the system */
unsigned _LMaxDegree; /*! @brief: maximal degree of the spherical harmonics basis*/
unsigned _nSystem; /*!< @brief total number of equations in the system */
unsigned _polyDegreeBasis; /*!< @brief maximal degree of the spherical harmonics basis*/
// System Matrix for x, y and z flux
// ==> not needed after computation of A+ and A- ==> maybe safe only temporarly and remove as member?
SymMatrix _Ax; /*! @brief: Flux Jacbioan in x direction */
SymMatrix _Ay; /*! @brief: Flux Jacbioan in x direction */
SymMatrix _Az; /*! @brief: Flux Jacbioan in x direction */
SymMatrix _Ax; /*!< @brief Flux Jacbioan in x direction */
SymMatrix _Ay; /*!< @brief Flux Jacbioan in x direction */
SymMatrix _Az; /*!< @brief Flux Jacbioan in x direction */
// Upwinding Matrices
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 */
Vector _scatterMatDiag; /*! @brief: diagonal of the scattering matrix (its a diagonal matrix by construction). Contains eigenvalues of the
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 */
Vector _scatterMatDiag; /*!< @brief diagonal of the scattering matrix (its a diagonal matrix by construction). Contains eigenvalues of the
scattering kernel. */
VectorVector _solDx; /*! @brief: temporary storage of x-derivatives of solution */
VectorVector _solDy; /*! @brief: temporary storage of y-derivatives of solution */
VectorVector _solDx; /*!< @brief Temporary storage of x-derivatives of solution */
VectorVector _solDy; /*!< @brief Temporary storage of y-derivatives of solution */
// ---- Member functions ----
// IO
/*! @brief Initializes the output groups and fields of this solver and names the fields */
void PrepareVolumeOutput() override;
/*! @brief Function that prepares VTK export and csv export of the current solver iteration
@returns: Mass of current iteration */
void WriteVolumeOutput( unsigned idx_pseudoTime ) override;
// Solver
void FVMUpdate( unsigned idx_energy ) override;
void FluxUpdate() override;
void IterPreprocessing( unsigned idx_pseudotime ) override;
void IterPostprocessing( unsigned idx_pseudotime );
void IterPostprocessing( unsigned idx_pseudotime ) override;
// Helper
void ComputeRadFlux();
// Initialization of the Solver
/*! @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
/*! @brief parameter functions for setting up system matrix
* @param l degree , it must hold: 0 <= l <=_nq
* @param k order , it must hold: -l <=k <= l
*/
double AParam( int l, int k ) const;
/*! @brief parameter functions for setting up system matrix
* @param l degree , it must hold: 0 <= l <=_nq
* @param k order , it must hold: -l <=k <= l
*/
double BParam( int l, int k ) const;
/*! @brief parameter functions for setting up system matrix
* @param l degree , it must hold: 0 <= l <=_nq
* @param k order , it must hold: -l <=k <= l
*/
double CParam( int l, int k ) const;
/*! @brief parameter functions for setting up system matrix
* @param l degree , it must hold: 0 <= l <=_nq
* @param k order , it must hold: -l <=k <= l
*/
double DParam( int l, int k ) const;
/*! @brief parameter functions for setting up system matrix
* @param l degree , it must hold: 0 <= l <=_nq
* @param k order , it must hold: -l <=k <= l
*/
double EParam( int l, int k ) const;
/*! @brief parameter functions for setting up system matrix
* @param l degree , it must hold: 0 <= l <=_nq
* @param k order , it must hold: -l <=k <= l
*/
double FParam( int l, int k ) const;
/*! @brief parameter functions for setting up system matrix
* @param l degree , it must hold: 0 <= l <=_nq
* @param k order , it must hold: -l <=k <= l
*/
double CTilde( int l, int k ) const;
/*! @brief parameter functions for setting up system matrix
* @param l degree , it must hold: 0 <= l <=_nq
* @param k order , it must hold: -l <=k <= l
*/
double DTilde( int l, int k ) const;
/*! @brief parameter functions for setting up system matrix
* @param l degree , it must hold: 0 <= l <=_nq
* @param k order , it must hold: -l <=k <= l
*/
double ETilde( int l, int k ) const;
/*! @brief parameter functions for setting up system matrix
* @param l degree , it must hold: 0 <= l <=_nq
* @param k order , it must hold: -l <=k <= l
*/
double FTilde( int l, int k ) const;
/*! @brief: mathematical + index functions. Helper functions for setting up system matrix.
* @param: k: arbitrary integer
*/
int Sgn( int k ) const;
/*! @brief mathematical + index functions. Helper functions for setting up system matrix.
* @param k arbitrary integer
*/
int kPlus( int k ) const;
/*! @brief mathematical + index functions. Helper functions for setting up system matrix.
* @param k arbitrary integer
*/
int kMinus( int k ) const;
/*! @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
/*! @brief computes the global index of the moment corresponding to basis function (l,k)
* @param l degree l, it must hold: 0 <= l <=_nq
* @param k order k, it must hold: -l <=k <= l
* @returns global index
*/
int GlobalIndex( int l, int k ) const;
/*! @brief : Checks, if index invariant for global index holds
* @returns : True, if invariant holds, else false
/*! @brief Checks, if index invariant for global index holds
* @returns True, if invariant holds, else false
*/
bool CheckIndex( int l, int k ) const;
/*! @brief: function for computing and setting up system matrices */
/*! @brief Function for computing and setting up system matrices */
void ComputeSystemMatrices();
/*! @brief: function for computing and setting up flux matrices for upwinding */
/*! @brief Function for computing and setting up flux matrices for upwinding */
void ComputeFluxComponents();
/*! @brief: fucntion for computing and setting up EV matrix for scattering kernel */
/*! @brief Function for computing and setting up EV matrix for scattering kernel */
void ComputeScatterMatrix();
/*! @brief: Computes Legedre polinomial of oder l at point x */
/*! @brief Computes Legedre polinomial of oder l at point x */
double LegendrePoly( double x, int l );
/*! @brief: Sets Entries of FluxMatrices to zero, if they are below double precision,
......
......@@ -69,17 +69,17 @@ class SolverBase
std::vector<double> _solverOutput; /*!< @brief LEGACY: Outputfield for solver ==> Will be replaced by _outputFields in the near future */
// Output related members
std::vector<std::vector<std::vector<double>>> _outputFields; /*!< @brief: Solver Output: dimensions (GroupID,FieldID,CellID).*/
std::vector<std::vector<std::string>> _outputFieldNames; /*!< @brief: Names of the outputFields: dimensions (GroupID,FieldID) */
std::vector<std::vector<std::vector<double>>> _outputFields; /*!< @brief Solver Output: dimensions (GroupID,FieldID,CellID).*/
std::vector<std::vector<std::string>> _outputFieldNames; /*!< @brief Names of the outputFields: dimensions (GroupID,FieldID) */
// we will have to add a further dimension for quadPoints and weights once we start with multilevel SN
// Output related members
std::vector<double> _screenOutputFields; /*!< @brief: Solver Output: dimensions (FieldID). */
std::vector<std::string> _screenOutputFieldNames; /*!< @brief: Names of the outputFields: dimensions (FieldID) */
std::vector<double> _screenOutputFields; /*!< @brief Solver Output: dimensions (FieldID). */
std::vector<std::string> _screenOutputFieldNames; /*!< @brief Names of the outputFields: dimensions (FieldID) */
// Output related members
std::vector<double> _historyOutputFields; /*!< @brief: Solver Output: dimensions (FieldID). */
std::vector<std::string> _historyOutputFieldNames; /*!< @brief: Names of the outputFields: dimensions (FieldID) */
std::vector<double> _historyOutputFields; /*!< @brief Solver Output: dimensions (FieldID). */
std::vector<std::string> _historyOutputFieldNames; /*!< @brief Names of the outputFields: dimensions (FieldID) */
// ---- Member functions ----
......@@ -87,19 +87,19 @@ class SolverBase
/*! @brief Performs preprocessing steps before the pseudo time iteration is started*/
virtual void SolverPreprocessing();
/*! @brief Performs preprocessing for the current solver iteration
@param idx_iter : current (peudo) time iteration */
@param idx_iter current (peudo) time iteration */
virtual void IterPreprocessing( unsigned idx_iter ) = 0;
/*! @brief Performs postprocessing for the current solver iteration */
virtual void IterPostprocessing( unsigned idx_pseudotime ) = 0;
/*! @brief Constructs the flux update for the current iteration and stores it in psiNew*/
virtual void FluxUpdate() = 0;
/*! @brief Computes the finite Volume update step for the current iteration
@param idx_iter : current (peudo) time iteration */
@param idx_iter current (peudo) time iteration */
virtual void FVMUpdate( unsigned idx_iter ) = 0;
// Helper
/*! @brief ComputeTimeStep calculates the maximal stable time step using the cfl number
@param used cfl number */
@param cfl Courant-Friedrichs-Levy condition number */
double ComputeTimeStep( double cfl ) const;
/*! @brief: Computes the flux of the solution to check conservation properties */
virtual void ComputeRadFlux() = 0;
......@@ -108,27 +108,27 @@ class SolverBase
/*! @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
@param idx_iter : current (pseudo) time iteration */
@param idx_iter current (pseudo) time iteration */
virtual void WriteVolumeOutput( unsigned idx_iter ) = 0;
/*! @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.
@param idx_iter : current (pseudo) time iteration */
@param idx_iter current (pseudo) time iteration */
void PrintVolumeOutput( int idx_iter ) const;
/*! @brief: Initialized the output fields and their Names for the screenoutput */
void PrepareScreenOutput();
/*! @brief Function that writes screen and history output fields
@param idx_iter : current (pseudo) time iteration */
@param idx_iter current (pseudo) time iteration */
void WriteScalarOutput( unsigned idx_iter );
/*! @brief Prints ScreenOutputFields to Screen and to logger. Write frequency is given by
option SCREEN_OUTPUT_FREQUENCY. Always prints last iteration.
@param idx_iter : current (pseudo) time iteration */
@param idx_iter current (pseudo) time iteration */
void PrintScreenOutput( unsigned idx_iter );
/*! @brief: Initialized the historyOutputFields and their Names for history output. Write frequency is given by
/*! @brief Initialized the historyOutputFields and their Names for history output. Write frequency is given by
option HISTORY_OUTPUT_FREQUENCY. Always prints last iteration. */
void PrepareHistoryOutput();
/*! @brief Prints HistoryOutputFields to logger
@param idx_iter : current (pseudo) time iteration */
@param idx_iter current (pseudo) time iteration */
void PrintHistoryOutput( unsigned idx_iter );
/*! @brief Pre Solver Screen and Logger Output */
void DrawPreSolverOutput();
......@@ -137,13 +137,13 @@ class SolverBase
public:
/*! @brief Solver constructor
* @param settings :config class that stores all needed config information */
* @param settings config class that stores all needed config information */
SolverBase( Config* settings );
virtual ~SolverBase();
/*! @brief Create constructor
* @param settings :config class that stores all needed config information
* @param settings config class that stores all needed config information
* @return pointer to SolverBase */
static SolverBase* Create( Config* settings );
......
......@@ -19,9 +19,9 @@
MNSolver::MNSolver( Config* settings ) : SolverBase( settings ) {
_LMaxDegree = settings->GetMaxMomentDegree();
_basis = SphericalBase::Create( _settings );
_nTotalEntries = _basis->GetBasisSize();
_polyDegreeBasis = settings->GetMaxMomentDegree();
_basis = SphericalBase::Create( _settings );
_nSystem = _basis->GetBasisSize();
// build quadrature object and store quadrature points and weights
_quadPoints = _quadrature->GetPoints();
......@@ -31,11 +31,11 @@ MNSolver::MNSolver( Config* settings ) : SolverBase( settings ) {
//_settings->SetNQuadPoints( _nq );
// Initialize temporary storages of alpha derivatives
_solDx = VectorVector( _nCells, Vector( _nTotalEntries, 0.0 ) );
_solDy = VectorVector( _nCells, Vector( _nTotalEntries, 0.0 ) );
_solDx = VectorVector( _nCells, Vector( _nSystem, 0.0 ) );
_solDy = VectorVector( _nCells, Vector( _nSystem, 0.0 ) );
// Initialize Scatter Matrix --
_scatterMatDiag = Vector( _nTotalEntries, 0.0 );
_scatterMatDiag = Vector( _nSystem, 0.0 );
ComputeScatterMatrix();
// Initialize Entropy
......@@ -45,10 +45,10 @@ MNSolver::MNSolver( Config* settings ) : SolverBase( settings ) {
_optimizer = OptimizerBase::Create( _settings );
// Initialize lagrange Multiplier
_alpha = VectorVector( _nCells, Vector( _nTotalEntries, 0.0 ) );
_alpha = VectorVector( _nCells, Vector( _nSystem, 0.0 ) );
// Initialize and Pre-Compute Moments at quadrature points
_moments = VectorVector( _nq, Vector( _nTotalEntries, 0.0 ) );
_moments = VectorVector( _nq, Vector( _nSystem, 0.0 ) );
ComputeMoments();
}
......@@ -63,7 +63,7 @@ void MNSolver::ComputeScatterMatrix() {
// --- Isotropic ---
_scatterMatDiag[0] = -1.0;
for( unsigned idx_diag = 1; idx_diag < _nTotalEntries; idx_diag++ ) {
for( unsigned idx_diag = 1; idx_diag < _nSystem; idx_diag++ ) {
_scatterMatDiag[idx_diag] = 0.0;
}
}
......@@ -84,11 +84,11 @@ Vector MNSolver::ConstructFlux( unsigned idx_cell ) {
//--- Integration of moments of flux ---
double entropyL, entropyR, entropyFlux;
Vector flux( _nTotalEntries, 0.0 );
Vector flux( _nSystem, 0.0 );
//--- Temporary storages of reconstructed alpha ---
Vector alphaL( _nTotalEntries, 0.0 );
Vector alphaR( _nTotalEntries, 0.0 );
Vector alphaL( _nSystem, 0.0 );
Vector alphaR( _nSystem, 0.0 );
for( unsigned idx_quad = 0; idx_quad < _nq; idx_quad++ ) {
entropyFlux = 0.0; // reset temorary flux
......@@ -152,7 +152,7 @@ void MNSolver::IterPreprocessing( unsigned /*idx_pseudotime*/ ) {
}
}
void MNSolver::IterPostprocessing( unsigned idx_pseudotime ) {
void MNSolver::IterPostprocessing( unsigned /*idx_iter*/ ) {
// --- Update Solution ---
_sol = _solNew;
......@@ -170,7 +170,7 @@ void MNSolver::ComputeRadFlux() {
void MNSolver::FluxUpdate() {
if( _reconsOrder > 1 ) {
_mesh->ReconstructSlopesU( _nTotalEntries, _solDx, _solDy, _alpha ); // unstructured reconstruction
_mesh->ReconstructSlopesU( _nSystem, _solDx, _solDy, _alpha ); // unstructured reconstruction
}
// Loop over the grid cells
......@@ -182,19 +182,19 @@ void MNSolver::FluxUpdate() {
}
}
void MNSolver::FVMUpdate( unsigned idx_energy ) {
void MNSolver::FVMUpdate( unsigned idx_iter ) {
// Loop over the grid cells
//#pragma omp parallel for
for( unsigned idx_cell = 0; idx_cell < _nCells; idx_cell++ ) {
// Dirichlet Boundaries stay
if( _boundaryCells[idx_cell] == BOUNDARY_TYPE::DIRICHLET ) continue;
// Flux update
for( unsigned idx_system = 0; idx_system < _nTotalEntries; idx_system++ ) {
for( unsigned idx_system = 0; idx_system < _nSystem; idx_system++ ) {
_solNew[idx_cell][idx_system] = _sol[idx_cell][idx_system] -
( _dE / _areas[idx_cell] ) * _solNew[idx_cell][idx_system] /* cell averaged flux */
- _dE * _sol[idx_cell][idx_system] *
( _sigmaT[idx_energy][idx_cell] /* absorbtion influence */
+ _sigmaS[idx_energy][idx_cell] * _scatterMatDiag[idx_system] ); /* scattering influence */
( _sigmaT[idx_iter][idx_cell] /* absorbtion influence */
+ _sigmaS[idx_iter][idx_cell] * _scatterMatDiag[idx_system] ); /* scattering influence */
}
// Source Term
_solNew[idx_cell][0] += _dE * _Q[0][idx_cell][0];
......@@ -224,11 +224,11 @@ void MNSolver::PrepareVolumeOutput() {
case MOMENTS:
// As many entries as there are moments in the system
_outputFields[idx_group].resize( _nTotalEntries );
_outputFieldNames[idx_group].resize( _nTotalEntries );
_outputFields[idx_group].resize( _nSystem );
_outputFieldNames[idx_group].resize( _nSystem );
if( _settings->GetSphericalBasisName() == SPHERICAL_HARMONICS ) {
for( int idx_l = 0; idx_l <= (int)_LMaxDegree; idx_l++ ) {
for( int idx_l = 0; idx_l <= (int)_polyDegreeBasis; idx_l++ ) {
for( int idx_k = -idx_l; idx_k <= idx_l; idx_k++ ) {
_outputFields[idx_group][_basis->GetGlobalIndexBasis( idx_l, idx_k )].resize( _nCells );
_outputFieldNames[idx_group][_basis->GetGlobalIndexBasis( idx_l, idx_k )] =
......@@ -237,7 +237,7 @@ void MNSolver::PrepareVolumeOutput() {
}
}
else {
for( unsigned idx_l = 0; idx_l <= _LMaxDegree; idx_l++ ) {
for( unsigned idx_l = 0; idx_l <= _polyDegreeBasis; idx_l++ ) {
unsigned maxOrder_k = _basis->GetCurrDegreeSize( idx_l );
for( unsigned idx_k = 0; idx_k < maxOrder_k; idx_k++ ) {
_outputFields[idx_group][_basis->GetGlobalIndexBasis( idx_l, idx_k )].resize( _nCells );
......@@ -250,11 +250,11 @@ void MNSolver::PrepareVolumeOutput() {
case DUAL_MOMENTS:
// As many entries as there are moments in the system
_outputFields[idx_group].resize( _nTotalEntries );
_outputFieldNames[idx_group].resize( _nTotalEntries );
_outputFields[idx_group].resize( _nSystem );
_outputFieldNames[idx_group].resize( _nSystem );
if( _settings->GetSphericalBasisName() == SPHERICAL_HARMONICS ) {
for( int idx_l = 0; idx_l <= (int)_LMaxDegree; idx_l++ ) {
for( int idx_l = 0; idx_l <= (int)_polyDegreeBasis; idx_l++ ) {
for( int idx_k = -idx_l; idx_k <= idx_l; idx_k++ ) {
_outputFields[idx_group][_basis->GetGlobalIndexBasis( idx_l, idx_k )].resize( _nCells );
_outputFieldNames[idx_group][_basis->GetGlobalIndexBasis( idx_l, idx_k )] =
......@@ -263,7 +263,7 @@ void MNSolver::PrepareVolumeOutput() {
}
}
else {
for( int idx_l = 0; idx_l <= (int)_LMaxDegree; idx_l++ ) {
for( int idx_l = 0; idx_l <= (int)_polyDegreeBasis; idx_l++ ) {
unsigned maxOrder_k = _basis->GetCurrDegreeSize( idx_l );
for( unsigned idx_k = 0; idx_k < maxOrder_k; idx_k++ ) {
_outputFields[idx_group][_basis->GetGlobalIndexBasis( idx_l, idx_k )].resize( _nCells );
......@@ -287,12 +287,12 @@ void MNSolver::PrepareVolumeOutput() {
}
}
void MNSolver::WriteVolumeOutput( unsigned idx_pseudoTime ) {
void MNSolver::WriteVolumeOutput( unsigned idx_iter ) {
unsigned nGroups = (unsigned)_settings->GetNVolumeOutput();
// Check if volume output fields are written to file this iteration
if( ( _settings->GetVolumeOutputFrequency() != 0 && idx_pseudoTime % (unsigned)_settings->GetVolumeOutputFrequency() == 0 ) ||
( idx_pseudoTime == _nEnergies - 1 ) /* need sol at last iteration */ ) {
if( ( _settings->GetVolumeOutputFrequency() != 0 && idx_iter % (unsigned)_settings->GetVolumeOutputFrequency() == 0 ) ||
( idx_iter == _nEnergies - 1 ) /* need sol at last iteration */ ) {
for( unsigned idx_group = 0; idx_group < nGroups; idx_group++ ) {
switch( _settings->GetVolumeOutput()[idx_group] ) {
......@@ -302,14 +302,14 @@ void MNSolver::WriteVolumeOutput( unsigned idx_pseudoTime ) {
}
break;
case MOMENTS:
for( unsigned idx_sys = 0; idx_sys < _nTotalEntries; idx_sys++ ) {
for( unsigned idx_sys = 0; idx_sys < _nSystem; idx_sys++ ) {
for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
_outputFields[idx_group][idx_sys][idx_cell] = _sol[idx_cell][idx_sys];
}
}
break;
case DUAL_MOMENTS:
for( unsigned idx_sys = 0; idx_sys < _nTotalEntries; idx_sys++ ) {
for( unsigned idx_sys = 0; idx_sys < _nSystem; idx_sys++ ) {
for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
_outputFields[idx_group][idx_sys][idx_cell] = _alpha[idx_cell][idx_sys];
}
......@@ -319,10 +319,10 @@ void MNSolver::WriteVolumeOutput( unsigned idx_pseudoTime ) {
// Compute total "mass" of the system ==> to check conservation properties
for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
double time = idx_pseudoTime * _dE;
double time = idx_iter * _dE;
_outputFields[idx_group][0][idx_cell] = _problem->GetAnalyticalSolution(
_mesh->GetCellMidPoints()[idx_cell][0], _mesh->GetCellMidPoints()[idx_cell][1], time, _sigmaS[idx_pseudoTime][idx_cell] );
_mesh->GetCellMidPoints()[idx_cell][0], _mesh->GetCellMidPoints()[idx_cell][1], time, _sigmaS[idx_iter][idx_cell] );
}
break;
......
......@@ -11,30 +11,30 @@
#include <mpi.h>
PNSolver::PNSolver( Config* settings ) : SolverBase( settings ) {
_LMaxDegree = settings->GetMaxMomentDegree();
_nTotalEntries = GlobalIndex( int( _LMaxDegree ), int( _LMaxDegree ) ) + 1;
_polyDegreeBasis = settings->GetMaxMomentDegree();
_nSystem = GlobalIndex( int( _polyDegreeBasis ), int( _polyDegreeBasis ) ) + 1;
// Initialize System Matrices
_Ax = SymMatrix( _nTotalEntries );