Commit 16e7ef5c authored by jonas.kusch's avatar jonas.kusch
Browse files

bugfix transport coefficients

parent e469a984
Pipeline #117323 failed with stages
in 26 minutes and 25 seconds
......@@ -8,11 +8,11 @@ CONTINUOUS_SLOWING_DOWN = YES
HYDROGEN_FILE = ENDL_H.txt
OXYGEN_FILE = ENDL_O.txt
KERNEL = ISOTROPIC_1D
CFL_NUMBER = 0.005
CFL_NUMBER = 0.01
TIME_FINAL = 1.0
CLEAN_FLUX_MATRICES = NO
BC_DIRICHLET = ( dirichlet )
BC_NEUMANN = ( wall_low, wall_up )
QUAD_TYPE = GAUSS_LEGENDRE_1D
QUAD_ORDER = 201
QUAD_ORDER = 11
......@@ -465,53 +465,24 @@ void ICRU::GetAngularScatteringXS( Matrix& angularXS, Vector& integratedXS ) {
}
void ICRU::GetTransportCoefficients( Matrix& xi ) {
//xi.resize(xi.rows(),ELAW.size());
Vector mu(_XMU.size());
for( unsigned n = 0; n < mu.size(); ++n ) mu[n] = fabs(_XMU[n]-1.0); // mu starts at 1 and goes to 0
for( unsigned n = 0; n < mu.size(); ++n ) mu[n] = 1.0-2.0*_XMU[n]; // mu starts at 1 and goes to 0
for( unsigned i = 0; i < _E.size(); ++i ) {
std::vector<double> dxse, dxsi, dxs;
//_E[i] = 10*1e6;
// get scattering cross sections for energy E[i]
angdcs( 3u, _E[i], dxse, dxsi, dxs );
/*
std::cout<<"mu = [ ";
for( unsigned k = 0; k < dxse.size(); ++k ){
std::cout<<_XMU[k]<<"; ";
}
std::cout<<"];"<<std::endl;
std::cout<<"sigmaE = [ ";
for( unsigned k = 0; k < dxse.size(); ++k ){
std::cout<<dxse[k]*H2OMolecularDensity*2*PI<<"; "; // TODO: fix to barn per ster
}
std::cout<<"];"<<std::endl;*/
//exit(EXIT_FAILURE);
//std::cout<<"size dxs -> "<<dxs.size()<<std::endl;
// compute moments with trapezoidal rule
for( unsigned n = 0; n < xi.rows(); ++n ){
xi(n,i) = 0.0;
// so far only integration on [1,0] -> angles from 180 to 0 degrees
// integration over mu
for( unsigned k = 0; k < dxs.size()-1; ++k ){
//std::cout<<n<<" on [0,1] : "<<pow(1.0-_XMU[k+1],n)<<" "<<( pow(1.0-_XMU[k+1],n)*dxs[k+1] + pow(1.0+_XMU[k],n)*dxs[k] )*(_XMU[k+1]-_XMU[k])<<std::endl;
//std::cout<<mu[k]<<" "<<dxs[k]<<std::endl;
xi(n,i) -= 0.5 * ( pow(1.0-mu[k+1],n)*dxs[k+1] + pow(1.0-mu[k],n)*dxs[k] )*(mu[k+1]-mu[k]);
//std::cout<<"XMU at "<<k<<":"<<_XMU[k]<<std::endl;
//std::cout<<"dxs at "<<k<<": "<<dxs[k]<<std::endl;
}
// integration over [-1,0] -> angles from -180 to 0 degrees
for( unsigned k = 0; k < dxs.size()-1; ++k ){
//std::cout<<n<<" on [-1,0] :"<<pow(1.0+_XMU[k+1],n)<<" "<<( pow(1.0+_XMU[k+1],n)*dxs[k+1] + pow(1.0+_XMU[k],n)*dxs[k] )*(_XMU[k+1]-_XMU[k])<<std::endl;
//xi(n,i) += 0.5 * ( pow(1.0+mu[k+1],n)*dxs[k+1] + pow(1.0+mu[k],n)*dxs[k] )*(-mu[k+1]+mu[k]);
//std::cout<<"XMU at "<<k<<":"<<_XMU[k]<<std::endl;
//std::cout<<"dxs at "<<k<<": "<<dxs[k]<<std::endl;
}
}
}
xi *= 2.0*PI*H2OMolecularDensity;
// account for 2 times integration domain over mu\in [0,1]
xi *= 2.0;
}
void ICRU::GetStoppingPower( Vector& stoppingPower ) {
......
......@@ -18,11 +18,13 @@ CSDSNSolverFP::CSDSNSolverFP( Config* settings ) : SNSolver( settings ) {
_angle = Vector( _settings->GetNQuadPoints(), 0.0 ); // my
_energies = Vector( _nEnergies, 0.0 ); // equidistant
_energyMin = 1e-4;
_energyMax = 10.0;
//_energyMax = 10.0;
_energyMax = 5.0;
// write equidistant energy grid
_dE = ComputeTimeStep( settings->GetCFL() );
_nEnergies = unsigned( ( _energyMax - _energyMin ) / _dE );
std::cout<<"nEnergies = "<<_nEnergies<<std::endl;
_energies.resize( _nEnergies );
for( unsigned n = 0; n < _nEnergies; ++n ) {
_energies[n] = _energyMin + ( _energyMax - _energyMin ) / ( _nEnergies - 1 ) * n;
......@@ -70,9 +72,31 @@ CSDSNSolverFP::CSDSNSolverFP( Config* settings ) : SNSolver( settings ) {
// read in medical data if radiation therapy option selected
if(_RT){
//_nEnergies = 1000;
//_energies.resize(_nEnergies);
//_xi = Matrix(6,_nEnergies);
//double minExp = -4.0;
//double maxExp = 1.5;
//for( int n = 0; n<_nEnergies; ++n){
// double exponent = minExp + ( maxExp - minExp ) / ( _nEnergies - 1 ) * n;
// _energies[n] = pow(10.0,exponent);
//}
ICRU database( abs(mu), _energies );
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
......@@ -88,7 +112,8 @@ CSDSNSolverFP::CSDSNSolverFP( Config* settings ) : SNSolver( settings ) {
}
void CSDSNSolverFP::Solve() {
std::cout << "Solve Fokker-Planck with Heney-Greenstein kernel" << std::endl;
std::cout << "Solve Fokker-Planck with Heney-Greenstein kernel using "<<_nEnergies<<" energies " << std::endl;
auto log = spdlog::get( "event" );
auto cellMids = _mesh->GetCellMidPoints();
......
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