added problem class; crosssections are now cell specific; added checkerboard...

added problem class; crosssections are now cell specific; added checkerboard testcase; moved mpi sn to separate file; scattering implementation missing

code/include/problems/checkerboard.h

+#ifndef CHECKERBOARD_H
+#define CHECKERBOARD_H
+
+#include "problembase.h"
+
+class Checkerboard : public ProblemBase
+{
+  private:
+    std::vector<Matrix> _scatteringXS;
+    std::vector<double> _totalXS;
+
+    Checkerboard() = delete;
+
+    bool isAbsorption( const Vector& pos );
+
+  public:
+    Checkerboard( Config* settings, Mesh* mesh );
+    virtual ~Checkerboard();
+
+    virtual std::vector<Matrix> GetScatteringXS( const double energy );
+    virtual std::vector<double> GetTotalXS( const double energy );
+    virtual std::vector<double> GetStoppingPower( const std::vector<double>& energies );
+    virtual VectorVector SetupIC();
+};
+
+#endif    // CHECKERBOARD_H

code/include/problems/electronrt.h

+#ifndef ELECTRONRT_H
+#define ELECTRONRT_H
+
+#include "problembase.h"
+
+class ElectronRT : public ProblemBase
+{
+  private:
+    /**
+     * @brief LoadXSH2O loads (total) scattering cross sections for water from file and saves them in _xsH2O and _totalxsH2O
+     * @param fileSigmaS is name of scattering cross section file
+     * @param fileSigmaT is name of total cross section file
+     */
+    void LoadXSH20( std::string fileSigmaS, std::string fileSigmaT );
+
+    ElectronRT() = delete;
+
+  public:
+    ElectronRT( Config* settings, Mesh* mesh );
+    virtual ~ElectronRT();
+
+    virtual std::vector<Matrix> GetScatteringXS( const double energy );
+    virtual std::vector<double> GetTotalXS( const double energy );
+    virtual std::vector<double> GetStoppingPower( const std::vector<double>& energies );
+    virtual VectorVector SetupIC();
+};
+
+#endif    // ELECTRONRT_H

code/include/problems/linesource.h

+#ifndef LINESOURCE_H
+#define LINESOURCE_H
+
+#include "problembase.h"
+
+class LineSource : public ProblemBase
+{
+  private:
+    LineSource() = delete;
+
+  public:
+    LineSource( Config* settings, Mesh* mesh );
+    virtual ~LineSource();
+
+    virtual std::vector<Matrix> GetScatteringXS( const double energy );
+    virtual std::vector<double> GetTotalXS( const double energy );
+    virtual std::vector<double> GetStoppingPower( const std::vector<double>& energies );
+    virtual VectorVector SetupIC();
+};
+
+#endif    // LINESOURCE_H

code/include/problems/problembase.h

+#ifndef PROBLEMBASE_H
+#define PROBLEMBASE_H
+
+#include "mesh.h"
+#include "physics.h"
+#include "settings/config.h"
+#include "typedef.h"
+
+class ProblemBase
+{
+
+  protected:
+    Config* _settings;
+    Mesh* _mesh;
+    Physics* _physics;
+
+    std::vector<double> _totalXS;
+    Matrix _scatteringXS;
+    std::vector<double> _density;
+    std::vector<double> _stoppingPower;
+
+    ProblemBase() = delete;
+
+  public:
+    /**
+     * @brief GetScatteringXS gives back vector of vectors of scattering cross sections for materials defined by density and energies in vector energy
+     * @param energies is vector with energies
+     * @param density is vector with patient densities (at different spatial cells)
+     * @param Omega are scattering angles
+     */
+    virtual std::vector<Matrix> GetScatteringXS( const double energy ) = 0;
+
+    /**
+     * @brief GetTotalXS gives back vector of vectors of total cross sections for materials defined by density and energies in vector energy
+     * @param energies is vector with energies
+     * @param density is vector with patient densities (at different spatial cells)
+     */
+    virtual std::vector<double> GetTotalXS( const double energy ) = 0;
+
+    /**
+     * @brief GetStoppingPower gives back vector of vectors of stopping powers for materials defined by density and energies in vector energy
+     * @param energies is vector with energies
+     * @param density is vector with patient densities (at different spatial cells)
+     * @param sH2O is vector of stopping powers in water
+     */
+    virtual std::vector<double> GetStoppingPower( const std::vector<double>& energies ) = 0;
+
+    /**
+     * @brief Setup the initial condition for the flux psi
+     */
+    virtual VectorVector SetupIC() = 0;
+
+    /**
+     * @brief Physics constructor
+     * @param settings stores all needed user information
+     */
+    ProblemBase( Config* settings, Mesh* mesh );
+    virtual ~ProblemBase();
+
+    /**
+     * @brief Create constructor
+     * @param settings stores all needed information
+     * @return pointer to Physics
+     */
+    static ProblemBase* Create( Config* settings, Mesh* mesh );
+};
+
+#endif

code/include/settings/config.h

@@ -44,9 +44,16 @@ class Config
     // Quadrature
     QUAD_NAME _quadName;       /*!< \brief Quadrature Name*/
     unsigned short _quadOrder; /*!< \brief Quadrature Order*/
+    unsigned _nQuadPoints;
+
+    // Mesh
+    unsigned _nCells;
+
     // Solver
     double _CFL;  /*!< \brief CFL Number for Solver*/
     double _tEnd; /*!< \brief Final Time for Simulation */
+    PROBLEM_NAME _problemName;
+
     // Boundary Conditions
     /*!< \brief List of all Pairs (marker, BOUNDARY_TYPE), e.g. (farfield,DIRICHLET).
          Each Boundary Conditions must have an entry in enum BOUNDARY_TYPE*/
@@ -176,10 +183,17 @@ class Config
     // Quadrature Structure
     QUAD_NAME inline GetQuadName() const { return _quadName; }
     unsigned short inline GetQuadOrder() const { return _quadOrder; }
+    void SetNQuadPoints( unsigned nq ) { _nQuadPoints = nq; }
+    unsigned GetNQuadPoints() { return _nQuadPoints; }
+
+    // Mesh Structure
+    void SetNCells( unsigned nCells ) { _nCells = nCells; }
+    unsigned GetNCells() { return _nCells; }
 
     // Solver Structure
     double inline GetCFL() const { return _CFL; }
     double inline GetTEnd() const { return _tEnd; }
+    PROBLEM_NAME inline GetProblemName() const { return _problemName; }
 
     // Boundary Conditions
     BOUNDARY_TYPE GetBoundaryType( std::string nameMarker ) const; /*! \brief Get Boundary Type of given marker */

code/include/settings/globalconstants.h

@@ -37,4 +37,9 @@ inline std::map<std::string, QUAD_NAME> Quadrature_Map{ { "MONTE_CARLO", QUAD_Mo
                                                         { "LEBEDEV", QUAD_Lebedev },
                                                         { "LDFESA", QUAD_LDFESA } };
 
+enum PROBLEM_NAME { PROBLEM_LineSource, PROBLEM_Checkerboard, PROBLEM_ElectronRT };
+
+inline std::map<std::string, PROBLEM_NAME> Problem_Map{
+    { "LINESOURCE", PROBLEM_LineSource }, { "CHECKERBOARD", PROBLEM_Checkerboard }, { "ELECTRONRT", PROBLEM_ElectronRT } };
+
 #endif    // GLOBAL_CONSTANTS_H

code/include/snsolver.h

@@ -24,7 +24,6 @@ class SNSolver : public Solver
      * @brief Output solution to VTK file
      */
     virtual void Save() const;
-    void SolveMPI();    // can be deleated later
 };
 
 #endif    // SNSOLVER_H

code/include/snsolverMPI.h

+#ifndef SNSOLVERMPI_H
+#define SNSOLVERMPI_H
+
+#include <mpi.h>
+
+#include "solver.h"
+#include "typedef.h"
+
+#include "settings/config.h"
+
+class SNSolverMPI : public Solver
+{
+  public:
+    /**
+     * @brief SNSolver constructor
+     * @param settings stores all needed information
+     */
+    SNSolverMPI( Config* settings );
+    /**
+     * @brief Solve functions runs main time loop
+     */
+    virtual void Solve();
+    /**
+     * @brief Output solution to VTK file
+     */
+    virtual void Save() const;
+};
+
+#endif    // SNSOLVERMPI_H

code/include/solver.h

@@ -6,6 +6,7 @@
 // include Matrix, Vector definitions
 #include "io.h"
 #include "numericalflux.h"
+#include "problems/problembase.h"
 #include "settings/config.h"
 #include "typedef.h"
 
@@ -14,16 +15,15 @@ class Solver
   protected:
     unsigned _nq;                                     // number of quadrature points
     unsigned _nCells;                                 // number of spatial cells
-    unsigned _nTimeSteps;                             // number of time steps, number of nodal energy values for CSD
-    double _dt;                                       // time step size
+    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> _sH20;                        // stopping power H2O, dim(_sH20) = _nTimeSteps
-    std::vector<double> _sigmaTH20;                   // total cross section, dim(_sigmaTH20) = _nTimeSteps
-    std::vector<Matrix> _sigmaSH20;                   // scattering cross section, dim(_sigmaSH20) = (_nTimeSteps,_nq,_nq)
+    std::vector<double> _s;                           // stopping power, dim(_s) = _nTimeSteps
     VectorVector _quadPoints;                         // quadrature points, dim(_quadPoints) = (_nTimeSteps,spatialDim)
     Vector _weights;                                  // quadrature weights, dim(_weights) = (_NCells)
     std::vector<bool> _boundaryCells;                 // boundary type for all cells, dim(_boundary) = (_NCells)
@@ -32,30 +32,7 @@ class Solver
     NumericalFlux* _g;    // numerical flux function
     Mesh* _mesh;          // mesh object for writing out information
     Config* _settings;
-
-    /**
-     * @brief LoadPatientDensity loads density of patient from MRT/CT scan and saves it in _density
-     * @param fileName is name of patient file
-     */
-    void LoadPatientDensity( std::string fileName );
-
-    /**
-     * @brief LoadStoppingPower loads stopping power of H20 and saves it in _sH20
-     * @param fileName is name of stopping power file
-     */
-    void LoadStoppingPower( std::string fileName );
-
-    /**
-     * @brief LoadSigmaS loads scattering cross section of H20 and saves it in _sigmaSH20
-     * @param fileName is name of scattering cross section file
-     */
-    void LoadSigmaS( std::string fileName );
-
-    /**
-     * @brief LoadSigmaT loads total cross section of H20 and saves it in _sigmaTH20
-     * @param fileName is name of scattering cross section file
-     */
-    void LoadSigmaT( std::string fileName );
+    ProblemBase* _problem;
 
     /**
      * @brief ComputeTimeStep calculates the maximal stable time step
@@ -63,11 +40,6 @@ class Solver
      */
     double ComputeTimeStep( double cfl ) const;
 
-    /**
-     * @brief SetupIC writes intial condition onto _psi
-     */
-    void SetupIC();
-
   public:
     /**
      * @brief Solver constructor

code/input/checkerboard.cfg

+% ---- File specifications ----
+OUTPUT_DIR = ../result
+OUTPUT_FILE = checkerboard
+LOG_DIR = ../result/logs
+MESH_FILE = checkerboard.su2
+
+% ---- Solver specifications ----
+CFL_NUMBER = 0.5
+TIME_FINAL = 2.0
+PROBLEM = CHECKERBOARD
+
+% ---- Boundary Conditions ----
+BC_DIRICHLET = ( void )
+QUAD_TYPE = LEBEDEV
+QUAD_ORDER = 19

code/input/checkerboard.py

@@ -12,24 +12,20 @@ def add_block(x0,y0,length,char_length,geom):
     ])
     return geom.add_polygon(coords, char_length)
 
-char_length = 0.05
+char_length = 0.06
 geom = pg.opencascade.Geometry()
-domain = add_block(0, 0, 1.2, char_length, geom)
-xpos = ypos = [0.1, 0.3, 0.5, 0.7, 0.9]
+domain = add_block(0, 0, 6, char_length, geom)
+xpos = ypos = [0.5, 1.5, 2.5, 3.5, 4.5]
 pattern = list(itertools.product(xpos, ypos))[::2]
 pattern.pop(len(pattern)-2)
 boxes = [domain]
 for pos in pattern:
-    boxes.append(add_block(pos[0], pos[1], 0.2, char_length, geom))
+    boxes.append(add_block(pos[0], pos[1], 1, char_length, geom))
 geom.boolean_fragments(boxes,[])
 geom.add_physical(domain.lines, label="void")
 
-#doesn't work on my machine for some reason, but should work in general
-#pg.generate_mesh(geom, msh_filename='checkerboard.msh', dim=2, extra_gmsh_arguments=['-save_all', '-optimize'])
-
-#workaround
 mesh_code = geom.get_code()
 with open("checkerboard.geo","w") as mesh_file:
     mesh_file.write(mesh_code)
-os.system('gmsh checkerboard.geo -2 -format su2 -save_all -optimize')
+os.system('gmsh checkerboard.geo -2 -format su2 -save_all')
 os.system('rm checkerboard.geo')

code/input/example.toml

-[io]
-outputDir = "../result"
-outputFile = "example"
-logDir = "../result/logs"
-meshFile = "checkerboard.su2"
-
-[solver]
-CFL = 0.9
-tEnd = 0.3
-boundaryConditions = [["void", "dirichlet"]]
-quadType = "MONTE_CARLO"
-quadOrder = 42 

code/input/linesource.cfg

+% ---- File specifications ----
+OUTPUT_DIR = ../result
+OUTPUT_FILE = linesource
+LOG_DIR = ../result/logs
+MESH_FILE = linesource.su2
+
+% ---- Solver specifications ----
+CFL_NUMBER = 0.5
+TIME_FINAL = 0.3
+PROBLEM = LINESOURCE
+
+% ---- Boundary Conditions ----
+BC_DIRICHLET = ( void )
+QUAD_TYPE = LEBEDEV
+QUAD_ORDER = 35

code/input/linesource.py

+import pygmsh as pg
+import numpy as np
+import itertools
+import os
+
+def add_block(x0,y0,length,char_length,geom):
+    coords = np.array([
+        [x0, y0, 0.0],
+        [x0+length, y0, 0.0],
+        [x0+length, y0+length, 0.0],
+        [x0, y0+length, 0.0]
+    ])
+    return geom.add_polygon(coords, char_length)
+
+char_length = 0.01
+geom = pg.opencascade.Geometry()
+domain = add_block(0, 0, 1, char_length, geom)
+geom.add_raw_code('psource = newp;\nPoint(psource) = {0.5, 0.5, 0.0, '+str(char_length)+'};\nPoint{psource} In Surface{'+domain.id+'};')
+geom.add_physical(domain.lines, label="void")
+
+mesh_code = geom.get_code()
+with open("linesource.geo","w") as mesh_file:
+    mesh_file.write(mesh_code)
+os.system('gmsh linesource.geo -2 -format su2 -save_all')
+os.system('rm linesource.geo')

code/src/problems/checkerboard.cpp

+#include "problems/checkerboard.h"
+
+Checkerboard::Checkerboard( Config* settings, Mesh* mesh ) : ProblemBase( settings, mesh ) {
+    _physics = nullptr;
+    _scatteringXS.resize( _mesh->GetNumCells(), Matrix( _settings->GetNQuadPoints(), _settings->GetNQuadPoints(), 0.0 ) );    // @TODO
+    _totalXS.resize( _mesh->GetNumCells(), 1.0 );
+    auto cellMids = _mesh->GetCellMidPoints();
+    for( unsigned j = 0; j < cellMids.size(); ++j ) {
+        if( isAbsorption( cellMids[j] ) ) {
+            _scatteringXS[j] = 0.0;
+            _totalXS[j]      = 10.0;
+        }
+    }
+}
+
+Checkerboard::~Checkerboard() {}
+
+std::vector<Matrix> Checkerboard::GetScatteringXS( const double energy ) { return _scatteringXS; }
+
+std::vector<double> Checkerboard::GetTotalXS( const double energy ) { return _totalXS; }
+
+std::vector<double> Checkerboard::GetStoppingPower( const std::vector<double>& energies ) {
+    // @TODO
+    return std::vector<double>( energies.size(), 0.0 );
+}
+
+VectorVector Checkerboard::SetupIC() {
+    VectorVector psi = std::vector( _mesh->GetNumCells(), Vector( _settings->GetNQuadPoints(), 1e-7 ) );
+    auto cellMids    = _mesh->GetCellMidPoints();
+    for( unsigned j = 0; j < cellMids.size(); ++j ) {
+        if( cellMids[j][0] >= 2.5 && cellMids[j][0] <= 3.5 && cellMids[j][1] >= 2.5 && cellMids[j][1] <= 3.5 ) {
+            psi[j] = 1.0;
+        }
+    }
+    return psi;
+}
+
+bool Checkerboard::isAbsorption( const Vector& pos ) {
+    std::vector<double> lbounds{ 0.5, 1.5, 2.5, 3.5, 4.5 };
+    std::vector<double> ubounds{ 1.5, 2.5, 3.5, 4.5, 5.5 };
+    for( unsigned k = 0; k < lbounds.size(); ++k ) {
+        for( unsigned l = 0; l < lbounds.size(); ++l ) {
+            if( ( l + k ) % 2 == 1 || ( k == 2 && l == 2 ) || ( k == 4 && l == 2 ) ) continue;
+            if( pos[0] >= lbounds[k] && pos[0] <= ubounds[k] && pos[1] >= lbounds[l] && pos[1] <= ubounds[l] ) {
+                return true;
+            }
+        }
+    }
+    return false;
+}

code/src/problems/electronrt.cpp

+#include "problems/electronrt.h"
+
+ElectronRT::ElectronRT( Config* settings, Mesh* mesh ) : ProblemBase( settings, mesh ) {
+    // @TODO
+    _physics = nullptr;
+}
+
+ElectronRT::~ElectronRT() {}
+
+std::vector<Matrix> ElectronRT::GetScatteringXS( const double energy ) {
+    // @TODO
+    return std::vector<Matrix>( _mesh->GetNumCells(), Matrix( _settings->GetNQuadPoints(), _settings->GetNQuadPoints(), 0.0 ) );
+}
+
+std::vector<double> ElectronRT::GetTotalXS( const double energy ) {
+    // @TODO
+    return std::vector<double>( _mesh->GetNumCells(), 0.0 );
+}
+
+std::vector<double> ElectronRT::GetStoppingPower( const std::vector<double>& energies ) {
+    // @TODO
+    return std::vector<double>( energies.size(), 0.0 );
+}
+
+VectorVector ElectronRT::SetupIC() {
+    // @TODO
+    return VectorVector( _mesh->GetNumCells(), Vector( _settings->GetNQuadPoints(), 1e-7 ) );
+}
+
+void ElectronRT::LoadXSH20( std::string fileSigmaS, std::string fileSigmaT ) {
+    // @TODO
+}

code/src/problems/linesource.cpp

+#include "problems/linesource.h"
+
+LineSource::LineSource( Config* settings, Mesh* mesh ) : ProblemBase( settings, mesh ) { _physics = nullptr; }
+
+LineSource::~LineSource(){};
+
+std::vector<Matrix> LineSource::GetScatteringXS( const double energy ) {
+    return std::vector<Matrix>( _mesh->GetNumCells(), Matrix( _settings->GetNQuadPoints(), _settings->GetNQuadPoints(), 0.0 ) );
+}
+
+std::vector<double> LineSource::GetTotalXS( const double energy ) { return std::vector<double>( _mesh->GetNumCells(), 0.0 ); }
+
+std::vector<double> LineSource::GetStoppingPower( const std::vector<double>& energies ) {
+    // @TODO
+    return std::vector<double>( energies.size(), 0.0 );
+}
+
+VectorVector LineSource::SetupIC() {
+    VectorVector psi = std::vector( _mesh->GetNumCells(), Vector( _settings->GetNQuadPoints(), 1e-7 ) );
+    auto cellMids    = _mesh->GetCellMidPoints();
+    Vector midPoint( 2, 0.5 );
+    for( unsigned j = 0; j < cellMids.size(); ++j ) {
+        if( norm( cellMids[j] - midPoint ) <= 0.1 ) {
+            psi[j] = 1.0;
+        }
+    }
+    return psi;
+}

code/src/problems/problembase.cpp

+#include "problems/problembase.h"
+#include "problems/checkerboard.h"
+#include "problems/electronrt.h"
+#include "problems/linesource.h"
+
+ProblemBase::ProblemBase( Config* settings, Mesh* mesh ) : _settings( settings ), _mesh( mesh ) {}
+
+ProblemBase::~ProblemBase() {}
+
+ProblemBase* ProblemBase::Create( Config* settings, Mesh* mesh ) {
+    auto name = settings->GetProblemName();
+    switch( name ) {
+        case PROBLEM_LineSource: return new LineSource( settings, mesh );
+        case PROBLEM_Checkerboard: return new Checkerboard( settings, mesh );
+        case PROBLEM_ElectronRT: return new ElectronRT( settings, mesh );
+        default: return new ElectronRT( settings, mesh );    // Use RadioTherapy as dummy
+    }
+}

code/src/settings/config.cpp

@@ -198,6 +198,8 @@ void Config::SetConfigOptions() {
     AddDoubleOption( "CFL_NUMBER", _CFL, 1.0 );
     /*!\brief TIME_FINAL \n DESCRIPTION: Final time for simulation \n DEFAULT 1.0 \ingroup Config.*/
     AddDoubleOption( "TIME_FINAL", _tEnd, 1.0 );
+    /*!\brief Problem \n DESCRIPTION: Type of problem setting \n DEFAULT PROBLEM_ElectronRT \ingroup Config.*/
+    AddEnumOption( "PROBLEM", _problemName, Problem_Map, PROBLEM_ElectronRT );
 
     // Mesh related options
     // Boundary Markers