tidy up of input files. Some minor documentation

code/include/problems/checkerboard.h

@@ -6,13 +6,13 @@
 class Checkerboard_SN : public ProblemBase
 {
   private:
-    Vector _scatteringXS;     /*! @brief Vector of scattering crosssections */
+    Vector _scatteringXS; /*! @brief Vector of scattering crosssections */
     Vector _totalXS;      /*! @brief Vector of total crosssections */
 
     Checkerboard_SN() = delete;
 
     bool isAbsorption( const Vector& pos ) const; /*! @return True if pos is in absorption region, False otherwise */
-    bool isSource( const Vector& pos ) const; /*! @return True if pos is in source region, False otherwise */
+    bool isSource( const Vector& pos ) const;     /*! @return True if pos is in source region, False otherwise */
 
   public:
     Checkerboard_SN( Config* settings, Mesh* mesh );
@@ -27,13 +27,13 @@ class Checkerboard_SN : public ProblemBase
 class Checkerboard_PN : public ProblemBase
 {
   private:
-    Vector _scatteringXS;     /*! @brief Vector of scattering crosssections */
-    Vector _totalXS;      /*! @brief Vector of total crosssections */
+    Vector _scatteringXS; /*! @brief Vector of scattering crosssections len: numCells  */
+    Vector _totalXS;      /*! @brief Vector of total crosssections. len: numCells*/
 
     Checkerboard_PN() = delete;
 
     bool isAbsorption( const Vector& pos ) const; /*! @return True if pos is in absorption region, False otherwise */
-    bool isSource( const Vector& pos ) const; /*! @return True if pos is in source region, False otherwise */
+    bool isSource( const Vector& pos ) const;     /*! @return True if pos is in source region, False otherwise */
 
     /**
      * @brief Gets the global index for given order l of Legendre polynomials and given

code/include/solvers/solverbase.h

@@ -32,11 +32,11 @@ class Solver
     double _dE;                      /*! @brief energy/time step size */
     Vector _energies;                // energy groups used in the simulation [keV]
     std::vector<double> _density;    // patient density, dim(_density) = _nCells
-    Vector _s;                       // stopping power, dim(_s) = _nTimeSteps
+    Vector _s;                       // stopping power, dim(_s) = _nTimeSteps (only for csdsolver)
     std::vector<VectorVector> _Q;    /*!  @brief  external source term */
 
-    VectorVector _sigmaS; /*!  @brief scattering cross section for all energies */
-    VectorVector _sigmaT; /*!  @brief total cross section for all energies */
+    VectorVector _sigmaS; /*!  @brief scattering cross section for all energies. len: _nEnergies x numCells */
+    VectorVector _sigmaT; /*!  @brief total cross section for all energies.  len: _nEnergies x numCells*/
 
     // quadrature related numbers
     QuadratureBase* _quadrature; /*! @brief quadrature to create members below */

code/include/toolboxes/textprocessingtoolbox.h

@@ -57,7 +57,7 @@ inline void PrintVectorVector( const VectorVector vectorIn ) {
             else
                 printf( "| %.2f ", vectorIn[idx_Outer][idx_Inner] );
         }
-        std::cout << " |\n";
+        printf( " |\n" );
     }
 }
 

code/input/1DMesh.geo

-// settings Kersting
-lc = 0.005;
-width = 1.0;
-length = 10.0;
-start = -2.5;
- //+
-Point(1) = {start, width, 0, lc};
-//+
-Point(2) = {start+length, width, 0, lc};
-//+
-Point(3) = {start, 0, 0, lc};
-//+
-Point(4) = {start+length, 0, 0, lc}; 
-//+
-Line(1) = {1, 2};
-//+
-Line(2) = {2, 4};
-//+
-Line(3) = {4, 3};
-//+
-Line(4) = {3, 1};
-//+
-Curve Loop(1) = {4, 1, 2, 3};
-//+
-Plane Surface(1) = {1};
-//+
-Physical Curve("dirichlet") = {4, 2};
-//+
-Physical Curve("wall_low") = {1};
-Physical Curve("wall_up") = {3};
-//+
-Transfinite Surface {1} = {1, 2, 4, 3};
-//+
-Transfinite Curve {4, 2} = 1 Using Progression 1;
-//+
-Transfinite Curve {1, 3} = 1000 Using Progression 1;
-//+
-Recombine Surface {1};

code/input/checkerboard.cfg

-% ---- File specifications ----
-OUTPUT_DIR = ../result
-OUTPUT_FILE = checkerboard
-LOG_DIR = ../result/logs
-MESH_FILE = checkerboard.su2
-
-% ---- Solver specifications ----
-CFL_NUMBER = 0.5
-TIME_FINAL = 3.2
-PROBLEM = CHECKERBOARD
-
-% ---- Boundary Conditions ----
-BC_NEUMANN = ( void )
-QUAD_TYPE = GAUSS_LEGENDRE_TENSORIZED
-QUAD_ORDER = 10
-
-% ----- Output ---- 
-%
-VOLUME_OUTPUT = (MINIMAL)
-VOLUME_OUTPUT_FREQUENCY = 1
-SCREEN_OUTPUT = (ITER, MASS, RMS_FLUX, VTK_OUTPUT, CSV_OUTPUT)
-SCREEN_OUTPUT_FREQUENCY = 1
-HISTORY_OUTPUT = (ITER, MASS, RMS_FLUX, VTK_OUTPUT, CSV_OUTPUT)
-HISTORY_OUTPUT_FREQUENCY = 1

code/input/checkerboard.py

-import pygmsh as pg
-import numpy as np
-import itertools
-import os
-
-def add_block(x0,y0,length,char_length,geom):
-    coords = np.array([
-        [x0, y0, 0.0],
-        [x0+length, y0, 0.0],
-        [x0+length, y0+length, 0.0],
-        [x0, y0+length, 0.0]
-    ])
-    return geom.add_polygon(coords, char_length)
-
-char_length = 0.075
-geom = pg.opencascade.Geometry()
-domain = add_block(0, 0, 7, char_length, geom)
-xpos = ypos = [1, 2, 3, 4, 5]
-pattern = list(itertools.product(xpos, ypos))[::2]
-pattern.pop(7)
-boxes = [domain]
-for pos in pattern:
-    boxes.append(add_block(pos[0], pos[1], 1, char_length, geom))
-geom.boolean_fragments(boxes,[])
-geom.add_physical(domain.lines, label="void")
-
-mesh_code = geom.get_code()
-with open("checkerboard.geo","w") as mesh_file:
-    mesh_file.write(mesh_code)
-os.system('gmsh checkerboard.geo -2 -format su2 -save_all')
-os.system('rm checkerboard.geo')

code/input/exampleCSDSN.cfg

@@ -12,7 +12,7 @@ OUTPUT_DIR = ../result
 OUTPUT_FILE = example_csd1
 LOG_DIR = ../result/logs
 % Mesh File
-MESH_FILE = 1DMesh.su2
+MESH_FILE = meshes/1DMesh.su2
 %
 PROBLEM = WATERPHANTOM
 %

code/input/exampleDataGen.cfg

-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%  Data Generator for                  % 
-%   Neural Entropy Closure             %
-%  Author <Steffen Schotthöfer>		   %
-%  Date   17.12.2020                   %
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%
-% ---- Global settings ----
-DATA_GENERATOR_MODE = YES
-TRAINING_SET_SIZE = 500
-MAX_VALUE_FIRST_MOMENT = 20
-REALIZABLE_SET_EPSILON_U0 = 0.05
-REALIZABLE_SET_EPSILON_U1 = 0.97
-%
-% ---- File specifications ----
-% 
-% Output directory
-OUTPUT_DIR = ../result
-% Output file
-OUTPUT_FILE = trainM0
-% Log directory
-LOG_DIR = ../result/logs
-%
-% --- Spherical Basis ----
-%
-% Choice of basis functions
-SPHERICAL_BASIS = SPHERICAL_MONOMIALS
-%
-% Maximal Moment degree
-MAX_MOMENT_SOLVER = 1
-%
-% --- Entropy settings ----
-ENTROPY_FUNCTIONAL = MAXWELL_BOLTZMANN
-ENTROPY_OPTIMIZER = NEWTON
-%
-% ----- Newton Solver Specifications ----
-%
-NEWTON_FAST_MODE = NO
-NEWTON_ITER = 1000000
-NEWTON_EPSILON = 0.01
-NEWTON_STEP_SIZE = 0.7
-NEWTON_LINE_SEARCH_ITER = 1000
-%
-% ----- Quadrature ----
-% 
-% Quadrature Rule
-QUAD_TYPE = GAUSS_LEGENDRE_TENSORIZED
-% Quadrature Order
-QUAD_ORDER = 8
-%

code/input/exampleMN.cfg

-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%  Linesource Benchmarking File MN     %
-%  Author <Steffen Schotthöfer>		   %
-%  Date   01.07.2020                   %
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%
-%
-% ---- File specifications ----
-% 
-% Output directory
-OUTPUT_DIR = ../result
-% Output file
-OUTPUT_FILE = exampleMN_mono
-% Log directory
-LOG_DIR = ../result/logs
-% Mesh File
-MESH_FILE = linesource.su2
-%MESH_FILE = linesource_debug.su2
-%
-% ---- Problem description ---
-%
-PROBLEM = LINESOURCE
-SCATTER_COEFF = 1
-%
-% ---- Solver specifications ----
-%
-% Choice of basis functions
-SPHERICAL_BASIS = SPHERICAL_MONOMIALS
-% Maximal Moment degree
-MAX_MOMENT_SOLVER = 1
-% Solver type
-SOLVER = MN_SOLVER
-% CFL number
-CFL_NUMBER = 0.7
-% Final time for simulation
-TIME_FINAL = 0.5
-
-% Reconstruction order (spatial flux)
-RECONS_ORDER = 1
-%
-% ---- Entropy settings ----
-%
-ENTROPY_FUNCTIONAL = MAXWELL_BOLTZMANN
-ENTROPY_OPTIMIZER = NEWTON
-NEURAL_MODEL = 4
-%
-% ----- Newton Solver Specifications ----
-%
-NEWTON_FAST_MODE = NO
-NEWTON_ITER = 1000
-NEWTON_EPSILON = 0.01
-NEWTON_STEP_SIZE = 0.7
-NEWTON_LINE_SEARCH_ITER = 1000
-%
-% ---- Boundary Conditions ----
-%
-% Example: BC_DIRICLET = (dummyMarker1, dummyMarker2)
-% Dirichlet Boundary
-BC_DIRICHLET = ( void )
-%
-% ----- Quadrature ----
-% 
-% Quadrature Rule
-QUAD_TYPE = GAUSS_LEGENDRE_TENSORIZED
-% Quadrature Order
-QUAD_ORDER = 8
-%
-% ----- Output ---- 
-%
-VOLUME_OUTPUT = (MINIMAL, MOMENTS, DUAL_MOMENTS)
-VOLUME_OUTPUT_FREQUENCY = 1
-SCREEN_OUTPUT = (ITER, MASS, RMS_FLUX, VTK_OUTPUT, CSV_OUTPUT)
-SCREEN_OUTPUT_FREQUENCY = 1
-HISTORY_OUTPUT = (ITER, MASS, RMS_FLUX, VTK_OUTPUT, CSV_OUTPUT)
-HISTORY_OUTPUT_FREQUENCY = 1

code/input/examplePN.cfg

-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%  Linesource Benchmarking File PN     %
-%  Author <Steffen Schotthöfer>		   %
-%  Date   01.07.2020                   %
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%
-%
-% ---- File specifications ----
-% 
-% Output directory
-OUTPUT_DIR = ../result
-% Output file
-OUTPUT_FILE = examplePN
-% Log directory
-LOG_DIR = ../result/logs
-% Mesh File
-MESH_FILE = linesource.su2
-%MESH_FILE = linesource_debug.su2
-
-%
-PROBLEM = LINESOURCE
-SCATTER_COEFF = 1
-%
-%
-% ---- Solver specifications ----
-%
-% Solver type
-SOLVER = PN_SOLVER
-% CFL number
-CFL_NUMBER = 0.7
-% Final time for simulation
-TIME_FINAL = 0.3
-% Maximal Moment degree
-MAX_MOMENT_SOLVER = 2
-% Reconstruction order
-RECONS_ORDER = 2
-
-% ---- Boundary Conditions ----
-% Example: BC_DIRICLET = (dummyMarker1, dummyMarker2)
-% Dirichlet Boundary
-BC_DIRICHLET = ( void )
-%
-% ----- Output ---- 
-%
-VOLUME_OUTPUT = (MINIMAL)
-

code/input/exampleSN.cfg

-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%  Linesource Benchmarking File SN     %
-%  Author <Steffen Schotthöfer>		   %
-%  Date   01.10.2020                   %
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%
-%
-% ---- File specifications ----
-% 
-% Output directory
-OUTPUT_DIR = ../result
-% Output file
-OUTPUT_FILE = example_SN
-% Log directory
-LOG_DIR = ../result/logs
-% Mesh File
-MESH_FILE = linesource.su2
-%
-PROBLEM = LINESOURCE
-%
-% ---- Solver specifications ----
-%
-% Solver type
-SOLVER = SN_SOLVER
-% CFL number
-CFL_NUMBER = 0.7
-% Final time for simulation
-TIME_FINAL = 0.5
-% Reconstruction order
-RECONS_ORDER = 2
-%
-
-% ---- Boundary Conditions ----
-% Example: BC_DIRICLET = (dummyMarker1, dummyMarker2)
-% Dirichlet Boundary
-BC_DIRICHLET = ( void )
-%
-%% Quadrature Specifications
-% Quadrature Type
-QUAD_TYPE = GAUSS_LEGENDRE_TENSORIZED
-% Quadrature Order
-QUAD_ORDER = 8
-%
-% ----- Output ---- 
-%
-VOLUME_OUTPUT = (ANALYTIC,  MINIMAL)
-VOLUME_OUTPUT_FREQUENCY = 1

code/input/exampleWaterphantomCSD.cfg

@@ -4,7 +4,7 @@
 OUTPUT_DIR = ../result
 OUTPUT_FILE = example_csdSN_1d_NOScatter
 LOG_DIR = ../result/logs
-MESH_FILE = 1DMesh.su2
+MESH_FILE = meshes/1DMesh.su2
 PROBLEM = WATERPHANTOM
 SOLVER = CSD_SN_SOLVER
 

code/input/linesource.py

-import pygmsh as pg
-import numpy as np
-import itertools
-import os
-
-def add_block(x0,y0,lengthX, lengthY,char_length,geom):
-    coords = np.array([
-        [x0, y0, 0.0],
-        [x0+lengthX, y0, 0.0],
-        [x0+lengthX, y0+lengthY, 0.0],
-        [x0, y0+lengthY, 0.0]
-    ])
-    return geom.add_polygon(coords, char_length)
-
-char_length = 0.015
-geom = pg.opencascade.Geometry()
-domain = add_block(-1, -1, 2,2, char_length, geom)
-geom.add_raw_code('psource = newp;\nPoint(psource) = {0.0, 0.0, 0.0, '+str(char_length)+'};\nPoint{psource} In Surface{'+domain.id+'};')
-geom.add_physical(domain.lines, label="void")
-
-mesh_code = geom.get_code()
-with open("linesource.geo","w") as mesh_file:
-    mesh_file.write(mesh_code)
-os.system('gmsh linesource.geo -2 -format su2 -save_all')
-os.system('rm linesource.geo')

code/src/problems/checkerboard.cpp

@@ -95,7 +95,7 @@ std::vector<VectorVector> Checkerboard_PN::GetExternalSource( const Vector& /*en
     VectorVector Q( _mesh->GetNumCells(), Vector( 1u, 0.0 ) );
     auto cellMids = _mesh->GetCellMidPoints();
     for( unsigned j = 0; j < cellMids.size(); ++j ) {
-        if( isSource( cellMids[j] ) ) Q[j] = 1.0 / std::sqrt( 4 * M_PI );    // isotropic source
+        if( isSource( cellMids[j] ) ) Q[j] = 10.0 / std::sqrt( 4 * M_PI );    // isotropic source
     }
     return std::vector<VectorVector>( 1u, Q );
 }

code/src/solvers/mnsolver.cpp

@@ -179,9 +179,6 @@ void MNSolver::FluxUpdate() {
         // Dirichlet Boundaries stayd
         if( _boundaryCells[idx_cell] == BOUNDARY_TYPE::DIRICHLET ) continue;
         _solNew[idx_cell] = ConstructFlux( idx_cell );
-        // if( norm( _solNew[idx_cell] ) > 0.001 ) {
-        //    std::cout << _solNew[idx_cell] << "\n";
-        //}
     }
 }
 
@@ -193,10 +190,9 @@ void MNSolver::FVMUpdate( unsigned idx_energy ) {
         if( _boundaryCells[idx_cell] == BOUNDARY_TYPE::DIRICHLET ) continue;
         // Flux update
         for( unsigned idx_sys = 0; idx_sys < _nTotalEntries; idx_sys++ ) {
-            _solNew[idx_cell][idx_sys] = _sol[idx_cell][idx_sys] - ( _dE / _areas[idx_cell] ) * _solNew[idx_cell][idx_sys] /* cell averaged flux */
-                                         - _dE * _sol[idx_cell][idx_sys] *
-                                               ( _sigmaT[idx_energy][idx_cell]                                 /* absorbtion influence */
-                                                 + _sigmaS[idx_energy][idx_cell] * _scatterMatDiag[idx_sys] ); /* scattering influence */
+            _solNew[idx_cell][idx_sys] = _sol[idx_cell][idx_sys] - ( _dE / _areas[idx_cell] ) * _solNew[idx_cell][idx_sys]; /* cell averaged flux */
+            -_dE* _sol[idx_cell][idx_sys] * ( _sigmaT[idx_energy][idx_cell]                                                 /* absorbtion influence */
+                                              + _sigmaS[idx_energy][idx_cell] * _scatterMatDiag[idx_sys] );                 /* scattering influence */
         }
         // Source Term
         _solNew[idx_cell][0] += _dE * _Q[0][idx_cell][0];

code/src/solvers/solverbase.cpp

@@ -10,6 +10,7 @@
 #include "solvers/csdsolvertrafofp.h"
 #include "solvers/csdsolvertrafofp2d.h"
 #include "solvers/csdsolvertrafofpsh2d.h"
+#include "toolboxes/textprocessingtoolbox.h"
 
 #include "solvers/mnsolver.h"
 #include "solvers/pnsolver.h"