Commit b788f931 authored by jannick.wolters's avatar jannick.wolters
Browse files

added endl config options

parent 3b1e3d37
Pipeline #100357 passed with stages
in 35 minutes and 17 seconds
......@@ -62,8 +62,9 @@ class Config
* to improve floating point accuracy */
bool _cleanFluxMat;
/*!< @brief If true, continuous slowing down approximation will be used */
bool _csd;
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 */
// Boundary Conditions
/*!< @brief List of all Pairs (marker, BOUNDARY_TYPE), e.g. (farfield,DIRICHLET).
......@@ -204,6 +205,8 @@ class Config
std::string inline GetOutputFile() const { return std::filesystem::path( _outputFile ).lexically_normal(); }
std::string inline GetLogDir() const { return std::filesystem::path( _logDir ).lexically_normal(); }
std::string inline GetCTFile() const { return std::filesystem::path( _ctFile ).lexically_normal(); }
std::string inline GetHydrogenFile() const { return std::filesystem::path( _hydrogenFile ).lexically_normal(); }
std::string inline GetOxygenFile() const { return std::filesystem::path( _oxygenFile ).lexically_normal(); }
// Quadrature Structure
QUAD_NAME inline GetQuadName() const { return _quadName; }
......
......@@ -18,15 +18,13 @@ class Physics
VectorVector _xsTransportH2O;
VectorVector _stpowH2O;
public:
// prototype data readers
std::tuple<std::vector<VectorVector>, std::vector<VectorVector>> ReadENDL( std::string filename );
VectorVector ReadStoppingPowers( std::string fileName );
// load and prepare data from database
void LoadDatabase( std::string fileName_H, std::string fileName_O, std::string fileName_stppower );
Physics() = 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
......@@ -35,6 +33,8 @@ class Physics
*/
VectorVector GetScatteringXS( Vector energies, Vector angle );
VectorVector GetScatteringXSE( Vector energies );
/**
* @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
......@@ -42,6 +42,8 @@ class Physics
*/
VectorVector GetTotalXS( Vector energies, Vector density );
VectorVector GetTotalXSE( Vector energies );
/**
* @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
......@@ -57,18 +59,17 @@ class Physics
*/
VectorVector GetTransportXS( Vector energies, Vector density );
VectorVector GetTransportXSE( Vector energies );
/**
* @brief Physics constructor
* @param settings stores all needed user information
*/
Physics();
Physics( std::string fileName_H, std::string fileName_O, std::string fileName_stppower );
/**
* @brief Create constructor
* @param settings stores all needed information
* @return pointer to Physics
* @brief Physics destructor
*/
static Physics* Create();
~Physics();
};
#endif
......@@ -18,9 +18,9 @@ class Checkerboard : public ProblemBase
Checkerboard( Config* settings, Mesh* mesh );
virtual ~Checkerboard();
virtual VectorVector GetScatteringXS( const std::vector<double>& energies );
virtual VectorVector GetTotalXS( const std::vector<double>& energies );
virtual std::vector<VectorVector> GetExternalSource( const std::vector<double>& energies );
virtual VectorVector GetScatteringXS( const Vector& energies );
virtual VectorVector GetTotalXS( const Vector& energies );
virtual std::vector<VectorVector> GetExternalSource( const Vector& energies );
virtual Vector GetStoppingPower( const Vector& energies );
virtual VectorVector SetupIC();
};
......
......@@ -20,11 +20,11 @@ class ElectronRT : public ProblemBase
ElectronRT( Config* settings, Mesh* mesh );
virtual ~ElectronRT();
virtual VectorVector GetScatteringXS( const std::vector<double>& energies );
virtual VectorVector GetTotalXS( const std::vector<double>& energies );
virtual VectorVector GetScatteringXS( const Vector& energies );
virtual VectorVector GetTotalXS( const Vector& energies );
virtual VectorVector GetScatteringXSE( const Vector& energies, const Vector& angles );
virtual Vector GetTotalXSE( const Vector& energies );
virtual std::vector<VectorVector> GetExternalSource( const std::vector<double>& energies );
virtual std::vector<VectorVector> GetExternalSource( const Vector& energies );
virtual Vector GetStoppingPower( const Vector& energies );
virtual VectorVector SetupIC();
std::vector<double> GetDensity( const VectorVector& cellMidPoints );
......
......@@ -12,9 +12,9 @@ class LineSource_SN : public ProblemBase
LineSource_SN( Config* settings, Mesh* mesh );
virtual ~LineSource_SN();
virtual VectorVector GetScatteringXS( const std::vector<double>& energies );
virtual VectorVector GetTotalXS( const std::vector<double>& energies );
virtual std::vector<VectorVector> GetExternalSource( const std::vector<double>& energies );
virtual VectorVector GetScatteringXS( const Vector& energies );
virtual VectorVector GetTotalXS( const Vector& energies );
virtual std::vector<VectorVector> GetExternalSource( const Vector& energies );
virtual Vector GetStoppingPower( const Vector& energies );
virtual VectorVector SetupIC();
};
......@@ -49,9 +49,9 @@ class LineSource_PN : public ProblemBase
LineSource_PN( Config* settings, Mesh* mesh );
virtual ~LineSource_PN();
virtual VectorVector GetScatteringXS( const std::vector<double>& energies ) override;
virtual VectorVector GetTotalXS( const std::vector<double>& energies ) override;
virtual std::vector<VectorVector> GetExternalSource( const std::vector<double>& energies ) override;
virtual VectorVector GetScatteringXS( const Vector& energies ) override;
virtual VectorVector GetTotalXS( const Vector& energies ) override;
virtual std::vector<VectorVector> GetExternalSource( const Vector& energies ) override;
virtual Vector GetStoppingPower( const Vector& energies ) override;
virtual VectorVector SetupIC() override;
};
......
......@@ -28,7 +28,7 @@ class ProblemBase
* in vector energy
* @param energy is the energy the cross section is queried for
*/
virtual VectorVector GetScatteringXS( const std::vector<double>& energies ) = 0;
virtual VectorVector GetScatteringXS( const Vector& energies ) = 0;
/**
* @brief GetTotalXS gives back vector of vectors of total cross sections for
......@@ -36,7 +36,7 @@ class ProblemBase
* @param energy is the energy the cross section is queried for
* @param density is vector with patient densities (at different spatial cells)
*/
virtual VectorVector GetTotalXS( const std::vector<double>& energies ) = 0;
virtual VectorVector GetTotalXS( const Vector& energies ) = 0;
/**
* @brief GetTotalXSE gives back vector of total cross sections for
* energies in vector energy
......@@ -56,7 +56,7 @@ class ProblemBase
* energy, cell and angle
* @param energies is vector with energies
*/
virtual std::vector<VectorVector> GetExternalSource( const std::vector<double>& energies ) = 0;
virtual std::vector<VectorVector> GetExternalSource( const Vector& energies ) = 0;
/**
* @brief GetStoppingPower gives back vector of vectors of stopping powers for
......
......@@ -12,7 +12,7 @@ class WaterPhantom : public ElectronRT
WaterPhantom( Config* settings, Mesh* mesh );
virtual ~WaterPhantom();
virtual std::vector<VectorVector> GetExternalSource( const std::vector<double>& energies );
virtual std::vector<VectorVector> GetExternalSource( const Vector& energies );
virtual Vector GetStoppingPower( const Vector& energies );
virtual VectorVector SetupIC();
std::vector<double> GetDensity( const VectorVector& cellMidPoints );
......
......@@ -11,10 +11,9 @@ class CSDSNSolver : public SNSolver
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 */
Physics* _physics; /*! @brief: class that hold physics information */
Vector _energies; /*! @brief: energy levels for CSD, lenght = _nEnergies */
Vector _angle; /*! @brief: angles for SN */
Vector _density; /*! @brief: patient density for each grid cell */
VectorVector _sigmaSE; /*! @brief scattering cross section for all energies*/
Vector _sigmaTE; /*! @brief total cross section for all energies*/
......
......@@ -22,12 +22,12 @@ class Solver
// --------- Often used variables of member classes for faster access ----
unsigned _nEnergies; /*! @brief number of energy/time steps, number of nodal energy values for CSD */
double _dE; /*! @brief energy/time step size */
std::vector<double> _energies; // energy groups used in the simulation [keV]
std::vector<double> _density; // patient density, dim(_density) = _nCells
Vector _s; // stopping power, dim(_s) = _nTimeSteps
std::vector<VectorVector> _Q; /*! @brief external source term */
unsigned _nEnergies; /*! @brief number of energy/time steps, number of nodal energy values for CSD */
double _dE; /*! @brief energy/time step size */
Vector _energies; // energy groups used in the simulation [keV]
std::vector<double> _density; // patient density, dim(_density) = _nCells
Vector _s; // stopping power, dim(_s) = _nTimeSteps
std::vector<VectorVector> _Q; /*! @brief external source term */
VectorVector _sigmaS; /*! @brief scattering cross section for all energies and spatial cells*/
VectorVector _sigmaT; /*! @brief total cross section for all energies and spatial cells*/
......
This diff is collapsed.
This diff is collapsed.
OUTPUT_DIR = ../result
OUTPUT_FILE = example_IC3_lCFL
LOG_DIR = ../result/logs
MESH_FILE = linesource.su2
PROBLEM = ELECTRONRT
SOLVER = SN_SOLVER
CONTINUOUS_SLOWING_DOWN = YES
HYDROGEN_FILE = ENDL_H.txt
OXYGEN_FILE = ENDL_O.txt
CFL_NUMBER = 0.9
TIME_FINAL = 0.5
CLEAN_FLUX_MATRICES = NO
BC_DIRICHLET = ( void )
QUAD_TYPE = MONTE_CARLO
QUAD_ORDER = 2
ESTAR: Stopping Powers and Range Tables for Electrons
WATER, LIQUID
Kinetic Total
Energy Stp. Pow.
MeV MeV cm2/g
1.000E-02 2.256E+01
1.250E-02 1.898E+01
1.500E-02 1.647E+01
1.750E-02 1.461E+01
2.000E-02 1.318E+01
2.500E-02 1.110E+01
3.000E-02 9.657E+00
3.500E-02 8.596E+00
4.000E-02 7.781E+00
4.500E-02 7.134E+00
5.000E-02 6.607E+00
5.500E-02 6.170E+00
6.000E-02 5.801E+00
7.000E-02 5.211E+00
8.000E-02 4.761E+00
9.000E-02 4.407E+00
1.000E-01 4.119E+00
1.250E-01 3.596E+00
1.500E-01 3.242E+00
1.750E-01 2.988E+00
2.000E-01 2.798E+00
2.500E-01 2.533E+00
3.000E-01 2.360E+00
3.500E-01 2.241E+00
4.000E-01 2.154E+00
4.500E-01 2.090E+00
5.000E-01 2.041E+00
5.500E-01 2.003E+00
6.000E-01 1.972E+00
7.000E-01 1.926E+00
8.000E-01 1.896E+00
9.000E-01 1.876E+00
1.000E+00 1.862E+00
1.250E+00 1.845E+00
1.500E+00 1.841E+00
1.750E+00 1.844E+00
2.000E+00 1.850E+00
2.500E+00 1.868E+00
3.000E+00 1.889E+00
3.500E+00 1.910E+00
4.000E+00 1.931E+00
4.500E+00 1.951E+00
5.000E+00 1.971E+00
5.500E+00 1.991E+00
6.000E+00 2.010E+00
7.000E+00 2.047E+00
8.000E+00 2.082E+00
9.000E+00 2.116E+00
1.000E+01 2.149E+00
1.250E+01 2.230E+00
1.500E+01 2.306E+00
1.750E+01 2.381E+00
2.000E+01 2.454E+00
2.500E+01 2.598E+00
3.000E+01 2.738E+00
3.500E+01 2.876E+00
4.000E+01 3.013E+00
4.500E+01 3.150E+00
5.000E+01 3.286E+00
5.500E+01 3.421E+00
6.000E+01 3.556E+00
7.000E+01 3.827E+00
8.000E+01 4.096E+00
9.000E+01 4.366E+00
1.000E+02 4.636E+00
1.250E+02 5.311E+00
1.500E+02 5.987E+00
1.750E+02 6.663E+00
2.000E+02 7.341E+00
2.500E+02 8.698E+00
3.000E+02 1.006E+01
3.500E+02 1.142E+01
4.000E+02 1.278E+01
4.500E+02 1.414E+01
5.000E+02 1.551E+01
5.500E+02 1.688E+01
6.000E+02 1.824E+01
7.000E+02 2.098E+01
8.000E+02 2.371E+01
9.000E+02 2.645E+01
1.000E+03 2.919E+01
\ No newline at end of file
......@@ -222,6 +222,12 @@ void Config::SetConfigOptions() {
/*! @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 );
/*! @brief HydogenFile \n DESCRIPTION: If the continuous slowing down approximation is used, this referes to the cross section file for hydrogen.
* . \n DEFAULT "h.dat" \ingroup Config */
AddStringOption( "HYDROGEN_FILE", _hydrogenFile, string( "ENDL_H.txt" ) );
/*! @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" ) );
// Entropy related options
/*! @brief Entropy Functional \n DESCRIPTION: Entropy functional used for the MN_Solver \n DEFAULT QUADRTATIC @ingroup Config. */
......@@ -371,6 +377,17 @@ void Config::SetPostprocessing() {
// if( !std::filesystem::exists( _meshFile ) ) {
// ErrorMessages::Error( "Path to mesh file <" + _meshFile + "> does not exist. Please check your config file.", CURRENT_FUNCTION );
//}
if( this->IsCSD() ) {
if( !std::filesystem::exists( this->GetHydrogenFile() ) ) {
ErrorMessages::Error( "Path to mesh file <" + this->GetHydrogenFile() + "> does not exist. Please check your config file.",
CURRENT_FUNCTION );
}
if( !std::filesystem::exists( this->GetOxygenFile() ) ) {
ErrorMessages::Error( "Path to mesh file <" + this->GetOxygenFile() + "> does not exist. Please check your config file.",
CURRENT_FUNCTION );
}
}
}
void Config::SetOutput() {
......
#include "physics.h"
Physics::Physics() {
// LoadDatabase
Physics::Physics( std::string fileName_H, std::string fileName_O, std::string fileName_stppower ) {
LoadDatabase( fileName_H, fileName_O, fileName_stppower );
}
Physics::~Physics() {}
void Physics::LoadDatabase( std::string fileName_H, std::string fileName_O, std::string fileName_stppower ) {
VectorVector transport_XS_H;
VectorVector transport_XS_O;
......@@ -165,6 +167,32 @@ VectorVector Physics::GetTotalXS( Vector energies, Vector density ) {
return total_XS;
}
VectorVector Physics::GetTotalXSE( Vector energies ) {
VectorVector total_XS;
double xsH2O_i;
Spline interp;
interp.set_points( _xsTotalH2O[0], _xsTotalH2O[1], false ); // false == linear interpolation
for( unsigned i = 0; i < energies.size(); i++ ) {
if( energies[i] < _xsTotalH2O[1][0] ) {
xsH2O_i = _xsTotalH2O[1][0];
}
else if( energies[i] > _xsTotalH2O[1][energies.size() - 1] ) {
xsH2O_i = _xsTotalH2O[1][energies.size() - 1];
}
else {
// Linear interpolation
xsH2O_i = interp( energies[i] );
}
total_XS[i] = xsH2O_i;
}
return total_XS;
}
Vector Physics::GetStoppingPower( Vector energies ) {
Vector stopping_power;
double stpw_H2O_i;
......@@ -215,7 +243,30 @@ VectorVector Physics::GetTransportXS( Vector energies, Vector density ) {
return transport_XS;
}
Physics* Physics::Create() { return new Physics(); }
VectorVector Physics::GetTransportXSE( Vector energies ) {
VectorVector transport_XS;
double xsH2O_i;
Spline interp;
interp.set_points( _xsTransportH2O[0], _xsTransportH2O[1], false ); // false == linear interpolation
for( unsigned i = 0; i < energies.size(); i++ ) {
if( energies[i] < _xsTransportH2O[1][0] ) {
xsH2O_i = _xsTransportH2O[1][0];
}
else if( energies[i] > _xsTransportH2O[1][energies.size() - 1] ) {
xsH2O_i = _xsTransportH2O[1][energies.size() - 1];
}
else {
// Linear interpolation
xsH2O_i = interp( energies[i] );
}
transport_XS[i] = xsH2O_i;
}
return transport_XS;
}
std::tuple<std::vector<VectorVector>, std::vector<VectorVector>> Physics::ReadENDL( std::string filename ) {
std::vector<VectorVector> _headers;
......
......@@ -17,11 +17,11 @@ Checkerboard::Checkerboard( Config* settings, Mesh* mesh ) : ProblemBase( settin
Checkerboard::~Checkerboard() {}
VectorVector Checkerboard::GetScatteringXS( const std::vector<double>& energies ) { return VectorVector( energies.size(), _scatteringXS ); }
VectorVector Checkerboard::GetScatteringXS( const Vector& energies ) { return VectorVector( energies.size(), _scatteringXS ); }
VectorVector Checkerboard::GetTotalXS( const std::vector<double>& energies ) { return VectorVector( energies.size(), _totalXS ); }
VectorVector Checkerboard::GetTotalXS( const Vector& energies ) { return VectorVector( energies.size(), _totalXS ); }
std::vector<VectorVector> Checkerboard::GetExternalSource( const std::vector<double>& energies ) {
std::vector<VectorVector> Checkerboard::GetExternalSource( const Vector& energies ) {
VectorVector Q( _mesh->GetNumCells(), Vector( 1u, 0.0 ) );
auto cellMids = _mesh->GetCellMidPoints();
for( unsigned j = 0; j < cellMids.size(); ++j ) {
......
......@@ -3,18 +3,17 @@
#include "common/mesh.h"
ElectronRT::ElectronRT( Config* settings, Mesh* mesh ) : ProblemBase( settings, mesh ) {
// @TODO
_physics = new Physics();
_physics = new Physics( settings->GetHydrogenFile(), settings->GetOxygenFile(), "../input/stopping_power.txt" ); // TODO
}
ElectronRT::~ElectronRT() {}
VectorVector ElectronRT::GetScatteringXS( const std::vector<double>& energies ) {
VectorVector ElectronRT::GetScatteringXS( const Vector& energies ) {
// @TODO
return VectorVector( energies.size(), Vector( _mesh->GetNumCells(), 0.0 ) );
}
VectorVector ElectronRT::GetTotalXS( const std::vector<double>& energies ) {
VectorVector ElectronRT::GetTotalXS( const Vector& energies ) {
// @TODO
return VectorVector( energies.size(), Vector( _mesh->GetNumCells(), 0.0 ) );
}
......@@ -32,7 +31,7 @@ Vector ElectronRT::GetTotalXSE( const Vector& energies ) {
return Vector( energies.size(), 0.0 );
}
std::vector<VectorVector> ElectronRT::GetExternalSource( const std::vector<double>& energies ) {
std::vector<VectorVector> ElectronRT::GetExternalSource( const Vector& energies ) {
// @TODO
return std::vector<VectorVector>( energies.size(), std::vector<Vector>( _mesh->GetNumCells(), Vector( _settings->GetNQuadPoints(), 0.0 ) ) );
}
......
......@@ -8,15 +8,11 @@ LineSource_SN::LineSource_SN( Config* settings, Mesh* mesh ) : ProblemBase( sett
LineSource_SN::~LineSource_SN(){};
VectorVector LineSource_SN::GetScatteringXS( const std::vector<double>& energies ) {
return VectorVector( energies.size(), Vector( _mesh->GetNumCells(), 1.0 ) );
}
VectorVector LineSource_SN::GetScatteringXS( const Vector& energies ) { return VectorVector( energies.size(), Vector( _mesh->GetNumCells(), 1.0 ) ); }
VectorVector LineSource_SN::GetTotalXS( const std::vector<double>& energies ) {
return VectorVector( energies.size(), Vector( _mesh->GetNumCells(), 1.0 ) );
}
VectorVector LineSource_SN::GetTotalXS( const Vector& energies ) { return VectorVector( energies.size(), Vector( _mesh->GetNumCells(), 1.0 ) ); }
std::vector<VectorVector> LineSource_SN::GetExternalSource( const std::vector<double>& energies ) {
std::vector<VectorVector> LineSource_SN::GetExternalSource( const Vector& energies ) {
return std::vector<VectorVector>( 1u, std::vector<Vector>( _mesh->GetNumCells(), Vector( 1u, 0.0 ) ) );
}
......@@ -64,15 +60,11 @@ LineSource_PN::LineSource_PN( Config* settings, Mesh* mesh ) : ProblemBase( sett
LineSource_PN::~LineSource_PN(){};
VectorVector LineSource_PN::GetScatteringXS( const std::vector<double>& energies ) {
return VectorVector( energies.size(), Vector( _mesh->GetNumCells(), 1.0 ) );
}
VectorVector LineSource_PN::GetScatteringXS( const Vector& energies ) { return VectorVector( energies.size(), Vector( _mesh->GetNumCells(), 1.0 ) ); }
VectorVector LineSource_PN::GetTotalXS( const std::vector<double>& energies ) {
return VectorVector( energies.size(), Vector( _mesh->GetNumCells(), 1.0 ) );
}
VectorVector LineSource_PN::GetTotalXS( const Vector& energies ) { return VectorVector( energies.size(), Vector( _mesh->GetNumCells(), 1.0 ) ); }
std::vector<VectorVector> LineSource_PN::GetExternalSource( const std::vector<double>& energies ) {
std::vector<VectorVector> LineSource_PN::GetExternalSource( const Vector& energies ) {
return std::vector<VectorVector>( 1u, std::vector<Vector>( _mesh->GetNumCells(), Vector( 1u, 0.0 ) ) );
}
......
......@@ -3,13 +3,12 @@
#include "common/mesh.h"
WaterPhantom::WaterPhantom( Config* settings, Mesh* mesh ) : ElectronRT( settings, mesh ) {
// @TODO get pointer to correct physics class
_physics = nullptr;
_physics = new Physics( settings->GetHydrogenFile(), settings->GetOxygenFile(), "../input/stopping_power.txt" ); // TODO
}
WaterPhantom::~WaterPhantom() {}
WaterPhantom::~WaterPhantom() { delete _physics; }
std::vector<VectorVector> WaterPhantom::GetExternalSource( const std::vector<double>& energies ) {
std::vector<VectorVector> WaterPhantom::GetExternalSource( const Vector& energies ) {
return std::vector<VectorVector>( energies.size(), std::vector<Vector>( _mesh->GetNumCells(), Vector( _settings->GetNQuadPoints(), 0.0 ) ) );
}
......
......@@ -30,7 +30,8 @@ Solver::Solver( Config* settings ) : _settings( settings ) {
// set time step
_dE = ComputeTimeStep( settings->GetCFL() );
_nEnergies = unsigned( settings->GetTEnd() / _dE );
for( unsigned i = 0; i < _nEnergies; ++i ) _energies.push_back( ( i + 1 ) * _dE );
_energies.resize( _nEnergies );
for( unsigned i = 0; i < _nEnergies; ++i ) _energies[i] = ( i + 1 ) * _dE;
// setup problem and store frequently used params
_problem = ProblemBase::Create( _settings, _mesh );
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment