Commit 1a11aef8 authored by Jonas Kusch's avatar Jonas Kusch
Browse files

rescale fix in scattering

parent e024fc99
Pipeline #107480 failed with stages
in 27 minutes and 59 seconds
......@@ -4,11 +4,11 @@ LOG_DIR = ../result/logs
MESH_FILE = linesource_pseudo_1D.su2
PROBLEM = WATERPHANTOM
SOLVER = CSD_SN_SOLVER
CONTINUOUS_SLOWING_DOWN = NO
CONTINUOUS_SLOWING_DOWN = YES
HYDROGEN_FILE = ENDL_H.txt
OXYGEN_FILE = ENDL_O.txt
KERNEL = ISOTROPIC_1D
CFL_NUMBER = 1.5
CFL_NUMBER = 0.1
TIME_FINAL = 1.0
CLEAN_FLUX_MATRICES = NO
......
......@@ -21,7 +21,7 @@ Matrix Isotropic1D::GetScatteringKernel() {
tmp += kernel( i, j );
}
for( unsigned j = 0; j < nq; ++j ) {
kernel( i, j ) /= tmp;
// kernel( i, j ) /= tmp;
}
}
return kernel;
......
......@@ -118,8 +118,9 @@ VectorVector Physics::GetScatteringXS( Vector energies, Vector angle ) {
for( unsigned a = 1; a < angle.size(); ++a ) {
integral_sXS += 0.5 * ( angle[a] - angle[a - 1] ) * ( sXS[a] + sXS[a - 1] );
}
intermediateGrid.back() /= integral_sXS; // re-normalize to yield a valid distribution in a statistical sense; behaves poorly due to great
// differences in order of magnitude
// intermediateGrid.back() /= integral_sXS; // re-normalize to yield a valid distribution in a statistical sense; behaves poorly due to
// great
// differences in order of magnitude
}
VectorVector xs( energies.size(), Vector( angle.size() ) );
for( unsigned j = 0; j < angle.size(); ++j ) {
......@@ -130,8 +131,8 @@ VectorVector Physics::GetScatteringXS( Vector energies, Vector angle ) {
Interpolation interp( dataEnergy, xsPerE, Interpolation::linear );
Vector eGrid = interp( energies );
for( unsigned i = 0; i < energies.size(); ++i ) {
xs[i][j] = eGrid[i] * integratedScatteringXS[i]; // multiply with the integrated scattering cross section to force the
// integration over the angle to yield integratedScatteringXS
// xs[i][j] = eGrid[i] * integratedScatteringXS[i]; // multiply with the integrated scattering cross section to force the
// integration over the angle to yield integratedScatteringXS
}
}
......
......@@ -15,7 +15,7 @@ CSDSNSolver::CSDSNSolver( Config* settings ) : SNSolver( settings ) {
// Set angle and energies
_angle = Vector( _settings->GetNQuadPoints(), 0.0 ); // my
_energies = Vector( _nEnergies, 0.0 ); // equidistant
double energyMin = 1e-1;
double energyMin = 1e-5;
double energyMax = 5e0;
// write equidistant energy grid
......@@ -39,14 +39,59 @@ CSDSNSolver::CSDSNSolver( Config* settings ) : SNSolver( settings ) {
_sigmaSE = _problem->GetScatteringXSE( _energies, muMatrix );
_sigmaTE = _problem->GetTotalXSE( _energies );
_s = _problem->GetStoppingPower( _energies );
_Q = _problem->GetExternalSource( _energies );
// compute scaling s.t. scattering kernel integrates to one for chosen quadrature
for( unsigned n = 0; n < _nEnergies; ++n ) {
for( unsigned p = 0; p < _nq; ++p ) {
double cp = 0.0;
for( unsigned q = 0; q < _nq; ++q ) {
cp += _weights[q] * _sigmaSE[n]( p, q );
}
for( unsigned q = 0; q < _nq; ++q ) {
_sigmaSE[n]( p, q ) = _sigmaTE[n] / cp * _sigmaSE[n]( p, q );
}
}
}
_s = _problem->GetStoppingPower( _energies );
_Q = _problem->GetExternalSource( _energies );
// 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, p ) = _weights[p];
}
}
/*
// test 1
for( unsigned n = 0; n < _nEnergies; ++n ) {
for( unsigned p = 0; p < _nq; ++p ) {
double tmp = 0.0;
for( unsigned q = 0; q < _nq; ++q ) {
tmp += _sigmaSE[n]( p, q ) * _scatteringKernel( q, p );
}
std::cout << tmp << " ";
}
std::cout << std::endl;
}
// exit( EXIT_FAILURE );
// test scattering
Vector ones( _nq, 1.0 );
for( unsigned n = 0; n < _nEnergies; ++n ) {
std::cout << ( _sigmaSE[_nEnergies - n - 1] * _scatteringKernel * ones )[0] << " " << _sigmaTE[_nEnergies - n - 1] << std::endl;
}*/
// exit( EXIT_FAILURE );
// Get patient density
_density = Vector( _nCells, 1.0 );
}
void CSDSNSolver::Solve() {
std::cout << "Solve" << std::endl;
auto log = spdlog::get( "event" );
// angular flux at next time step (maybe store angular flux at all time steps, since time becomes energy?)
......@@ -69,9 +114,10 @@ void CSDSNSolver::Solve() {
}
// 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];
double tmp = 0.0;
_energies[0] = 0.0;
for( unsigned n = 1; n < _nEnergies; ++n ) {
tmp = tmp + _dE * 0.5 * ( 1.0 / _s[n] + 1.0 / _s[n - 1] );
_energies[n] = tmp;
}
......@@ -145,6 +191,7 @@ void CSDSNSolver::Solve() {
// Save( n );
dFlux = blaze::l2Norm( fluxNew - fluxOld );
fluxOld = fluxNew;
// std::cout << _energies[n] << std::endl;
if( rank == 0 ) log->info( "{:03.8f} {:01.5e} {:01.5e}", _energies[n], _dE / densityMin, dFlux );
if( std::isinf( dFlux ) || std::isnan( dFlux ) ) break;
}
......
......@@ -49,6 +49,7 @@ TEST_CASE( "checking angular integral of scattering cross sections to be unity",
for( unsigned a = 1; a < angles.size(); ++a ) {
integral_sXS[e] += 0.5 * ( angles[a] - angles[a - 1] ) * ( sXS[e][a] + sXS[e][a - 1] );
}
std::cout << integral_sXS[e] << " " << tXS[e] << std::endl;
REQUIRE( std::fabs( integral_sXS[e] - tXS[e] ) / std::fabs( tXS[e] ) < 0.05 );
}
}
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