Commit 4a4101d1 authored by steffen.schotthoefer's avatar steffen.schotthoefer
Browse files

adapted mn and pn solver to scatterXS and totalXS notation as snsolver


Former-commit-id: 66aa4556
parent aef02006
......@@ -27,9 +27,6 @@ class MNSolver : public Solver
unsigned _nTotalEntries; /*! @brief: Total number of equations in the system */
unsigned short _LMaxDegree; /*! @brief: Max Order of Moments */
// Solver specific physics
VectorVector _sigmaA; /*! @brief: Absorption coefficient for all energies*/
// Moment basis
SphericalHarmonics* _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 */
......@@ -82,6 +79,9 @@ class MNSolver : public Solver
/*! @brief : Pre-Compute Moments at all quadrature points. */
void ComputeMoments();
/*! @brief: fucntion for computing and setting up EV matrix for scattering kernel */
void ComputeScatterMatrix();
/*! @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*/
void ComputeRealizableSolution( unsigned idx_cell );
......
......@@ -19,7 +19,7 @@ class PNSolver : public Solver
unsigned _nTotalEntries; /*! @brief: total number of equations in the system */
unsigned _LMaxDegree; /*! @brief: maximal degree of the spherical harmonics basis*/
VectorVector _sigmaA; /*! @brief: Absorption coefficient for all energies*/
// VectorVector _sigmaA; /*! @brief: Absorption coefficient for all energies*/
// 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 */
......@@ -39,7 +39,8 @@ class PNSolver : public Solver
// 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) */
Vector _scatterMatDiag; /*! @brief: diagonal of the scattering matrix (its a diagonal matrix by construction). Contains eigenvalues of the
scattering kernel. */
// Output related members
std::vector<std::vector<std::vector<double>>> _outputFields; /*! @brief: Solver Output: dimensions (GroupID,FieldID,CellID). !Protoype output for
......@@ -97,7 +98,7 @@ class PNSolver : public Solver
void ComputeSystemMatrices();
/*! @brief: function for computing and setting up flux matrices for upwinding */
void ComputeFluxComponents();
/*! @brief: fucntion for computing and setting up diagonal matrix for scattering kernel */
/*! @brief: fucntion for computing and setting up EV matrix for scattering kernel */
void ComputeScatterMatrix();
/*! @brief: Computes Legedre polinomial of oder l at point x */
double LegendrePoly( double x, int l );
......
......@@ -31,19 +31,17 @@ MNSolver::MNSolver( Config* settings ) : Solver( settings ) {
_quadPointsSphere = _quadrature->GetPointsSphere();
_settings->SetNQuadPoints( _nq );
// transform sigmaT and sigmaS in sigmaA.
_sigmaA = VectorVector( _nEnergies, Vector( _nCells, 0 ) ); // Get rid of this extra vektor!
// This must be shifted to Problem Class
for( unsigned n = 0; n < _nEnergies; n++ ) {
for( unsigned j = 0; j < _nCells; j++ ) {
_sigmaA[n][j] = 0; //_sigmaT[n][j] - _sigmaS[n][j];
_sigmaS[n][j] = 0;
_sigmaT[n][j] = 1; //_sigmaT[n][j] - _sigmaS[n][j];
_sigmaS[n][j] = 1;
}
}
// Initialize Scatter Matrix
_scatterMatDiag = Vector( _nTotalEntries, 1.0 );
_scatterMatDiag[0] = 0.0; // First entry is zero by construction.
// Initialize Scatter Matrix --
_scatterMatDiag = Vector( _nTotalEntries, 0.0 );
ComputeScatterMatrix();
// Initialize Entropy
_entropy = EntropyBase::Create( _settings );
......@@ -70,6 +68,15 @@ MNSolver::~MNSolver() {
delete _basis;
}
void MNSolver::ComputeScatterMatrix() {
// --- Isotropic ---
_scatterMatDiag[0] = -1.0;
for( unsigned idx_diag = 1; idx_diag < _nTotalEntries; idx_diag++ ) {
_scatterMatDiag[idx_diag] = 0.0;
}
}
int MNSolver::GlobalIndex( int l, int k ) const {
int numIndicesPrevLevel = l * l; // number of previous indices untill level l-1
int prevIndicesThisLevel = k + l; // number of previous indices in current level
......@@ -173,7 +180,7 @@ void MNSolver::Solve() {
psiNew[idx_cell][idx_system] = _sol[idx_cell][idx_system] -
( _dE / _areas[idx_cell] ) * psiNew[idx_cell][idx_system] /* cell averaged flux */
- _dE * _sol[idx_cell][idx_system] *
( _sigmaA[idx_energy][idx_cell] /* absorbtion influence */
( _sigmaT[idx_energy][idx_cell] /* absorbtion influence */
+ _sigmaS[idx_energy][idx_cell] * _scatterMatDiag[idx_system] ); /* scattering influence */
}
}
......
......@@ -13,12 +13,10 @@ PNSolver::PNSolver( Config* settings ) : Solver( settings ) {
_LMaxDegree = settings->GetMaxMomentDegree();
_nTotalEntries = GlobalIndex( int( _LMaxDegree ), int( _LMaxDegree ) ) + 1;
// transform sigmaT and sigmaS in sigmaA.
_sigmaA = VectorVector( _nEnergies, Vector( _nCells, 0 ) ); // Get rid of this extra vektor!
// This must be shifted to Problem Class
for( unsigned n = 0; n < _nEnergies; n++ ) {
for( unsigned j = 0; j < _nCells; j++ ) {
_sigmaA[n][j] = 0; //_sigmaT[n][j] - _sigmaS[n][j];
_sigmaT[n][j] = 1;
_sigmaS[n][j] = 1;
}
}
......@@ -136,7 +134,7 @@ void PNSolver::Solve() {
psiNew[idx_cell][idx_system] = _sol[idx_cell][idx_system] -
( _dE / _areas[idx_cell] ) * psiNew[idx_cell][idx_system] /* cell averaged flux */
- _dE * _sol[idx_cell][idx_system] *
( _sigmaA[idx_energy][idx_cell] /* absorbtion influence */
( _sigmaT[idx_energy][idx_cell] /* absorbtion influence */
+ _sigmaS[idx_energy][idx_cell] * _scatterMatDiag[idx_system] ); /* scattering influence */
}
}
......@@ -329,9 +327,9 @@ void PNSolver::ComputeFluxComponents() {
void PNSolver::ComputeScatterMatrix() {
// --- Isotropic ---
_scatterMatDiag[0] = 0.0;
_scatterMatDiag[0] = -1.0;
for( unsigned idx_diag = 1; idx_diag < _nTotalEntries; idx_diag++ ) {
_scatterMatDiag[idx_diag] = 1.0;
_scatterMatDiag[idx_diag] = 0.0;
}
}
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment