Commit b8ff24fb authored by steffen.schotthoefer's avatar steffen.schotthoefer
Browse files

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

parents 6304a7ee 180366b9
Pipeline #91042 passed with stages
in 18 minutes and 31 seconds
......@@ -74,9 +74,6 @@ class PNSolver : public Solver
void ComputeScatterMatrix();
// Computes Legedre polinomial of oder l at point x
double Legendre( double x, int l );
// Computes Moments for given zeroth moment initial conditions
void ComputeICMoments();
};
#endif // PNSOLVER_H
......@@ -3,14 +3,14 @@
#include "problembase.h"
class LineSource : public ProblemBase
class LineSource_SN : public ProblemBase
{
private:
LineSource() = delete;
LineSource_SN() = delete;
public:
LineSource( Config* settings, Mesh* mesh );
virtual ~LineSource();
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 );
......@@ -19,4 +19,30 @@ class LineSource : public ProblemBase
virtual VectorVector SetupIC();
};
class LineSource_PN : public ProblemBase
{
private:
LineSource_PN() = delete;
/**
* @brief Gets the global index for given order l of Legendre polynomials and given
* order k of Legendre functions.
* Note: This is code doubling from PNSolver::GlobalIndex
* @param l : order of Legendre polynomial
* @param k : order of Legendre function
* @returns global index
*/
int GlobalIndex( int l, int k ) const;
public:
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 std::vector<double> GetStoppingPower( const std::vector<double>& energies ) override;
virtual VectorVector SetupIC() override;
};
#endif // LINESOURCE_H
......@@ -21,26 +21,31 @@ class ProblemBase
public:
/**
* @brief GetScatteringXS gives back vector of vectors of scattering cross sections for materials defined by density and energies in vector energy
* @brief GetScatteringXS gives back vector (each energy) of vectors (each grid cell)
* of scattering cross sections for materials defined by density and energies
* in vector energy
* @param energy is the energy the cross section is queried for
*/
virtual VectorVector GetScatteringXS( const std::vector<double>& energies ) = 0;
/**
* @brief GetTotalXS gives back vector of vectors of total cross sections for materials defined by density and energies in vector energy
* @brief GetTotalXS gives back vector of vectors of total cross sections for
* materials defined by density and energies in vector energy
* @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;
/**
* @brief GetStoppingPower gives back vector of vectors of source terms for each energy, cell and angle
* @brief GetExternalSource gives back vector of vectors of source terms for each
* energy, cell and angle
* @param energies is vector with energies
*/
virtual std::vector<VectorVector> GetExternalSource( const std::vector<double>& energies ) = 0;
/**
* @brief GetStoppingPower gives back vector of vectors of stopping powers for materials defined by density and energies in vector energy
* @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
*/
virtual std::vector<double> GetStoppingPower( const std::vector<double>& energies ) = 0;
......
......@@ -11,7 +11,9 @@ OUTPUT_FILE = example
% Log directory
LOG_DIR = ../result/logs
% Mesh File
MESH_FILE = checkerboard.su2
MESH_FILE = linesource.su2
%
PROBLEM = LINESOURCE
%
% ---- Solver specifications ----
%
......
......@@ -247,6 +247,7 @@ void PNSolver::ComputeFluxComponents() {
}
void PNSolver::ComputeScatterMatrix() {
// right now independent of spatial coordinates
// G[i][i] = 1- g_l
// g_l = 2*pi*\int_-1 ^1 k(x,my)P_l(my)dmy
......@@ -260,9 +261,7 @@ void PNSolver::ComputeScatterMatrix() {
unsigned idx_global = 0;
for( int idx_lOrder = 0; idx_lOrder <= int( _nq ); idx_lOrder++ ) {
std::cout << "l: " << idx_lOrder << std::endl;
for( int idx_kOrder = -idx_lOrder; idx_kOrder <= idx_lOrder; idx_kOrder++ ) {
std::cout << "k :" << idx_kOrder << std::endl;
idx_global = unsigned( GlobalIndex( idx_lOrder, idx_kOrder ) );
_scatterMatDiag[idx_global] = 0;
......@@ -281,8 +280,6 @@ void PNSolver::ComputeScatterMatrix() {
}
}
void PNSolver::ComputeICMoments() {}
double PNSolver::Legendre( double x, int l ) {
// Pre computed low order polynomials for faster computation
switch( l ) {
......@@ -369,31 +366,30 @@ void PNSolver::Solve() {
psiNew[idx_cell][idx_system] = _psi[idx_cell][idx_system] /* value at last time */
- ( _dE / _areas[idx_cell] ) * psiNew[idx_cell][idx_system] /* cell averaged flux */
- _dE * _psi[idx_cell][idx_system] *
( _sigmaA[idx_energy][idx_cell] /* absorbtion influence */
+ _sigmaS[idx_energy][idx_cell] * _scatterMatDiag[idx_system] /* scattering influence */
);
( _sigmaA[idx_energy][idx_cell] /* absorbtion influence */
+ _sigmaS[idx_energy][idx_cell] * _scatterMatDiag[idx_system] ); /* scattering influence */
}
}
// TODO: Adapt to PN
// TODO: Add Source Term!!! (not neccessary in LineSource
// TODO: figure out a more elegant way
// add external source contribution
if( _Q.size() == 1u ) { // constant source for all energies
if( _Q[0][idx_cell].size() == 1u ) // isotropic source
psiNew[idx_cell] += _dE * _Q[0][idx_cell][0];
else
psiNew[idx_cell] += _dE * _Q[0][idx_cell];
}
else {
if( _Q[0][idx_cell].size() == 1u ) // isotropic source
psiNew[idx_cell] += _dE * _Q[idx_energy][idx_cell][0];
else
psiNew[idx_cell] += _dE * _Q[idx_energy][idx_cell];
}
// if( _Q.size() == 1u ) { // constant source for all energies
// if( _Q[0][idx_cell].size() == 1u ) // isotropic source
// psiNew[idx_cell] += _dE * _Q[0][idx_cell][0];
// else
// psiNew[idx_cell] += _dE * _Q[0][idx_cell];
//}
// else {
// if( _Q[0][idx_cell].size() == 1u ) // isotropic source
// psiNew[idx_cell] += _dE * _Q[idx_energy][idx_cell][0];
// else
// psiNew[idx_cell] += _dE * _Q[idx_energy][idx_cell];
//}
}
_psi = psiNew;
for( unsigned i = 0; i < _nCells; ++i ) {
fluxNew[i] = dot( _psi[i], _weights );
fluxNew[i] = _psi[i][0]; // zeroth moment is raditation densitiy we are interested in
}
dFlux = blaze::l2Norm( fluxNew - fluxOld );
fluxOld = fluxNew;
......@@ -405,7 +401,7 @@ void PNSolver::Save() const {
std::vector<std::string> fieldNames{ "flux" };
std::vector<double> flux( _nCells, 0.0 );
for( unsigned i = 0; i < _nCells; ++i ) {
flux[i] = dot( _psi[i], _weights );
flux[i] = _psi[i][0];
}
std::vector<std::vector<double>> scalarField( 1, flux );
std::vector<std::vector<std::vector<double>>> results{ scalarField };
......
#include "problems/linesource.h"
LineSource::LineSource( Config* settings, Mesh* mesh ) : ProblemBase( settings, mesh ) { _physics = nullptr; }
// ---- LineSource_SN ----
LineSource::~LineSource(){};
LineSource_SN::LineSource_SN( Config* settings, Mesh* mesh ) : ProblemBase( settings, mesh ) { _physics = nullptr; }
VectorVector LineSource::GetScatteringXS( const std::vector<double>& energies ) {
LineSource_SN::~LineSource_SN(){};
VectorVector LineSource_SN::GetScatteringXS( const std::vector<double>& energies ) {
return VectorVector( energies.size(), Vector( _mesh->GetNumCells(), 1.0 ) );
}
VectorVector LineSource::GetTotalXS( const std::vector<double>& energies ) {
VectorVector LineSource_SN::GetTotalXS( const std::vector<double>& energies ) {
return VectorVector( energies.size(), Vector( _mesh->GetNumCells(), 1.0 ) );
}
std::vector<VectorVector> LineSource::GetExternalSource( const std::vector<double>& energies ) {
std::vector<VectorVector> LineSource_SN::GetExternalSource( const std::vector<double>& energies ) {
return std::vector<VectorVector>( 1u, std::vector<Vector>( _mesh->GetNumCells(), Vector( 1u, 0.0 ) ) );
}
std::vector<double> LineSource::GetStoppingPower( const std::vector<double>& energies ) {
std::vector<double> LineSource_SN::GetStoppingPower( const std::vector<double>& energies ) {
// @TODO
return std::vector<double>( energies.size(), 0.0 );
}
VectorVector LineSource::SetupIC() {
VectorVector LineSource_SN::SetupIC() {
VectorVector psi( _mesh->GetNumCells(), Vector( _settings->GetNQuadPoints(), 1e-10 ) );
auto cellMids = _mesh->GetCellMidPoints();
double t = 3.2e-4; // pseudo time for gaussian smoothing
......@@ -32,3 +34,49 @@ VectorVector LineSource::SetupIC() {
}
return psi;
}
// ---- LineSource_PN ----
int LineSource_PN::GlobalIndex( int l, int k ) const {
int numIndicesPrevLevel = l * l; // number of previous indices untill level l-1
int prevIndicesThisLevel = k + l; // number of previous indices in current level
return numIndicesPrevLevel + prevIndicesThisLevel;
}
LineSource_PN::LineSource_PN( Config* settings, Mesh* mesh ) : ProblemBase( settings, mesh ) { _physics = nullptr; }
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::GetTotalXS( const std::vector<double>& energies ) {
return VectorVector( energies.size(), Vector( _mesh->GetNumCells(), 1.0 ) );
}
std::vector<VectorVector> LineSource_PN::GetExternalSource( const std::vector<double>& energies ) {
return std::vector<VectorVector>( 1u, std::vector<Vector>( _mesh->GetNumCells(), Vector( 1u, 0.0 ) ) );
}
std::vector<double> LineSource_PN::GetStoppingPower( const std::vector<double>& energies ) {
// @TODO
return std::vector<double>( energies.size(), 0.0 );
}
VectorVector LineSource_PN::SetupIC() {
// Compute number of equations in the system
int ntotalEquations = GlobalIndex( _settings->GetNQuadPoints(), _settings->GetNQuadPoints() ) + 1;
VectorVector psi( _mesh->GetNumCells(), Vector( ntotalEquations, 0 ) ); // zero could lead to problems?
VectorVector cellMids = _mesh->GetCellMidPoints();
// Initial condition is dirac impulse at (x,y) = (0,0) ==> constant in angle ==> all moments are zero.
double t = 3.2e-4; // pseudo time for gaussian smoothing (Approx to dirac impulse)
for( unsigned j = 0; j < cellMids.size(); ++j ) {
double x = cellMids[j][0];
double y = cellMids[j][1];
psi[j][0] = 1 / ( 4 * M_PI * t ) * std::exp( -( x * x + y * y ) / ( 4 * t ) ) / ( 4 * M_PI );
}
return psi;
}
......@@ -10,7 +10,9 @@ ProblemBase::~ProblemBase() {}
ProblemBase* ProblemBase::Create( Config* settings, Mesh* mesh ) {
auto name = settings->GetProblemName();
switch( name ) {
case PROBLEM_LineSource: return new LineSource( settings, mesh );
case PROBLEM_LineSource:
if( settings->GetSolverName() == PN_SOLVER ) return new LineSource_PN( settings, mesh );
return new LineSource_SN( settings, mesh ); // default
case PROBLEM_Checkerboard: return new Checkerboard( settings, mesh );
case PROBLEM_ElectronRT: return new ElectronRT( settings, mesh );
default: return new ElectronRT( settings, mesh ); // Use RadioTherapy as dummy
......
Supports Markdown
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