code docu for PN solver, change option MAX_MOMENT to specify PN_SOLVER

code/include/settings/config.h

@@ -51,11 +51,11 @@ class Config
     unsigned _nCells;
 
     // Solver
-    double _CFL;                    /*!< @brief CFL Number for Solver*/
-    double _tEnd;                   /*!< @brief Final Time for Simulation */
-    PROBLEM_NAME _problemName;      /*!< @brief Name of predefined Problem   */
-    SOLVER_NAME _solverName;        /*!< @brief Name of the used Solver */
-    unsigned short _maxMomentOrder; /*!< @brief Maximal Order of Moments for PN and MN Solver */
+    double _CFL;                     /*!< @brief CFL Number for Solver*/
+    double _tEnd;                    /*!< @brief Final Time for Simulation */
+    PROBLEM_NAME _problemName;       /*!< @brief Name of predefined Problem   */
+    SOLVER_NAME _solverName;         /*!< @brief Name of the used Solver */
+    unsigned short _maxMomentDegree; /*!< @brief Maximal Order of Moments for PN and MN Solver */
 
     /*!< @brief If true, very low entries (10^-10 or smaller) of the flux matrices will be set to zero,
      * to improve floating point accuracy */
@@ -203,7 +203,7 @@ class Config
     unsigned GetNCells() { return _nCells; }
 
     // Solver Structure
-    unsigned short inline GetMaxMomentOrder() const { return _maxMomentOrder; }
+    unsigned short inline GetMaxMomentDegree() const { return _maxMomentDegree; }
     double inline GetCFL() const { return _CFL; }
     double inline GetTEnd() const { return _tEnd; }
     PROBLEM_NAME inline GetProblemName() const { return _problemName; }

code/include/solvers/pnsolver.h

@@ -8,48 +8,45 @@
 class PNSolver : public Solver
 {
   public:
-    /**
-     * @brief PNSolver constructor
-     * @param settings stores all needed information
+    /*! @brief PNSolver constructor
+     *  @param settings stores all needed information
      */
     PNSolver( Config* settings );
-    /**
-     * @brief Solve functions runs main time loop
-     */
-    virtual void Solve() override;
-    /**
-     * @brief Output solution to VTK file
-     */
-    void Save() const override;
 
-    void Save( int currEnergy ) const override;
+    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 */
 
   protected:
-    // moment orders for P_N
-    unsigned _nTotalEntries;    // total number of equations in the system
+    unsigned _nTotalEntries; /*! @brief: total number of equations in the system */
+    unsigned _LMaxDegree;    /*! @brief: maximal degree of the spherical harmonics basis*/
 
-    VectorVector _sigmaA;    // absorbtion 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;
-    SymMatrix _Ay;
-    SymMatrix _Az;
+    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;
-    Matrix _AxMinus;
-    Matrix _AxAbs;
-    Matrix _AyPlus;
-    Matrix _AyMinus;
-    Matrix _AyAbs;
-    Matrix _AzPlus;
-    Matrix _AzMinus;
-    Matrix _AzAbs;
-    double _combinedSpectralRadius;
-
-    Vector _scatterMatDiag;    // diagonal of the scattering matrix (its a diagonal matrix by construction)
-
-    // parameter functions for setting up system matrix
+    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) */
+
+    /*! @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
+     */
     double AParam( int l, int k ) const;
     double BParam( int l, int k ) const;
     double CParam( int l, int k ) const;
@@ -62,7 +59,9 @@ class PNSolver : public Solver
     double ETilde( int l, int k ) const;
     double FTilde( int l, int k ) const;
 
-    // mathematical + index functions
+    /*! @brief: mathematical + index functions. Helper functions for setting up system matrix.
+     *  @param: k: arbitrary integer
+     */
     int Sgn( int k ) const;
     int kPlus( int k ) const;
     int kMinus( int k ) const;
@@ -79,16 +78,18 @@ class PNSolver : public Solver
      */
     bool CheckIndex( int l, int k ) const;
 
-    // function for computing and setting up system matrices
+    /*! @brief: function for computing and setting up system matrices */
     void ComputeSystemMatrices();
-    // function for computing and setting up flux matrices for upwinding
+    /*! @brief:  function for computing and setting up flux matrices for upwinding */
     void ComputeFluxComponents();
-    // fucntion for computing and setting up diagonal matrix for scattering kernel
+    /*! @brief:  fucntion for computing and setting up diagonal matrix for scattering kernel */
     void ComputeScatterMatrix();
-    // 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 );
-    // Sets Entries of FluxMatrices to zero, if they are below double accuracy, to prevent floating point
-    // inaccuracies lateron
+
+    /*! @brief: Sets Entries of FluxMatrices to zero, if they are below double precision,
+     *          to prevent floating point inaccuracies later in the solver
+     */
     void CleanFluxMatrices();
 };
 

code/src/settings/config.cpp

@@ -196,7 +196,7 @@ void Config::SetConfigOptions() {
 
     // Solver related options
     /*! @brief MAX_MOMENT_ORDER \n: DESCRIPTON: Specifies the maximal order of Moments for PN and SN Solver */
-    AddUnsignedShortOption( "MAX_MOMENT_ORDER", _maxMomentOrder, 1 );
+    AddUnsignedShortOption( "MAX_MOMENT_SOLVER", _maxMomentDegree, 1 );
     /*! @brief CFL \n DESCRIPTION: CFL number \n DEFAULT 1.0 @ingroup Config.*/
     AddDoubleOption( "CFL_NUMBER", _CFL, 1.0 );
     /*! @brief TIME_FINAL \n DESCRIPTION: Final time for simulation \n DEFAULT 1.0 @ingroup Config.*/

code/src/solvers/mnsolver.cpp

 #include "solvers/mnsolver.h"
 
-MNSolver::MNSolver( Config* settings ) : Solver( settings ), _nMaxMomentsOrder( settings->GetMaxMomentOrder() ), _basis( _nMaxMomentsOrder ) {
+MNSolver::MNSolver( Config* settings ) : Solver( settings ), _nMaxMomentsOrder( settings->GetMaxMomentDegree() ), _basis( _nMaxMomentsOrder ) {
     // Is this good (fast) code using a constructor list?
 
     _nTotalEntries = GlobalIndex( _nMaxMomentsOrder, int( _nMaxMomentsOrder ) ) + 1;

code/src/solvers/pnsolver.cpp

@@ -5,7 +5,8 @@
 #include <mpi.h>
 
 PNSolver::PNSolver( Config* settings ) : Solver( settings ) {
-    _nTotalEntries = GlobalIndex( _nq, int( _nq ) ) + 1;
+    _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!
@@ -332,7 +333,7 @@ void PNSolver::ComputeFluxComponents() {
     // std::cout << "Spectral Radius Y " << blaze::max( blaze::abs( eigenValuesY ) ) << "\n";
     // std::cout << "Spectral Radius Z " << blaze::max( blaze::abs( eigenValues ) ) << "\n";
 
-    _combinedSpectralRadius = blaze::max( blaze::abs( eigenValues + eigenValuesX + eigenValuesY ) );
+    // _combinedSpectralRadius = blaze::max( blaze::abs( eigenValues + eigenValuesX + eigenValuesY ) );
     // std::cout << "Spectral Radius combined " << _combinedSpectralRadius << "\n";
 }
 
@@ -463,7 +464,7 @@ int PNSolver::GlobalIndex( int l, int k ) const {
 }
 
 bool PNSolver::CheckIndex( int l, int k ) const {
-    if( l >= 0 && l <= int( _nq ) ) {
+    if( l >= 0 && l <= int( _LMaxDegree ) ) {
         if( k >= -l && k <= l ) return true;
     }
     return false;

code/tests/input/linesource_PN.cfg

@@ -9,7 +9,8 @@ CFL_NUMBER = 0.5
 TIME_FINAL = 0.5
 PROBLEM = LINESOURCE
 SOLVER = PN_SOLVER
+MAX_MOMENT_SOLVER = 2
 
 % ---- Boundary Conditions ----
 BC_NEUMANN = ( void )
-QUAD_ORDER = 2
+