Commit da3b6414 authored by jonas.kusch's avatar jonas.kusch
Browse files

Cleaning up code

parent 4e5ed7a0
Pipeline #117719 failed with stages
in 26 minutes and 43 seconds
......@@ -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
......
#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"
ProblemBase::ProblemBase( Config* settings, Mesh* mesh ) {
_settings = settings;
......
......@@ -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;
}
......
......@@ -81,29 +81,19 @@ 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);
}
void CSDSolverTrafoFP::Solve() {
......@@ -113,9 +103,7 @@ void CSDSolverTrafoFP::Solve() {
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
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 ) );
}
......@@ -209,12 +197,6 @@ void CSDSolverTrafoFP::Solve() {
}
}
for( unsigned j = 0; j < _nCells; ++j ) {
if( _boundaryCells[j] == BOUNDARY_TYPE::DIRICHLET ) continue;
psi1[j] = blaze::solve( _IL, _sol[j] );
psi1[j] = _alpha * _L * psi1[j] + _alpha2 * _L * _sol[j];
}
// add FP scattering term implicitly
for( unsigned j = 0; j < _nCells; ++j ) {
if( _boundaryCells[j] == BOUNDARY_TYPE::DIRICHLET ) continue;
......@@ -240,7 +222,7 @@ void CSDSolverTrafoFP::Solve() {
_normals[j][idx_neighbor] ) /
_areas[j];
}
// tteamime update angular flux with numerical flux and total scattering cross section
// 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];
}
}
......
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