Commit 233152ab 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 1f67c0e9 a1509ffa
Pipeline #117971 failed with stages
in 28 minutes and 39 seconds
......@@ -55,12 +55,13 @@ enum BOUNDARY_TYPE { DIRICHLET, NEUMANN, NONE, INVALID };
/*! @brief Enum for all currently available quadratures in rtsn.
* Option enums are written in capital letters with underscores as spaces (e.g option "time integration" has option enum "TIME_INTEGRATION")
*/
enum QUAD_NAME { QUAD_MonteCarlo, QUAD_GaussLegendreTensorized, QUAD_GaussLegendre1D, QUAD_LevelSymmetric, QUAD_Lebedev, QUAD_LDFESA };
enum QUAD_NAME { QUAD_MonteCarlo, QUAD_GaussLegendreTensorized, QUAD_GaussLegendre1D, QUAD_LevelSymmetric, QUAD_Lebedev, QUAD_LDFESA, QUAD_Product };
/*! @brief Conversion Map String to enum
*/
inline std::map<std::string, QUAD_NAME> Quadrature_Map{ { "MONTE_CARLO", QUAD_MonteCarlo },
{ "GAUSS_LEGENDRE_TENSORIZED", QUAD_GaussLegendreTensorized },
{ "PRODUCT", QUAD_Product },
{ "GAUSS_LEGENDRE_1D", QUAD_GaussLegendre1D },
{ "LEVEL_SYMMETRIC", QUAD_LevelSymmetric },
{ "LEBEDEV", QUAD_Lebedev },
......@@ -85,7 +86,7 @@ inline std::map<std::string, PROBLEM_NAME> Problem_Map{ { "LINESOURCE", PROBLEM_
{ "WATERPHANTOM", PROBLEM_WaterPhantom },
{ "AIRCAVITY", PROBLEM_AirCavity },
{ "MUSCLEBONELUNG", PROBLEM_MuscleBoneLung },
{ "PHANTOM2D", PROBLEM_Phantom2D},
{ "PHANTOM2D", PROBLEM_Phantom2D },
{ "LINESOURCE_PSEUDO_1D", PROBLEM_LineSource_Pseudo_1D },
{ "LINESOURCE_PSEUDO_1D_PHYSICS", PROBLEM_LineSource_Pseudo_1D_Physics } };
......@@ -95,10 +96,25 @@ 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, CSD_SN_FOKKERPLANCK_TRAFO_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,
CSD_SN_FOKKERPLANCK_TRAFO_SOLVER_2D,
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 }, { "CSD_SN_FOKKERPLANCK_TRAFO_SOLVER", CSD_SN_FOKKERPLANCK_TRAFO_SOLVER }, { "PN_SOLVER", PN_SOLVER }, { "MN_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 },
{ "CSD_SN_FOKKERPLANCK_TRAFO_SOLVER", CSD_SN_FOKKERPLANCK_TRAFO_SOLVER },
{ "CSD_SN_FOKKERPLANCK_TRAFO_SOLVER_2D", CSD_SN_FOKKERPLANCK_TRAFO_SOLVER_2D },
{ "PN_SOLVER", PN_SOLVER },
{ "MN_SOLVER", MN_SOLVER } };
// Entropy functional
enum ENTROPY_NAME { QUADRATIC, MAXWELL_BOLZMANN, BOSE_EINSTEIN, FERMI_DIRAC };
......
#ifndef PRODUCTQUADRATURE_H
#define PRODUCTQUADRATURE_H
#include "quadraturebase.h"
class ProductQuadrature : public QuadratureBase
{
// Implementation is done accordingly to Kendall Atkinson 1981, Australian Matematical Society.
private:
double Pythag( const double a, const double b );
std::pair<Vector, Matrix> ComputeEigenValTriDiagMatrix( const Matrix& mat );
bool CheckOrder();
public:
ProductQuadrature( unsigned order );
virtual ~ProductQuadrature() {}
inline void SetName() override { _name = "Product quadrature"; }
inline void SetNq() override { _nq = 4 * pow( GetOrder(), 2 ); }
void SetPointsAndWeights() override;
void SetConnectivity() override;
inline VectorVector GetPointsSphere() const override { return _pointsSphere; }
};
#endif // PRODUCTQUADRATURE_H
#ifndef CSDSOLVERTRAFOFP2D_H
#define CSDSOLVERTRAFOFP2D_H
#include "icru.h"
#include "solvers/snsolver.h"
class Physics;
class CSDSolverTrafoFP2D : 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 */
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 */
VectorVector _quadPoints;
VectorVector _quadPointsSphere;
Vector _weights;
Vector _mu;
Vector _phi;
Vector _wp;
Vector _wa;
double _alpha;
double _alpha2;
double _beta;
Vector _xi1;
Vector _xi2;
Matrix _xi;
unsigned _FPMethod;
bool _RT;
double _energyMin;
double _energyMax;
void GenerateEnergyGrid( bool refinement );
public:
/**
* @brief CSDSolverTrafoFP2D constructor
* @param settings stores all needed information
*/
CSDSolverTrafoFP2D( 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 // CSDSOLVERTRAFOFP2D_H
#include "quadratures/productquadrature.h"
#include "toolboxes/errormessages.h"
ProductQuadrature::ProductQuadrature( unsigned order ) : QuadratureBase( order ) {
SetName();
SetNq();
SetPointsAndWeights();
SetConnectivity();
}
void ProductQuadrature::SetPointsAndWeights() {
Vector nodes1D( 2 * _order ), weights1D( 2 * _order );
// construct companion matrix
Matrix CM( 2 * _order, 2 * _order, 0.0 );
for( unsigned i = 0; i < 2 * _order - 1; ++i ) {
CM( i + 1, i ) = std::sqrt( 1 / ( 4 - 1 / std::pow( static_cast<double>( i + 1 ), 2 ) ) );
CM( i, i + 1 ) = std::sqrt( 1 / ( 4 - 1 / std::pow( static_cast<double>( i + 1 ), 2 ) ) );
}
// compute eigenvalues and -vectors of the companion matrix
auto evSys = ComputeEigenValTriDiagMatrix( CM );
for( unsigned i = 0; i < 2 * _order; ++i ) {
if( std::fabs( evSys.first[i] ) < 1e-15 ) // avoid rounding errors
nodes1D[i] = 0;
else
nodes1D[i] = evSys.first[i];
weights1D[i] = 2 * std::pow( evSys.second( 0, i ), 2 );
}
// sort nodes increasingly and also reorder weigths for consistency
std::vector<unsigned> sortOrder( nodes1D.size() );
std::iota( sortOrder.begin(), sortOrder.end(), 0 );
std::sort( sortOrder.begin(), sortOrder.end(), [&]( unsigned i, unsigned j ) { return nodes1D[i] < nodes1D[j]; } );
Vector sorted_nodes( static_cast<unsigned>( sortOrder.size() ) ), sorted_weights( static_cast<unsigned>( sortOrder.size() ) );
std::transform( sortOrder.begin(), sortOrder.end(), sorted_nodes.begin(), [&]( unsigned i ) { return nodes1D[i]; } );
std::transform( sortOrder.begin(), sortOrder.end(), sorted_weights.begin(), [&]( unsigned i ) { return weights1D[i]; } );
nodes1D = sorted_nodes;
weights1D = sorted_weights;
// setup equidistant angle phi around z axis
Vector phi( 2 * _order );
for( unsigned i = 0; i < 2 * _order; ++i ) {
phi[i] = ( i + 0.5 ) * M_PI / _order;
}
// unsigned range = std::floor( _order / 2.0 ); // comment (steffen): why do we only need half of the points:
//=> In 2D we would count everything twice. (not wrong with scaling
// resize points and weights
_points.resize( _nq );
_pointsSphere.resize( _nq );
for( auto& p : _points ) {
p.resize( 3 );
}
for( auto& p : _pointsSphere ) {
p.resize( 2 );
}
_weights.resize( _nq );
// transform tensorized (x,y,z)-grid to spherical grid points
for( unsigned j = 0; j < 2 * _order; ++j ) {
for( unsigned i = 0; i < 2 * _order; ++i ) {
_points[j * ( 2 * _order ) + i][0] = sqrt( 1 - nodes1D[j] * nodes1D[j] ) * std::cos( phi[i] );
_points[j * ( 2 * _order ) + i][1] = sqrt( 1 - nodes1D[j] * nodes1D[j] ) * std::sin( phi[i] );
_points[j * ( 2 * _order ) + i][2] = nodes1D[j];
_pointsSphere[j * ( 2 * _order ) + i][0] = nodes1D[j]; // my
_pointsSphere[j * ( 2 * _order ) + i][1] = phi[i]; // phi
_weights[j * ( 2 * _order ) + i] = M_PI / _order * weights1D[j];
}
}
}
void ProductQuadrature::SetConnectivity() { // TODO
// Not initialized for this quadrature.
VectorVectorU connectivity;
_connectivity = connectivity;
}
std::pair<Vector, Matrix> ProductQuadrature::ComputeEigenValTriDiagMatrix( const Matrix& mat ) {
// copied from 'Numerical Recipes' and updated + modified to work with blaze
unsigned n = mat.rows();
Vector d( n, 0.0 ), e( n, 0.0 );
Matrix z( n, n, 0.0 );
for( unsigned i = 0; i < n; ++i ) {
d[i] = mat( i, i );
z( i, i ) = 1.0;
i == 0 ? e[i] = 0.0 : e[i] = mat( i, i - 1 );
}
int m, l, iter, i, k;
m = l = iter = i = k = 0;
double s, r, p, g, f, dd, c, b;
s = r = p = g = f = dd = c = b = 0.0;
const double eps = std::numeric_limits<double>::epsilon();
for( i = 1; i < static_cast<int>( n ); i++ ) e[i - 1] = e[i];
e[n - 1] = 0.0;
for( l = 0; l < static_cast<int>( n ); l++ ) {
iter = 0;
do {
for( m = l; m < static_cast<int>( n ) - 1; m++ ) {
dd = std::fabs( d[m] ) + std::fabs( d[m + 1] );
if( std::fabs( e[m] ) <= eps * dd ) break;
}
if( m != l ) {
if( iter++ == 30 ) ErrorMessages::Error( "Solving the tridiagonal matrix took too many iterations!", CURRENT_FUNCTION );
g = ( d[l + 1] - d[l] ) / ( 2.0 * e[l] );
r = Pythag( g, 1.0 );
g = d[m] - d[l] + e[l] / ( g + std::copysign( r, g ) );
s = c = 1.0;
p = 0.0;
for( i = m - 1; i >= l; i-- ) {
f = s * e[i];
b = c * e[i];
e[i + 1] = ( r = Pythag( f, g ) );
if( r == 0.0 ) {
d[i + 1] -= p;
e[m] = 0.0;
break;
}
s = f / r;
c = g / r;
g = d[i + 1] - p;
r = ( d[i] - g ) * s + 2.0 * c * b;
d[i + 1] = g + ( p = s * r );
g = c * r - b;
for( k = 0; k < static_cast<int>( n ); k++ ) {
f = z( static_cast<unsigned>( k ), static_cast<unsigned>( i ) + 1 );
z( static_cast<unsigned>( k ), static_cast<unsigned>( i ) + 1 ) =
s * z( static_cast<unsigned>( k ), static_cast<unsigned>( i ) ) + c * f;
z( static_cast<unsigned>( k ), static_cast<unsigned>( i ) ) =
c * z( static_cast<unsigned>( k ), static_cast<unsigned>( i ) ) - s * f;
}
}
if( r == 0.0 && i >= l ) continue;
d[l] -= p;
e[l] = g;
e[m] = 0.0;
}
} while( m != l );
}
return std::make_pair( d, z );
}
double ProductQuadrature::Pythag( const double a, const double b ) {
// copied from 'Numerical Recipes'
double absa = std::fabs( a ), absb = std::fabs( b );
return ( absa > absb ? absa * std::sqrt( 1.0 + ( absb / absa ) * ( absb / absa ) )
: ( absb == 0.0 ? 0.0 : absb * std::sqrt( 1.0 + ( absa / absb ) * ( absa / absb ) ) ) );
}
bool ProductQuadrature::CheckOrder() {
if( _order % 2 == 1 ) { // order needs to be even
ErrorMessages::Error( "ERROR! Order " + std::to_string( _order ) + " for " + GetName() + " not available. \n Order must be an even number. ",
CURRENT_FUNCTION );
}
return true;
}
#include "quadratures/quadraturebase.h"
#include "quadratures/productquadrature.h"
#include "quadratures/qgausslegendre1D.h"
#include "quadratures/qgausslegendretensorized.h"
#include "quadratures/qldfesa.h"
......@@ -14,6 +15,7 @@ QuadratureBase* QuadratureBase::CreateQuadrature( QUAD_NAME name, unsigned order
switch( name ) {
case QUAD_MonteCarlo: return new QMonteCarlo( order );
case QUAD_GaussLegendreTensorized: return new QGaussLegendreTensorized( order );
case QUAD_Product: return new ProductQuadrature( order );
case QUAD_GaussLegendre1D: return new QGaussLegendre1D( order );
case QUAD_LevelSymmetric: return new QLevelSymmetric( order );
case QUAD_LDFESA: return new QLDFESA( order );
......
#include "solvers/csdsolvertrafofp2d.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>
CSDSolverTrafoFP2D::CSDSolverTrafoFP2D( Config* settings ) : SNSolver( settings ) {
_dose = std::vector<double>( _settings->GetNCells(), 0.0 );
// Set angle and energies
_energies = Vector( _nEnergies, 0.0 ); // equidistant
_energyMin = 1e-4 * 0.511;
_energyMax = 10e0;
// write equidistant energy grid (false) or refined grid (true)
GenerateEnergyGrid( false );
// create quadrature
unsigned order = _quadrature->GetOrder();
unsigned nq = _settings->GetNQuadPoints();
_quadPoints = _quadrature->GetPoints();
_weights = _quadrature->GetWeights();
_quadPointsSphere = _quadrature->GetPointsSphere();
// transform structured quadrature
_mu = Vector( 2 * order );
_phi = Vector( 2 * order );
_wp = Vector( 2 * order );
_wa = Vector( 2 * order );
// create quadrature 1D to compute mu grid
QuadratureBase* quad1D = QuadratureBase::CreateQuadrature( QUAD_GaussLegendre1D, 2 * order );
Vector w = quad1D->GetWeights();
VectorVector muVec = quad1D->GetPoints();
for( unsigned k = 0; k < 2 * order; ++k ) {
_mu[k] = muVec[k][0];
_wp[k] = w[k];
}
for( unsigned i = 0; i < 2 * order; ++i ) {
_phi[i] = ( i + 0.5 ) * M_PI / order;
_wa[i] = M_PI / order;
}
// setup Laplace Beltrami matrix L in slab geometry
_L = Matrix( nq, nq, 0.0 );
_FPMethod = 2;
double DMinus = 0.0;
double DPlus = 0.0;
Vector gamma( order, 0.0 );
double dPlus;
double c, K;
double dMinus = 0.0;
DPlus = DMinus - 2 * _mu[0] * w[0];
for( unsigned j = 0; j < 2 * order - 1; ++j ) {
DMinus = DPlus;
DPlus = DMinus - 2 * _mu[j] * w[j];
dPlus = ( sqrt( 1 - _mu[j + 1] * _mu[j + 1] ) - sqrt( 1 - _mu[j] * _mu[j] ) ) / ( _mu[j + 1] - _mu[j] );
c = ( DPlus * dPlus - DMinus * dMinus ) / _wp[j];
K = 2 * ( 1 - _mu[j] * _mu[j] ) + c * sqrt( 1 - _mu[j] * _mu[j] );
gamma[j] = M_PI * M_PI * K / ( 2 * order * ( 1 - std::cos( M_PI / order ) ) );
dMinus = dPlus;
}
DPlus = 0.0;
// implementation of 2D spherical Laplacian according to SN book, equation (1.136)
for( unsigned j = 0; j < 2 * order; ++j ) {
DMinus = DPlus;
DPlus = DMinus - 2 * _mu[j] * _wp[j];
for( unsigned i = 0; i < 2 * order; ++i ) {
if( j > 0 ) {
_L( j * 2 * order + i, ( j - 1 ) * 2 * order + i ) = DMinus / ( _mu[j] - _mu[j - 1] ) / _wp[j];
_L( j * 2 * order + i, j * 2 * order + i ) = -DMinus / ( _mu[j] - _mu[j - 1] ) / _wp[j];
}
if( i > 0 ) {
_L( j * 2 * order + i, j * 2 * order + i - 1 ) = 1.0 / ( 1 - _mu[j] * _mu[j] ) * gamma[j] / ( _phi[i] - _phi[i - 1] ) / _wa[i];
_L( j * 2 * order + i, j * 2 * order + i ) += -1.0 / ( 1 - _mu[j] * _mu[j] ) * gamma[j] / ( _phi[i] - _phi[i - 1] ) / _wa[i];
}
if( j < 2 * order - 1 ) {
_L( j * 2 * order + i, ( j + 1 ) * 2 * order + i ) = DPlus / ( _mu[j + 1] - _mu[j] ) / _wp[j];
_L( j * 2 * order + i, j * 2 * order + i ) += -DPlus / ( _mu[j + 1] - _mu[j] ) / _wp[j];
}
if( i < 2 * order - 1 ) {
_L( j * 2 * order + i, j * 2 * order + i + 1 ) = 1.0 / ( 1 - _mu[j] * _mu[j] ) * gamma[j] / ( _phi[i + 1] - _phi[i] ) / _wa[i];
_L( j * 2 * order + i, j * 2 * order + i ) += -1.0 / ( 1 - _mu[j] * _mu[j] ) * gamma[j] / ( _phi[i + 1] - _phi[i] ) / _wa[i];
}
}
}
// 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( 4, _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;
}
// initialize stopping power vector
_s = Vector( _nEnergies, 1.0 );
_RT = true;
// 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;
*/
}
_density = std::vector<double>( _nCells, 1.0 );
// exit(EXIT_SUCCESS);
}
void CSDSolverTrafoFP2D::Solve() {
auto log = spdlog::get( "event" );
// save original energy field for boundary conditions
auto energiesOrig = _energies;
// 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
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
}
}
// 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?
double xi1 = _xi( 1, _nEnergies - n - 1 );
double xi2 = _xi( 2, _nEnergies - n - 1 );
double xi3 = _xi( 3, _nEnergies - n - 1 );
// 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;
}
_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.0 * pow( _energyMax - energiesOrig[_nEnergies - n - 1], 2 ) ) * _density[0] * _s[_nEnergies - n - 1];
}
}
}