added scattering kernel and scatteringXS

code/include/problems/checkerboard.h

@@ -6,8 +6,8 @@
 class Checkerboard : public ProblemBase
 {
   private:
-    std::vector<Matrix> _scatteringXS;
-    std::vector<double> _totalXS;
+    Vector _scatteringXS;
+    Vector _totalXS;
 
     Checkerboard() = delete;
 
@@ -17,8 +17,8 @@ class Checkerboard : public ProblemBase
     Checkerboard( Config* settings, Mesh* mesh );
     virtual ~Checkerboard();
 
-    virtual std::vector<Matrix> GetScatteringXS( const double energy );
-    virtual std::vector<double> GetTotalXS( const double energy );
+    virtual VectorVector GetScatteringXS( const std::vector<double>& energies );
+    virtual VectorVector GetTotalXS( const std::vector<double>& energies );
     virtual std::vector<double> GetStoppingPower( const std::vector<double>& energies );
     virtual VectorVector SetupIC();
 };

code/include/problems/electronrt.h

@@ -19,8 +19,8 @@ class ElectronRT : public ProblemBase
     ElectronRT( Config* settings, Mesh* mesh );
     virtual ~ElectronRT();
 
-    virtual std::vector<Matrix> GetScatteringXS( const double energy );
-    virtual std::vector<double> GetTotalXS( const double energy );
+    virtual VectorVector GetScatteringXS( const std::vector<double>& energies );
+    virtual VectorVector GetTotalXS( const std::vector<double>& energies );
     virtual std::vector<double> GetStoppingPower( const std::vector<double>& energies );
     virtual VectorVector SetupIC();
 };

code/include/problems/linesource.h

@@ -12,8 +12,8 @@ class LineSource : public ProblemBase
     LineSource( Config* settings, Mesh* mesh );
     virtual ~LineSource();
 
-    virtual std::vector<Matrix> GetScatteringXS( const double energy );
-    virtual std::vector<double> GetTotalXS( const double energy );
+    virtual VectorVector GetScatteringXS( const std::vector<double>& energies );
+    virtual VectorVector GetTotalXS( const std::vector<double>& energies );
     virtual std::vector<double> GetStoppingPower( const std::vector<double>& energies );
     virtual VectorVector SetupIC();
 };

code/include/problems/problembase.h

@@ -14,8 +14,6 @@ class ProblemBase
     Mesh* _mesh;
     Physics* _physics;
 
-    std::vector<double> _totalXS;
-    Matrix _scatteringXS;
     std::vector<double> _density;
     std::vector<double> _stoppingPower;
 
@@ -24,24 +22,20 @@ class ProblemBase
   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
+     * @param energy is the energy the cross section is queried for
      */
-    virtual std::vector<Matrix> GetScatteringXS( const double energy ) = 0;
+    virtual VectorVector GetScatteringXS( const std::vector<double>& energies ) = 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 energy is the energy the cross section is queried for
      * @param density is vector with patient densities (at different spatial cells)
      */
-    virtual std::vector<double> GetTotalXS( const double energy ) = 0;
+    virtual VectorVector GetTotalXS( const std::vector<double>& energies ) = 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;
 

code/include/snsolver.h

@@ -10,6 +10,7 @@
 
 class SNSolver : public Solver
 {
+  private:
   public:
     /**
      * @brief SNSolver constructor

code/include/solver.h

@@ -24,6 +24,9 @@ class Solver
     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
+    Matrix _scatteringKernel;                         // scattering kernel for the quadrature
     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)
@@ -39,6 +42,7 @@ class Solver
      * @param cfl is cfl number
      */
     double ComputeTimeStep( double cfl ) const;
+    Matrix ComputeScatteringKernel() const;
 
   public:
     /**

code/input/checkerboard.cfg

@@ -12,4 +12,4 @@ PROBLEM = CHECKERBOARD
 % ---- Boundary Conditions ----
 BC_DIRICHLET = ( void )
 QUAD_TYPE = LEBEDEV
-QUAD_ORDER = 19
+QUAD_ORDER = 35

code/src/problems/checkerboard.cpp

@@ -2,7 +2,7 @@
 
 Checkerboard::Checkerboard( Config* settings, Mesh* mesh ) : ProblemBase( settings, mesh ) {
     _physics = nullptr;
-    _scatteringXS.resize( _mesh->GetNumCells(), Matrix( _settings->GetNQuadPoints(), _settings->GetNQuadPoints(), 0.0 ) );    // @TODO
+    _scatteringXS.resize( _mesh->GetNumCells(), 1.0 );
     _totalXS.resize( _mesh->GetNumCells(), 1.0 );
     auto cellMids = _mesh->GetCellMidPoints();
     for( unsigned j = 0; j < cellMids.size(); ++j ) {
@@ -15,9 +15,9 @@ Checkerboard::Checkerboard( Config* settings, Mesh* mesh ) : ProblemBase( settin
 
 Checkerboard::~Checkerboard() {}
 
-std::vector<Matrix> Checkerboard::GetScatteringXS( const double energy ) { return _scatteringXS; }
+VectorVector Checkerboard::GetScatteringXS( const std::vector<double>& energies ) { return VectorVector( energies.size(), _scatteringXS ); }
 
-std::vector<double> Checkerboard::GetTotalXS( const double energy ) { return _totalXS; }
+VectorVector Checkerboard::GetTotalXS( const std::vector<double>& energies ) { return VectorVector( energies.size(), _totalXS ); }
 
 std::vector<double> Checkerboard::GetStoppingPower( const std::vector<double>& energies ) {
     // @TODO

code/src/problems/electronrt.cpp

@@ -7,14 +7,14 @@ ElectronRT::ElectronRT( Config* settings, Mesh* mesh ) : ProblemBase( settings,
 
 ElectronRT::~ElectronRT() {}
 
-std::vector<Matrix> ElectronRT::GetScatteringXS( const double energy ) {
+VectorVector ElectronRT::GetScatteringXS( const std::vector<double>& energies ) {
     // @TODO
-    return std::vector<Matrix>( _mesh->GetNumCells(), Matrix( _settings->GetNQuadPoints(), _settings->GetNQuadPoints(), 0.0 ) );
+    return VectorVector( energies.size(), Vector( _mesh->GetNumCells(), 0.0 ) );
 }
 
-std::vector<double> ElectronRT::GetTotalXS( const double energy ) {
+VectorVector ElectronRT::GetTotalXS( const std::vector<double>& energies ) {
     // @TODO
-    return std::vector<double>( _mesh->GetNumCells(), 0.0 );
+    return VectorVector( energies.size(), Vector( _mesh->GetNumCells(), 0.0 ) );
 }
 
 std::vector<double> ElectronRT::GetStoppingPower( const std::vector<double>& energies ) {

code/src/problems/linesource.cpp

@@ -4,11 +4,13 @@ LineSource::LineSource( Config* settings, Mesh* mesh ) : ProblemBase( settings,
 
 LineSource::~LineSource(){};
 
-std::vector<Matrix> LineSource::GetScatteringXS( const double energy ) {
-    return std::vector<Matrix>( _mesh->GetNumCells(), Matrix( _settings->GetNQuadPoints(), _settings->GetNQuadPoints(), 0.0 ) );
+VectorVector LineSource::GetScatteringXS( const std::vector<double>& energies ) {
+    return VectorVector( energies.size(), Vector( _mesh->GetNumCells(), 0.0 ) );
 }
 
-std::vector<double> LineSource::GetTotalXS( const double energy ) { return std::vector<double>( _mesh->GetNumCells(), 0.0 ); }
+VectorVector LineSource::GetTotalXS( const std::vector<double>& energies ) {
+    return VectorVector( energies.size(), Vector( _mesh->GetNumCells(), 0.0 ) );
+}
 
 std::vector<double> LineSource::GetStoppingPower( const std::vector<double>& energies ) {
     // @TODO

code/src/snsolver.cpp

@@ -14,8 +14,6 @@ void SNSolver::Solve() {
 
     // loop over energies (pseudo-time)
     for( unsigned n = 0; n < _nEnergies; ++n ) {
-        auto sigmaT = _problem->GetTotalXS( _energies[n] );
-        auto sigmaS = _problem->GetScatteringXS( _energies[n] );
         // loop over all spatial cells
         for( unsigned j = 0; j < _nCells; ++j ) {
             if( _boundaryCells[j] ) continue;
@@ -28,10 +26,10 @@ void SNSolver::Solve() {
                     psiNew[j][k] -= _g->Flux( _quadPoints[k], _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[j] * _psi[j][k];
+                psiNew[j][k] = _psi[j][k] + ( _dE / _areas[j] ) * psiNew[j][k] - _dE * _sigmaT[n][j] * _psi[j][k];
             }
             // compute scattering effects
-            psiNew[j] += sigmaS[j] * _psi[j] * _weights;    // multiply scattering matrix with psi
+            psiNew[j] += _scatteringKernel * _sigmaS[n][j] * _psi[j] * _weights;    // multiply scattering matrix with psi
         }
         _psi = psiNew;
         if( rank == 0 ) log->info( "{:03.8f}   {:01.5e}", _energies[n], 0.0 );

code/src/solver.cpp

@@ -21,16 +21,19 @@ Solver::Solver( Config* settings ) : _settings( settings ) {
     _weights          = q->GetWeights();
     _nq               = q->GetNq();
     _settings->SetNQuadPoints( _nq );
+    _scatteringKernel = ComputeScatteringKernel();
 
     // set time step
     _dE        = ComputeTimeStep( settings->GetCFL() );
     _nEnergies = unsigned( settings->GetTEnd() / _dE );
-    for( unsigned i = 0; i < _nEnergies; ++i ) _energies.push_back( i * _dE );
+    for( unsigned i = 0; i < _nEnergies; ++i ) _energies.push_back( ( i + 1 ) * _dE );
 
     // setup problem
     _problem = ProblemBase::Create( _settings, _mesh );
     _psi     = _problem->SetupIC();
     _s       = _problem->GetStoppingPower( _energies );
+    _sigmaT  = _problem->GetTotalXS( _energies );
+    _sigmaS  = _problem->GetScatteringXS( _energies );
 
     // setup numerical flux
     _g = NumericalFlux::Create( settings );
@@ -50,4 +53,25 @@ double Solver::ComputeTimeStep( double cfl ) const {
     return cfl * maxEdge;
 }
 
+Matrix Solver::ComputeScatteringKernel() const {
+    Matrix kernel( _nq, _nq, 0.0 );
+    for( unsigned i = 0; i < _nq; ++i ) {
+        auto omega = _quadPoints[i];
+        for( unsigned j = 0; j < _nq; ++j ) {
+            auto omegaprime = _quadPoints[j];
+            double eps      = 4.0 / _nq;
+            auto x          = ( 1 - dot( omega, omegaprime ) );
+            double val      = ( _weights[j] * 1.0 / eps * exp( -( x * x ) / ( eps * eps ) ) );
+            if( val < 1e-8 )
+                kernel( i, j ) = 0.0;
+            else
+                kernel( i, j ) = val;
+        }
+
+        column( kernel, i ) /= sum( column( kernel, i ) );
+        column( kernel, i ) /= _weights;
+    }
+    return kernel;
+}
+
 Solver* Solver::Create( Config* settings ) { return new SNSolver( settings ); }