dry run MN Solver works. Flux Jacobians must be reworked

code/include/quadratures/qgausslegendretensorized.h

@@ -14,7 +14,7 @@ class QGaussLegendreTensorized : public QuadratureBase
     QGaussLegendreTensorized( unsigned order );
     virtual ~QGaussLegendreTensorized() {}
 
-    inline void SetName() override { _name = "Tensorized Gauss-Legendre quadrature."; }
+    inline void SetName() override { _name = "Tensorized Gauss-Legendre quadrature"; }
     inline void SetNq() override { _nq = pow( GetOrder(), 2 ); }
     void SetPointsAndWeights() override;
     void SetConnectivity() override;

code/include/settings/globalconstants.h

@@ -49,8 +49,8 @@ enum KERNEL_NAME { KERNEL_Isotropic };
 inline std::map<std::string, KERNEL_NAME> Kernel_Map{ { "ISOTROPIC", KERNEL_Isotropic } };
 
 // Solver name
-enum SOLVER_NAME { SN_SOLVER, PN_SOLVER };
+enum SOLVER_NAME { SN_SOLVER, PN_SOLVER, MN_SOLVER };
 
-inline std::map<std::string, SOLVER_NAME> Solver_Map{ { "SN_SOLVER", SN_SOLVER }, { "PN_SOLVER", PN_SOLVER } };
+inline std::map<std::string, SOLVER_NAME> Solver_Map{ { "SN_SOLVER", SN_SOLVER }, { "PN_SOLVER", PN_SOLVER }, { "MN_SOLVER", MN_SOLVER } };
 
 #endif    // GLOBAL_CONSTANTS_H

code/include/solvers/mnsolver.h

@@ -4,7 +4,7 @@
 #include "solverbase.h"
 #include "sphericalharmonics.h"
 
-class MNSolver : Solver
+class MNSolver : public Solver
 {
   public:
     /**
@@ -20,6 +20,8 @@ class MNSolver : Solver
      * @brief Solve functions runs main time loop
      */
     void Solve() override;
+    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 */
 
   private:
     unsigned _nTotalEntries;          /*! @brief: Total number of equations in the system */

code/include/solvers/snsolver.h

 #ifndef SNSOLVER_H
 #define SNSOLVER_H
 
-#include <mpi.h>
-
 #include "solvers/solverbase.h"
 
 class SNSolver : public Solver

code/src/io.cpp

@@ -279,7 +279,7 @@ std::string ParseArguments( int argc, char* argv[] ) {
             inputFile = std::string( argv[i] );
             std::ifstream f( inputFile );
             if( !f.is_open() ) {
-                ErrorMessages::OptionNotSetError( "Unable to open inputfile '" + inputFile + "' !", CURRENT_FUNCTION );
+                ErrorMessages::Error( "Unable to open inputfile '" + inputFile + "' !", CURRENT_FUNCTION );
             }
         }
     }

code/src/quadratures/qgausslegendretensorized.cpp

@@ -151,7 +151,8 @@ double QGaussLegendreTensorized::Pythag( const double a, const double b ) {
 
 bool QGaussLegendreTensorized::CheckOrder() {
     if( _order % 2 == 1 ) {    // order needs to be even
-        ErrorMessages::Error( "ERROR! Order " + std::to_string( _order ) + " for " + GetName() + " not available. ", CURRENT_FUNCTION );
+        ErrorMessages::Error( "ERROR! Order " + std::to_string( _order ) + " for " + GetName() + " not available. \n Order must be an even number. ",
+                              CURRENT_FUNCTION );
     }
     return true;
 }

code/src/solvers/mnsolver.cpp

 #include "solvers/mnsolver.h"
+#include "toolboxes/textprocessingtoolbox.h"
+#include <mpi.h>
 
 MNSolver::MNSolver( Config* settings ) : Solver( settings ), _nMaxMomentsOrder( settings->GetMaxMomentDegree() ), _basis( _nMaxMomentsOrder ) {
     // Is this good (fast) code using a constructor list?
@@ -24,6 +26,7 @@ MNSolver::MNSolver( Config* settings ) : Solver( settings ), _nMaxMomentsOrder(
 
     // Fill System Matrices
     ComputeSystemMatrices();
+    TextProcessingToolbox::PrintVectorVector( _A );
 }
 
 MNSolver::~MNSolver() { delete _quadrature; }
@@ -144,3 +147,23 @@ void MNSolver::Solve() {
         if( rank == 0 ) log->info( "{:03.8f}   {:01.5e}", _energies[idx_energy], dFlux );
     }
 }
+
+void MNSolver::Save() const {
+    std::vector<std::string> fieldNames{ "flux" };
+    std::vector<double> flux;
+    flux.resize( _nCells );
+
+    for( unsigned i = 0; i < _nCells; ++i ) {
+        flux[i] = _psi[i][0];
+    }
+    std::vector<std::vector<double>> scalarField( 1, flux );
+    std::vector<std::vector<std::vector<double>>> results{ scalarField };
+    ExportVTK( _settings->GetOutputFile(), results, fieldNames, _mesh );
+}
+
+void MNSolver::Save( int currEnergy ) const {
+    std::vector<std::string> fieldNames{ "flux" };
+    std::vector<std::vector<double>> scalarField( 1, _solverOutput );
+    std::vector<std::vector<std::vector<double>>> results{ scalarField };
+    ExportVTK( _settings->GetOutputFile() + "_" + std::to_string( currEnergy ), results, fieldNames, _mesh );
+}

code/src/solvers/snsolver.cpp

 #include "solvers/snsolver.h"
+#include <mpi.h>
 
 SNSolver::SNSolver( Config* settings ) : Solver( settings ) {}
 

code/src/solvers/solverbase.cpp

@@ -3,6 +3,7 @@
 #include "mesh.h"
 #include "quadratures/quadraturebase.h"
 #include "settings/globalconstants.h"
+#include "solvers/mnsolver.h"
 #include "solvers/pnsolver.h"
 #include "solvers/snsolver.h"
 
@@ -71,6 +72,7 @@ Solver* Solver::Create( Config* settings ) {
     switch( settings->GetSolverName() ) {
         case SN_SOLVER: return new SNSolver( settings );
         case PN_SOLVER: return new PNSolver( settings );
+        case MN_SOLVER: return new MNSolver( settings );
         default: return new SNSolver( settings );
     }
 }