Merge branch 'CSDEquation' into 'sprint1'

csd solver added

See merge request !14

code/include/problems/electronrt.h

@@ -24,6 +24,7 @@ class ElectronRT : public ProblemBase
     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    // ELECTRONRT_H

code/include/problems/problembase.h

@@ -50,6 +50,12 @@ class ProblemBase
      */
     virtual std::vector<double> GetStoppingPower( const std::vector<double>& energies ) = 0;
 
+    /**
+     * @brief GetDensity gives back vector of densities for every spatial cell
+     * @param cellMidPoints is vector with cell mid points
+     */
+    virtual std::vector<double> GetDensity( const VectorVector& cellMidPoints );
+
     /**
      * @brief Setup the initial condition for the flux psi
      */

code/include/settings/config.h

@@ -61,6 +61,9 @@ class Config
      * to improve floating point accuracy */
     bool _cleanFluxMat;
 
+    /*!< @brief If true, continuous slowing down approximation will be used */
+    bool _csd;
+
     // 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*/
@@ -209,6 +212,7 @@ class Config
     PROBLEM_NAME inline GetProblemName() const { return _problemName; }
     SOLVER_NAME inline GetSolverName() const { return _solverName; }
     bool inline GetCleanFluxMat() const { return _cleanFluxMat; }
+    bool inline IsCSD() const { return _csd; }
     // Boundary Conditions
     BOUNDARY_TYPE GetBoundaryType( std::string nameMarker ) const; /*! @brief Get Boundary Type of given marker */
 

code/include/solvers/csdsnsolver.h

+#ifndef CSDSNSOLVER_H
+#define CSDSNSOLVER_H
+
+#include <mpi.h>
+
+#include "solvers/solverbase.h"
+
+class CSDSNSolver : public Solver
+{
+  private:
+  public:
+    /**
+     * @brief CSDSNSolver constructor
+     * @param settings stores all needed information
+     */
+    CSDSNSolver( Config* settings );
+    /**
+     * @brief Solve functions runs main time loop
+     */
+    virtual void Solve();
+    /**
+     * @brief Output solution to VTK file
+     */
+    virtual void Save() const;
+    virtual void Save( int currEnergy ) const;
+};
+
+#endif    // CSDSNSOLVER_H

code/input/example.cfg

@@ -13,12 +13,14 @@ LOG_DIR = ../result/logs
 % Mesh File
 MESH_FILE = linesource.su2
 %
-PROBLEM = LINESOURCE
+PROBLEM = ELECTRONRT
 %
 % ---- Solver specifications ----
 %
 % Solver type
-SOLVER = PN_SOLVER
+SOLVER = SN_SOLVER
+% Slowing down approximation
+CONTINUOUS_SLOWING_DOWN = YES
 % CFL number
 CFL_NUMBER = 0.8
 % Final time for simulation

code/src/problems/electronrt.cpp

@@ -24,7 +24,7 @@ std::vector<VectorVector> ElectronRT::GetExternalSource( const std::vector<doubl
 
 std::vector<double> ElectronRT::GetStoppingPower( const std::vector<double>& energies ) {
     // @TODO
-    return std::vector<double>( energies.size(), 0.0 );
+    return std::vector<double>( energies.size(), 1.0 );
 }
 
 VectorVector ElectronRT::SetupIC() {
@@ -35,3 +35,5 @@ VectorVector ElectronRT::SetupIC() {
 void ElectronRT::LoadXSH20( std::string fileSigmaS, std::string fileSigmaT ) {
     // @TODO
 }
+
+std::vector<double> ElectronRT::GetDensity( const VectorVector& cellMidPoints ) { return std::vector<double>( cellMidPoints.size(), 1.0 ); }

code/src/problems/problembase.cpp

@@ -21,3 +21,5 @@ ProblemBase* ProblemBase::Create( Config* settings, Mesh* mesh ) {
         default: return new ElectronRT( settings, mesh );    // Use RadioTherapy as dummy
     }
 }
+
+std::vector<double> ProblemBase::GetDensity( const VectorVector& cellMidPoints ) { return std::vector<double>( cellMidPoints.size(), 1.0 ); }

code/src/settings/config.cpp

@@ -208,6 +208,9 @@ void Config::SetConfigOptions() {
     /*! @brief CleanFluxMatrices \n DESCRIPTION:  If true, very low entries (10^-10 or smaller) of the flux matrices will be set to zero,
      * to improve floating point accuracy \n DEFAULT false \ingroup Config */
     AddBoolOption( "CLEAN_FLUX_MATRICES", _cleanFluxMat, false );
+    /*! @brief ContinuousSlowingDown \n DESCRIPTION: If true, the program uses the continuous slowing down approximation to treat energy dependent
+     * problems. \n DEFAULT false \ingroup Config */
+    AddBoolOption( "CONTINUOUS_SLOWING_DOWN", _csd, false );
 
     // Mesh related options
     // Boundary Markers

code/src/solvers/csdsnsolver.cpp

+#include "solvers/csdsnsolver.h"
+
+CSDSNSolver::CSDSNSolver( Config* settings ) : Solver( settings ) {}
+
+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;
+    double dFlux        = 1e10;
+    Vector fluxNew( _nCells, 0.0 );
+    Vector fluxOld( _nCells, 0.0 );
+    int rank;
+    MPI_Comm_rank( MPI_COMM_WORLD, &rank );
+    if( rank == 0 ) log->info( "{:10}   {:10}", "E", "dFlux" );
+
+    // 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
+        }
+    }
+
+    // store transformed energies ETilde instead of E in _energies vector (cf. Dissertation Kerstion Kuepper, Eq. 1.25)
+    double tmp = 0.0;
+    for( unsigned n = 0; n < _nEnergies; ++n ) {
+        tmp          = tmp + _dE / _s[n];
+        _energies[n] = tmp;
+    }
+
+    // store transformed energies ETildeTilde instead of ETilde in _energies vector (cf. Dissertation Kerstion Kuepper, Eq. 1.25)
+    for( unsigned n = 0; n < _nEnergies; ++n ) {
+        _energies[n] = _energies[_nEnergies - 1] - _energies[n];
+    }
+
+    // loop over energies (pseudo-time)
+    for( unsigned n = 1; n < _nEnergies; ++n ) {
+        _dE = _energies[n] - _energies[n - 1];
+        // loop over all spatial cells
+        for( unsigned j = 0; j < _nCells; ++j ) {
+            if( _boundaryCells[j] == BOUNDARY_TYPE::DIRICHLET ) continue;
+            // loop over all ordinates
+            for( unsigned k = 0; k < _nq; ++k ) {
+                psiNew[j][k] = 0.0;
+                // loop over all neighbor cells (edges) of cell j and compute numerical fluxes
+                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] / _density[j], _psi[j][k], _psi[j][k], _normals[j][l] );
+                    else
+                        psiNew[j][k] += _g->Flux( _quadPoints[k] / _density[j], _psi[j][k], _psi[_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];
+            }
+            // compute scattering effects
+            psiNew[j] += _dE * _sigmaS[n][j] * _scatteringKernel * _psi[j];    // multiply scattering matrix with psi
+
+            // TODO: figure out a more elegant way
+            // 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];
+                else
+                    psiNew[j] += _dE * _Q[0][j] * _s[n];
+            }
+            else {
+                if( _Q[0][j].size() == 1u )    // isotropic source
+                    psiNew[j] += _dE * _Q[n][j][0] * _s[n];
+                else
+                    psiNew[j] += _dE * _Q[n][j] * _s[n];
+            }
+        }
+        _psi = psiNew;
+
+        // 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
+            }
+        }
+
+        for( unsigned i = 0; i < _nCells; ++i ) {
+            fluxNew[i]       = dot( psiNew[i], _weights );
+            _solverOutput[i] = fluxNew[i];
+        }
+        Save( n );
+        dFlux   = blaze::l2Norm( fluxNew - fluxOld );
+        fluxOld = fluxNew;
+        if( rank == 0 ) log->info( "{:03.8f}   {:01.5e}", _energies[n], dFlux );
+    }
+}
+
+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::vector<std::vector<double>>> results{ scalarField };
+    ExportVTK( _settings->GetOutputFile(), results, fieldNames, _mesh );
+}
+
+void CSDSNSolver::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/solverbase.cpp

@@ -3,6 +3,7 @@
 #include "mesh.h"
 #include "quadratures/quadraturebase.h"
 #include "settings/globalconstants.h"
+#include "solvers/csdsnsolver.h"
 #include "solvers/pnsolver.h"
 #include "solvers/snsolver.h"
 
@@ -38,6 +39,7 @@ Solver::Solver( Config* settings ) : _settings( settings ) {
     _sigmaT  = _problem->GetTotalXS( _energies );
     _sigmaS  = _problem->GetScatteringXS( _energies );
     _Q       = _problem->GetExternalSource( _energies );
+    _density = _problem->GetDensity( _mesh->GetCellMidPoints() );    // double check: do we need to give the cell midpoints to fct?
 
     // setup numerical flux
     _g = NumericalFlux::Create( settings );
@@ -62,21 +64,20 @@ double Solver::ComputeTimeStep( double cfl ) const {
 
 Solver* Solver::Create( Config* settings ) {
     switch( settings->GetSolverName() ) {
-        case SN_SOLVER: return new SNSolver( settings );
-        case PN_SOLVER: return new PNSolver( settings );
+        case SN_SOLVER: {
+            if( settings->IsCSD() )
+                return new CSDSNSolver( settings );
+            else
+                return new SNSolver( settings );
+        }
+        case PN_SOLVER: {
+            if( settings->IsCSD() ) {
+                std::cerr << "[Solver]: Continuous slowing down option not available for PN Solver!" << std::endl;
+                return NULL;
+            }
+            else
+                return new PNSolver( settings );
+        }
         default: return new SNSolver( settings );
     }
 }
-
-void Solver::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() + "_" + std::to_string( _nEnergies ), results, fieldNames, _mesh );
-}