code clean up, Changed linesource_PN - now unit tests fail!

code/include/entropies/entropybase.h

@@ -6,7 +6,9 @@
 class EntropyBase
 {
   public:
-    inline EntropyBase(){};
+    inline EntropyBase() {}
+
+    virtual inline ~EntropyBase() {}
 
     static EntropyBase* Create( Config* settings );
 

code/include/entropies/maxwellboltzmannentropy.h

@@ -7,17 +7,19 @@
 class MaxwellBoltzmannEntropy : public EntropyBase
 {
   public:
-    inline MaxwellBoltzmannEntropy() : EntropyBase(){};
+    inline MaxwellBoltzmannEntropy() : EntropyBase() {}
 
-    inline double Entropy( double z ) override { return z * log( z ) - z; };
+    inline ~MaxwellBoltzmannEntropy() {}
 
-    inline double EntropyPrime( double z ) override { return log( z ); };
+    inline double Entropy( double z ) override { return z * log( z ) - z; }
 
-    inline double EntropyDual( double y ) override { return exp( y ); };
+    inline double EntropyPrime( double z ) override { return log( z ); }
 
-    inline double EntropyPrimeDual( double y ) override { return exp( y ); };
+    inline double EntropyDual( double y ) override { return exp( y ); }
 
-    inline bool CheckDomain( double z ) override { return ( z >= 0 && std::isnan( z ) ); };
+    inline double EntropyPrimeDual( double y ) override { return exp( y ); }
+
+    inline bool CheckDomain( double z ) override { return ( z >= 0 && std::isnan( z ) ); }
 };
 
 #endif    // MAXWELLBOLTZMANNENTROPY_H

code/include/entropies/quadraticentropy.h

@@ -7,17 +7,19 @@
 class QuadraticEntropy : public EntropyBase
 {
   public:
-    inline QuadraticEntropy() : EntropyBase(){};
+    inline QuadraticEntropy() : EntropyBase() {}
 
-    inline double Entropy( double z ) override { return 0.5 * z * z; };
+    inline ~QuadraticEntropy() {}
 
-    inline double EntropyPrime( double z ) override { return z; };
+    inline double Entropy( double z ) override { return 0.5 * z * z; }
 
-    inline double EntropyDual( double y ) override { return 0.5 * y * y; };
+    inline double EntropyPrime( double z ) override { return z; }
 
-    inline double EntropyPrimeDual( double y ) override { return y; };
+    inline double EntropyDual( double y ) override { return 0.5 * y * y; }
 
-    inline bool CheckDomain( double z ) override { return std::isnan( z ); };
+    inline double EntropyPrimeDual( double y ) override { return y; }
+
+    inline bool CheckDomain( double z ) override { return std::isnan( z ); }
 };
 
 #endif    // QUADRATICENTROPY_H

code/include/optimizers/newtonoptimizer.h

@@ -6,7 +6,9 @@
 class NewtonOptimizer : public OptimizerBase
 {
   public:
-    inline NewtonOptimizer( Config* settings ) : OptimizerBase( settings ){};
+    inline NewtonOptimizer( Config* settings ) : OptimizerBase( settings ) {}
+
+    inline ~NewtonOptimizer() {}
 
     virtual Vector Solve( Vector u ) override;
 };

code/include/optimizers/optimizerbase.h

@@ -10,6 +10,8 @@ class OptimizerBase
   public:
     OptimizerBase( Config* settings );
 
+    virtual inline ~OptimizerBase() { delete _entropy; }
+
     static OptimizerBase* Create( Config* settings );
 
     /*! @brief  : Computes the optimal Lagrange multilpiers for the dual entropy minimization problem

code/include/quadratures/qmontecarlo.h

@@ -7,7 +7,7 @@ class QMonteCarlo : public QuadratureBase
 {
   public:
     QMonteCarlo( unsigned order );
-    virtual ~QMonteCarlo() {}
+    inline ~QMonteCarlo() {}
 
     inline void SetName() override { _name = "Monte Carlo Quadrature."; }
     inline void SetNq() override { _nq = GetOrder(); }

code/include/solvers/mnsolver.h

@@ -32,8 +32,9 @@ class MNSolver : public Solver
     VectorVector _sigmaA;      /*! @brief: Absorbtion coefficient for all energies */
     SphericalHarmonics _basis; /*! @brief: Class to compute and store current spherical harmonics basis */
     Vector _scatterMatDiag;    /*! @brief: Diagonal of the scattering matrix (its a diagonal matrix by construction) */
-    VectorVector _A;             /*! @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: Diagonals of the system matrices  (each sysmatric is a diagonal matrix by construction)
+    //                                          layout: Nx2 (in 2d), N = size of system */
+
     QuadratureBase* _quadrature; /*! @brief: Quadrature rule to compute flux Jacobian */
 
     EntropyBase* _entropy; /*! @brief: Class to handle entropy functionals */
@@ -48,9 +49,9 @@ class MNSolver : public Solver
      */
     int GlobalIndex( int l, int k ) const;
 
-    /*! @brief : function for computing and setting up the System Matrices.
-                 The System Matrices are diagonal, so there is no need to diagonalize*/
-    void ComputeSystemMatrices();
+    // /*! @brief : function for computing and setting up the System Matrices.
+    //              The System Matrices are diagonal, so there is no need to diagonalize*/
+    // void ComputeSystemMatrices();
 
     /*! @brief : Construct flux by computing the Moment of the FVM discretization at the interface of cell
      *  @param : idx_cell = current cell id

code/src/problems/linesource.cpp

@@ -66,7 +66,7 @@ std::vector<double> LineSource_PN::GetStoppingPower( const std::vector<double>&
 
 VectorVector LineSource_PN::SetupIC() {
     // Compute number of equations in the system
-    int ntotalEquations = GlobalIndex( _settings->GetNQuadPoints(), _settings->GetNQuadPoints() ) + 1;
+    int ntotalEquations = GlobalIndex( _settings->GetMaxMomentDegree(), _settings->GetMaxMomentDegree() ) + 1;
 
     VectorVector psi( _mesh->GetNumCells(), Vector( ntotalEquations, 0 ) );    // zero could lead to problems?
     VectorVector cellMids = _mesh->GetCellMidPoints();
@@ -76,7 +76,7 @@ VectorVector LineSource_PN::SetupIC() {
     for( unsigned j = 0; j < cellMids.size(); ++j ) {
         double x  = cellMids[j][0];
         double y  = cellMids[j][1];
-        psi[j][0] = 1.0 / ( 4.0 * M_PI * t ) * std::exp( -( x * x + y * y ) / ( 4 * t ) );    // / ( 4 * M_PI );
+        psi[j][0] = 1.0 / ( 4.0 * M_PI * t ) * std::exp( -( x * x + y * y ) / ( 4 * t ) );
     }
 
     // for( unsigned j = 0; j < cellMids.size(); ++j ) {

code/src/problems/problembase.cpp

@@ -11,7 +11,7 @@ ProblemBase* ProblemBase::Create( Config* settings, Mesh* mesh ) {
     auto name = settings->GetProblemName();
     switch( name ) {
         case PROBLEM_LineSource: {
-            if( settings->GetSolverName() == PN_SOLVER )
+            if( settings->GetSolverName() == PN_SOLVER || settings->GetSolverName() == MN_SOLVER )
                 return new LineSource_PN( settings, mesh );
             else
                 return new LineSource_SN( settings, mesh );    // default

code/src/settings/config.cpp

@@ -192,7 +192,7 @@ void Config::SetConfigOptions() {
      * Config*/
     AddEnumOption( "QUAD_TYPE", _quadName, Quadrature_Map, QUAD_MonteCarlo );
     /*!\brief QUAD_ORDER \n DESCRIPTION: Order of Quadrature rule \n DEFAULT 2 \ingroup Config.*/
-    AddUnsignedShortOption( "QUAD_ORDER", _quadOrder, 2 );
+    AddUnsignedShortOption( "QUAD_ORDER", _quadOrder, 1 );
 
     // Solver related options
     /*! @brief MAX_MOMENT_ORDER \n: DESCRIPTON: Specifies the maximal order of Moments for PN and SN Solver */

code/src/solvers/mnsolver.cpp

@@ -6,7 +6,7 @@ MNSolver::MNSolver( Config* settings ) : Solver( settings ), _nMaxMomentsOrder(
     // Is this good (fast) code using a constructor list?
 
     _nTotalEntries = GlobalIndex( _nMaxMomentsOrder, int( _nMaxMomentsOrder ) ) + 1;
-    _quadrature    = QuadratureBase::CreateQuadrature( settings->GetQuadName(), settings->GetQuadOrder() );
+    _quadrature    = QuadratureBase::CreateQuadrature( _settings->GetQuadName(), settings->GetQuadOrder() );
     // transform sigmaT and sigmaS in sigmaA.
     _sigmaA = VectorVector( _nEnergies, Vector( _nCells, 0 ) );    // Get rid of this extra vektor!
 
@@ -21,10 +21,12 @@ MNSolver::MNSolver( Config* settings ) : Solver( settings ), _nMaxMomentsOrder(
     _scatterMatDiag    = Vector( _nTotalEntries, 1.0 );
     _scatterMatDiag[0] = 0.0;    // First entry is zero by construction.
 
-    // Initialize System Matrices
-    _A = VectorVector( _nTotalEntries );
+    // Initialize Entropy
+    _entropy = EntropyBase::Create( _settings );
+
+    // Initialize Optimizer
+    _optimizer = OptimizerBase::Create( _settings );
 
-    _entropy = EntropyBase::Create( settings );
     // Initialize lagrange Multipliers
     _alpha = VectorVector( _nCells );
     for( unsigned idx_cells = 0; idx_cells < _nTotalEntries; idx_cells++ ) {
@@ -32,7 +34,11 @@ MNSolver::MNSolver( Config* settings ) : Solver( settings ), _nMaxMomentsOrder(
     }
 }
 
-MNSolver::~MNSolver() { delete _quadrature; }
+MNSolver::~MNSolver() {
+    delete _quadrature;
+    delete _entropy;
+    delete _optimizer;
+}
 
 int MNSolver::GlobalIndex( int l, int k ) const {
     int numIndicesPrevLevel  = l * l;    // number of previous indices untill level l-1
@@ -40,29 +46,6 @@ int MNSolver::GlobalIndex( int l, int k ) const {
     return numIndicesPrevLevel + prevIndicesThisLevel;
 }
 
-void MNSolver::ComputeSystemMatrices() {
-
-    // Initialize entries of _A
-    for( unsigned idx_sys = 0; idx_sys < _nTotalEntries; idx_sys++ ) {
-        Vector entry( 2, 0.0 );
-        _A[idx_sys] = entry;
-    }
-
-    // Compute Integrals _A[idx_system][i]=<Omega_i*m^(l,k)>
-    Vector funcEval( _nTotalEntries, 0.0 );
-    for( unsigned i = 0; i < _nq; i++ ) {
-        double x = _quadrature->GetPoints()[i][0];
-        double y = _quadrature->GetPoints()[i][1];
-        double z = _quadrature->GetPoints()[i][2];
-        double w = _quadrature->GetWeights()[i];
-        funcEval = _basis.ComputeSphericalBasis( x, y, z );    // = m_l^k
-        for( unsigned idx_system = 0; idx_system < _nTotalEntries; idx_system++ ) {
-            _A[idx_system][0] += w * x * funcEval[idx_system];
-            _A[idx_system][1] += w * y * funcEval[idx_system];
-        }
-    }
-}
-
 Vector MNSolver::ConstructFlux( unsigned idx_cell, unsigned idx_neigh ) {
 
     // ---- Integration of Moment of flux ----
@@ -116,16 +99,21 @@ void MNSolver::Solve() {
 
     // Loop over energies (pseudo-time of continuous slowing down approach)
 
-    for( unsigned idx_energy = 0; idx_energy < _nEnergies; ++idx_energy ) {
+    for( unsigned idx_energy = 0; idx_energy < _nEnergies; idx_energy++ ) {
 
         // Loop over the grid cells
         for( unsigned idx_cell = 0; idx_cell < _nCells; idx_cell++ ) {
 
+            // Reset temporary variable
+            for( unsigned idx_sys = 0; idx_sys < _nTotalEntries; idx_sys++ ) {
+                psiNew[idx_cell][idx_sys] = 0.0;
+            }
+
             // ------- Reconstruction Step -------
 
             _alpha[idx_cell] = _optimizer->Solve( _psi[idx_cell] );
 
-            // ------- Compute Fluxes ---------
+            // ------- Flux Computation Step ---------
 
             // Loop over neighbors
             for( unsigned idx_neigh = 0; idx_neigh < _neighbors[idx_cell].size(); idx_neigh++ ) {
@@ -134,6 +122,7 @@ void MNSolver::Solve() {
             }
 
             // ------ FVM Step ------
+
             // NEED TO VECTORIZE
             for( unsigned idx_system = 0; idx_system < _nTotalEntries; idx_system++ ) {
 
@@ -157,7 +146,7 @@ void MNSolver::Solve() {
         Save( idx_energy );
         dFlux   = blaze::l2Norm( fluxNew - fluxOld );
         fluxOld = fluxNew;
-        if( rank == 0 ) log->info( "{:03.8f}   {:01.5e}", _energies[idx_energy], dFlux );
+        if( rank == 0 ) log->info( "{:03.8f}   {:01.5e} {:01.5e}", _energies[idx_energy], dFlux, mass );
     }
 }
 

code/src/solvers/pnsolver.cpp

@@ -88,8 +88,6 @@ void PNSolver::Solve() {
 
     Save( -1 );    // Save initial condition
 
-    VectorVector cellMids = _mesh->GetCellMidPoints();
-
     int rank;
 
     unsigned idx_system = 0;