merged current master

.gitlab-ci.yml

@@ -8,7 +8,7 @@ unit_tests:
   script:
     - git submodule update --init --recursive
     - cd code/build/release
-    - cmake -G Ninja -DCMAKE_BUILD_TYPE=Debug ../../
+    - cmake -G Ninja -DCMAKE_BUILD_TYPE=Debug -Dbuild_tests=ON ../../
     - ninja
     - OMP_NUM_THREADS=1 ./unit_tests -r junit > unit_tests_report.xml
     - gcovr -e ../../ext/ -e ../../tests/ -r ../../

code/CMakeLists.txt

 cmake_minimum_required( VERSION 3.12.4 )
-project( KiT-RT VERSION 0.0.2 LANGUAGES CXX )
+project( KiT-RT VERSION 0.1.0 LANGUAGES CXX )
 
 ### OPTIONS #####################################
 option( build_tests "enables compiling of unit_tests" OFF )

code/include/common/optionstructure.h

@@ -22,7 +22,7 @@ class OptionBase
     std::vector<std::string> _value; /*! @brief: String name of the option */
 
   public:
-    OptionBase(){};
+    OptionBase() {}
     virtual ~OptionBase() = 0;
 
     virtual std::string SetValue( std::vector<std::string> value ); /*! @brief: Set string name of the option */
@@ -45,7 +45,7 @@ class OptionDouble : public OptionBase
   public:
     OptionDouble( std::string option_field_name, double& option_field, double default_value );
 
-    ~OptionDouble() override{};
+    ~OptionDouble() override {}
 
     std::string SetValue( std::vector<std::string> option_value ) override;
 
@@ -61,7 +61,7 @@ class OptionString : public OptionBase
   public:
     OptionString( std::string option_field_name, std::string& option_field, std::string default_value );
 
-    ~OptionString() override{};
+    ~OptionString() override {}
 
     std::string SetValue( std::vector<std::string> option_value ) override;
 
@@ -77,7 +77,7 @@ class OptionInt : public OptionBase
   public:
     OptionInt( std::string option_field_name, int& option_field, int default_value );
 
-    ~OptionInt() override{};
+    ~OptionInt() override {}
 
     std::string SetValue( std::vector<std::string> option_value ) override;
 
@@ -93,7 +93,7 @@ class OptionULong : public OptionBase
   public:
     OptionULong( std::string option_field_name, unsigned long& option_field, unsigned long default_value );
 
-    ~OptionULong() override{};
+    ~OptionULong() override {}
 
     std::string SetValue( std::vector<std::string> option_value ) override;
 
@@ -109,7 +109,7 @@ class OptionUShort : public OptionBase
   public:
     OptionUShort( std::string option_field_name, unsigned short& option_field, unsigned short default_value );
 
-    ~OptionUShort() override{};
+    ~OptionUShort() override {}
 
     std::string SetValue( std::vector<std::string> option_value ) override;
 
@@ -125,7 +125,7 @@ class OptionLong : public OptionBase
   public:
     OptionLong( std::string option_field_name, long& option_field, long default_value );
 
-    ~OptionLong() override{};
+    ~OptionLong() override {}
 
     std::string SetValue( std::vector<std::string> option_value ) override;
 
@@ -141,7 +141,7 @@ class OptionBool : public OptionBase
   public:
     OptionBool( std::string option_field_name, bool& option_field, bool default_value );
 
-    ~OptionBool() override{};
+    ~OptionBool() override {}
 
     std::string SetValue( std::vector<std::string> option_value ) override;
 
@@ -157,7 +157,7 @@ class OptionStringList : public OptionBase
   public:
     OptionStringList( std::string option_field_name, unsigned short& list_size, std::vector<std::string>& option_field );
 
-    ~OptionStringList() override{};
+    ~OptionStringList() override {}
 
     std::string SetValue( std::vector<std::string> option_value ) override;
 

code/include/icru.h → code/include/problems/icru.h

@@ -5,6 +5,7 @@
  * Disclaimer: This is copied from the code distributed with the ICRU 77 report
  * Link: https://icru.org/home/reports/elastic-scattering-of-electrons-and-positrons-icru-report-77
  */
+
 #ifndef ICRU_H
 #define ICRU_H
 

code/include/quadratures/productquadrature.h → code/include/quadratures/qproduct.h

+/*! @file: qproduct.h
+ *  @brief: Product quadrature implementation. Implementation is done accordingly to Kendall Atkinson 1981, Australian Matematical Society.
+ *  @author: J. Kusch
+ */
+
 #ifndef PRODUCTQUADRATURE_H
 #define PRODUCTQUADRATURE_H
 
 #include "quadraturebase.h"
 
-class ProductQuadrature : public QuadratureBase
+class QProduct : public QuadratureBase
 {
-    // Implementation is done accordingly to Kendall Atkinson 1981, Australian Matematical Society.
   private:
     double Pythag( const double a, const double b );
     std::pair<Vector, Matrix> ComputeEigenValTriDiagMatrix( const Matrix& mat );
     bool CheckOrder();
 
   public:
-    ProductQuadrature( Config* settings );
-    ProductQuadrature( unsigned order );
+    QProduct( Config* settings );
+    QProduct( unsigned order );
 
-    virtual ~ProductQuadrature() {}
+    virtual ~QProduct() {}
 
     inline void SetName() override { _name = "Product quadrature"; }
     inline void SetNq() override { _nq = 4 * pow( GetOrder(), 2 ); }

code/include/quadratures/quadraturebase.h

+/*! @file: quadraturebase.h
+ *  @brief: Base class for all quadrature rules in KiT-RT
+ *  @author: S. Schotthöfer
+ */
+
 #ifndef QUADRATURE_H
 #define QUADRATURE_H
 
@@ -9,8 +14,13 @@ class Config;
 class QuadratureBase
 {
   public:
+    /*! @brief: Constructor using settings class. This is the recommended constructor.
+     *  @param: Config* settings: Settings class storing all important options.
+     */
     QuadratureBase( Config* settings );
-    QuadratureBase( unsigned order );    // Only for productquadrature... refactor this!
+    /*! @brief: Constructor using directly the order of the quadrature. Not applicable for GaussLegendre, that need additional options.
+                It sets member _settings = nulltpr.*/
+    QuadratureBase( unsigned order );
     virtual ~QuadratureBase() {}
 
     // Aux functions
@@ -37,13 +47,13 @@ class QuadratureBase
     /*! @brief Creates a quadrature rule with a given name and a given order.
      *  @param: Config* settings: Settings to handle quadrature options
      *  @returns Quadrature* quadrature: returns pointer to instance of the given derived quadrature class */
-    static QuadratureBase* CreateQuadrature( Config* settings );
+    static QuadratureBase* Create( Config* settings );
 
     /*! @brief Creates a quadrature rule with a given name and a given order.
      *  @param: name: name of quadrature as enum
      *  @param: quadOrder: order of quadrature
      *  @returns Quadrature* quadrature: returns pointer to instance of the given derived quadrature class */
-    static QuadratureBase* CreateQuadrature( QUAD_NAME name, unsigned quadOrder );
+    static QuadratureBase* Create( QUAD_NAME name, unsigned quadOrder );
 
     // Getter
     inline std::string GetName() const { return _name; }      /*! @returns std::string _name:  name of the quadrature */

code/include/solvers/csdsnsolver.h

 #ifndef CSDSNSOLVER_H
 #define CSDSNSOLVER_H
 
-#include "icru.h"
 #include "solvers/snsolver.h"
 
 class CSDSNSolver : public SNSolver

code/include/solvers/csdsnsolverfp.h

 #ifndef CSDSNSOLVERFP_H
 #define CSDSNSOLVERFP_H
 
-#include "icru.h"
 #include "solvers/snsolver.h"
 
 class Physics;
@@ -22,9 +21,9 @@ class CSDSNSolverFP : public SNSolver
     Matrix _IL; /*!  @brief Laplace Beltrami Matrix */
 
     double _alpha; /*!  @brief  Coefficient of GFP operators (see Olbrant 2010, eq. (8)*/
-    double _beta; /*!  @brief  Coefficient of GFP operators (see Olbrant 2010, eq. (8)*/
+    double _beta;  /*!  @brief  Coefficient of GFP operators (see Olbrant 2010, eq. (8)*/
 
-    Matrix _xi;  /*!  @brief matrix of transport coefficients */
+    Matrix _xi; /*!  @brief matrix of transport coefficients */
 
     bool _RT; /*!  @brief radiotherapy application (on/off), if true use crosssections + stopping powers from database  */
 

code/include/solvers/csdsnsolvernotrafo.h

 #ifndef CSDSNSOLVERNOTRAFO_H
 #define CSDSNSOLVERNOTRAFO_H
 
-#include "icru.h"
 #include "solvers/snsolver.h"
 
 class Physics;

code/include/solvers/csdsolvertrafofp.h

 #ifndef CSDSOLVERTRAFOFP_H
 #define CSDSOLVERTRAFOFP_H
 
-#include "icru.h"
 #include "solvers/snsolver.h"
 
 class Physics;
@@ -21,15 +20,16 @@ class CSDSolverTrafoFP : public SNSolver
     Matrix _L;  /*!  @brief Laplace Beltrami Matrix */
     Matrix _IL; /*!  @brief Laplace Beltrami Matrix */
 
-    double _alpha; /*!  @brief  Coefficient of GFP operators (see Olbrant 2010, Appendix B)*/
+    double _alpha;  /*!  @brief  Coefficient of GFP operators (see Olbrant 2010, Appendix B)*/
     double _alpha2; /*!  @brief  Coefficient of GFP operators (see Olbrant 2010, Appendix B)*/
-    double _beta; /*!  @brief  Coefficient of GFP operators (see Olbrant 2010, Appendix B)*/
+    double _beta;   /*!  @brief  Coefficient of GFP operators (see Olbrant 2010, Appendix B)*/
 
-    Matrix _xi;  /*!  @brief matrix of transport coefficients */
+    Matrix _xi; /*!  @brief matrix of transport coefficients */
     Vector _xi1;
     Vector _xi2;
-    
-    unsigned _FPMethod; /*!  @brief Encodes different ways of computing coefficients alpha, alpha2 & beta, _FPMethod == 1, 2 ,3 stand for methods with increasing accuracy (see Olbrant 2010, Appendix B)*/
+
+    unsigned _FPMethod; /*!  @brief Encodes different ways of computing coefficients alpha, alpha2 & beta, _FPMethod == 1, 2 ,3 stand for methods with
+                           increasing accuracy (see Olbrant 2010, Appendix B)*/
 
     bool _RT; /*!  @brief radiotherapy application (on/off), if true use crosssections + stopping powers from database  */
 

code/include/solvers/csdsolvertrafofp2d.h

 #ifndef CSDSOLVERTRAFOFP2D_H
 #define CSDSOLVERTRAFOFP2D_H
 
-#include "icru.h"
 #include "solvers/snsolver.h"
 
 class Physics;

code/include/solvers/csdsolvertrafofpsh2d.h

@@ -8,12 +8,11 @@
 #include "common/config.h"
 #include "common/io.h"
 #include "fluxes/numericalflux.h"
-#include "icru.h"
 #include "kernels/scatteringkernelbase.h"
 #include "problems/problembase.h"
 #include "quadratures/quadraturebase.h"
 #include "solvers/snsolver.h"
-#include "sphericalharmonics.h"
+#include "toolboxes/sphericalharmonics.h"
 
 class Physics;
 

code/include/solvers/mnsolver.h

@@ -42,6 +42,9 @@ class MNSolver : public Solver
                                            Layout: _nCells x _nTotalEntries*/
     OptimizerBase* _optimizer; /*! @brief: Class to solve minimal entropy problem */
 
+    VectorVector _solDx; /*! @brief:  temporary storage of x-derivatives of alpha */
+    VectorVector _solDy; /*! @brief:  temporary storage of y-derivatives of alpha */
+
     // ---- Private Member functions ---
 
     // IO

code/include/solvers/pnsolver.h

 #ifndef PNSOLVER_H
 #define PNSOLVER_H
 
-#include "solverbase.h"
+#include "solvers/solverbase.h"
 
 class PNSolver : public Solver
 {
@@ -38,6 +38,9 @@ class PNSolver : public Solver
     Vector _scatterMatDiag; /*! @brief: diagonal of the scattering matrix (its a diagonal matrix by construction). Contains eigenvalues of the
                                scattering kernel.  */
 
+    VectorVector _solDx; /*! @brief:  temporary storage of x-derivatives of solution */
+    VectorVector _solDy; /*! @brief:  temporary storage of y-derivatives of solution */
+
     // ---- Member functions ----
 
     // IO

code/include/solvers/solverbase.h

@@ -14,6 +14,7 @@ class Mesh;
 class Config;
 class ProblemBase;
 class QuadratureBase;
+class Reconstructor;
 
 class Solver
 {
@@ -51,6 +52,14 @@ class Solver
     /*! @brief edge neighbor cell ids, dim(_neighbors) = (_NCells,nEdgesPerCell) */
     std::vector<std::vector<unsigned>> _neighbors;
 
+    // slope related params
+    Reconstructor* _reconstructor; /*! @brief reconstructor object for high-order scheme */
+    unsigned _reconsOrder; /*! @brief reconstruction order (current: 1 & 2) */
+    VectorVector _psiDx; /*! @brief slope of solutions in X direction */
+    VectorVector _psiDy; /*! @brief slope of solutions in Y direction */
+    VectorVector _cellMidPoints; /*! @brief middle point locations of elements */
+    std::vector<std::vector<Vector>> _interfaceMidPoints; /*! @brief middle point locations of edges */
+
     // Solution related members
     VectorVector _sol;                 /*! @brief solution of the PDE, e.g. angular flux or moments */
     std::vector<double> _solverOutput; /*! @brief LEGACY: Outputfield for solver ==> Will be replaced by _outputFields in the near future */

code/include/toolboxes/reconstructor.h

@@ -15,12 +15,18 @@ class Config;
 
 class Reconstructor
 {
+  protected:
+    unsigned _reconsOrder;
+
   public:
     /**
      * @brief Reconstruction
      * @param settings
      */
     Reconstructor( Config* settings );
+    static Reconstructor* Create( Config* settings );
+
+    unsigned inline GetReconsOrder() { return _reconsOrder; }
 
     /*! Method 1: structured developing
      * @brief Slope of angular flux psi inside a given cell

code/input/exampleMN.cfg

@@ -36,6 +36,8 @@ MAX_MOMENT_SOLVER = 0
 %% Entropy settings
 ENTROPY_FUNCTIONAL = MAXWELL_BOLTZMANN
 ENTROPY_OPTIMIZER = NEWTON
+% Reconstruction order
+RECONS_ORDER = 2
 %
 % ----- Newton Solver Specifications ----
 %