phantom 2D running

code/input/phantom2d_csd.cfg

@@ -3,16 +3,15 @@ OUTPUT_FILE = phantom2d_csd_10MeV_fine
 LOG_DIR = ../result/logs
 MESH_FILE = phantom2D.su2
 PROBLEM = PHANTOM2D
-SOLVER = CSD_SN_FOKKERPLANCK_SOLVER
+SOLVER = CSD_SN_FOKKERPLANCK_TRAFO_SOLVER_2D
 CONTINUOUS_SLOWING_DOWN = YES
 HYDROGEN_FILE = ENDL_H.txt
 OXYGEN_FILE = ENDL_O.txt
 KERNEL = ISOTROPIC
-CFL_NUMBER = 0.01
+CFL_NUMBER = 0.15
 TIME_FINAL = 1.0
 CLEAN_FLUX_MATRICES = NO
 
-BC_DIRICHLET = ( dirichlet )
-BC_NEUMANN = ( wall_low, wall_up )
-QUAD_TYPE = GAUSS_LEGENDRE_TENSORIZED
-QUAD_ORDER = 10
\ No newline at end of file
+BC_DIRICHLET = ( void )
+QUAD_TYPE = PRODUCT
+QUAD_ORDER = 6

code/python/mesh_from_image.py

@@ -21,7 +21,7 @@ def generate(image_name, mesh_name):
     width, height = gsImage.shape
     xRes = dimensions[0]/width
     yRes = dimensions[1]/height
-    char_length = min(xRes,yRes)*5 #arbitrary
+    char_length = min(xRes,yRes)*10 #arbitrary
     geom = pg.opencascade.Geometry()
     domain = add_rectangle(0.0, 0.0, dimensions[0], dimensions[1], char_length, geom)
     geom.add_physical(domain.lines, label="void")

code/src/problems/phantom2d.cpp

@@ -12,11 +12,14 @@ std::vector<VectorVector> Phantom2D::GetExternalSource( const Vector& energies )
 
 VectorVector Phantom2D::SetupIC() {
     VectorVector psi( _mesh->GetNumCells(), Vector( _settings->GetNQuadPoints(), 1e-10 ) );
-    auto cellMids = _mesh->GetCellMidPoints();
-    double s      = 0.1;
+    auto cellMids         = _mesh->GetCellMidPoints();
+    double s              = 0.1;
+    double enterPositionX = 0.0;
+    double enterPositionY = 0.5;
     for( unsigned j = 0; j < cellMids.size(); ++j ) {
-        double x                                = cellMids[j][0];
-        psi[j][_settings->GetNQuadPoints() - 1] = 1.0 / ( s * sqrt( 2 * M_PI ) ) * std::exp( -x * x / ( 2 * s * s ) );
+        double x = cellMids[j][0] - enterPositionX;
+        double y = cellMids[j][1] - enterPositionY;
+        psi[j]   = 1.0 / ( s * sqrt( 2 * M_PI ) ) * std::exp( -( x * x + y * y ) / ( 2 * s * s ) );
     }
     return psi;
 }
@@ -35,16 +38,16 @@ std::vector<double> Phantom2D::GetDensity( const VectorVector& cellMidPoints ) {
     unsigned m = gsImage.rows();
     unsigned n = gsImage.columns();
 
-    Vector x( m + 1 ), y( n + 1 );
-    for( unsigned i = 0; i < m + 1; ++i ) {
-        x[i] = static_cast<double>( i ) / static_cast<double>( m ) * ( xMax - xMin );
+    Vector x( m ), y( n );
+    for( unsigned i = 0; i < m; ++i ) {
+        x[i] = xMin + static_cast<double>( i ) / static_cast<double>( m - 1 ) * ( xMax - xMin );
     }
-    for( unsigned i = 0; i < n + 1; ++i ) y[i] = static_cast<double>( i ) / static_cast<double>( n ) * ( yMax - yMin );
+    for( unsigned i = 0; i < n; ++i ) y[i] = yMin + static_cast<double>( i ) / static_cast<double>( n - 1 ) * ( yMax - yMin );
 
     Interpolation interp( x, y, gsImage );
     std::vector<double> result( _mesh->GetNumCells(), 0.0 );
     for( unsigned i = 0; i < _mesh->GetNumCells(); ++i ) {
-        result[i] = std::clamp( interp( cellMidPoints[i][0], cellMidPoints[i][1] ), 0.0, 1.0 );
+        result[i] = std::clamp( interp( cellMidPoints[i][0], cellMidPoints[i][1] ), 0.0, 1.85 );
     }
     return result;
 }

code/src/solvers/csdsolvertrafofp2d.cpp

@@ -16,7 +16,7 @@ CSDSolverTrafoFP2D::CSDSolverTrafoFP2D( Config* settings ) : SNSolver( settings
     // Set angle and energies
     _energies  = Vector( _nEnergies, 0.0 );    // equidistant
     _energyMin = 1e-4 * 0.511;
-    _energyMax = 10e0;
+    _energyMax = 50e0;    // 50e0;
 
     // write equidistant energy grid (false) or refined grid (true)
     GenerateEnergyGrid( false );
@@ -145,25 +145,28 @@ CSDSolverTrafoFP2D::CSDSolverTrafoFP2D( Config* settings ) : SNSolver( settings
         */
     }
 
-    _density = std::vector<double>( _nCells, 1.0 );
+    //_density = std::vector<double>( _nCells, 1.0 );
     // exit(EXIT_SUCCESS);
+    double densityMin = 0.1;
+    for( unsigned j = 0; j < _nCells; ++j ) {
+        if( _density[j] < densityMin ) _density[j] = densityMin;
+    }
 }
 
 void CSDSolverTrafoFP2D::Solve() {
-    std::cout << "Solve" << std::endl;
     auto log = spdlog::get( "event" );
 
     // save original energy field for boundary conditions
     auto energiesOrig = _energies;
 
     // setup incoming BC on left
-    _sol = VectorVector( _density.size(), Vector( _settings->GetNQuadPoints(), 0.0 ) );    // hard coded IC, needs to be changed
-    for( unsigned k = 0; k < _nq; ++k ) {
-        if( _quadPoints[k][0] > 0 && !_RT ) _sol[0][k] = 1e5 * exp( -10.0 * pow( 1.0 - _quadPoints[k][0], 2 ) );
-    }
+    //_sol = VectorVector( _density.size(), Vector( _settings->GetNQuadPoints(), 0.0 ) );    // hard coded IC, needs to be changed
+    // for( unsigned k = 0; k < _nq; ++k ) {
+    //    if( _quadPoints[k][0] > 0 && !_RT ) _sol[0][k] = 1e5 * exp( -10.0 * pow( 1.0 - _quadPoints[k][0], 2 ) );
+    //}
     // hard coded boundary type for 1D testcases (otherwise cells will be NEUMANN)
-    _boundaryCells[0]           = BOUNDARY_TYPE::DIRICHLET;
-    _boundaryCells[_nCells - 1] = BOUNDARY_TYPE::DIRICHLET;
+    //_boundaryCells[0]           = BOUNDARY_TYPE::DIRICHLET;
+    //_boundaryCells[_nCells - 1] = BOUNDARY_TYPE::DIRICHLET;
 
     // setup identity matrix for FP scattering
     Matrix identity( _nq, _nq, 0.0 );
@@ -238,7 +241,7 @@ void CSDSolverTrafoFP2D::Solve() {
         _IL = identity - _beta * _L;
 
         // write BC for water phantom
-        if( _RT ) {
+        if( _RT && false ) {
             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 ) ) *
@@ -301,6 +304,7 @@ void CSDSolverTrafoFP2D::Solve() {
             log->info( "{:03.8f}  {:03.8f}  {:01.5e}  {:01.5e}", energiesOrig[_nEnergies - n - 1], _energies[n], _dE / densityMin, dFlux );
         if( std::isinf( dFlux ) || std::isnan( dFlux ) ) break;
     }
+    for( unsigned j = 0; j < _nCells; ++j ) _solverOutput[j] = _density[j];
     Save( 1 );
 }