maxMoment option added, started rebuilding mnsolver

parmetis @ e68c8e39

-Subproject commit 54affdf3f9e1d55c0047c0401401da308aea8878
+Subproject commit e68c8e3947f9383e28cb8623d116ca996b3651a9

code/include/settings/config.h

@@ -51,10 +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 */
+    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 */
 
     /*!< @brief If true, very low entries (10^-10 or smaller) of the flux matrices will be set to zero,
      * to improve floating point accuracy */
@@ -202,6 +203,7 @@ class Config
     unsigned GetNCells() { return _nCells; }
 
     // Solver Structure
+    unsigned short inline GetMaxMomentOrder() const { return _maxMomentOrder; }
     double inline GetCFL() const { return _CFL; }
     double inline GetTEnd() const { return _tEnd; }
     PROBLEM_NAME inline GetProblemName() const { return _problemName; }

code/include/solvers/mnsolver.h

@@ -21,6 +21,8 @@ class MNSolver : Solver
   private:
     /*! @brief: Total number of equations in the system */
     unsigned _nTotalEntries;
+    /*! @brief: Max Order of Moments */
+    unsigned short _nMaxMomentsOrder;
 
     /*! @brief: Absorbtion coefficient for all energies */
     VectorVector _sigmaA;
@@ -31,14 +33,17 @@ class MNSolver : Solver
     /*! @brief: Diagonal of the scattering matrix (its a diagonal matrix by construction) */
     Vector _scatterMatDiag;
 
-    /*! @brief: Diagonal of the system matrix in x direction (its a diagonal matrix by construction) */
-    Vector _Ax;
-
-    /*! @brief: Diagonal of the system matrix in y direction (its a diagonal matrix by construction) */
-    Vector _Ay;
-
-    /*! @brief: Diagonal of the system matrix in z direction (its a diagonal matrix by construction) */
-    Vector _Az;
+    /*! @brief: Diagonals of the system matrices  (each sysmatric is a diagonal matrix by construction)
+                layout: Nx2 (in 2d), N = size of system */
+    VectorVector _A;
+    ///*! @brief: Diagonal of the system matrix in x direction (its a diagonal matrix by construction) */
+    // Vector _Ax;
+    //
+    ///*! @brief: Diagonal of the system matrix in y direction (its a diagonal matrix by construction) */
+    // Vector _Ay;
+    //
+    ///*! @brief: Diagonal of the system matrix in z direction (its a diagonal matrix by construction) */
+    // Vector _Az;
 
     /*! @brief : computes the global index of the moment corresponding to basis function (l,k)
      *  @param : degree l, it must hold: 0 <= l <=_nq

code/src/settings/config.cpp

@@ -195,6 +195,8 @@ void Config::SetConfigOptions() {
     AddUnsignedShortOption( "QUAD_ORDER", _quadOrder, 2 );
 
     // 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 );
     /*! @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 <mnsolver.h>
 
-MNSolver::MNSolver( Config* settings )
-    : Solver( settings ), _basis( _nq ) /* We need max degree +1, since we need the N+1st moment to reconstruct phi */
-{
+MNSolver::MNSolver( Config* settings ) : Solver( settings ), _nMaxMomentsOrder( settings->GetMaxMomentOrder() ), _basis( _nMaxMomentsOrder ) {
     // Is this good (fast) code using a constructor list?
 
-    _nTotalEntries = GlobalIndex( _nq, int( _nq ) ) + 1;
+    _nTotalEntries = GlobalIndex( _nMaxMomentsOrder, int( _nMaxMomentsOrder ) ) + 1;
 
     // transform sigmaT and sigmaS in sigmaA.
     _sigmaA = VectorVector( _nEnergies, Vector( _nCells, 0 ) );    // Get rid of this extra vektor!
@@ -22,9 +20,7 @@ MNSolver::MNSolver( Config* settings )
     _scatterMatDiag[0] = 0.0;    // First entry is zero by construction.
 
     // Initialize System Matrices
-    _Ax = Vector( _nTotalEntries, 0.0 );
-    _Ay = Vector( _nTotalEntries, 0.0 );
-    _Az = Vector( _nTotalEntries, 0.0 );
+    _A = VectorVector( _nTotalEntries );
 
     // Fill System Matrices
     ComputeSystemMatrices();
@@ -36,8 +32,14 @@ int MNSolver::GlobalIndex( int l, int k ) const {
     return numIndicesPrevLevel + prevIndicesThisLevel;
 }
 
-void MNSolver::ComputeSystemMatrices() {}
-
+void MNSolver::ComputeSystemMatrices() {
+    for( int l_idx = 0; l_idx <= _nMaxMomentsOrder; l_idx++ ) {
+        for( unsigned k_idx = -l_idx; k_idx <= l_idx; k_idx++ ) {
+            _A[idx_sys] = Vector( 2, 0 );
+            // compute entries. =<Omega_i*m^(l,k)>
+        }
+    }
+}
 void MNSolver::Solve() {
 
     int rank;
@@ -73,7 +75,7 @@ void MNSolver::Solve() {
             if( _boundaryCells[idx_cell] == BOUNDARY_TYPE::DIRICHLET ) continue;    // Dirichlet cells stay at IC, farfield assumption
 
             // Loop over all equations of the system.
-            for( int idx_lDegree = 0; idx_lDegree <= int( _nq ); idx_lDegree++ ) {
+            for( int idx_lDegree = 0; idx_lDegree <= int( _nMaxMomentsOrder ); idx_lDegree++ ) {
                 for( int idx_kOrder = -idx_lDegree; idx_kOrder <= idx_lDegree; idx_kOrder++ ) {
                     idx_system                   = unsigned( GlobalIndex( idx_lDegree, idx_kOrder ) );
                     psiNew[idx_cell][idx_system] = 0.0;