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

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

parents 0ad69b14 233152ab
Pipeline #117984 failed with stages
in 28 minutes and 48 seconds
......@@ -54,9 +54,12 @@ class ICRU
double _BETA2;
double _R1;
Vector _E, _QMU; // input energy and mu vectors
Vector _E, /*! @brief: User queried Energy */
_QMU; /*! @brief:User queried mu */
std::vector<double> _ET, _ETL;
std::vector<double> _XMU, _XMUL;
std::vector<double> _XMU, /* angular variable mu of dataset */
_XMUL;
std::vector<double> _ECS, _PCS;
std::vector<double> _ETCS1, _ETCS2;
std::vector<double> _PTCS1, _PTCS2;
......@@ -142,8 +145,20 @@ class ICRU
}
~ICRU(){};
/*! @brief: Computes the angular scattering cross sections and integrated XS by interpolating tabular data
* @arg: Matrix& angularXS = Matrix where the angular XS will be saved
* @arg: Vector& integratedXS = Vector, where the integratedXS will be saved
* @returns: void */
void GetAngularScatteringXS( Matrix& angularXS, Vector& integratedXS );
/*! @brief: Computes the stopping power by interpolating tabular data
* @arg: Vector& stoppingPower = Vector, where the stopping power will be saved
* @returns: void */
void GetStoppingPower( Vector& stoppingPower );
/*! @brief: Computes the transport coefficients by interpolating tabular data
* @arg: Matrix& xi = Vector, where the stopping power will be saved
* @returns: void */
void GetTransportCoefficients( Matrix& xi );
};
......
......@@ -540,19 +540,24 @@ double ICRU::DODCSC( double RMU, double EL ) {
void ICRU::GetAngularScatteringXS( Matrix& angularXS, Vector& integratedXS ) {
angularXS.resize( _QMU.size(), _E.size() );
std::vector<double> mu( _XMU.size() );
for( unsigned n = 0; n < mu.size(); ++n ) mu[n] = 1.0 - 2.0 * _XMU[n];
Vector QMU_trafo( _QMU.size() );
// for( unsigned n = 0; n < mu.size(); ++n ) mu[n] = 1.0 - 2.0 * _XMU[n]; // trafo
std::sort( mu.begin(), mu.end() );
// Ruecktrafo für _QM
for( unsigned n = 0; n < _QMU.size(); ++n ) {
QMU_trafo[n] = ( 1 - _QMU[n] ) / 2.0; // trafo
}
// std::sort( mu.begin(), mu.end() );
integratedXS.resize( _E.size() );
for( unsigned i = 0; i < _E.size(); ++i ) {
std::vector<double> dxse, dxsi, dxs;
angdcs( 3u, _E[i], dxse, dxsi, dxs );
Interpolation interpDXS( mu, dxs );
blaze::column( angularXS, i ) = interpDXS( _QMU );
for( unsigned j = 0; j < mu.size() - 1; ++j ) {
integratedXS[i] += 0.5 * ( mu[j + 1] - mu[j] ) * ( dxs[j + 1] + dxs[j] );
angdcs( 3u, _E[i], dxse, dxsi, dxs ); // speedup here possible
Interpolation interpDXS( _XMU, dxs );
blaze::column( angularXS, i ) = interpDXS( QMU_trafo );
for( unsigned j = 0; j < _XMU.size() - 1; ++j ) {
integratedXS[i] += 0.5 * ( _XMU[j + 1] - _XMU[j] ) * ( dxs[j + 1] + dxs[j] );
}
}
angularXS *= H2OMolecularDensity;
......
......@@ -13,13 +13,13 @@ CSDSNSolver::CSDSNSolver( Config* settings ) : SNSolver( settings ) {
std::cout << "Start CSDN Constructor\n";
_dose = std::vector<double>( _settings->GetNCells(), 0.0 );
// Set angle and energies
// --- Set angle and energies
_angle = Vector( _settings->GetNQuadPoints(), 0.0 ); // my
_energies = Vector( _nEnergies, 0.0 ); // equidistant
double energyMin = 2;
double energyMin = 1e-1;
double energyMax = 5e0;
// write equidistant energy grid
// --- Write equidistant energy grid
_dE = ComputeTimeStep( settings->GetCFL() );
_nEnergies = unsigned( ( energyMax - energyMin ) / _dE );
_energies.resize( _nEnergies );
......@@ -27,11 +27,16 @@ CSDSNSolver::CSDSNSolver( Config* settings ) : SNSolver( settings ) {
_energies[n] = energyMin + ( energyMax - energyMin ) / ( _nEnergies - 1 ) * n;
}
// --- Create Quadrature Points from Quad Rule
/* Legacy Code -- unused: Will be deleted
std::vector<double> buffer;
for( auto q : _quadPoints )
if( q[0] > 0.0 ) buffer.push_back( q[0] );
std::sort( buffer.begin(), buffer.end() );
Vector posMu( buffer.size(), buffer.data() );
*/
// create 1D quadrature
// write mu grid
Matrix muMatrix( _settings->GetNQuadPoints(), _settings->GetNQuadPoints() );
......@@ -50,10 +55,12 @@ CSDSNSolver::CSDSNSolver( Config* settings ) : SNSolver( settings ) {
// store Matrix with mu values in vector format to call GetScatteringXS
for( unsigned i = 0; i < muMatrix.rows(); ++i ) {
for( unsigned j = 0; j < muMatrix.columns(); ++j ) {
angleVec[i * muMatrix.columns() + j] = std::fabs( muMatrix( i, j ) ); // icru xs go from 0 to 1 due to symmetry
angleVec[i * muMatrix.columns() + j] = muMatrix( i, j );
}
}
std::cout << "Here 3\n";
ICRU database( angleVec, _energies );
Matrix total;
database.GetAngularScatteringXS( total, _sigmaTE );
......@@ -84,6 +91,7 @@ CSDSNSolver::CSDSNSolver( Config* settings ) : SNSolver( settings ) {
// _s = Vector( _nEnergies, 1.0 );
//_s = _problem->GetStoppingPower( _energies );
database.GetStoppingPower( _s );
_Q = _problem->GetExternalSource( _energies );
// recompute scattering kernel. TODO: add this to kernel function
......@@ -184,6 +192,9 @@ void CSDSNSolver::Solve() {
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];
if( _sol[0][k] < 0 ) {
ErrorMessages::Error( "boundary negative", CURRENT_FUNCTION );
}
}
}
......@@ -230,9 +241,8 @@ void CSDSNSolver::Solve() {
}
*/
}
std::cout << psiNew[0] << " | " << _sol[0] << std::endl;
for( unsigned j = 0; j < _nCells; ++j ) {
for( unsigned j = 1; j < _nCells; ++j ) {
if( _boundaryCells[j] == BOUNDARY_TYPE::DIRICHLET ) continue;
_sol[j] = psiNew[j];
}
......@@ -244,6 +254,7 @@ void CSDSNSolver::Solve() {
_density[j]; // update dose with trapezoidal rule
_solverOutput[j] = fluxNew[j];
}
// std::cout << "dose at boundary: " << _dose[0] << " \n| " << _sol[0] << std::endl;
Save( n );
dFlux = blaze::l2Norm( fluxNew - fluxOld );
......
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