Skip to content
GitLab
Commit e469a984 authored by jonas.kusch's avatar jonas.kusch
Browse files

fix in ICRU read, CSF with FP works

parent 2e8704d9
Pipeline #111962 failed with stages
in 26 minutes and 55 seconds
Hide whitespace changes
Inline Side-by-side
code/input/example_csd_1d.cfg
View file @ e469a984
OUTPUT_DIR = ../result
OUTPUT_FILE = example_csd_1d
OUTPUT_FILE = example_csd_1d_10MeV_fine
LOG_DIR = ../result/logs
MESH_FILE = 1DMesh.su2
PROBLEM = WATERPHANTOM
......@@ -8,11 +8,11 @@ CONTINUOUS_SLOWING_DOWN = YES
HYDROGEN_FILE = ENDL_H.txt
OXYGEN_FILE = ENDL_O.txt
KERNEL = ISOTROPIC_1D
CFL_NUMBER = 0.001
CFL_NUMBER = 0.005
TIME_FINAL = 1.0
CLEAN_FLUX_MATRICES = NO
BC_DIRICHLET = ( dirichlet )
BC_NEUMANN = ( wall_low, wall_up )
QUAD_TYPE = GAUSS_LEGENDRE_1D
QUAD_ORDER = 21
QUAD_ORDER = 201
code/src/icru.cpp
View file @ e469a984
......@@ -467,7 +467,7 @@ 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);
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 i = 0; i < _E.size(); ++i ) {
std::vector<double> dxse, dxsi, dxs;
//_E[i] = 10*1e6;
......@@ -491,7 +491,7 @@ void ICRU::GetTransportCoefficients( Matrix& xi ) {
// compute moments with trapezoidal rule
for( unsigned n = 0; n < xi.rows(); ++n ){
xi(n,i) = 0.0;
// so far only integration on [0,1]
// so far only integration on [1,0] -> angles from 180 to 0 degrees
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;
......@@ -499,16 +499,19 @@ void ICRU::GetTransportCoefficients( Matrix& xi ) {
//std::cout<<"XMU at "<<k<<":"<<_XMU[k]<<std::endl;
//std::cout<<"dxs at "<<k<<": "<<dxs[k]<<std::endl;
}
// integration over [-1,0]
// 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]);
//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 ) {
......
code/src/solvers/csdsnsolverfp.cpp
View file @ e469a984
......@@ -17,8 +17,8 @@ CSDSNSolverFP::CSDSNSolverFP( Config* settings ) : SNSolver( settings ) {
// Set angle and energies
_angle = Vector( _settings->GetNQuadPoints(), 0.0 ); // my
_energies = Vector( _nEnergies, 0.0 ); // equidistant
_energyMin = 50;
_energyMax = 51;
_energyMin = 1e-4;
_energyMax = 10.0;
// write equidistant energy grid
_dE = ComputeTimeStep( settings->GetCFL() );
......@@ -96,24 +96,8 @@ void CSDSNSolverFP::Solve() {
// 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 ) {
if(_RT){
_sol[0][k] = 1e5 * exp( -200.0 * pow( 1.0 - _quadPoints[k][0], 2 ) ) * exp(-50*(_energyMin-_energyMax));
}else
_sol[0][k] = 1e5 * exp( -10.0 * pow( 1.0 - _quadPoints[k][0], 2 ) );
}
if( _quadPoints[k][0] > 0 && !_RT) _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;
......@@ -165,18 +149,11 @@ void CSDSNSolverFP::Solve() {
_IL = identity - _beta * _L;
std::cout << "xi0 = " << _xi(0,_nEnergies-n-1) << std::endl;
std::cout << "xi1 = " << _xi(1,_nEnergies-n-1) << std::endl;
std::cout << "xi2 = " << _xi(2,_nEnergies-n-1) << std::endl;
std::cout << "alpha = " << _alpha << std::endl;
std::cout << "beta = " << _beta << std::endl;
// write BC
if(_RT){
for( unsigned k = 0; k < _nq; ++k ) {
if( _quadPoints[k][0] > 0 ) {
_sol[0][k] = 1e5 * exp( -200.0 * pow( 1.0 - _quadPoints[k][0], 2 ) ) * exp(-50*(_energyMin-_energies[n]));
_sol[0][k] = 1e5 * exp( -200.0 * pow( 1.0 - _quadPoints[k][0], 2 ) ) * exp(-50*pow(_energyMax-_energies[n],2));
}
}
}
......@@ -187,15 +164,6 @@ void CSDSNSolverFP::Solve() {
for( unsigned j = 0; j < _nCells; ++j ) {
if( _boundaryCells[j] == BOUNDARY_TYPE::DIRICHLET ) continue;
psi1[j] = blaze::solve( _IL, _sol[j] );
/*
if( norm( _IL * psi1[j] - _sol[j] ) > 1e-3 ) {
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_FAILURE );
}*/
psi1[j] = _alpha * _L * psi1[j];
}
// advection step
......
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

Imprint