start to include waterphantom test case as validation test

Former-commit-id: f9a09199

code/include/common/config.h

@@ -67,9 +67,10 @@ class Config
     bool _cleanFluxMat;
     bool _allGaussPts; /*!< @brief If true, the SN Solver uses all Gauss pts in the quadrature */
 
-    bool _csd;                 /*!< @brief If true, continuous slowing down approximation will be used */
-    std::string _hydrogenFile; /*!< @brief Name of hydrogen cross section file */
-    std::string _oxygenFile;   /*!< @brief Name of oxygen cross section file */
+    bool _csd;                      /*!< @brief If true, continuous slowing down approximation will be used */
+    std::string _hydrogenFile;      /*!< @brief Name of hydrogen cross section file path*/
+    std::string _oxygenFile;        /*!< @brief Name of oxygen cross section file path */
+    std::string _stoppingPowerFile; /*!< @brief Name of stopping power file path */
 
     // Boundary Conditions
     /*!< @brief List of all Pairs (marker, BOUNDARY_TYPE), e.g. (farfield,DIRICHLET).
@@ -230,7 +231,7 @@ class Config
     std::string inline GetOutputDir() const { return std::filesystem::path( _outputDir ).lexically_normal(); }
     std::string inline GetOutputFile() const { return std::filesystem::path( _outputFile ).lexically_normal(); }
     std::string inline GetOxygenFile() const { return std::filesystem::path( _oxygenFile ).lexically_normal(); }
-
+    std::string inline GetStoppingPowerFile() const { return std::filesystem::path( _stoppingPowerFile ).lexically_normal(); }
     // Quadrature Structure
     unsigned GetNQuadPoints() { return _nQuadPoints; }
     QUAD_NAME inline GetQuadName() const { return _quadName; }

code/src/common/config.cpp

@@ -247,6 +247,8 @@ void Config::SetConfigOptions() {
     /*! @brief OxygenFile \n DESCRIPTION: If the continuous slowing down approximation is used, this referes to the cross section file for oxygen.
      * . \n DEFAULT "o.dat" \ingroup Config */
     AddStringOption( "OXYGEN_FILE", _oxygenFile, string( "ENDL_O.txt" ) );
+    /*! @brief StoppingPowerFile \n DESCRIPTION: Only temporary added. \ingroup Config */
+    AddStringOption( "STOPPING_POWER_FILE", _stoppingPowerFile, string( "stopping_power.txt" ) );
     /*! @brief SN_ALL_GAUSS_PTS \n DESCRIPTION: If true, the SN Solver uses all Gauss Quadrature Points for 2d. \n DEFAULT false \ingroup Config */
     AddBoolOption( "SN_ALL_GAUSS_PTS", _allGaussPts, false );
 
@@ -393,13 +395,14 @@ void Config::SetPostprocessing() {
     if( _inputDir[_inputDir.size() - 1] != '/' ) _inputDir.append( "/" );
 
     // setup relative paths
-    _logDir       = _inputDir + _logDir;
-    _outputDir    = _inputDir + _outputDir;
-    _meshFile     = _inputDir + _meshFile;
-    _outputFile   = _outputDir + _outputFile;
-    _ctFile       = _inputDir + _ctFile;
-    _hydrogenFile = _inputDir + _hydrogenFile;
-    _oxygenFile   = _inputDir + _oxygenFile;
+    _logDir            = _inputDir + _logDir;
+    _outputDir         = _inputDir + _outputDir;
+    _meshFile          = _inputDir + _meshFile;
+    _outputFile        = _outputDir + _outputFile;
+    _ctFile            = _inputDir + _ctFile;
+    _hydrogenFile      = _inputDir + _hydrogenFile;
+    _oxygenFile        = _inputDir + _oxygenFile;
+    _stoppingPowerFile = _inputDir + _stoppingPowerFile;
 
     // create directories if they dont exist
     if( !std::filesystem::exists( _outputDir ) ) std::filesystem::create_directory( _outputDir );

code/src/problems/electronrt.cpp

@@ -4,44 +4,44 @@
 
 // Constructor: Legacy code, physics class is no longer used
 ElectronRT::ElectronRT( Config* settings, Mesh* mesh ) : ProblemBase( settings, mesh ) {
-    _physics = new Physics( settings->GetHydrogenFile(), settings->GetOxygenFile(), "../input/stopping_power.txt" );
+    _physics = new Physics( settings->GetHydrogenFile(), settings->GetOxygenFile(), settings->GetStoppingPowerFile() );
 }
 
 ElectronRT::~ElectronRT() { delete _physics; }
 
 VectorVector ElectronRT::GetScatteringXS( const Vector& energies ) {
-    // @TODO 
-    // Specified in subclasses 
+    // @TODO
+    // Specified in subclasses
     return VectorVector( energies.size(), Vector( _mesh->GetNumCells(), 0.0 ) );
 }
 
 VectorVector ElectronRT::GetTotalXS( const Vector& energies ) {
-    // @TODO 
-    // Specified in subclasses 
+    // @TODO
+    // Specified in subclasses
     return VectorVector( energies.size(), Vector( _mesh->GetNumCells(), 0.0 ) );
 }
 
 std::vector<Matrix> ElectronRT::GetScatteringXSE( const Vector& energies, const Matrix& angles ) {
-    // @TODO 
-    // Specified in subclasses 
+    // @TODO
+    // Specified in subclasses
     return _physics->GetScatteringXS( energies, angles );
 }
 
 Vector ElectronRT::GetTotalXSE( const Vector& energies ) {
-    // @TODO 
-    // Specified in subclasses 
+    // @TODO
+    // Specified in subclasses
     return _physics->GetTotalXSE( energies );
 }
 
 std::vector<VectorVector> ElectronRT::GetExternalSource( const Vector& energies ) {
-    // @TODO 
-    // Specified in subclasses 
+    // @TODO
+    // Specified in subclasses
     return std::vector<VectorVector>( energies.size(), std::vector<Vector>( _mesh->GetNumCells(), Vector( _settings->GetNQuadPoints(), 0.0 ) ) );
 }
 
 VectorVector ElectronRT::SetupIC() {
-    // @TODO 
-    // Specified in subclasses 
+    // @TODO
+    // Specified in subclasses
     return VectorVector( _mesh->GetNumCells(), Vector( _settings->GetNQuadPoints(), 1e-10 ) );
 }
 
@@ -50,6 +50,5 @@ void ElectronRT::LoadXSH20( std::string /*fileSigmaS*/, std::string /*fileSigmaT
     // Specified in subclasses
 }
 
-
 // Default densities = 1, for patient files this is overwritten with "real" densities
 std::vector<double> ElectronRT::GetDensity( const VectorVector& cellMidPoints ) { return std::vector<double>( cellMidPoints.size(), 1.0 ); }

code/src/quadratures/quadraturebase.cpp

@@ -27,8 +27,10 @@ QuadratureBase* QuadratureBase::Create( Config* settings ) {
         case QUAD_LDFESA: return new QLDFESA( settings );
         case QUAD_Lebedev: return new QLebedev( settings );
         case QUAD_Product: return new QProduct( settings );
-        default: ErrorMessages::Error( "Creator for the chose quadrature does not yet exist. This is is the fault of the coder!", CURRENT_FUNCTION );
+        default: ErrorMessages::Error( "Creator for the chosen quadrature does not yet exist. This is is the fault of the coder!", CURRENT_FUNCTION );
     }
+
+    return nullptr;
 }
 
 QuadratureBase* QuadratureBase::Create( QUAD_NAME name, unsigned quadOrder ) {
@@ -45,6 +47,7 @@ QuadratureBase* QuadratureBase::Create( QUAD_NAME name, unsigned quadOrder ) {
         case QUAD_Product: return new QProduct( quadOrder );
         default: ErrorMessages::Error( "Creator for the chose quadrature does not yet exist. This is is the fault of the coder!", CURRENT_FUNCTION );
     }
+    return nullptr;
 }
 
 double QuadratureBase::Integrate( double( f )( double x0, double x1, double x2 ) ) {

code/src/solvers/solverbase.cpp

@@ -7,6 +7,8 @@
 #include "problems/problembase.h"
 #include "quadratures/quadraturebase.h"
 #include "solvers/csdsnsolver.h"
+#include "solvers/csdsolvertrafofp.h"
+
 #include "solvers/mnsolver.h"
 #include "solvers/pnsolver.h"
 #include "solvers/snsolver.h"
@@ -70,11 +72,13 @@ Solver::~Solver() {
 Solver* Solver::Create( Config* settings ) {
     switch( settings->GetSolverName() ) {
         case SN_SOLVER: return new SNSolver( settings );
-        case CSD_SN_SOLVER: return new CSDSNSolver( settings );
         case PN_SOLVER: return new PNSolver( settings );
         case MN_SOLVER: return new MNSolver( settings );
-        default: return new SNSolver( settings );
+        case CSD_SN_SOLVER: return new CSDSNSolver( settings );
+        case CSD_SN_FOKKERPLANCK_TRAFO_SOLVER: return new CSDSolverTrafoFP( settings );
+        default: ErrorMessages::Error( "Creator for the chosen solver does not yet exist. This is is the fault of the coder!", CURRENT_FUNCTION );
     }
+    return nullptr;
 }
 
 void Solver::Solve() {

code/tests/test_cases.cpp

@@ -189,6 +189,30 @@ TEST_CASE( "MN_SOLVER", "[validation_tests]" ) {
     }
 }
 
+TEST_CASE( "CSD_SN_FP_SOLVER", "[validation_tests]" ) {
+    std::string sn_fileDir = "input/validation_tests/CSD_SN_FP_solver/";
+    SECTION( "waterphantom 1D" ) {
+        std::string config_file_name = std::string( TESTS_PATH ) + sn_fileDir + "waterphantom_1D.cfg";
+
+        Config* config = new Config( config_file_name );
+        Solver* solver = Solver::Create( config );
+
+        solver->Solve();
+        solver->PrintVolumeOutput();
+
+        auto test      = readVTKFile( std::string( TESTS_PATH ) + "result/rtsn_test_waterphantom_1D_CSD_FP_1.vtk" );
+        auto reference = readVTKFile( std::string( TESTS_PATH ) + sn_fileDir + "waterphantom_1D_CSD_SN_FP_reference.vtk" );
+
+        double eps = 1e-3;
+        REQUIRE( test.size() == reference.size() );
+        bool errorWithinBounds = true;
+        for( unsigned i = 0; i < test.size(); ++i ) {
+            if( std::fabs( test[i] - reference[i] ) > eps ) errorWithinBounds = false;
+        }
+        REQUIRE( errorWithinBounds );
+    }
+}
+
 // --- Validation Tests Output ---
 void tokenize( std::string const& str, const char delim, std::vector<std::string>& out ) {
     size_t start;