Commit 14f7860e authored by jonas.kusch's avatar jonas.kusch
Browse files

FP with errors

parent 1ed01498
Pipeline #111607 failed with stages
in 26 minutes and 51 seconds
......@@ -89,10 +89,10 @@ 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, PN_SOLVER, MN_SOLVER };
enum SOLVER_NAME { SN_SOLVER, CSD_SN_SOLVER, CSD_SN_NOTRAFO_SOLVER, CSD_SN_FOKKERPLANCK_SOLVER, 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 }, { "PN_SOLVER", PN_SOLVER }, { "MN_SOLVER", MN_SOLVER } };
{ "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 }, { "PN_SOLVER", PN_SOLVER }, { "MN_SOLVER", MN_SOLVER } };
// Entropy functional
enum ENTROPY_NAME { QUADRATIC, MAXWELL_BOLZMANN, BOSE_EINSTEIN, FERMI_DIRAC };
......
#ifndef SLABGEOHG_H
#define SLABGEOHG_H
#include "problembase.h"
class SlabGeoHG : public ProblemBase
{
private:
SlabGeoHG() = delete;
public:
SlabGeoHG( Config* settings, Mesh* mesh );
virtual ~SlabGeoHG();
virtual VectorVector GetScatteringXS( const Vector& energies );
virtual VectorVector GetTotalXS( const Vector& energies );
virtual std::vector<VectorVector> GetExternalSource( const Vector& energies );
virtual VectorVector SetupIC();
};
#endif // SLABGEOHG_H
#ifndef CSDSNSOLVERFP_H
#define CSDSNSOLVERFP_H
#include "icru.h"
#include "solvers/snsolver.h"
class Physics;
class CSDSNSolverFP : 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 */
Vector _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 _IL; /*! @brief Laplace Beltrami Matrix */
double _alpha;
double _beta;
public:
/**
* @brief CSDSNSolverFP constructor
* @param settings stores all needed information
*/
CSDSNSolverFP( 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 // CSDSNSOLVERFP_H
......@@ -3,16 +3,16 @@ OUTPUT_FILE = example_csd_1d
LOG_DIR = ../result/logs
MESH_FILE = 1DMesh.su2
PROBLEM = WATERPHANTOM
SOLVER = CSD_SN_NOTRAFO_SOLVER
SOLVER = CSD_SN_FOKKERPLANCK_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.001
TIME_FINAL = 1.0
CLEAN_FLUX_MATRICES = NO
BC_DIRICHLET = ( dirichlet )
BC_NEUMANN = ( wall_low, wall_up )
QUAD_TYPE = GAUSS_LEGENDRE_1D
QUAD_ORDER = 101
QUAD_ORDER = 21
#include "problems/slabgeohg.h"
#include "common/config.h"
#include "common/mesh.h"
#include "physics.h"
// ---- SlabGeoHG ----
SlabGeoHG::SlabGeoHG( Config* settings, Mesh* mesh ) : ProblemBase( settings, mesh ) { _physics = nullptr; }
SlabGeoHG::~SlabGeoHG(){};
VectorVector SlabGeoHG::GetScatteringXS( const Vector& energies ) { return VectorVector( energies.size(), Vector( _mesh->GetNumCells(), 0.0 ) ); }
VectorVector SlabGeoHG::GetTotalXS( const Vector& energies ) { return VectorVector( energies.size(), Vector( _mesh->GetNumCells(), 0.0 ) ); }
std::vector<VectorVector> SlabGeoHG::GetExternalSource( const Vector& energies ) {
return std::vector<VectorVector>( 1u, std::vector<Vector>( _mesh->GetNumCells(), Vector( 1u, 0.0 ) ) );
}
VectorVector SlabGeoHG::SetupIC() {
VectorVector psi( _mesh->GetNumCells(), Vector( _settings->GetNQuadPoints(), 0.0 ) );
auto boundaryCells = _mesh->GetBoundaryTypes();
auto cellMids = _mesh->GetCellMidPoints();
double t = 3.2e-4; // pseudo time for gaussian smoothing
return psi;
}
#include "solvers/csdsnsolverfp.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"
#include "common/mesh.h"
// externals
#include "spdlog/spdlog.h"
#include <mpi.h>
CSDSNSolverFP::CSDSNSolverFP( Config* settings ) : SNSolver( settings ) {
_dose = std::vector<double>( _settings->GetNCells(), 0.0 );
// Set angle and energies
_angle = Vector( _settings->GetNQuadPoints(), 0.0 ); // my
_energies = Vector( _nEnergies, 0.0 ); // equidistant
double energyMin = 1e-1;
double energyMax = 1e0;
// write equidistant energy grid
_dE = ComputeTimeStep( settings->GetCFL() );
_nEnergies = unsigned( ( energyMax - energyMin ) / _dE );
_energies.resize( _nEnergies );
for( unsigned n = 0; n < _nEnergies; ++n ) {
_energies[n] = energyMin + ( energyMax - energyMin ) / ( _nEnergies - 1 ) * n;
}
// create 1D quadrature
unsigned nq = _settings->GetNQuadPoints();
QuadratureBase* quad1D = QuadratureBase::CreateQuadrature(QUAD_GaussLegendre1D,nq);
Vector w = quad1D->GetWeights();
VectorVector muVec = quad1D->GetPoints();
Vector mu(nq);
for(unsigned k=0;k<nq;++k){
mu[k] = muVec[k][0];
}
// setup Laplace Beltrami matrix L in slab geometry
_L = Matrix(nq,nq,0.0);
double DMinus = 0.0;
double DPlus = 0.0;
for( unsigned k = 0; k < nq; ++k ){
DMinus = DPlus;
DPlus = DMinus-2*mu[k]*w[k];
if( k > 0 ){
_L(k,k-1) = DMinus/(mu[k]-mu[k-1])/w[k];
_L(k,k) = -DMinus/(mu[k]-mu[k-1])/w[k];
}
if( k < nq-1 ){
_L(k,k+1) = DPlus/(mu[k+1]-mu[k])/w[k];
_L(k,k) += -DPlus/(mu[k+1]-mu[k])/w[k];
}
}
// Heney-Greenstein parameter
double g = 0.8;
// determine momente of Heney-Greenstein
double xi1 = 1.0-g; // paper Olbrant, Frank (11)
double xi2 = 4.0/3.0 - 2.0 * g + 2.0/3.0 * g*g;
// determine alpha and beta
_alpha = xi1/2.0 + xi2/8.0;
_beta = xi2/8.0/xi1;
std::cout<<"alpha = "<<_alpha<<std::endl;
std::cout<<"beta = "<<_beta<<std::endl;
Matrix identity(nq,nq,0.0);
for( unsigned k = 0; k < nq; ++k ) identity(k,k) = 1.0;
_IL = identity - _beta*_L;
std::cout<<_IL<<std::endl;
_s = Vector(_nEnergies,1.0);
// 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;
}
_scatteringKernel( p, p ) = _weights[p];
}
_density = Vector( _nCells, 1.0 );
//exit(EXIT_SUCCESS);
}
void CSDSNSolverFP::Solve() {
std::cout << "Solve Fokker-Planck with Heney-Greenstein kernel" << std::endl;
auto log = spdlog::get( "event" );
//int* ipiv = new int[_nq];
//for( unsigned i = 0; i < _nq ;++i ) ipiv[i] = i+1;
blaze::DynamicVector<int> ipiv = blaze::linspace( _nq, 0, static_cast<int>( _nq )-1 );
//blaze::gesv( A, b, ipiv.data() );
auto cellMids = _mesh->GetCellMidPoints();
// setup IC and incoming BC on left
//auto cellMids = _settings->GetCellMidPoints();
_sol = std::vector<Vector>( _nCells, Vector( _nq, 0.0 ) );
for( unsigned j = 0; j < _nCells; ++j ) {
//if(_boundaryCells[j] == BOUNDARY_TYPE::NEUMANN) std::cout<<"BOUNDARY CELL DETECTED. x = "<<cellMids[j][0]<<std::endl;
if(_boundaryCells[j] == BOUNDARY_TYPE::DIRICHLET && cellMids[j][0] < 1.0){
std::cout<<"BOUNDARY CELL is left!"<<std::endl;
for( unsigned k = 0; k<_nq; ++k ){
if(_quadPoints[k][0] > 0){
_sol[j][k] = 1e5*exp(-10.0*pow(1.0-_quadPoints[k][0],2));
}
}
}
}
for( unsigned k = 0; k<_nq; ++k ){
if(_quadPoints[k][0] > 0){
_sol[0][k] = 1e5*exp(-10.0*pow(1.0-_quadPoints[k][0],2));
}
}
_boundaryCells[0] = BOUNDARY_TYPE::DIRICHLET;
_boundaryCells[_nCells-1] = BOUNDARY_TYPE::DIRICHLET;
// angular flux at next time step (maybe store angular flux at all time steps, since time becomes energy?)
VectorVector psiNew( _nCells, Vector( _nq, 0.0 ) );
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" );
// revert energies and stopping power
Vector sSave = _s;
for( unsigned n = 0; n < _nEnergies; ++n ) {
_energies[n] = _energies[_nEnergies - 1] - _energies[n];
_s[n] = sSave[_nEnergies - 1 - n];
}
for( unsigned j = 0; j < _nCells; ++j ) {
for( unsigned k = 0; k < _nq; ++k ) {
fluxOld[j] += _weights[k] * _sol[j][k] * _s[0];
}
}
// determine minimal density for CFL computation
double densityMin = _density[0];
for( unsigned j = 1; j < _nCells; ++j ) {
if( densityMin > _density[j] ) densityMin = _density[j];
}
VectorVector psi1 = _sol;
// 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?
// loop over all spatial cells
// scattering step
for( unsigned j = 0; j < _nCells; ++j ) {
if( _boundaryCells[j] == BOUNDARY_TYPE::DIRICHLET ) continue;
psi1[j] = _sol[j];
Vector psiVec = _sol[j];
blaze::gesv(_IL,psi1[j],ipiv.data());
if(norm(_IL*psi1[j]-_sol[j])> 0.1){
std::cout<<_sol[j]<<std::endl;
std::cout<<"----------------------------"<<std::endl;
std::cout<<psi1[j]<<std::endl;
std::cout<<"res = "<<norm(_IL*psi1[j]-_sol[j])<<std::endl;
exit(EXIT_SUCCESS);
}
psi1[j] = _alpha*_L*psi1[j];
}
// advection step
for( unsigned j = 0; j < _nCells; ++j ) {
if( _boundaryCells[j] == BOUNDARY_TYPE::DIRICHLET ) continue;
// loop over all ordinates
for( unsigned i = 0; i < _nq; ++i ) {
psiNew[j][i] = 0.0;
// loop over all neighbor cells (edges) of cell j and compute numerical fluxes
for( unsigned idx_neighbor = 0; idx_neighbor < _neighbors[j].size(); ++idx_neighbor ) {
// store flux contribution on psiNew_sigmaS to save memory
if( _boundaryCells[j] == BOUNDARY_TYPE::NEUMANN && _neighbors[j][idx_neighbor] == _nCells )
continue; // adiabatic wall, add nothing
else
psiNew[j][i] += _g->Flux( _quadPoints[i], _sol[j][i], _sol[_neighbors[j][idx_neighbor]][i], _normals[j][idx_neighbor] ) /
_areas[j] / ( _density[j] * _s[n + 1] );
}
// time update angular flux with numerical flux and total scattering cross section
psiNew[j][i] = _sol[j][i] - _dE * psiNew[j][i] - _dE * psi1[j][i];
}
}
_sol = psiNew;
for( unsigned j = 0; j < _nCells; ++j ) {
fluxNew[j] = 0.0;
for( unsigned k = 0; k < _nq; ++k ) {
fluxNew[j] += _weights[k] * _sol[j][k] * _s[n + 1];
}
_dose[j] += 0.5 * _dE * ( fluxNew[j] + fluxOld[j] ); // update dose with trapezoidal rule
_solverOutput[j] = fluxNew[j];
}
// Save( n );
dFlux = blaze::l2Norm( fluxNew - fluxOld );
fluxOld = fluxNew;
if( rank == 0 ) log->info( "{:03.8f} {:01.5e} {:01.5e}", _energies[n], _dE / densityMin, dFlux );
if( std::isinf( dFlux ) || std::isnan( dFlux ) ) break;
}
Save( 1 );
}
void CSDSNSolverFP::Save() const {
std::vector<std::string> fieldNames{ "dose", "normalized dose" };
std::vector<std::vector<double>> dose( 1, _dose );
std::vector<std::vector<double>> normalizedDose( 1, _dose );
double maxDose = *std::max_element( _dose.begin(), _dose.end() );
for( unsigned i = 0; i < _dose.size(); ++i ) normalizedDose[0][i] /= maxDose;
std::vector<std::vector<std::vector<double>>> results{ dose, normalizedDose };
ExportVTK( _settings->GetOutputFile(), results, fieldNames, _mesh );
}
void CSDSNSolverFP::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 );
}
......@@ -11,6 +11,7 @@
#include "solvers/pnsolver.h"
#include "solvers/snsolver.h"
#include "solvers/csdsnsolvernotrafo.h"
#include "solvers/csdsnsolverfp.h"
Solver::Solver( Config* settings ) : _settings( settings ) {
// @TODO save parameters from settings class
......@@ -74,6 +75,7 @@ Solver* Solver::Create( Config* settings ) {
case SN_SOLVER: return new SNSolver( settings );
case CSD_SN_SOLVER: return new CSDSNSolver( settings );
case CSD_SN_NOTRAFO_SOLVER: return new CSDSNSolverNoTrafo( settings );
case CSD_SN_FOKKERPLANCK_SOLVER: return new CSDSNSolverFP( settings );
case PN_SOLVER: return new PNSolver( settings );
case MN_SOLVER: return new MNSolver( settings );
default: return new SNSolver( settings );
......
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