solver trafo FP added

code/include/common/globalconstants.h

@@ -89,10 +89,10 @@ enum KERNEL_NAME { KERNEL_Isotropic, KERNEL_Isotropic1D };
 inline std::map<std::string, KERNEL_NAME> Kernel_Map{ { "ISOTROPIC", KERNEL_Isotropic }, { "ISOTROPIC_1D", KERNEL_Isotropic1D } };
 
 // Solver name
-enum SOLVER_NAME { SN_SOLVER, CSD_SN_SOLVER, CSD_SN_NOTRAFO_SOLVER, CSD_SN_FOKKERPLANCK_SOLVER, PN_SOLVER, MN_SOLVER };
+enum SOLVER_NAME { SN_SOLVER, CSD_SN_SOLVER, CSD_SN_NOTRAFO_SOLVER, CSD_SN_FOKKERPLANCK_SOLVER, CSD_SN_FOKKERPLANCK_TRAFO_SOLVER, PN_SOLVER, MN_SOLVER };
 
 inline std::map<std::string, SOLVER_NAME> Solver_Map{
-    { "SN_SOLVER", SN_SOLVER }, { "CSD_SN_SOLVER", CSD_SN_SOLVER },{ "CSD_SN_NOTRAFO_SOLVER", CSD_SN_NOTRAFO_SOLVER }, { "CSD_SN_FOKKERPLANCK_SOLVER", CSD_SN_FOKKERPLANCK_SOLVER }, { "PN_SOLVER", PN_SOLVER }, { "MN_SOLVER", MN_SOLVER } };
+    { "SN_SOLVER", SN_SOLVER }, { "CSD_SN_SOLVER", CSD_SN_SOLVER },{ "CSD_SN_NOTRAFO_SOLVER", CSD_SN_NOTRAFO_SOLVER }, { "CSD_SN_FOKKERPLANCK_SOLVER", CSD_SN_FOKKERPLANCK_SOLVER }, { "CSD_SN_FOKKERPLANCK_TRAFO_SOLVER", CSD_SN_FOKKERPLANCK_TRAFO_SOLVER }, { "PN_SOLVER", PN_SOLVER }, { "MN_SOLVER", MN_SOLVER } };
 
 // Entropy functional
 enum ENTROPY_NAME { QUADRATIC, MAXWELL_BOLZMANN, BOSE_EINSTEIN, FERMI_DIRAC };

code/include/solvers/csdsolvertrafofp.h

+#ifndef CSDSOLVERTRAFOFP_H
+#define CSDSOLVERTRAFOFP_H
+
+#include "icru.h"
+#include "solvers/snsolver.h"
+
+class Physics;
+
+class CSDSolverTrafoFP : public SNSolver
+{
+  private:
+    std::vector<double> _dose; /*! @brief: TODO */
+
+    // Physics acess
+    Vector _energies; /*! @brief: energy levels for CSD, lenght = _nEnergies */
+    Vector _angle;    /*! @brief: angles for SN */
+    Vector _density;  /*! @brief: patient density for each grid cell */
+
+    std::vector<Matrix> _sigmaSE; /*!  @brief scattering cross section for all energies*/
+    Vector _sigmaTE;              /*!  @brief total cross section for all energies*/
+
+    Matrix _L; /*!  @brief Laplace Beltrami Matrix */
+    Matrix _IL; /*!  @brief Laplace Beltrami Matrix */
+
+    double _alpha;
+    double _beta;
+
+    Vector _xi1;
+    Vector _xi2;
+    Matrix _xi;
+
+    bool _RT;
+
+    double _energyMin;
+    double _energyMax;
+
+  public:
+    /**
+     * @brief CSDSolverTrafoFP constructor
+     * @param settings stores all needed information
+     */
+    CSDSolverTrafoFP( 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 // CSDSOLVERTRAFOFP_H

code/input/example_csd_1d.cfg

@@ -3,16 +3,16 @@ OUTPUT_FILE = example_csd_1d_10MeV_fine
 LOG_DIR = ../result/logs
 MESH_FILE = 1DMesh.su2
 PROBLEM = WATERPHANTOM
-SOLVER = CSD_SN_FOKKERPLANCK_SOLVER
+SOLVER = CSD_SN_FOKKERPLANCK_TRAFO_SOLVER
 CONTINUOUS_SLOWING_DOWN = YES
 HYDROGEN_FILE = ENDL_H.txt
 OXYGEN_FILE = ENDL_O.txt
 KERNEL = ISOTROPIC_1D
-CFL_NUMBER = 0.01
+CFL_NUMBER = 0.008
 TIME_FINAL = 1.0
 CLEAN_FLUX_MATRICES = NO
 
 BC_DIRICHLET = ( dirichlet )
 BC_NEUMANN = ( wall_low, wall_up )
 QUAD_TYPE = GAUSS_LEGENDRE_1D
-QUAD_ORDER = 11
+QUAD_ORDER = 26

code/src/solvers/csdsolvertrafofp.cpp

+#include "solvers/csdsolvertrafofp.h"
+#include "common/config.h"
+#include "common/io.h"
+#include "fluxes/numericalflux.h"
+#include "kernels/scatteringkernelbase.h"
+#include "problems/problembase.h"
+#include "quadratures/quadraturebase.h"
+
+// externals
+#include "spdlog/spdlog.h"
+#include <mpi.h>
+
+CSDSolverTrafoFP::CSDSolverTrafoFP( Config* settings ) : SNSolver( settings ) {
+    _dose = std::vector<double>( _settings->GetNCells(), 0.0 );
+
+    // Set angle and energies
+    _angle           = Vector( _settings->GetNQuadPoints(), 0.0 );    // my
+    _energies        = Vector( _nEnergies, 0.0 );                     // equidistant
+     _energyMin = 1e-4;
+     //_energyMax = 10.0;
+     _energyMax = 5.0;
+    // write equidistant energy grid
+
+    _dE        = ComputeTimeStep( settings->GetCFL() );
+    _nEnergies = unsigned( ( _energyMax - _energyMin ) / _dE );
+    std::cout<<"nEnergies = "<<_nEnergies<<std::endl;
+    _energies.resize( _nEnergies );
+    for( unsigned n = 0; n < _nEnergies; ++n ) {
+        _energies[n] = _energyMin + ( _energyMax - _energyMin ) / ( _nEnergies - 1 ) * n;
+    }
+
+    // create 1D quadrature
+    unsigned nq            = _settings->GetNQuadPoints();
+    QuadratureBase* quad1D = QuadratureBase::CreateQuadrature( QUAD_GaussLegendre1D, nq );
+    Vector w               = quad1D->GetWeights();
+    VectorVector muVec     = quad1D->GetPoints();
+    Vector mu( nq );
+    for( unsigned k = 0; k < nq; ++k ) {
+        mu[k] = muVec[k][0];
+    }
+
+    // setup Laplace Beltrami matrix L in slab geometry
+    _L = Matrix( nq, nq, 0.0 );
+
+    double DMinus = 0.0;
+    double DPlus  = 0.0;
+    for( unsigned k = 0; k < nq; ++k ) {
+        DMinus = DPlus;
+        DPlus  = DMinus - 2 * mu[k] * w[k];
+        if( k > 0 ) {
+            _L( k, k - 1 ) = DMinus / ( mu[k] - mu[k - 1] ) / w[k];
+            _L( k, k )     = -DMinus / ( mu[k] - mu[k - 1] ) / w[k];
+        }
+        if( k < nq - 1 ) {
+            _L( k, k + 1 ) = DPlus / ( mu[k + 1] - mu[k] ) / w[k];
+            _L( k, k ) += -DPlus / ( mu[k + 1] - mu[k] ) / w[k];
+        }
+    }
+
+    // Heney-Greenstein parameter
+    double g = 0.8;
+
+    // determine momente of Heney-Greenstein
+    _xi1 = Vector(_nEnergies, 1.0 - g); // paper Olbrant, Frank (11)
+    _xi2 = Vector(_nEnergies, 4.0 / 3.0 - 2.0 * g + 2.0 / 3.0 * g * g);
+    _xi = Matrix(3,_nEnergies);
+    for(unsigned n = 0; n<_nEnergies; ++n ){
+        _xi(1,n) = 1.0 - g;
+        _xi(2,n) = 4.0 / 3.0 - 2.0 * g + 2.0 / 3.0 * g * g;
+    }
+
+    _s = Vector( _nEnergies, 1.0 );
+
+    _RT = false;
+
+    // read in medical data if radiation therapy option selected
+    if(_RT){
+        /*
+        _nEnergies = 100;
+        _energies.resize(_nEnergies);
+        _xi = Matrix(6,_nEnergies);
+        double minExp = -4.0;
+        double  maxExp = 1.5;
+        for( int n = 0; n<_nEnergies; ++n){
+            double exponent = minExp + ( maxExp - minExp ) / ( _nEnergies - 1 ) * n;
+            _energies[n] = pow(10.0,exponent);
+        }*/
+        ICRU database( mu, _energies );
+        database.GetTransportCoefficients( _xi );
+        database.GetStoppingPower( _s );
+/*
+        // print coefficients
+        std::cout<<"E = [";
+        for( unsigned n = 0; n<_nEnergies; ++n){
+            std::cout<<_energies[n]<<"; ";
+        }
+        std::cout<<"];"<<std::endl;
+        std::cout<<"xi = [";
+        for( unsigned n = 0; n<_nEnergies; ++n){
+            std::cout<<_xi(0,n)<<" "<<_xi(1,n)<<" "<<_xi(2,n)<<" "<<_xi(3,n)<<" "<<_xi(4,n)<<" "<<_xi(5,n)<<"; ";
+        }
+        std::cout<<"];"<<std::endl;*/
+
+    }
+
+    // recompute scattering kernel. TODO: add this to kernel function
+    for( unsigned p = 0; p < _nq; ++p ) {
+        for( unsigned q = 0; q < _nq; ++q ) {
+            _scatteringKernel( p, q ) = 0.0;
+        }
+        _scatteringKernel( p, p ) = _weights[p];
+    }
+
+    _density = Vector( _nCells, 1.0 );
+    //exit(EXIT_SUCCESS);
+}
+
+void CSDSolverTrafoFP::Solve() {
+    std::cout << "Solve" << std::endl;
+    auto log = spdlog::get( "event" );
+
+    auto energiesOrig = _energies;
+
+    // setup IC and incoming BC on left
+    // auto cellMids = _settings->GetCellMidPoints();
+    _sol = std::vector<Vector>( _nCells, Vector( _nq, 0.0 ) );
+    for( unsigned k = 0; k < _nq; ++k ) {
+        if( _quadPoints[k][0] > 0 && !_RT) _sol[0][k] = 1e5 * exp( -10.0 * pow( 1.0 - _quadPoints[k][0], 2 ) );
+    }
+    _boundaryCells[0]           = BOUNDARY_TYPE::DIRICHLET;
+    _boundaryCells[_nCells - 1] = BOUNDARY_TYPE::DIRICHLET;
+
+    Matrix identity( _nq, _nq, 0.0 );
+    for( unsigned k = 0; k < _nq; ++k ) identity( k, k ) = 1.0;
+
+    // angular flux at next time step (maybe store angular flux at all time steps, since time becomes energy?)
+    VectorVector psiNew( _nCells, Vector( _nq, 0.0 ) );
+    double dFlux = 1e10;
+    Vector fluxNew( _nCells, 0.0 );
+    Vector fluxOld( _nCells, 0.0 );
+    //for( unsigned j = 0; j < _nCells; ++j ) {
+    //    fluxOld[j] = dot( _sol[j], _weights );
+    //}
+
+    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 ) {
+            _sol[j][k] = _sol[j][k] * _density[j] * _s[_nEnergies - 1];    // note that _s[_nEnergies - 1] is stopping power at highest energy
+        }
+    }
+
+    VectorVector psi1 = _sol;
+
+    // store transformed energies ETilde instead of E in _energies vector (cf. Dissertation Kerstion Kuepper, Eq. 1.25)
+    double tmp   = 0.0;
+    _energies[0] = 0.0;
+    for( unsigned n = 1; n < _nEnergies; ++n ) {
+        tmp          = tmp + _dE * 0.5 * ( 1.0 / _s[n] + 1.0 / _s[n - 1] );
+        _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];
+    }
+
+    // determine minimal density for CFL computation
+    double densityMin = _density[0];
+    for( unsigned j = 1; j < _nCells; ++j ) {
+        if( densityMin > _density[j] ) densityMin = _density[j];
+    }
+
+    // cross sections do not need to be transformed to ETilde energy grid since e.g. TildeSigmaT(ETilde) = SigmaT(E(ETilde))
+
+    // loop over energies (pseudo-time)
+    for( unsigned n = 0; n < _nEnergies - 1; ++n ) {
+        _dE = fabs( _energies[n + 1] - _energies[n] );    // is the sign correct here?
+
+        // set alpha and beta
+        _alpha = _xi(1,_nEnergies-n-1) / 2.0 + _xi(2,_nEnergies-n-1) / 8.0;
+        _beta  = _xi(2,_nEnergies-n-1) / 8.0 / _xi(1,_nEnergies-n-1);
+
+        _IL = identity - _beta * _L;
+
+        // write BC for water phantom
+        if(_RT){
+            for( unsigned k = 0; k < _nq; ++k ) {
+                if( _quadPoints[k][0] > 0 ) {
+                        //_sol[0][k] = 1e5 * exp( -200.0 * pow( 1.0 - _quadPoints[k][0], 2 ) ) * exp(-50*pow(_energies[0]-_energies[n],2))* _density[0] * _s[_nEnergies - n- 1];
+                        _sol[0][k] = 1e5 * exp( -200.0 * pow( 1.0 - _quadPoints[k][0], 2 ) ) * exp(-50.0*pow(_energyMax-energiesOrig[_nEnergies - n - 1],2))* _density[0] * _s[_nEnergies - n - 1];
+                        //psiNew[0][k] = _sol[0][k];
+                }
+            }
+        }
+
+        for( unsigned j = 0; j < _nCells; ++j ) {
+            if( _boundaryCells[j] == BOUNDARY_TYPE::DIRICHLET ) continue;
+            psi1[j] = blaze::solve( _IL, _sol[j] );
+            psi1[j] = _alpha * _L * psi1[j];
+        }
+
+        // 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 i = 0; i < _nq; ++i ) {
+                psiNew[j][i] = 0.0;
+                // loop over all neighbor cells (edges) of cell j and compute numerical fluxes
+                for( unsigned idx_neighbor = 0; idx_neighbor < _neighbors[j].size(); ++idx_neighbor ) {
+                    // store flux contribution on psiNew_sigmaS to save memory
+                    if( _boundaryCells[j] == BOUNDARY_TYPE::NEUMANN && _neighbors[j][idx_neighbor] == _nCells )
+                        continue;    // adiabatic wall, add nothing
+                    else
+                        psiNew[j][i] += _g->Flux( _quadPoints[i],
+                                                  _sol[j][i] / _density[j],
+                                                  _sol[_neighbors[j][idx_neighbor]][i] / _density[_neighbors[j][idx_neighbor]],
+                                                  _normals[j][idx_neighbor] ) / _areas[j];
+                }
+                // time update angular flux with numerical flux and total scattering cross section
+                psiNew[j][i] = _sol[j][i] - _dE * psiNew[j][i] + _dE * psi1[j][i];
+            }
+        }
+        _sol = psiNew;
+
+        for( unsigned j = 0; j < _nCells; ++j ) {
+            fluxNew[j] = dot( psiNew[j], _weights );
+            if(n > 0){
+                _dose[j] += 0.5 * _dE * ( fluxNew[j] * _s[_nEnergies - n - 1] + fluxOld[j] * _s[_nEnergies - n] ) / _density[j];    // update dose with trapezoidal rule
+            }else{
+                _dose[j] += _dE * fluxNew[j] * _s[_nEnergies - n - 1]/ _density[j];
+            }
+            _solverOutput[j] = fluxNew[j];
+        }
+
+        // Save( n );
+        dFlux   = blaze::l2Norm( fluxNew - fluxOld );
+        fluxOld = fluxNew;
+        if( rank == 0 ) log->info( "{:03.8f}  {:01.5e}  {:01.5e}", _energies[n], _dE / densityMin, dFlux );
+        if( std::isinf( dFlux ) || std::isnan( dFlux ) ) break;
+    }
+}
+
+void CSDSolverTrafoFP::Save() const {
+    std::vector<std::string> fieldNames{ "dose", "normalized dose" };
+    std::vector<std::vector<double>> dose( 1, _dose );
+    std::vector<std::vector<double>> normalizedDose( 1, _dose );
+    double maxDose = *std::max_element( _dose.begin(), _dose.end() );
+    for( unsigned i = 0; i < _dose.size(); ++i ) normalizedDose[0][i] /= maxDose;
+    std::vector<std::vector<std::vector<double>>> results{ dose, normalizedDose };
+    ExportVTK( _settings->GetOutputFile(), results, fieldNames, _mesh );
+}
+
+void CSDSolverTrafoFP::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

@@ -12,6 +12,7 @@
 #include "solvers/snsolver.h"
 #include "solvers/csdsnsolvernotrafo.h"
 #include "solvers/csdsnsolverfp.h"
+#include "solvers/csdsolvertrafofp.h"
 
 Solver::Solver( Config* settings ) : _settings( settings ) {
     // @TODO save parameters from settings class
@@ -76,6 +77,7 @@ Solver* Solver::Create( Config* settings ) {
         case CSD_SN_SOLVER: return new CSDSNSolver( settings );
         case CSD_SN_NOTRAFO_SOLVER: return new CSDSNSolverNoTrafo( settings );
         case CSD_SN_FOKKERPLANCK_SOLVER: return new CSDSNSolverFP( settings );
+        case CSD_SN_FOKKERPLANCK_TRAFO_SOLVER: return new CSDSolverTrafoFP( settings );
         case PN_SOLVER: return new PNSolver( settings );
         case MN_SOLVER: return new MNSolver( settings );
         default: return new SNSolver( settings );