Commit 4e554169 authored by Steffen Schotthöfer's avatar Steffen Schotthöfer
Browse files

config addon for csdpnsolver

parent 2bfe8941
Pipeline #142954 failed with stage
OUTPUT_DIR = ../result OUTPUT_DIR = ../result
OUTPUT_FILE = IsotropicPointSource_sph_nq6 OUTPUT_FILE = IsotropicPointSource_sph_nq6
LOG_DIR = ../result/logs LOG_DIR = ../result/logs
......
...@@ -563,6 +563,7 @@ void Config::SetPostprocessing() { ...@@ -563,6 +563,7 @@ void Config::SetPostprocessing() {
case CSD_SN_FOKKERPLANCK_TRAFO_SOLVER_2D: // Fallthrough case CSD_SN_FOKKERPLANCK_TRAFO_SOLVER_2D: // Fallthrough
case CSD_SN_FOKKERPLANCK_TRAFO_SH_SOLVER_2D: // Fallthrough case CSD_SN_FOKKERPLANCK_TRAFO_SH_SOLVER_2D: // Fallthrough
case CSD_SN_SOLVER: // Fallthrough case CSD_SN_SOLVER: // Fallthrough
case CSD_PN_SOLVER: // Fallthrough
supportedGroups = { MINIMAL, MEDICAL }; supportedGroups = { MINIMAL, MEDICAL };
if( supportedGroups.end() == std::find( supportedGroups.begin(), supportedGroups.end(), _volumeOutput[idx_volOutput] ) ) { if( supportedGroups.end() == std::find( supportedGroups.begin(), supportedGroups.end(), _volumeOutput[idx_volOutput] ) ) {
...@@ -694,7 +695,7 @@ bool Config::TokenizeString( string& str, string& option_name, vector<string>& o ...@@ -694,7 +695,7 @@ bool Config::TokenizeString( string& str, string& option_name, vector<string>& o
pos = str.find( "=" ); pos = str.find( "=" );
if( pos == string::npos ) { if( pos == string::npos ) {
string errmsg = "Error in Config::TokenizeString(): line in the configuration file with no \"=\" sign. "; string errmsg = "Error in Config::TokenizeString(): line in the configuration file with no \"=\" sign. ";
errmsg += "\nLook for: \n str.length() = " + std::to_string(str.length()); errmsg += "\nLook for: \n str.length() = " + std::to_string( str.length() );
spdlog::error( errmsg ); spdlog::error( errmsg );
throw( -1 ); throw( -1 );
} }
......
...@@ -12,28 +12,49 @@ std::vector<VectorVector> IsotropicSource2D::GetExternalSource( const Vector& en ...@@ -12,28 +12,49 @@ std::vector<VectorVector> IsotropicSource2D::GetExternalSource( const Vector& en
VectorVector IsotropicSource2D::SetupIC() { VectorVector IsotropicSource2D::SetupIC() {
VectorVector psi( _mesh->GetNumCells(), Vector( _settings->GetNQuadPoints(), 1e-10 ) ); VectorVector psi( _mesh->GetNumCells(), Vector( _settings->GetNQuadPoints(), 1e-10 ) );
auto cellMids = _mesh->GetCellMidPoints(); auto cellMids = _mesh->GetCellMidPoints();
double enterPositionX = 0.0; // double enterPositionX = 0.5; // 0.0;
double enterPositionY = 0.5; // double enterPositionY = 0.5;
auto boundaryCells = _mesh->GetBoundaryTypes(); // auto boundaryCells = _mesh->GetBoundaryTypes();
// Case 1: Ingoing radiation in just one cell
// find cell that best matches enter position // find cell that best matches enter position
double dist = 1000.0; // double dist = 1000.0;
unsigned indexSource = 0; // unsigned indexSource = 0;
for( unsigned j = 0; j < cellMids.size(); ++j ) { for( unsigned j = 0; j < cellMids.size(); ++j ) {
if( boundaryCells[j] == BOUNDARY_TYPE::DIRICHLET ) { // if( boundaryCells[j] == BOUNDARY_TYPE::DIRICHLET ) {
double x = cellMids[j][0];
double y = cellMids[j][1];
if( x >= 0.49 && x <= 0.5 && y >= 0.49 && y <= 0.5 ) {
psi[j] = Vector( _settings->GetNQuadPoints(), 1.0 );
}
//}
}
// psi[indexSource] = Vector( _settings->GetNQuadPoints(), 1.0 );
/*
// Case 2: Ingoing radiation as Gauss curve
double t = 1e-5; // pseudo time for gaussian smoothing
for( unsigned j = 0; j < cellMids.size(); ++j ) {
double x = cellMids[j][0] - enterPositionX; double x = cellMids[j][0] - enterPositionX;
double y = cellMids[j][1] - enterPositionY; double y = cellMids[j][1] - enterPositionY;
if( sqrt( x * x + y * y ) < dist ) { psi[j] = 1.0 / ( 4.0 * M_PI * t ) * std::exp( -( x * x + y * y ) / ( 4 * t ) );
dist = sqrt( x * x + y * y );
indexSource = j;
}
} }
} */
psi[indexSource] = Vector( _settings->GetNQuadPoints(), 1.0 );
return psi; return psi;
} }
std::vector<double> IsotropicSource2D::GetDensity( const VectorVector& /*cellMidPoints*/ ) { std::vector<double> IsotropicSource2D::GetDensity( const VectorVector& /*cellMidPoints*/ ) {
return std::vector<double>( _settings->GetNCells(), 1.0 ); double rhoL = 1.0;
double rhoR = 5.0;
std::vector<double> rho( _settings->GetNCells(), rhoL );
// use values rhoR on right third of domain
auto cellMids = _mesh->GetCellMidPoints();
for( unsigned j = 0; j < cellMids.size(); ++j ) {
double x = cellMids[j][0];
if( x >= 0.56 ) {
rho[j] = rhoR;
}
}
return rho;
} }
...@@ -49,8 +49,12 @@ CSDPNSolver::CSDPNSolver( Config* settings ) : PNSolver( settings ) { ...@@ -49,8 +49,12 @@ CSDPNSolver::CSDPNSolver( Config* settings ) : PNSolver( settings ) {
double y = _cellMidPoints[i][1]; double y = _cellMidPoints[i][1];
const double stddev = .005; const double stddev = .005;
f[i] = normpdf( x, pos_beam[0], stddev ) * normpdf( y, pos_beam[1], stddev ); f[i] = normpdf( x, pos_beam[0], stddev ) * normpdf( y, pos_beam[1], stddev );
for( unsigned j = 0; j < _nSystem; j++ ) {
// Vector IC = f * SMMoments; // must be VectorVector
}
} }
Vector IC = f * StarMAPmoments; Vector SMMoments = StarMAPmoments;
Vector IC = f * SMMoments; // must be VectorVector
Matrix sigma_t( _energies.size(), sigma_t.rows() ); Matrix sigma_t( _energies.size(), sigma_t.rows() );
for( unsigned i = 0; i < _nSystem; ++i ) { for( unsigned i = 0; i < _nSystem; ++i ) {
......
Markdown is supported
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