started cleanup of solvers

code/include/kernels/scatteringkernelbase.h

@@ -15,7 +15,7 @@ class ScatteringKernel
 
   public:
     ScatteringKernel( QuadratureBase* quad );
-    ~ScatteringKernel();
+    virtual ~ScatteringKernel();
 
     virtual Matrix GetScatteringKernel() = 0;
 

code/include/solvers/csdsnsolver.h

@@ -7,6 +7,7 @@ class CSDSNSolver : public Solver
 {
   private:
     std::vector<double> _dose;
+    Matrix _scatteringKernel; /*!  @brief scattering kernel for the quadrature */
 
   public:
     /**

code/include/solvers/mnsolver.h

@@ -31,8 +31,7 @@ class MNSolver : public Solver
     unsigned _nTotalEntries;          /*! @brief: Total number of equations in the system */
     unsigned short _nMaxMomentsOrder; /*! @brief: Max Order of Moments */
 
-    VectorVector _sigmaA; /*! @brief: Absorbtion coefficient for all energies */
-
+    VectorVector _sigmaA;       /*!  @brief: Absorption coefficient for all energies*/
     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 */
 

code/include/solvers/pnsolver.h

@@ -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 */

code/include/solvers/snsolver.h

@@ -6,6 +6,8 @@
 class SNSolver : public Solver
 {
   private:
+    Matrix _scatteringKernel; /*!  @brief scattering kernel for the quadrature */
+
   public:
     /**
      * @brief SNSolver constructor

code/include/solvers/solverbase.h

@@ -15,32 +15,42 @@ class QuadratureBase;
 class Solver
 {
   protected:
-    unsigned _nq;                                     // number of quadrature points
-    unsigned _nCells;                                 // number of spatial cells
-    unsigned _nEnergies;                              // number of energy/time steps, number of nodal energy values for CSD
-    double _dE;                                       // energy/time step size
-    std::vector<double> _energies;                    // energy groups used in the simulation [keV]
-    VectorVector _psi;                                // angular flux vector, dim(_psi) = (_NCells,_nq)
-    std::vector<double> _areas;                       // surface area of all spatial cells, dim(_areas) = _NCells
-    std::vector<std::vector<Vector>> _normals;        // edge normals multiplied by edge length, dim(_normals) = (_NCells,nEdgesPerCell,spatialDim)
-    std::vector<std::vector<unsigned>> _neighbors;    // edge normals multiplied by edge length, dim(_neighbors) = (_NCells,nEdgesPerCell)
-    std::vector<double> _density;                     // patient density, dim(_density) = _nCells
-    std::vector<double> _s;                           // stopping power, dim(_s) = _nTimeSteps
-    VectorVector _sigmaS;                             // scattering cross section for all energies
-    VectorVector _sigmaT;                             // total cross section for all energies
-    std::vector<VectorVector> _Q;                     // external source term
-    Matrix _scatteringKernel;                         // scattering kernel for the quadrature
-    VectorVector _quadPoints;                         // quadrature points, dim(_quadPoints) = (_nSystem,spatialDim)
-    Vector _weights;                                  // quadrature weights, dim(_weights) = (_NCells)
-    std::vector<BOUNDARY_TYPE> _boundaryCells;        // boundary type for all cells, dim(_boundary) = (_NCells)
-    std::vector<double> _solverOutput;                // PROTOTYPE: Outputfield for solver
+    Mesh* _mesh;           /*! @brief mesh object for writing out information */
+    NumericalFlux* _g;     /*! @brief class for numerical flux */
+    Config* _settings;     /*! @brief config class for global information */
+    ProblemBase* _problem; /*! @brief problem class for initial conditions */
 
-    // we will have to add a further dimension for quadPoints and weights once we start with multilevel SN
+    // --------- Often used variables of member classes for faster access ----
+
+    unsigned _nEnergies;              /*! @brief number of energy/time steps, number of nodal energy values for CSD */
+    double _dE;                       /*! @brief energy/time step size */
+    std::vector<double> _energies;    // energy groups used in the simulation [keV]
+    std::vector<double> _density;     // patient density, dim(_density) = _nCells
+    std::vector<double> _s;           // stopping power, dim(_s) = _nTimeSteps
+    std::vector<VectorVector> _Q;     /*!  @brief  external source term */
+
+    VectorVector _sigmaS; /*!  @brief scattering cross section for all energies */
+    VectorVector _sigmaT; /*!  @brief total cross section for all energies */
+
+    // quadrature related numbers
+    unsigned _nq;             /*! @brief number of quadrature points */
+    VectorVector _quadPoints; /*!  @brief quadrature points, dim(_quadPoints) = (_nSystem,spatialDim) */
+    Vector _weights;          /*!  @brief quadrature weights, dim(_weights) = (_NCells) */
 
-    NumericalFlux* _g;    // numerical flux function
-    Mesh* _mesh;          // mesh object for writing out information
-    Config* _settings;
-    ProblemBase* _problem;
+    // Mesh related members
+    unsigned _nCells;                          /*! @brief number of spatial cells */
+    std::vector<BOUNDARY_TYPE> _boundaryCells; /*! boundary type for all cells, dim(_boundary) = (_NCells) */
+    std::vector<double> _areas;                /*! @brief surface area of all spatial cells, dim(_areas) = _NCells */
+    /*! @brief edge normals multiplied by edge length, dim(_normals) = (_NCells,nEdgesPerCell,spatialDim) */
+    std::vector<std::vector<Vector>> _normals;
+    /*! @brief edge neighbor cell ids, dim(_neighbors) = (_NCells,nEdgesPerCell) */
+    std::vector<std::vector<unsigned>> _neighbors;
+
+    // Solution related members
+    VectorVector _sol;                 /*! @brief solution of the PDE, e.g. angular flux or moments */
+    std::vector<double> _solverOutput; /*! @brief PROTOTYPE: Outputfield for solver */
+
+    // we will have to add a further dimension for quadPoints and weights once we start with multilevel SN
 
     /**
      * @brief ComputeTimeStep calculates the maximal stable time step

code/src/solvers/csdsnsolver.cpp

 #include "fluxes/numericalflux.h"
 #include "io.h"
+#include "kernels/scatteringkernelbase.h"
 #include "settings/config.h"
 
 #include "solvers/csdsnsolver.h"
 
-CSDSNSolver::CSDSNSolver( Config* settings ) : Solver( settings ) { _dose = std::vector<double>( _settings->GetNCells(), 0.0 ); }
+CSDSNSolver::CSDSNSolver( Config* settings ) : Solver( settings ) {
+    _dose = std::vector<double>( _settings->GetNCells(), 0.0 );
+
+    QuadratureBase* quad = QuadratureBase::CreateQuadrature( settings->GetQuadName(), settings->GetQuadOrder() );
+    ScatteringKernel* k  = ScatteringKernel::CreateScatteringKernel( settings->GetKernelName(), quad );
+    _scatteringKernel    = k->GetScatteringKernel();
+    delete quad;
+    delete k;
+}
 
 void CSDSNSolver::Solve() {
     auto log = spdlog::get( "event" );
 
     // angular flux at next time step (maybe store angular flux at all time steps, since time becomes energy?)
-    VectorVector psiNew = _psi;
+    VectorVector psiNew = _sol;
     double dFlux        = 1e10;
     Vector fluxNew( _nCells, 0.0 );
     Vector fluxOld( _nCells, 0.0 );
     for( unsigned j = 0; j < _nCells; ++j ) {
-        fluxOld[j] = dot( _psi[j], _weights );
+        fluxOld[j] = dot( _sol[j], _weights );
     }
     int rank;
     MPI_Comm_rank( MPI_COMM_WORLD, &rank );
@@ -24,7 +33,7 @@ void CSDSNSolver::Solve() {
     // do substitution from psi to psiTildeHat (cf. Dissertation Kerstion Kuepper, Eq. 1.23)
     for( unsigned j = 0; j < _nCells; ++j ) {
         for( unsigned k = 0; k < _nq; ++k ) {
-            _psi[j][k] = _psi[j][k] * _density[j] * _s[_nEnergies - 1];    // note that _s[_nEnergies - 1] is stopping power at highest energy
+            _sol[j][k] = _sol[j][k] * _density[j] * _s[_nEnergies - 1];    // note that _s[_nEnergies - 1] is stopping power at highest energy
         }
     }
 
@@ -60,16 +69,16 @@ void CSDSNSolver::Solve() {
                     //                            _normals[idx_cell][idx_neighbor] );
                     else
                         psiNew[idx_cell][idx_ord] += _g->Flux( _quadPoints[idx_ord] / _density[idx_cell],
-                                                               _psi[idx_cell][idx_ord],
-                                                               _psi[_neighbors[idx_cell][idx_neighbor]][idx_ord],
+                                                               _sol[idx_cell][idx_ord],
+                                                               _sol[_neighbors[idx_cell][idx_neighbor]][idx_ord],
                                                                _normals[idx_cell][idx_neighbor] );
                 }
                 // time update angular flux with numerical flux and total scattering cross section
-                psiNew[idx_cell][idx_ord] = _psi[idx_cell][idx_ord] - ( _dE / _areas[idx_cell] ) * psiNew[idx_cell][idx_ord] -
-                                            _dE * _sigmaT[idx_energy][idx_cell] * _psi[idx_cell][idx_ord];
+                psiNew[idx_cell][idx_ord] = _sol[idx_cell][idx_ord] - ( _dE / _areas[idx_cell] ) * psiNew[idx_cell][idx_ord] -
+                                            _dE * _sigmaT[idx_energy][idx_cell] * _sol[idx_cell][idx_ord];
             }
             // compute scattering effects
-            psiNew[idx_cell] += _dE * _sigmaS[idx_energy][idx_cell] * _scatteringKernel * _psi[idx_cell];    // multiply scattering matrix with psi
+            psiNew[idx_cell] += _dE * _sigmaS[idx_energy][idx_cell] * _scatteringKernel * _sol[idx_cell];    // multiply scattering matrix with psi
 
             // TODO: figure out a more elegant way
             // add external source contribution
@@ -86,12 +95,12 @@ void CSDSNSolver::Solve() {
                     psiNew[idx_cell] += _dE * _Q[idx_energy][idx_cell] * _s[_nEnergies - idx_energy - 1];
             }
         }
-        _psi = psiNew;
+        _sol = psiNew;
 
         // do backsubstitution from psiTildeHat to psi (cf. Dissertation Kerstion Kuepper, Eq. 1.23)
         for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
             for( unsigned idx_ord = 0; idx_ord < _nq; ++idx_ord ) {
-                psiNew[idx_cell][idx_ord] = _psi[idx_cell][idx_ord] * _density[idx_cell] *
+                psiNew[idx_cell][idx_ord] = _sol[idx_cell][idx_ord] * _density[idx_cell] *
                                             _s[_nEnergies - idx_energy - 1];    // note that _s[0] is stopping power at lowest energy
             }
         }

code/src/solvers/mnsolver.cpp

@@ -113,15 +113,15 @@ void MNSolver::Solve() {
     auto log = spdlog::get( "event" );
 
     // angular flux at next time step (maybe store angular flux at all time steps, since time becomes energy?)
-    VectorVector psiNew = _psi;
+    VectorVector psiNew = _sol;
     double dFlux        = 1e10;
     Vector fluxNew( _nCells, 0.0 );
     Vector fluxOld( _nCells, 0.0 );
 
     double mass1 = 0;
     for( unsigned i = 0; i < _nCells; ++i ) {
-        _solverOutput[i] = _psi[i][0];
-        mass1 += _psi[i][0] * _areas[i];
+        _solverOutput[i] = _sol[i][0];
+        mass1 += _sol[i][0] * _areas[i];
     }
 
     dFlux   = blaze::l2Norm( fluxNew - fluxOld );
@@ -145,7 +145,7 @@ void MNSolver::Solve() {
 
             // ------- Reconstruction Step -------
 
-            _alpha[idx_cell] = _psi[idx_cell];    // _optimizer->Solve( _psi[idx_cell] );
+            _alpha[idx_cell] = _sol[idx_cell];    // _optimizer->Solve( _psi[idx_cell] );
 
             // ------- Flux Computation Step ---------
 
@@ -159,21 +159,21 @@ void MNSolver::Solve() {
             // NEED TO VECTORIZE
             for( unsigned idx_system = 0; idx_system < _nTotalEntries; idx_system++ ) {
 
-                psiNew[idx_cell][idx_system] = _psi[idx_cell][idx_system] -
+                psiNew[idx_cell][idx_system] = _sol[idx_cell][idx_system] -
                                                ( _dE / _areas[idx_cell] ) * psiNew[idx_cell][idx_system] /* cell averaged flux */
-                                               - _dE * _psi[idx_cell][idx_system] *
+                                               - _dE * _sol[idx_cell][idx_system] *
                                                      ( _sigmaA[idx_energy][idx_cell]                                    /* absorbtion influence */
                                                        + _sigmaS[idx_energy][idx_cell] * _scatterMatDiag[idx_system] ); /* scattering influence */
             }
         }
-        _psi = psiNew;
+        _sol = psiNew;
 
         // pseudo time iteration output
         double mass = 0.0;
         for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
-            fluxNew[idx_cell]       = _psi[idx_cell][0];    // zeroth moment is raditation densitiy we are interested in
-            _solverOutput[idx_cell] = _psi[idx_cell][0];
-            mass += _psi[idx_cell][0] * _areas[idx_cell];
+            fluxNew[idx_cell]       = _sol[idx_cell][0];    // zeroth moment is raditation densitiy we are interested in
+            _solverOutput[idx_cell] = _sol[idx_cell][0];
+            mass += _sol[idx_cell][0] * _areas[idx_cell];
         }
 
         dFlux   = blaze::l2Norm( fluxNew - fluxOld );
@@ -189,7 +189,7 @@ void MNSolver::Save() const {
     flux.resize( _nCells );
 
     for( unsigned i = 0; i < _nCells; ++i ) {
-        flux[i] = _psi[i][0];
+        flux[i] = _sol[i][0];
     }
     std::vector<std::vector<double>> scalarField( 1, flux );
     std::vector<std::vector<std::vector<double>>> results{ scalarField };

code/src/solvers/pnsolver.cpp

@@ -79,15 +79,15 @@ void PNSolver::Solve() {
     auto log = spdlog::get( "event" );
 
     // angular flux at next time step (maybe store angular flux at all time steps, since time becomes energy?)
-    VectorVector psiNew = _psi;
+    VectorVector psiNew = _sol;
     double dFlux        = 1e10;
     Vector fluxNew( _nCells, 0.0 );
     Vector fluxOld( _nCells, 0.0 );
 
     double mass1 = 0;
     for( unsigned i = 0; i < _nCells; ++i ) {
-        _solverOutput[i] = _psi[i][0];
-        mass1 += _psi[i][0];
+        _solverOutput[i] = _sol[i][0];
+        mass1 += _sol[i][0];
     }
 
     dFlux   = blaze::l2Norm( fluxNew - fluxOld );
@@ -121,7 +121,7 @@ void PNSolver::Solve() {
                 // Compute flux contribution and store in psiNew to save memory
                 if( _boundaryCells[idx_cell] == BOUNDARY_TYPE::NEUMANN && _neighbors[idx_cell][idx_neighbor] == _nCells )
                     psiNew[idx_cell] += _g->Flux(
-                        _AxPlus, _AxMinus, _AyPlus, _AyMinus, _AzPlus, _AzMinus, _psi[idx_cell], _psi[idx_cell], _normals[idx_cell][idx_neighbor] );
+                        _AxPlus, _AxMinus, _AyPlus, _AyMinus, _AzPlus, _AzMinus, _sol[idx_cell], _sol[idx_cell], _normals[idx_cell][idx_neighbor] );
                 else
                     psiNew[idx_cell] += _g->Flux( _AxPlus,
                                                   _AxMinus,
@@ -129,8 +129,8 @@ void PNSolver::Solve() {
                                                   _AyMinus,
                                                   _AzPlus,
                                                   _AzMinus,
-                                                  _psi[idx_cell],
-                                                  _psi[_neighbors[idx_cell][idx_neighbor]],
+                                                  _sol[idx_cell],
+                                                  _sol[_neighbors[idx_cell][idx_neighbor]],
                                                   _normals[idx_cell][idx_neighbor] );
             }
 
@@ -139,21 +139,21 @@ void PNSolver::Solve() {
                 for( int idx_kOrder = -idx_lOrder; idx_kOrder <= idx_lOrder; idx_kOrder++ ) {
                     idx_system = unsigned( GlobalIndex( idx_lOrder, idx_kOrder ) );
 
-                    psiNew[idx_cell][idx_system] = _psi[idx_cell][idx_system] -
+                    psiNew[idx_cell][idx_system] = _sol[idx_cell][idx_system] -
                                                    ( _dE / _areas[idx_cell] ) * psiNew[idx_cell][idx_system] /* cell averaged flux */
-                                                   - _dE * _psi[idx_cell][idx_system] *
+                                                   - _dE * _sol[idx_cell][idx_system] *
                                                          ( _sigmaA[idx_energy][idx_cell]                                    /* absorbtion influence */
                                                            + _sigmaS[idx_energy][idx_cell] * _scatterMatDiag[idx_system] ); /* scattering influence */
                 }
             }
         }
-        _psi = psiNew;
+        _sol = psiNew;
 
         double mass = 0.0;
         for( unsigned i = 0; i < _nCells; ++i ) {
-            fluxNew[i]       = _psi[i][0];    // zeroth moment is raditation densitiy we are interested in
-            _solverOutput[i] = _psi[i][0];
-            mass += _psi[i][0] * _areas[i];
+            fluxNew[i]       = _sol[i][0];    // zeroth moment is raditation densitiy we are interested in
+            _solverOutput[i] = _sol[i][0];
+            mass += _sol[i][0] * _areas[i];
         }
 
         dFlux   = blaze::l2Norm( fluxNew - fluxOld );
@@ -371,7 +371,7 @@ void PNSolver::Save() const {
     flux.resize( _nCells );
 
     for( unsigned i = 0; i < _nCells; ++i ) {
-        flux[i] = _psi[i][0];
+        flux[i] = _sol[i][0];
     }
     std::vector<std::vector<double>> scalarField( 1, flux );
     std::vector<std::vector<std::vector<double>>> results{ scalarField };

code/src/solvers/snsolver.cpp

 #include "solvers/snsolver.h"
 #include "fluxes/numericalflux.h"
 #include "io.h"
+#include "kernels/scatteringkernelbase.h"
 #include "mesh.h"
 #include "settings/config.h"
 
 #include <mpi.h>
 
-SNSolver::SNSolver( Config* settings ) : Solver( settings ) {}
+SNSolver::SNSolver( Config* settings ) : Solver( settings ) {
+
+    QuadratureBase* quad = QuadratureBase::CreateQuadrature( settings->GetQuadName(), settings->GetQuadOrder() );
+    ScatteringKernel* k  = ScatteringKernel::CreateScatteringKernel( settings->GetKernelName(), quad );
+    _scatteringKernel    = k->GetScatteringKernel();
+    delete k;
+    delete quad;
+}
 
 void SNSolver::Solve() {
     auto log = spdlog::get( "event" );
 
     // angular flux at next time step (maybe store angular flux at all time steps, since time becomes energy?)
-    VectorVector psiNew = _psi;
+    VectorVector psiNew = _sol;
 
     // reconstruction order
     unsigned reconsOrder = _settings->GetReconsOrder();
@@ -63,7 +71,7 @@ void SNSolver::Solve() {
 
         // reconstruct slopes for higher order solver
         if( reconsOrder > 1 ) {
-            _mesh->ReconstructSlopesU( _nq, psiDx, psiDy, _psi );    // unstructured reconstruction
+            _mesh->ReconstructSlopesU( _nq, psiDx, psiDy, _sol );    // unstructured reconstruction
             //_mesh->ReconstructSlopesS( _nq, psiDx, psiDy, _psi );    // structured reconstruction (not stable currently)
         }
 
@@ -77,24 +85,24 @@ void SNSolver::Solve() {
                 for( unsigned l = 0; l < _neighbors[j].size(); ++l ) {
                     // store flux contribution on psiNew_sigmaS to save memory
                     if( _boundaryCells[j] == BOUNDARY_TYPE::NEUMANN && _neighbors[j][l] == _nCells )
-                        psiNew[j][k] += _g->Flux( _quadPoints[k], _psi[j][k], _psi[j][k], _normals[j][l] );
+                        psiNew[j][k] += _g->Flux( _quadPoints[k], _sol[j][k], _sol[j][k], _normals[j][l] );
                     else {
                         switch( reconsOrder ) {
                             // first order solver
-                            case 1: psiNew[j][k] += _g->Flux( _quadPoints[k], _psi[j][k], _psi[_neighbors[j][l]][k], _normals[j][l] ); break;
+                            case 1: psiNew[j][k] += _g->Flux( _quadPoints[k], _sol[j][k], _sol[_neighbors[j][l]][k], _normals[j][l] ); break;
                             // second order solver
                             case 2:
                                 // left status of interface
-                                psiL = _psi[j][k] + psiDx[j][k] * ( interfaceMidPoints[j][l][0] - cellMidPoints[j][0] ) +
+                                psiL = _sol[j][k] + psiDx[j][k] * ( interfaceMidPoints[j][l][0] - cellMidPoints[j][0] ) +
                                        psiDy[j][k] * ( interfaceMidPoints[j][l][1] - cellMidPoints[j][1] );
                                 // right status of interface
-                                psiR = _psi[_neighbors[j][l]][k] +
+                                psiR = _sol[_neighbors[j][l]][k] +
                                        psiDx[_neighbors[j][l]][k] * ( interfaceMidPoints[j][l][0] - cellMidPoints[_neighbors[j][l]][0] ) +
                                        psiDy[_neighbors[j][l]][k] * ( interfaceMidPoints[j][l][1] - cellMidPoints[_neighbors[j][l]][1] );
                                 // positivity check (if not satisfied, deduce to first order)
                                 if( psiL < 0.0 || psiR < 0.0 ) {
-                                    psiL = _psi[j][k];
-                                    psiR = _psi[_neighbors[j][l]][k];
+                                    psiL = _sol[j][k];
+                                    psiR = _sol[_neighbors[j][l]][k];
                                 }
                                 // flux evaluation
                                 psiNew[j][k] += _g->Flux( _quadPoints[k], psiL, psiR, _normals[j][l] );
@@ -102,15 +110,15 @@ void SNSolver::Solve() {
                             // higher order solver
                             case 3: std::cout << "higher order is WIP" << std::endl; break;
                             // default: first order solver
-                            default: psiNew[j][k] += _g->Flux( _quadPoints[k], _psi[j][k], _psi[_neighbors[j][l]][k], _normals[j][l] );
+                            default: psiNew[j][k] += _g->Flux( _quadPoints[k], _sol[j][k], _sol[_neighbors[j][l]][k], _normals[j][l] );
                         }
                     }
                 }
                 // time update angular flux with numerical flux and total scattering cross section
-                psiNew[j][k] = _psi[j][k] - ( _dE / _areas[j] ) * psiNew[j][k] - _dE * _sigmaT[n][j] * _psi[j][k];
+                psiNew[j][k] = _sol[j][k] - ( _dE / _areas[j] ) * psiNew[j][k] - _dE * _sigmaT[n][j] * _sol[j][k];
             }
             // compute scattering effects
-            psiNew[j] += _dE * _sigmaS[n][j] * _scatteringKernel * _psi[j];    // multiply scattering matrix with psi
+            psiNew[j] += _dE * _sigmaS[n][j] * _scatteringKernel * _sol[j];    // multiply scattering matrix with psi
 
             // TODO: figure out a more elegant way
             // add external source contribution
@@ -127,9 +135,9 @@ void SNSolver::Solve() {
                     psiNew[j] += _dE * _Q[n][j];
             }
         }
-        _psi = psiNew;
+        _sol = psiNew;
         for( unsigned i = 0; i < _nCells; ++i ) {
-            fluxNew[i]       = dot( _psi[i], _weights );
+            fluxNew[i]       = dot( _sol[i], _weights );
             _solverOutput[i] = fluxNew[i];
         }
         // Save( n );
@@ -143,7 +151,7 @@ void SNSolver::Save() const {
     std::vector<std::string> fieldNames{ "flux" };
     std::vector<double> flux( _nCells, 0.0 );
     for( unsigned i = 0; i < _nCells; ++i ) {
-        flux[i] = dot( _psi[i], _weights );
+        flux[i] = dot( _sol[i], _weights );
     }
     std::vector<std::vector<double>> scalarField( 1, flux );
     std::vector<std::vector<std::vector<double>>> results{ scalarField };

code/src/solvers/snsolverMPI.cpp

@@ -69,7 +69,7 @@ void SNSolverMPI::Save() const {
     std::vector<std::string> fieldNames{ "flux" };
     std::vector<double> flux( _nCells, 0.0 );
     for( unsigned i = 0; i < _nCells; ++i ) {
-        flux[i] = dot( _psi[i], _weights );
+        flux[i] = dot( _sol[i], _weights );
     }
     std::vector<std::vector<double>> scalarField( 1, flux );
     std::vector<std::vector<std::vector<double>>> results{ scalarField };

code/src/solvers/solverbase.cpp

 #include "solvers/solverbase.h"
 #include "fluxes/numericalflux.h"
 #include "io.h"
-#include "kernels/scatteringkernelbase.h"
 #include "mesh.h"
 #include "problems/problembase.h"
 #include "quadratures/quadraturebase.h"
@@ -27,8 +26,7 @@ Solver::Solver( Config* settings ) : _settings( settings ) {
     _weights             = quad->GetWeights();
     _nq                  = quad->GetNq();
     _settings->SetNQuadPoints( _nq );
-    ScatteringKernel* k = ScatteringKernel::CreateScatteringKernel( settings->GetKernelName(), quad );
-    _scatteringKernel   = k->GetScatteringKernel();
+    delete quad;
 
     // set time step
     _dE        = ComputeTimeStep( settings->GetCFL() );
@@ -37,7 +35,7 @@ Solver::Solver( Config* settings ) : _settings( settings ) {
 
     // setup problem
     _problem = ProblemBase::Create( _settings, _mesh );
-    _psi     = _problem->SetupIC();
+    _sol     = _problem->SetupIC();
     _s       = _problem->GetStoppingPower( _energies );
     _sigmaT  = _problem->GetTotalXS( _energies );
     _sigmaS  = _problem->GetScatteringXS( _energies );
@@ -84,7 +82,7 @@ void Solver::Save() const {
     flux.resize( _nCells );
 
     for( unsigned i = 0; i < _nCells; ++i ) {
-        flux[i] = _psi[i][0];
+        flux[i] = _sol[i][0];
     }
     std::vector<std::vector<double>> scalarField( 1, flux );
     std::vector<std::vector<std::vector<double>>> results{ scalarField };

code/tests/test_cases.cpp

@@ -2,6 +2,7 @@
 #include <vtkUnstructuredGridReader.h>
 
 #include "catch.hpp"
+#include "io.h"    // THIS SUCKS, ITS ONLY THERE BC I DONT KNOW WHAT IS EREALLY NEEDED FOR vtk
 #include "mesh.h"
 #include "settings/config.h"
 #include "solvers/solverbase.h"