Commit 8a93b37b authored by Jonas Kusch's avatar Jonas Kusch
Browse files

csd solver added

parent 76aea3fa
Pipeline #93474 passed with stages
in 28 minutes and 57 seconds
......@@ -24,6 +24,7 @@ class ElectronRT : public ProblemBase
virtual std::vector<VectorVector> GetExternalSource( const std::vector<double>& energies );
virtual std::vector<double> GetStoppingPower( const std::vector<double>& energies );
virtual VectorVector SetupIC();
std::vector<double> GetDensity( const VectorVector& cellMidPoints );
};
#endif // ELECTRONRT_H
......@@ -50,6 +50,12 @@ class ProblemBase
*/
virtual std::vector<double> GetStoppingPower( const std::vector<double>& energies ) = 0;
/**
* @brief GetDensity gives back vector of densities for every spatial cell
* @param cellMidPoints is vector with cell mid points
*/
virtual std::vector<double> GetDensity( const VectorVector& cellMidPoints );
/**
* @brief Setup the initial condition for the flux psi
*/
......
......@@ -60,6 +60,9 @@ class Config
* to improve floating point accuracy */
bool _cleanFluxMat;
/*!< @brief If true, continuous slowing down approximation will be used */
bool _csd;
// Boundary Conditions
/*!< @brief List of all Pairs (marker, BOUNDARY_TYPE), e.g. (farfield,DIRICHLET).
Each Boundary Conditions must have an entry in enum BOUNDARY_TYPE*/
......@@ -207,6 +210,7 @@ class Config
PROBLEM_NAME inline GetProblemName() const { return _problemName; }
SOLVER_NAME inline GetSolverName() const { return _solverName; }
bool inline GetCleanFluxMat() const { return _cleanFluxMat; }
bool inline IsCSD() const { return _csd; }
// Boundary Conditions
BOUNDARY_TYPE GetBoundaryType( std::string nameMarker ) const; /*! @brief Get Boundary Type of given marker */
......
#ifndef CSDSNSOLVER_H
#define CSDSNSOLVER_H
#include <mpi.h>
#include "solvers/solverbase.h"
class CSDSNSolver : public Solver
{
private:
public:
/**
* @brief CSDSNSolver constructor
* @param settings stores all needed information
*/
CSDSNSolver( 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 // CSDSNSOLVER_H
......@@ -13,12 +13,14 @@ LOG_DIR = ../result/logs
% Mesh File
MESH_FILE = linesource.su2
%
PROBLEM = LINESOURCE
PROBLEM = ELECTRONRT
%
% ---- Solver specifications ----
%
% Solver type
SOLVER = PN_SOLVER
SOLVER = SN_SOLVER
% Slowing down approximation
CONTINUOUS_SLOWING_DOWN = YES
% CFL number
CFL_NUMBER = 0.8
% Final time for simulation
......
......@@ -24,7 +24,7 @@ std::vector<VectorVector> ElectronRT::GetExternalSource( const std::vector<doubl
std::vector<double> ElectronRT::GetStoppingPower( const std::vector<double>& energies ) {
// @TODO
return std::vector<double>( energies.size(), 0.0 );
return std::vector<double>( energies.size(), 1.0 );
}
VectorVector ElectronRT::SetupIC() {
......@@ -35,3 +35,5 @@ VectorVector ElectronRT::SetupIC() {
void ElectronRT::LoadXSH20( std::string fileSigmaS, std::string fileSigmaT ) {
// @TODO
}
std::vector<double> ElectronRT::GetDensity( const VectorVector& cellMidPoints ) { return std::vector<double>( cellMidPoints.size(), 1.0 ); }
......@@ -21,3 +21,5 @@ ProblemBase* ProblemBase::Create( Config* settings, Mesh* mesh ) {
default: return new ElectronRT( settings, mesh ); // Use RadioTherapy as dummy
}
}
std::vector<double> ProblemBase::GetDensity( const VectorVector& cellMidPoints ) { return std::vector<double>( cellMidPoints.size(), 1.0 ); }
......@@ -206,6 +206,9 @@ void Config::SetConfigOptions() {
/*! @brief CleanFluxMatrices \n DESCRIPTION: If true, very low entries (10^-10 or smaller) of the flux matrices will be set to zero,
* to improve floating point accuracy \n DEFAULT false \ingroup Config */
AddBoolOption( "CLEAN_FLUX_MATRICES", _cleanFluxMat, false );
/*! @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 );
// Mesh related options
// Boundary Markers
......
#include "solvers/csdsnsolver.h"
CSDSNSolver::CSDSNSolver( Config* settings ) : Solver( settings ) {}
void CSDSNSolver::Solve() {
auto log = spdlog::get( "event" );
// angular flux at next time step (maybe store angular flux at all time steps, since time becomes energy?)
VectorVector psiNew = _psi;
double dFlux = 1e10;
Vector fluxNew( _nCells, 0.0 );
Vector fluxOld( _nCells, 0.0 );
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 ) {
_psi[j][k] = _psi[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;
for( unsigned n = 0; n < _nEnergies; ++n ) {
tmp = tmp + _dE / _s[n];
_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];
}
// loop over energies (pseudo-time)
for( unsigned n = 1; n < _nEnergies; ++n ) {
_dE = _energies[n] - _energies[n - 1];
// loop over all spatial cells
for( unsigned j = 0; j < _nCells; ++j ) {
if( _boundaryCells[j] == BOUNDARY_TYPE::DIRICHLET ) continue;
// loop over all ordinates
for( unsigned k = 0; k < _nq; ++k ) {
psiNew[j][k] = 0.0;
// loop over all neighbor cells (edges) of cell j and compute numerical fluxes
for( unsigned l = 0; l < _neighbors[j].size(); ++l ) {
// store flux contribution on psiNew_sigmaS to save memory
if( _boundaryCells[j] == BOUNDARY_TYPE::NEUMANN && _neighbors[j][l] == _nCells )
psiNew[j][k] += _g->Flux( _quadPoints[k] / _density[j], _psi[j][k], _psi[j][k], _normals[j][l] );
else
psiNew[j][k] += _g->Flux( _quadPoints[k] / _density[j], _psi[j][k], _psi[_neighbors[j][l]][k], _normals[j][l] );
}
// time update angular flux with numerical flux and total scattering cross section
psiNew[j][k] = _psi[j][k] - ( _dE / _areas[j] ) * psiNew[j][k] - _dE * _sigmaT[n][j] * _psi[j][k];
}
// compute scattering effects
psiNew[j] += _dE * _sigmaS[n][j] * _scatteringKernel * _psi[j]; // multiply scattering matrix with psi
// TODO: figure out a more elegant way
// add external source contribution
if( _Q.size() == 1u ) { // constant source for all energies
if( _Q[0][j].size() == 1u ) // isotropic source
psiNew[j] += _dE * _Q[0][j][0] * _s[n];
else
psiNew[j] += _dE * _Q[0][j] * _s[n];
}
else {
if( _Q[0][j].size() == 1u ) // isotropic source
psiNew[j] += _dE * _Q[n][j][0] * _s[n];
else
psiNew[j] += _dE * _Q[n][j] * _s[n];
}
}
_psi = psiNew;
// do backsubstitution from psiTildeHat to psi (cf. Dissertation Kerstion Kuepper, Eq. 1.23)
for( unsigned j = 0; j < _nCells; ++j ) {
for( unsigned k = 0; k < _nq; ++k ) {
psiNew[j][k] = _psi[j][k] * _density[j] * _s[0]; // note that _s[0] is stopping power at highest energy
}
}
for( unsigned i = 0; i < _nCells; ++i ) {
fluxNew[i] = dot( psiNew[i], _weights );
_solverOutput[i] = fluxNew[i];
}
Save( n );
dFlux = blaze::l2Norm( fluxNew - fluxOld );
fluxOld = fluxNew;
if( rank == 0 ) log->info( "{:03.8f} {:01.5e}", _energies[n], dFlux );
}
}
void CSDSNSolver::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 );
}
std::vector<std::vector<double>> scalarField( 1, flux );
std::vector<std::vector<std::vector<double>>> results{ scalarField };
ExportVTK( _settings->GetOutputFile(), results, fieldNames, _mesh );
}
void CSDSNSolver::Save( int currEnergy ) const {
std::vector<std::string> fieldNames{ "flux" };
std::vector<std::vector<double>> scalarField( 1, _solverOutput );
std::vector<std::vector<std::vector<double>>> results{ scalarField };
ExportVTK( _settings->GetOutputFile() + "_" + std::to_string( currEnergy ), results, fieldNames, _mesh );
}
......@@ -3,6 +3,7 @@
#include "mesh.h"
#include "quadratures/quadraturebase.h"
#include "settings/globalconstants.h"
#include "solvers/csdsnsolver.h"
#include "solvers/pnsolver.h"
#include "solvers/snsolver.h"
......@@ -38,6 +39,7 @@ Solver::Solver( Config* settings ) : _settings( settings ) {
_sigmaT = _problem->GetTotalXS( _energies );
_sigmaS = _problem->GetScatteringXS( _energies );
_Q = _problem->GetExternalSource( _energies );
_density = _problem->GetDensity( _mesh->GetCellMidPoints() ); // double check: do we need to give the cell midpoints to fct?
// setup numerical flux
_g = NumericalFlux::Create( settings );
......@@ -62,21 +64,20 @@ double Solver::ComputeTimeStep( double cfl ) const {
Solver* Solver::Create( Config* settings ) {
switch( settings->GetSolverName() ) {
case SN_SOLVER: return new SNSolver( settings );
case PN_SOLVER: return new PNSolver( settings );
case SN_SOLVER: {
if( settings->IsCSD() )
return new CSDSNSolver( settings );
else
return new SNSolver( settings );
}
case PN_SOLVER: {
if( settings->IsCSD() ) {
std::cerr << "[Solver]: Continuous slowing down option not available for PN Solver!" << std::endl;
return NULL;
}
else
return new PNSolver( settings );
}
default: return new SNSolver( settings );
}
}
void Solver::Save() const {
std::vector<std::string> fieldNames{ "flux" };
std::vector<double> flux;
flux.resize( _nCells );
for( unsigned i = 0; i < _nCells; ++i ) {
flux[i] = _psi[i][0];
}
std::vector<std::vector<double>> scalarField( 1, flux );
std::vector<std::vector<std::vector<double>>> results{ scalarField };
ExportVTK( _settings->GetOutputFile() + "_" + std::to_string( _nEnergies ), results, fieldNames, _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