config addon for csdpnsolver

Former-commit-id: 4e554169

code/input/pointSource_sph_csd.cfg

+
 OUTPUT_DIR = ../result
 OUTPUT_FILE = IsotropicPointSource_sph_nq6
 LOG_DIR = ../result/logs

code/src/common/config.cpp

@@ -563,6 +563,7 @@ void Config::SetPostprocessing() {
                 case CSD_SN_FOKKERPLANCK_TRAFO_SOLVER_2D:       // Fallthrough
                 case CSD_SN_FOKKERPLANCK_TRAFO_SH_SOLVER_2D:    // Fallthrough
                 case CSD_SN_SOLVER:                             // Fallthrough
+                case CSD_PN_SOLVER:                             // Fallthrough
                     supportedGroups = { MINIMAL, MEDICAL };
                     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
     pos = str.find( "=" );
     if( pos == string::npos ) {
         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 );
         throw( -1 );
     }

code/src/problems/isotropicsource2d.cpp

@@ -12,28 +12,49 @@ std::vector<VectorVector> IsotropicSource2D::GetExternalSource( const Vector& en
 
 VectorVector IsotropicSource2D::SetupIC() {
     VectorVector psi( _mesh->GetNumCells(), Vector( _settings->GetNQuadPoints(), 1e-10 ) );
-    auto cellMids         = _mesh->GetCellMidPoints();
-    double enterPositionX = 0.0;
-    double enterPositionY = 0.5;
-    auto boundaryCells    = _mesh->GetBoundaryTypes();
-
+    auto cellMids = _mesh->GetCellMidPoints();
+    // double enterPositionX = 0.5;    // 0.0;
+    // double enterPositionY = 0.5;
+    // auto boundaryCells    = _mesh->GetBoundaryTypes();
+    // Case 1: Ingoing radiation in just one cell
     // find cell that best matches enter position
-    double dist          = 1000.0;
-    unsigned indexSource = 0;
+    // double dist = 1000.0;
+    // unsigned indexSource = 0;
     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 y = cellMids[j][1] - enterPositionY;
-            if( sqrt( x * x + y * y ) < dist ) {
-                dist        = sqrt( x * x + y * y );
-                indexSource = j;
-            }
+            psi[j]   = 1.0 / ( 4.0 * M_PI * t ) * std::exp( -( x * x + y * y ) / ( 4 * t ) );
         }
-    }
-    psi[indexSource] = Vector( _settings->GetNQuadPoints(), 1.0 );
+        */
+
     return psi;
 }
 
 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;
 }

code/src/solvers/csdpnsolver.cpp

@@ -49,8 +49,12 @@ CSDPNSolver::CSDPNSolver( Config* settings ) : PNSolver( settings ) {
         double y            = _cellMidPoints[i][1];
         const double stddev = .005;
         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() );
     for( unsigned i = 0; i < _nSystem; ++i ) {