Commit a0c042d8 authored by Steffen Schotthöfer's avatar Steffen Schotthöfer
Browse files

Merge branch 'sprint1' of https://git.scc.kit.edu/rtsn/rtsn into sprint1

parents 961255cb dc93fa76
......@@ -75,7 +75,8 @@ enum PROBLEM_NAME {
PROBLEM_LineSource_Pseudo_1D,
PROBLEM_LineSource_Pseudo_1D_Physics,
PROBLEM_AirCavity,
PROBLEM_MuscleBoneLung
PROBLEM_MuscleBoneLung,
PROBLEM_Phantom2D
};
inline std::map<std::string, PROBLEM_NAME> Problem_Map{ { "LINESOURCE", PROBLEM_LineSource },
......@@ -84,6 +85,7 @@ inline std::map<std::string, PROBLEM_NAME> Problem_Map{ { "LINESOURCE", PROBLEM_
{ "WATERPHANTOM", PROBLEM_WaterPhantom },
{ "AIRCAVITY", PROBLEM_AirCavity },
{ "MUSCLEBONELUNG", PROBLEM_MuscleBoneLung },
{ "PHANTOM2D", PROBLEM_Phantom2D},
{ "LINESOURCE_PSEUDO_1D", PROBLEM_LineSource_Pseudo_1D },
{ "LINESOURCE_PSEUDO_1D_PHYSICS", PROBLEM_LineSource_Pseudo_1D_Physics } };
......
......@@ -54,7 +54,7 @@ class ICRU
double _BETA2;
double _R1;
Vector _E, _QMU; // input energy and mu vectors
Vector _E, _QMU; // input energy and mu vectors
std::vector<double> _ET, _ETL;
std::vector<double> _XMU, _XMUL;
std::vector<double> _ECS, _PCS;
......@@ -85,6 +85,7 @@ class ICRU
double DODCSC( double RMU, double EL );
void DOICSS( double WR, double E, double RMU, double& DCSD, double& DCSC );
void angdcs( unsigned ncor, double e, std::vector<double>& dxse, std::vector<double>& dxsi, std::vector<double>& dxs );
void angdcs( unsigned ncor, std::vector<std::vector<double>>& dxs );
ICRU() = delete;
......@@ -143,7 +144,7 @@ class ICRU
void GetAngularScatteringXS( Matrix& angularXS, Vector& integratedXS );
void GetStoppingPower( Vector& stoppingPower );
void GetTransportCoefficients(Matrix& xi);
void GetTransportCoefficients( Matrix& xi );
};
#endif // ICRU_H
#ifndef PHANTOM2D
#define PHANTOM2D
#include "electronrt.h"
#include <fstream>
#include <numeric>
#include "common/config.h"
#include "common/io.h"
#include "common/mesh.h"
#include "interpolation.h"
class Phantom2D : public ElectronRT
{
private:
Phantom2D() = delete;
public:
Phantom2D( Config* settings, Mesh* mesh );
virtual ~Phantom2D();
virtual std::vector<VectorVector> GetExternalSource( const Vector& energies );
virtual VectorVector SetupIC();
std::vector<double> GetDensity( const VectorVector& cellMidPoints );
};
#endif // PHANTOM2D
\ No newline at end of file
......@@ -14,7 +14,6 @@ class CSDSNSolver : public SNSolver
// Physics acess
Vector _energies; /*! @brief: energy levels for CSD, lenght = _nEnergies */
Vector _angle; /*! @brief: angles for SN */
std::vector<double> _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*/
......
......@@ -14,12 +14,11 @@ class CSDSNSolverFP : public SNSolver
// Physics acess
Vector _energies; /*! @brief: energy levels for CSD, lenght = _nEnergies */
Vector _angle; /*! @brief: angles for SN */
std::vector<double> _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 _L; /*! @brief Laplace Beltrami Matrix */
Matrix _IL; /*! @brief Laplace Beltrami Matrix */
double _alpha;
......
......@@ -14,7 +14,6 @@ class CSDSNSolverNoTrafo : public SNSolver
// Physics acess
Vector _energies; /*! @brief: energy levels for CSD, lenght = _nEnergies */
Vector _angle; /*! @brief: angles for SN */
std::vector<double> _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*/
......
......@@ -14,7 +14,6 @@ class CSDSolverTrafoFP : public SNSolver
// Physics acess
Vector _energies; /*! @brief: energy levels for CSD, lenght = _nEnergies */
Vector _angle; /*! @brief: angles for SN */
std::vector<double> _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*/
......@@ -30,6 +29,8 @@ class CSDSolverTrafoFP : public SNSolver
Vector _xi2;
Matrix _xi;
unsigned _FPMethod;
bool _RT;
double _energyMin;
......
OUTPUT_DIR = ../result
OUTPUT_FILE = aircavity_csd_1d_10MeV_fine
OUTPUT_FILE = example_csd_1d_10MeV_fine
LOG_DIR = ../result/logs
MESH_FILE = 1DMesh.su2
PROBLEM = AIRCAVITY
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.0001
TIME_FINAL = 1.0
CLEAN_FLUX_MATRICES = NO
BC_DIRICHLET = ( dirichlet )
BC_NEUMANN = ( wall_low, wall_up )
QUAD_TYPE = GAUSS_LEGENDRE_1D
QUAD_ORDER = 32
QUAD_ORDER = 26
......@@ -80,6 +80,95 @@ void ICRU::angdcs( unsigned ncor, double e, std::vector<double>& dxse, std::vect
dxs[IA] = dxse[IA] + dxsi[IA];
}
}
void ICRU::angdcs( unsigned ncor, std::vector<std::vector<double>>& dxs ) {
unsigned NM = materialID;
if( NM < 1 || NM > 280 ) {
ErrorMessages::Error( "The material number " + std::to_string( NM ) + " is not defined.", CURRENT_FUNCTION );
}
std::string PATHMT = "../data/material/";
std::stringstream ss;
ss << std::setw( 3 ) << std::setfill( '0' ) << materialID;
std::string FILE1 = PATHMT + "mater" + ss.str() + ".stp";
std::ifstream f( FILE1 );
std::string line;
for( unsigned i = 0; i < 5; ++i ) std::getline( f, line );
// double RHO = std::stod( line.substr( line.find_first_of( "=" ) + 1, 16 ) );
std::getline( f, line );
_NELEM = std::stoul( line.substr( line.find_first_of( "=" ) + 1, 3 ) );
double ATW = 0.0, ZTOT = 0.0;
_IZ.resize( _NELEM );
_STF.resize( _NELEM );
for( unsigned IEL = 0; IEL < _NELEM; ++IEL ) {
std::getline( f, line );
_IZ[IEL] = std::stoul( line.substr( line.find_first_of( "=" ) + 1, 3 ) );
_STF[IEL] = std::stod( line.substr( line.find_last_of( "=" ) + 1, 3 ) );
ZTOT += _STF[IEL] * _IZ[IEL];
ATW += ELAW[_IZ[IEL]] * _STF[IEL];
}
std::getline( f, line );
std::getline( f, line );
double EXPOT = std::stod( line.substr( line.find_first_of( "=" ) + 1, 16 ) );
std::getline( f, line );
_NOSC = std::stoul( line.substr( line.find_first_of( "=" ) + 1, line.back() ) );
std::getline( f, line );
std::getline( f, line );
std::getline( f, line );
double EXPOTR = 0.0;
_F.resize( _NOSC );
_WRI.resize( _NOSC );
double buffer;
for( unsigned I = 0; I < _NOSC; ++I ) {
f >> buffer >> _F[I] >> buffer >> _WRI[I];
EXPOTR += _F[I] * std::log( _WRI[I] );
}
f.close();
EXPOTR = std::exp( EXPOTR / ZTOT );
if( std::fabs( EXPOT - EXPOTR ) > 1.0e-5 * EXPOT ) {
ErrorMessages::Error( "Inconsistent oscillator data", CURRENT_FUNCTION );
}
unsigned NDAT = _NA;
ELINIT( _IZ, _STF );
std::vector<double> dxse;
std::vector<double> dxsi;
dxs.resize( _E.size() );
for( unsigned n = 0; n < _E.size(); ++n ) {
DCSEL0( _E[n] );
dxse.resize( NDAT );
dxsi.resize( NDAT );
dxs[n].resize( NDAT );
for( unsigned IA = 0; IA < NDAT; ++IA ) {
dxse[IA] = DCSEL( _XMU[IA] );
}
if( ncor == 1u || ncor == 2u ) { // 1: Inelastic DCS = elastic DCS/Z | 2: Morse's incoherent scattering function.
ErrorMessages::Error( "Unsupported scattering model", CURRENT_FUNCTION );
}
else if( ncor == 3u ) { // Sternheimer -Liljequist GOS model.
for( unsigned IA = 0; IA < NDAT; ++IA ) {
dxsi[IA] = DCSIN3( _E[n], _XMU[IA] );
}
}
else { // No inelastic scattering correction
for( unsigned IA = 0; IA < NDAT; ++IA ) {
dxsi[IA] = 0.0e0;
}
}
for( unsigned IA = 0; IA < _NA; ++IA ) {
dxs[n][IA] = dxse[IA] + dxsi[IA];
}
}
}
void ICRU::ELINIT( std::vector<double> IZ, std::vector<double> STF ) {
unsigned IE = 0;
double FGRID = 10.0e0;
......@@ -168,7 +257,7 @@ void ICRU::ELINIT( std::vector<double> IZ, std::vector<double> STF ) {
void ICRU::DCSEL0( double E ) {
if( E < 49.9999e0 || E > 1.0001e8 ) {
std::cout<<"Energy "<<E<<std::endl;
std::cout << "Energy " << E << std::endl;
ErrorMessages::Error( "Outside of supported energy range", CURRENT_FUNCTION );
}
double EL = std::log( E );
......@@ -450,6 +539,8 @@ double ICRU::DODCSC( double RMU, double EL ) {
void ICRU::GetAngularScatteringXS( Matrix& angularXS, Vector& integratedXS ) {
angularXS.resize( _QMU.size(), _E.size() );
Vector mu( _XMU.size() );
for( unsigned n = 0; n < mu.size(); ++n ) mu[n] = 1.0 - 2.0 * _XMU[n];
integratedXS.resize( _E.size() );
for( unsigned i = 0; i < _E.size(); ++i ) {
std::vector<double> dxse, dxsi, dxs;
......@@ -465,24 +556,23 @@ void ICRU::GetAngularScatteringXS( Matrix& angularXS, Vector& integratedXS ) {
}
void ICRU::GetTransportCoefficients( Matrix& xi ) {
Vector mu(_XMU.size());
for( unsigned n = 0; n < mu.size(); ++n ) mu[n] = 1.0-2.0*_XMU[n]; // mu starts at 1 and goes to 0
Vector mu( _XMU.size() );
std::vector<std::vector<double>> dxs;
// get scattering cross sections for all energies E
angdcs( 3u, dxs );
for( unsigned n = 0; n < mu.size(); ++n ) mu[n] = 1.0 - 2.0 * _XMU[n]; // mu starts at 1 and goes to 0
#pragma omp parallel for
for( unsigned i = 0; i < _E.size(); ++i ) {
std::vector<double> dxse, dxsi, dxs;
// get scattering cross sections for energy E[i]
angdcs( 3u, _E[i], dxse, dxsi, dxs );
// compute moments with trapezoidal rule
for( unsigned n = 0; n < xi.rows(); ++n ){
xi(n,i) = 0.0;
for( unsigned n = 0; n < xi.rows(); ++n ) {
xi( n, i ) = 0.0;
// integration over mu
for( unsigned k = 0; k < dxs.size()-1; ++k ){
xi(n,i) -= 0.5 * ( pow(1.0-mu[k+1],n)*dxs[k+1] + pow(1.0-mu[k],n)*dxs[k] )*(mu[k+1]-mu[k]);
for( unsigned k = 0; k < dxs[i].size() - 1; ++k ) {
xi( n, i ) -= 0.5 * ( pow( 1.0 - mu[k + 1], n ) * dxs[i][k + 1] + pow( 1.0 - mu[k], n ) * dxs[i][k] ) * ( mu[k + 1] - mu[k] );
}
}
}
xi *= 2.0*PI*H2OMolecularDensity;
xi *= 2.0 * PI * H2OMolecularDensity;
}
void ICRU::GetStoppingPower( Vector& stoppingPower ) {
......
......@@ -21,10 +21,16 @@ VectorVector AirCavity1D::SetupIC() {
return psi;
}
std::vector<double> AirCavity1D::GetDensity( const VectorVector& cellMidPoints ) {
std::vector<double> densities ( 4*cellMidPoints.size()/9, 1.0 );
std::vector<double> air( 2*cellMidPoints.size()/9, 0.00125 );
std::vector<double> water ( 3*cellMidPoints.size()/9, 1.0 );
densities.insert(densities.end(),air.begin(),air.end());
densities.insert(densities.end(),water.begin(),water.end());
return densities; }
std::vector<double> AirCavity1D::GetDensity( const VectorVector& cellMidPoints ) {
std::vector<double> densities( cellMidPoints.size(), 1.0 );
/*std::vector<double> densities( 4 * cellMidPoints.size() / 9, 1.0 );
std::vector<double> air( 2 * cellMidPoints.size() / 9, 0.00125 );
std::vector<double> water( 3 * cellMidPoints.size() / 9, 1.0 );
densities.insert( densities.end(), air.begin(), air.end() );
densities.insert( densities.end(), water.begin(), water.end() );
return densities;*/
for( unsigned j = 0; j < cellMidPoints.size(); ++j ) {
if( cellMidPoints[j][0] > 1.5 - 2.5 && cellMidPoints[j][0] < 2.0 - 2.5 ) densities[j] = 0.01;
}
return densities;
}
......@@ -6,7 +6,7 @@ ElectronRT::ElectronRT( Config* settings, Mesh* mesh ) : ProblemBase( settings,
_physics = new Physics( settings->GetHydrogenFile(), settings->GetOxygenFile(), "../input/stopping_power.txt" );
}
ElectronRT::~ElectronRT() {}
ElectronRT::~ElectronRT() { delete _physics; }
VectorVector ElectronRT::GetScatteringXS( const Vector& energies ) {
// @TODO
......
......@@ -23,9 +23,9 @@ VectorVector MuscleBoneLung::SetupIC() {
std::vector<double> MuscleBoneLung::GetDensity( const VectorVector& cellMidPoints ) {
std::cout<<"Length of mesh "<< cellMidPoints.size() << std::endl;
std::vector<double> densities ( 167, 1.05 ); //muscle layer
std::vector<double> bone( 167, 1.92 );
std::vector<double> lung ( 665, 0.26 );
std::vector<double> densities ( 167, 1.04 ); //muscle layer
std::vector<double> bone( 167, 1.85 );
std::vector<double> lung ( 665, 0.3 ); //maybe this is just the lung tissue and after that there should be air?
densities.insert(densities.end(),bone.begin(),bone.end());
densities.insert(densities.end(),lung.begin(),lung.end());
std::cout<<"Length of densities "<< densities.size() << std::endl;
......
#include "problems/phantom2d.h"
#include "common/config.h"
#include "common/mesh.h"
Phantom2D::Phantom2D( Config* settings, Mesh* mesh ) : ElectronRT( settings, mesh ) {}
Phantom2D::~Phantom2D() { delete _physics; }
std::vector<VectorVector> Phantom2D::GetExternalSource( const Vector& energies ) {
return std::vector<VectorVector>( energies.size(), std::vector<Vector>( _mesh->GetNumCells(), Vector( _settings->GetNQuadPoints(), 0.0 ) ) );
}
VectorVector Phantom2D::SetupIC() {
VectorVector psi( _mesh->GetNumCells(), Vector( _settings->GetNQuadPoints(), 1e-10 ) );
auto cellMids = _mesh->GetCellMidPoints();
double s = 0.1;
for( unsigned j = 0; j < cellMids.size(); ++j ) {
double x = cellMids[j][0];
psi[j][_settings->GetNQuadPoints() - 1] = 1.0 / ( s * sqrt( 2 * M_PI ) ) * std::exp( -x * x / ( 2 * s * s ) );
}
return psi;
}
std::vector<double> Phantom2D::GetDensity( const VectorVector& cellMidPoints ) {
std::string imageFile = "../tests/input/phantom.png";
std::string meshFile = _settings->GetMeshFile();
Matrix gsImage = createSU2MeshFromImage( imageFile, meshFile );
auto bounds = _mesh->GetBounds();
double xMin = bounds[0].first;
double xMax = bounds[0].second;
double yMin = bounds[1].first;
double yMax = bounds[1].second;
unsigned m = gsImage.rows();
unsigned n = gsImage.columns();
Vector x( m + 1 ), y( n + 1 );
for( unsigned i = 0; i < m + 1; ++i ) {
x[i] = static_cast<double>( i ) / static_cast<double>( m ) * ( xMax - xMin );
}
for( unsigned i = 0; i < n + 1; ++i ) y[i] = static_cast<double>( i ) / static_cast<double>( n ) * ( yMax - yMin );
Interpolation interp( x, y, gsImage );
std::vector<double> result( _mesh->GetNumCells(), 0.0 );
for( unsigned i = 0; i < _mesh->GetNumCells(); ++i ) {
result[i] = std::clamp( interp( cellMidPoints[i][0], cellMidPoints[i][1] ), 0.0, 1.0 );
}
}
\ No newline at end of file
#include "common/config.h"
#include "problems/aircavity1d.h"
#include "problems/checkerboard.h"
#include "problems/electronrt.h"
#include "problems/linesource.h"
#include "problems/musclebonelung.h"
#include "problems/problembase.h"
#include "problems/waterphantom.h"
#include "problems/aircavity1d.h"
#include "problems/musclebonelung.h"
#include "problems/phantom2d.h"
ProblemBase::ProblemBase( Config* settings, Mesh* mesh ) {
_settings = settings;
......@@ -34,6 +35,7 @@ ProblemBase* ProblemBase::Create( Config* settings, Mesh* mesh ) {
case PROBLEM_AirCavity: return new AirCavity1D( settings, mesh );
case PROBLEM_MuscleBoneLung: return new MuscleBoneLung( settings, mesh );
case PROBLEM_WaterPhantom: return new WaterPhantom( settings, mesh );
case PROBLEM_Phantom2D: return new Phantom2D( settings, mesh );
case PROBLEM_LineSource_Pseudo_1D: return new LineSource_SN_Pseudo1D( settings, mesh );
case PROBLEM_LineSource_Pseudo_1D_Physics: return new LineSource_SN_Pseudo1D_Physics( settings, mesh );
default: return new ElectronRT( settings, mesh ); // Use RadioTherapy as dummy
......
......@@ -4,20 +4,14 @@
WaterPhantom::WaterPhantom( Config* settings, Mesh* mesh ) : ElectronRT( settings, mesh ) {}
WaterPhantom::~WaterPhantom() { delete _physics; }
WaterPhantom::~WaterPhantom() {}
std::vector<VectorVector> WaterPhantom::GetExternalSource( const Vector& energies ) {
return std::vector<VectorVector>( energies.size(), std::vector<Vector>( _mesh->GetNumCells(), Vector( _settings->GetNQuadPoints(), 0.0 ) ) );
}
VectorVector WaterPhantom::SetupIC() {
VectorVector psi( _mesh->GetNumCells(), Vector( _settings->GetNQuadPoints(), 1e-10 ) );
auto cellMids = _mesh->GetCellMidPoints();
double s = 0.1;
for( unsigned j = 0; j < cellMids.size(); ++j ) {
double x = cellMids[j][0];
psi[j][_settings->GetNQuadPoints() - 1] = 1.0 / ( s * sqrt( 2 * M_PI ) ) * std::exp( -x * x / ( 2 * s * s ) );
}
VectorVector psi( _mesh->GetNumCells(), Vector( _settings->GetNQuadPoints(), 0.0 ) );
return psi;
}
......
......@@ -16,10 +16,10 @@ CSDSolverTrafoFP::CSDSolverTrafoFP( Config* settings ) : SNSolver( settings ) {
// Set angle and energies
_energies = Vector( _nEnergies, 0.0 ); // equidistant
_energyMin = 1e-4 * 0.511;
_energyMax = 5e0;
_energyMax = 10e0;
// write equidistant energy grid (false) or refined grid (true)
GenerateEnergyGrid( true );
GenerateEnergyGrid( false );
// create 1D quadrature
unsigned nq = _settings->GetNQuadPoints();
......@@ -34,6 +34,8 @@ CSDSolverTrafoFP::CSDSolverTrafoFP( Config* settings ) : SNSolver( settings ) {
// setup Laplace Beltrami matrix L in slab geometry
_L = Matrix( nq, nq, 0.0 );
_FPMethod = 2;
double DMinus = 0.0;
double DPlus = 0.0;
for( unsigned k = 0; k < nq; ++k ) {
......@@ -61,6 +63,7 @@ CSDSolverTrafoFP::CSDSolverTrafoFP( Config* settings ) : SNSolver( settings ) {
_xi( 2, n ) = 4.0 / 3.0 - 2.0 * g + 2.0 / 3.0 * g * g;
}
// initialize stopping power vector
_s = Vector( _nEnergies, 1.0 );
_RT = true;
......@@ -81,50 +84,44 @@ CSDSolverTrafoFP::CSDSolverTrafoFP( Config* settings ) : SNSolver( settings ) {
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;
// 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)<<"; ";
}
_scatteringKernel( p, p ) = _weights[p];
std::cout<<"];"<<std::endl;
*/
}
_density = std::vector<double>( _nCells, 1.0 );
// exit(EXIT_SUCCESS);
_density = std::vector<double> ( _nCells, 1.0 );
//exit(EXIT_SUCCESS);
}
void CSDSolverTrafoFP::Solve() {
std::cout << "Solve" << std::endl;
auto log = spdlog::get( "event" );
// save original energy field for boundary conditions
auto energiesOrig = _energies;
// setup IC and incoming BC on left
// auto cellMids = _settings->GetCellMidPoints();
_sol = std::vector<Vector>( _nCells, Vector( _nq, 0.0 ) );
// setup incoming BC on left
_sol = VectorVector( _density.size(), Vector( _settings->GetNQuadPoints(), 0.0 ) ); // hard coded IC, needs to be changed
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 ) );
}
// hard coded boundary type for 1D testcases (otherwise cells will be NEUMANN)
_boundaryCells[0] = BOUNDARY_TYPE::DIRICHLET;
_boundaryCells[_nCells - 1] = BOUNDARY_TYPE::DIRICHLET;
// setup identity matrix for FP scattering
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?)
// angular flux at next time step
VectorVector psiNew( _nCells, Vector( _nq, 0.0 ) );
double dFlux = 1e10;
Vector fluxNew( _nCells, 0.0 );
......@@ -144,8 +141,6 @@ void CSDSolverTrafoFP::Solve() {
}
}
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;
......@@ -175,23 +170,22 @@ void CSDSolverTrafoFP::Solve() {
double xi2 = _xi( 2, _nEnergies - n - 1 );
double xi3 = _xi( 3, _nEnergies - n - 1 );
// FP1
/*
_alpha = 0.0;
_alpha2 = xi1 / 2.0;
_beta = 0.0;
*/
// FP2
_alpha = xi1 / 2.0 + xi2 / 8.0;
_alpha2 = 0.0;
_beta = xi2 / 8.0 / xi1;
/*
// FP3
_alpha = xi2 * ( 27.0 * xi2 * xi2 + 5.0 * xi3 * xi3 - 24.0 * xi2 * xi3 ) / ( 8.0 * xi3 * ( 3.0 * xi2 - 2.0 * xi3 ) );
_beta = xi3 / ( 6.0 * ( 3.0 * xi2 - 2.0 * xi3 ) );
_alpha2 = xi1 / 2.0 - 9.0 / 8.0 * xi2 * xi2 / xi3 + 3.0 / 8.0 * xi2;
*/
// setup coefficients in FP step
if( _FPMethod == 1 ) {
_alpha = 0.0;
_alpha2 = xi1 / 2.0;
_beta = 0.0;
}
else if( _FPMethod == 2 ) {
_alpha = xi1 / 2.0 + xi2 / 8.0;
_alpha2 = 0.0;
_beta = xi2 / 8.0 / xi1;
}
else if( _FPMethod == 3 ) {
_alpha = xi2 * ( 27.0 * xi2 * xi2 + 5.0 * xi3 * xi3 - 24.0 * xi2 * xi3 ) / ( 8.0 * xi3 * ( 3.0 * xi2 - 2.0 * xi3 ) );
_beta = xi3 / ( 6.0 * ( 3.0 * xi2 - 2.0 * xi3 ) );
_alpha2 = xi1 / 2.0 - 9.0 / 8.0 * xi2 * xi2 / xi3 + 3.0 / 8.0 * xi2;
}