waterphantom first functionality added, CFL fixed, small changes csd solver

code/include/problems/waterphantom.h

+#ifndef WATERPHANTOM_H
+#define WATERPHANTOM_H
+
+#include "electronrt.h"
+
+class WaterPhantom : public ElectronRT
+{
+  private:
+    WaterPhantom() = delete;
+
+  public:
+    WaterPhantom( Config* settings, Mesh* mesh );
+    virtual ~WaterPhantom();
+
+    virtual std::vector<VectorVector> GetExternalSource( const std::vector<double>& energies );
+    virtual std::vector<double> GetStoppingPower( const std::vector<double>& energies );
+    virtual VectorVector SetupIC();
+    std::vector<double> GetDensity( const VectorVector& cellMidPoints );
+};
+
+#endif    // WATERPHANTOM_H

code/include/settings/globalconstants.h

@@ -38,10 +38,12 @@ inline std::map<std::string, QUAD_NAME> Quadrature_Map{ { "MONTE_CARLO", QUAD_Mo
                                                         { "LDFESA", QUAD_LDFESA } };
 
 // Problem name
-enum PROBLEM_NAME { PROBLEM_LineSource, PROBLEM_Checkerboard, PROBLEM_ElectronRT };
+enum PROBLEM_NAME { PROBLEM_LineSource, PROBLEM_Checkerboard, PROBLEM_ElectronRT, PROBLEM_WaterPhantom };
 
-inline std::map<std::string, PROBLEM_NAME> Problem_Map{
-    { "LINESOURCE", PROBLEM_LineSource }, { "CHECKERBOARD", PROBLEM_Checkerboard }, { "ELECTRONRT", PROBLEM_ElectronRT } };
+inline std::map<std::string, PROBLEM_NAME> Problem_Map{ { "LINESOURCE", PROBLEM_LineSource },
+                                                        { "CHECKERBOARD", PROBLEM_Checkerboard },
+                                                        { "ELECTRONRT", PROBLEM_ElectronRT },
+                                                        { "WATERPHANTOM", PROBLEM_WaterPhantom } };
 
 // Kernel name
 enum KERNEL_NAME { KERNEL_Isotropic };

code/include/solvers/csdsnsolver.h

@@ -8,6 +8,8 @@
 class CSDSNSolver : public Solver
 {
   private:
+    std::vector<double> _dose;
+
   public:
     /**
      * @brief CSDSNSolver constructor

code/src/problems/linesource.cpp

@@ -20,7 +20,7 @@ std::vector<VectorVector> LineSource_SN::GetExternalSource( const std::vector<do
 
 std::vector<double> LineSource_SN::GetStoppingPower( const std::vector<double>& energies ) {
     // @TODO
-    return std::vector<double>( energies.size(), 0.0 );
+    return std::vector<double>( energies.size(), 1.0 );
 }
 
 VectorVector LineSource_SN::SetupIC() {

code/src/problems/problembase.cpp

@@ -2,6 +2,7 @@
 #include "problems/checkerboard.h"
 #include "problems/electronrt.h"
 #include "problems/linesource.h"
+#include "problems/waterphantom.h"
 
 ProblemBase::ProblemBase( Config* settings, Mesh* mesh ) : _settings( settings ), _mesh( mesh ) {}
 
@@ -18,6 +19,7 @@ ProblemBase* ProblemBase::Create( Config* settings, Mesh* mesh ) {
         }
         case PROBLEM_Checkerboard: return new Checkerboard( settings, mesh );
         case PROBLEM_ElectronRT: return new ElectronRT( settings, mesh );
+        case PROBLEM_WaterPhantom: return new WaterPhantom( settings, mesh );
         default: return new ElectronRT( settings, mesh );    // Use RadioTherapy as dummy
     }
 }

code/src/problems/waterphantom.cpp

+#include "problems/waterphantom.h"
+
+WaterPhantom::WaterPhantom( Config* settings, Mesh* mesh ) : ElectronRT( settings, mesh ) {
+    // @TODO get pointer to correct physics class
+    _physics = nullptr;
+}
+
+WaterPhantom::~WaterPhantom() {}
+
+std::vector<VectorVector> WaterPhantom::GetExternalSource( const std::vector<double>& energies ) {
+    return std::vector<VectorVector>( energies.size(), std::vector<Vector>( _mesh->GetNumCells(), Vector( _settings->GetNQuadPoints(), 0.0 ) ) );
+}
+
+std::vector<double> WaterPhantom::GetStoppingPower( const std::vector<double>& energies ) {
+    // @TODO get correct stopping power
+    return std::vector<double>( energies.size(), 1.0 );
+}
+
+VectorVector WaterPhantom::SetupIC() {
+    VectorVector psi( _mesh->GetNumCells(), Vector( _settings->GetNQuadPoints(), 1e-10 ) );
+    auto cellMids = _mesh->GetCellMidPoints();
+    double t      = 3.2e-4;    // pseudo time for gaussian smoothing
+    double s      = 0.1;
+    for( unsigned j = 0; j < cellMids.size(); ++j ) {
+        double x = cellMids[j][0];
+        psi[j]   = 1.0 / ( s * sqrt( 2 * M_PI ) ) * std::exp( -x * x / ( 2 * s * s ) );
+    }
+    return psi;
+}
+
+std::vector<double> WaterPhantom::GetDensity( const VectorVector& cellMidPoints ) { return std::vector<double>( cellMidPoints.size(), 1.0 ); }

code/src/solvers/csdsnsolver.cpp

 #include "solvers/csdsnsolver.h"
 
-CSDSNSolver::CSDSNSolver( Config* settings ) : Solver( settings ) {}
+CSDSNSolver::CSDSNSolver( Config* settings ) : Solver( settings ) { _dose = std::vector<double>( _settings->GetNCells(), 0.0 ); }
 
 void CSDSNSolver::Solve() {
     auto log = spdlog::get( "event" );
@@ -10,6 +10,9 @@ void CSDSNSolver::Solve() {
     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 );
+    }
     int rank;
     MPI_Comm_rank( MPI_COMM_WORLD, &rank );
     if( rank == 0 ) log->info( "{:10}   {:10}", "E", "dFlux" );
@@ -35,7 +38,7 @@ void CSDSNSolver::Solve() {
 
     // loop over energies (pseudo-time)
     for( unsigned n = 1; n < _nEnergies; ++n ) {
-        _dE = _energies[n] - _energies[n - 1];
+        _dE = fabs( _energies[n] - _energies[n - 1] );    // is the sign correct here?
         // loop over all spatial cells
         for( unsigned j = 0; j < _nCells; ++j ) {
             if( _boundaryCells[j] == BOUNDARY_TYPE::DIRICHLET ) continue;
@@ -60,15 +63,15 @@ void CSDSNSolver::Solve() {
             // add external source contribution
             if( _Q.size() == 1u ) {            // constant source for all energies
                 if( _Q[0][j].size() == 1u )    // isotropic source
-                    psiNew[j] += _dE * _Q[0][j][0] * _s[n];
+                    psiNew[j] += _dE * _Q[0][j][0] * _s[_nEnergies - n - 1];
                 else
-                    psiNew[j] += _dE * _Q[0][j] * _s[n];
+                    psiNew[j] += _dE * _Q[0][j] * _s[_nEnergies - n - 1];
             }
             else {
                 if( _Q[0][j].size() == 1u )    // isotropic source
-                    psiNew[j] += _dE * _Q[n][j][0] * _s[n];
+                    psiNew[j] += _dE * _Q[n][j][0] * _s[_nEnergies - n - 1];
                 else
-                    psiNew[j] += _dE * _Q[n][j] * _s[n];
+                    psiNew[j] += _dE * _Q[n][j] * _s[_nEnergies - n - 1];
             }
         }
         _psi = psiNew;
@@ -76,13 +79,14 @@ void CSDSNSolver::Solve() {
         // do backsubstitution from psiTildeHat to psi (cf. Dissertation Kerstion Kuepper, Eq. 1.23)
         for( unsigned j = 0; j < _nCells; ++j ) {
             for( unsigned k = 0; k < _nq; ++k ) {
-                psiNew[j][k] = _psi[j][k] * _density[j] * _s[0];    // note that _s[0] is stopping power at highest energy
+                psiNew[j][k] = _psi[j][k] * _density[j] * _s[_nEnergies - n - 1];    // note that _s[0] is stopping power at lowest energy
             }
         }
 
-        for( unsigned i = 0; i < _nCells; ++i ) {
-            fluxNew[i]       = dot( psiNew[i], _weights );
-            _solverOutput[i] = fluxNew[i];
+        for( unsigned j = 0; j < _nCells; ++j ) {
+            fluxNew[j] = dot( psiNew[j], _weights );
+            _dose[j] += 0.5 * _dE * ( fluxNew[j] + fluxOld[j] ) / _density[j];    // update dose with trapezoidal rule
+            _solverOutput[j] = fluxNew[j];
         }
         Save( n );
         dFlux   = blaze::l2Norm( fluxNew - fluxOld );
@@ -92,12 +96,8 @@ void CSDSNSolver::Solve() {
 }
 
 void CSDSNSolver::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 );
-    }
-    std::vector<std::vector<double>> scalarField( 1, flux );
+    std::vector<std::string> fieldNames{ "dose" };
+    std::vector<std::vector<double>> scalarField( 1, _dose );
     std::vector<std::vector<std::vector<double>>> results{ scalarField };
     ExportVTK( _settings->GetOutputFile(), results, fieldNames, _mesh );
 }

code/src/solvers/solverbase.cpp

@@ -52,14 +52,14 @@ Solver::Solver( Config* settings ) : _settings( settings ) {
 }
 
 double Solver::ComputeTimeStep( double cfl ) const {
-    double maxEdge = -1.0;
+    double minEdge = 100000.0;
     for( unsigned j = 0; j < _nCells; ++j ) {
         for( unsigned l = 0; l < _normals[j].size(); ++l ) {
             double currentEdge = _areas[j] / norm( _normals[j][l] );
-            if( currentEdge > maxEdge ) maxEdge = currentEdge;
+            if( currentEdge < minEdge ) minEdge = currentEdge;
         }
     }
-    return cfl * maxEdge;
+    return cfl * minEdge;
 }
 
 Solver* Solver::Create( Config* settings ) {