Merge branch 'master' into 'doc'

Pull updates from master into doc

See merge request !38

README.md

@@ -98,7 +98,7 @@ Append `HINTS VTK_INSTALL_DIR` to the `find_package( VTK ... )` line in the CMak
 find_package( VTK REQUIRED COMPONENTS vtkIOGeometry vtkFiltersCore HINTS ~/VTK-install )
 ```
 
-Compile it
+Compile it (make sure you use GCC as compiler: module load compiler/gnu/10.2) 
 ```bash
 module load devel/cmake/3.16
 module load compiler/gnu/9.2
@@ -164,6 +164,8 @@ Please stick to the following coding style for easier code readability:
  - class variables start with an underscore and lowercase letters e.g. `_foo`
  - functions start with a capital letter e.g. `GetSettings()`
  - any variable/function names have capital letters at each individual word e.g. `GetAllCellsAdjacentTo(Cell i)`
+ - Abstract base classes for inheritance structures end with the suffix "Base" e.g. QuadratureBase, whith child classes 
+   QGaussLegendre or QProduct
 
 Please also use the provided `code/.clang-format` style format to format your code before pushing your latest commits.
 Some editors offer to automatically apply the style format upon saving a file (e.g. `Qtcreator`).

Catch2 @ 0e9bae1c

-Subproject commit ff349a50bfc6214b4081f4ca63c7de35e2162f60
+Subproject commit 0e9bae1cdb5d7bcd967bcab10539805cc3bccb32

blaze @ 69e2fe55

-Subproject commit 0380370f0626de2ad0ec7ea815803e22bbf6b42e
+Subproject commit 69e2fe554e72d90efc39ceec72d8cbddf88802a8

spdlog @ 83b91499

-Subproject commit cbe9448650176797739dbab13961ef4c07f4290f
+Subproject commit 83b9149930f392d7797b54fe97a66ab3f2120671

code/include/common/config.h

 /*!
  * @file config.h
- * @brief Class to handle all options and their pre and postprocessing.
+ * @brief: Class to handle all options and their pre and postprocessing.
+ *         DO NOT CREATE SETTERS FOR THIS CLASS! ALL OPTIONS ARE CONSTANT (after SetPostprocessing).
+ *
  * @author S. Schotthöfer
  *
- * Disclaimer: This class structure was copied and modifed with open source permission from SU2 v7.0.3 https://su2code.github.io/
+ * Disclaimer: This class structure was copied and (heavily) modifed with open source permission from SU2 v7.0.3 https://su2code.github.io/
  */
 
 #ifndef CONFIG_H
@@ -67,7 +69,7 @@ class Config
     bool _cleanFluxMat;
     bool _allGaussPts; /*!< @brief If true, the SN Solver uses all Gauss pts in the quadrature */
 
-    bool _csd;                 /*!< @brief If true, continuous slowing down approximation will be used */
+    bool _csd;                 // LEGACY! /*!< @brief If true, continuous slowing down approximation will be used */
     std::string _hydrogenFile; /*!< @brief Name of hydrogen cross section file */
     std::string _oxygenFile;   /*!< @brief Name of oxygen cross section file */
 
@@ -83,12 +85,15 @@ class Config
     // Scattering Kernel
     KERNEL_NAME _kernelName; /*!< @brief Scattering Kernel Name*/
 
+    // Spherical Basis
+    SPHERICAL_BASIS_NAME _sphericalBasisName; /*!< @brief: Name of the basis on the unit sphere */
+
     // Optimizer
     OPTIMIZER_NAME _entropyOptimizerName; /*!< @brief Choice of optimizer */
     double _optimizerEpsilon;             /*!< @brief termination criterion epsilon for Newton Optmizer */
-    unsigned short _newtonIter;           /*!< @brief Maximal Number of newton iterations */
+    unsigned long _newtonIter;            /*!< @brief Maximal Number of newton iterations */
     double _newtonStepSize;               /*!< @brief Stepsize factor for newton optimizer */
-    unsigned short _newtonLineSearchIter; /*!< @brief Maximal Number of line search iterations for newton optimizer */
+    unsigned long _newtonLineSearchIter;  /*!< @brief Maximal Number of line search iterations for newton optimizer */
     bool _newtonFastMode;                 /*!< @brief If true, we skip the NewtonOptimizer for quadratic entropy and assign alpha = u */
 
     // Output Options
@@ -104,6 +109,13 @@ class Config
     std::vector<SCALAR_OUTPUT> _historyOutput; /*!< @brief Output groups for screen output*/
     unsigned short _historyOutputFrequency;    /*!< @brief Frequency of screen output*/
 
+    // Data Generator Settings
+    /*!< @brief Check, if data generator mode is active. If yes, no solver is called, but instead the data generator is executed */
+    bool _dataGeneratorMode;
+    unsigned long _tainingSetSize;         /*!< @brief Size of training data set for data generator */
+    unsigned long _maxValFirstMoment;      /*!< @brief Size of training data set for data generator */
+    double _boundaryDistanceRealizableSet; /*! @brief Distance of the sampled moments to the boundary of the realizable set */
+
     // --- Parsing Functionality and Initializing of Options ---
     /*!
      * @brief Set default values for all options not yet set.
@@ -256,9 +268,9 @@ class Config
 
     //  Optimizer
     double inline GetNewtonOptimizerEpsilon() const { return _optimizerEpsilon; }
-    unsigned inline GetNewtonIter() const { return _newtonIter; }
+    unsigned long inline GetNewtonIter() const { return _newtonIter; }
     double inline GetNewtonStepSize() const { return _newtonStepSize; }
-    unsigned inline GetNewtonMaxLineSearches() const { return _newtonLineSearchIter; }
+    unsigned long inline GetNewtonMaxLineSearches() const { return _newtonLineSearchIter; }
     bool inline GetNewtonFastMode() const { return _newtonFastMode; }
     OPTIMIZER_NAME inline GetOptimizerName() const { return _entropyOptimizerName; }
 
@@ -268,6 +280,8 @@ class Config
     // Scattering Kernel
     KERNEL_NAME inline GetKernelName() const { return _kernelName; }
 
+    // Basis name
+    SPHERICAL_BASIS_NAME inline GetSphericalBasisName() const { return _sphericalBasisName; }
     // Output Structure
     std::vector<VOLUME_OUTPUT> inline GetVolumeOutput() { return _volumeOutput; }
     unsigned short inline GetNVolumeOutput() { return _nVolumeOutput; }
@@ -280,15 +294,20 @@ class Config
     std::vector<SCALAR_OUTPUT> inline GetHistoryOutput() { return _historyOutput; }
     unsigned short inline GetNHistoryOutput() { return _nHistoryOutput; }
     unsigned short inline GetHistoryOutputFrequency() { return _historyOutputFrequency; }
-    // ---- Setters for option structure
 
+    // Data generator
+    bool inline GetDataGeneratorMode() { return _dataGeneratorMode; }
+    unsigned long inline GetTrainingDataSetSize() { return _tainingSetSize; }
+    unsigned long inline GetMaxValFirstMoment() { return _maxValFirstMoment; }
+    double GetBoundaryDistanceRealizableSet() { return _boundaryDistanceRealizableSet; }
+
+    // ---- Setters for option structure
+    // This section is dangerous
     // Quadrature Structure
-    void SetNQuadPoints( unsigned nq ) { _nQuadPoints = nq; }
+    void SetNQuadPoints( unsigned nq ) { _nQuadPoints = nq; }           /*! @brief Never change the nq! This is only for the test framework. */
     void SetQuadName( QUAD_NAME quadName ) { _quadName = quadName; }    /*! @brief Never change the quadName! This is only for the test framework. */
     void SetQuadOrder( unsigned quadOrder ) { _quadOrder = quadOrder; } /*! @brief Never change the quadOrder! This is only for the test framework. */
     void SetSNAllGaussPts( bool useall ) { _allGaussPts = useall; }     /*! @brief Never change the this! This is only for the test framework. */
-    // Mesh Structure
-    void SetNCells( unsigned nCells ) { _nCells = nCells; }
 };
 
 #endif    // CONFIG_H

code/include/common/globalconstants.h

@@ -141,4 +141,10 @@ enum SCALAR_OUTPUT { ITER, MASS, RMS_FLUX, VTK_OUTPUT, CSV_OUTPUT };
 inline std::map<std::string, SCALAR_OUTPUT> ScalarOutput_Map{
     { "ITER", ITER }, { "MASS", MASS }, { "RMS_FLUX", RMS_FLUX }, { "VTK_OUTPUT", VTK_OUTPUT }, { "CSV_OUTPUT", CSV_OUTPUT } };
 
+// Spherical Basis Name
+enum SPHERICAL_BASIS_NAME { SPHERICAL_HARMONICS, SPHERICAL_MONOMIALS };
+
+inline std::map<std::string, SPHERICAL_BASIS_NAME> SphericalBasis_Map{ { "SPHERICAL_HARMONICS", SPHERICAL_HARMONICS },
+                                                                       { "SPHERICAL_MONOMIALS", SPHERICAL_MONOMIALS } };
+
 #endif    // GLOBAL_CONSTANTS_H

code/include/optimizers/mloptimizer.h

@@ -16,8 +16,8 @@ class MLOptimizer : public OptimizerBase
 
     inline ~MLOptimizer();
 
-    void Solve( Vector& lambda, Vector& u, VectorVector& moments, unsigned idx_cell = 0 ) override;
-    void SolveMultiCell( VectorVector& lambda, VectorVector& u, VectorVector& moments ) override;
+    void Solve( Vector& lambda, Vector& u, const VectorVector& moments, unsigned idx_cell = 0 ) override;
+    void SolveMultiCell( VectorVector& lambda, VectorVector& u, const VectorVector& moments ) override;
 
   private:
     double* callNetwork( const unsigned input_size, double* input );

code/include/optimizers/newtonoptimizer.h

@@ -18,21 +18,21 @@ class NewtonOptimizer : public OptimizerBase
 
     ~NewtonOptimizer();
 
-    void Solve( Vector& lambda, Vector& sol, VectorVector& moments, unsigned idx_cell = 0 ) override;
-    void SolveMultiCell( VectorVector& lambda, VectorVector& sol, VectorVector& moments ) override;
+    void Solve( Vector& lambda, Vector& sol, const VectorVector& moments, unsigned idx_cell = 0 ) override;
+    void SolveMultiCell( VectorVector& lambda, VectorVector& sol, const VectorVector& moments ) override;
 
-  private:
     /*! @brief: Computes the objective function
                 grad = <eta(alpha*m)> - alpha*sol */
-    double ComputeObjFunc( Vector& alpha, Vector& sol, VectorVector& moments );
+    double ComputeObjFunc( Vector& alpha, Vector& sol, const VectorVector& moments );
 
+  private:
     /*! @brief: Computes gradient of objective function and stores it in grad
                 grad = <m*eta*'(alpha*m)> - sol */
-    void ComputeGradient( Vector& alpha, Vector& sol, VectorVector& moments, Vector& grad );
+    void ComputeGradient( Vector& alpha, Vector& sol, const VectorVector& moments, Vector& grad );
 
     /*! @brief: Computes hessian of objective function and stores it in hessian
                 grad = <mXm*eta*'(alpha*m)> */
-    void ComputeHessian( Vector& alpha, VectorVector& moments, Matrix& hessian );
+    void ComputeHessian( Vector& alpha, const VectorVector& moments, Matrix& hessian );
 
     QuadratureBase* _quadrature; /*! @brief: used quadrature */    // THis is memory doubling! Try to use a pointer.
     unsigned _nq;                                                  /*! @brief: number of quadrature points */

code/include/optimizers/optimizerbase.h

@@ -26,9 +26,9 @@ class OptimizerBase
     /*! @brief  : Computes the optimal Lagrange multilpiers for the dual entropy minimization problem
      *  @param  : Vector  u = pointer to vector of given moments. // Maybe use pointer for performance?
      *  @return : Vector  alpha = optimal lagrange multipliers. Has the same length as Vector u. */
-    virtual void Solve( Vector& lambda, Vector& u, VectorVector& moments, unsigned idx_cell = 0 ) = 0;
+    virtual void Solve( Vector& lambda, Vector& u, const VectorVector& moments, unsigned idx_cell = 0 ) = 0;
 
-    virtual void SolveMultiCell( VectorVector& lambda, VectorVector& u, VectorVector& moments ) = 0;
+    virtual void SolveMultiCell( VectorVector& lambda, VectorVector& u, const VectorVector& moments ) = 0;
 
   protected:
     EntropyBase* _entropy; /*! @brief: Class to handle entropy functional evaluations */

code/include/problems/linesource.h

@@ -82,16 +82,6 @@ class LineSource_PN : public LineSource
   private:
     LineSource_PN() = delete;
 
-    /**
-     * @brief Gets the global index for given order l of Legendre polynomials and given
-     *        order k of Legendre functions.
-     *        Note: This is code doubling from PNSolver::GlobalIndex
-     * @param l : order of Legendre polynomial
-     * @param k : order of Legendre function
-     * @returns global index
-     */
-    int GlobalIndex( int l, int k ) const;
-
   public:
     LineSource_PN( Config* settings, Mesh* mesh );
     ~LineSource_PN();

code/include/quadratures/qgausslegendre1D.h

+/*!
+ * @file config.h
+ * @brief Class to compute 1D Gauss Legendre Quadrature
+ * @author J. Kusch
+ *
+ * Disclaimer: This class structure was copied and modifed with open source permission from SU2 v7.0.3 https://su2code.github.io/
+ */
+
 #ifndef QGAUSSLEGENDRE1D_H
 #define QGAUSSLEGENDRE1D_H
 
@@ -15,10 +23,26 @@ class QGaussLegendre1D : public QuadratureBase
     QGaussLegendre1D( unsigned quadOrder );
     virtual ~QGaussLegendre1D() {}
 
-    inline void SetName() override { _name = "Tensorized Gauss-Legendre quadrature."; }
+    inline void SetName() override { _name = "Tensorized Gauss-Legendre quadrature 1D."; }
     inline void SetNq() override { _nq = _order; }
     void SetPointsAndWeights() override;
     void SetConnectivity() override;
+
+    /*! @brief Integrates f(x,y,z) with the quadrature.
+     *  @param double(f)( double x0, double x1, double x2 ) : density function that depends on a three spatial dimensions.
+     *  @returns double result: result of the quadrature rule */
+    double Integrate( double( f )( double x0, double x1, double x2 ) ) override;
+
+    /*! @brief Integrates f(x,y,z) with the quadrature.
+     *  @param double(f)( double my, double phi ) : density function that depends on a spherical coordinates.
+     *  @returns double result: result of the quadrature rule */
+    double IntegrateSpherical( double( f )( double my, double phi ) ) override;
+
+    /*! @brief Integrates vector valued f(x,y,z) with the quadrature. Each dimension is integrated by itself.
+     *  @param : double(f)( double x0, double x1, double x2 ) : density function that depends on a three spatial dimensions.
+     *  @param :  len : lenght of vector
+     *  @returns double result: result of the quadrature rule (vector valued) */
+    std::vector<double> Integrate( std::vector<double>( f )( double x0, double x1, double x2 ), unsigned /* len */ ) override;
 };
 
 #endif    // QGAUSSLEGENDRE1D_H

code/include/quadratures/qgausslegendretensorized.h

@@ -20,8 +20,6 @@ class QGaussLegendreTensorized : public QuadratureBase
     void SetNq() override;
     void SetPointsAndWeights() override;
     void SetConnectivity() override;
-
-    inline VectorVector GetPointsSphere() const override { return _pointsSphere; }
 };
 
 #endif    // QGAUSSLEGENDRETENSORIZED_H

code/include/quadratures/qmontecarlo.h

@@ -14,7 +14,6 @@ class QMonteCarlo : public QuadratureBase
     inline void SetNq() override { _nq = GetOrder(); }
     void SetPointsAndWeights() override;
     void SetConnectivity() override;
-    VectorVector GetPointsSphere() const override;
 };
 
 #endif    // QMONTECARLO_H

code/include/quadratures/qproduct.h

@@ -25,8 +25,6 @@ class QProduct : public QuadratureBase
     inline void SetNq() override { _nq = 4 * pow( GetOrder(), 2 ); }
     void SetPointsAndWeights() override;
     void SetConnectivity() override;
-
-    inline VectorVector GetPointsSphere() const override { return _pointsSphere; }
 };
 
 #endif    // PRODUCTQUADRATURE_H

code/include/quadratures/quadraturebase.h

@@ -35,13 +35,18 @@ class QuadratureBase
     /*! @brief Integrates f(x,y,z) with the quadrature.
      *  @param double(f)( double x0, double x1, double x2 ) : density function that depends on a three spatial dimensions.
      *  @returns double result: result of the quadrature rule */
-    double Integrate( double( f )( double x0, double x1, double x2 ) );
+    virtual double Integrate( double( f )( double x0, double x1, double x2 ) );
+
+    /*! @brief Integrates f(x,y,z) with the quadrature.
+     *  @param double(f)( double my, double phi ) : density function that depends on a spherical coordinates.
+     *  @returns double result: result of the quadrature rule */
+    virtual double IntegrateSpherical( double( f )( double my, double phi ) );
 
     /*! @brief Integrates vector valued f(x,y,z) with the quadrature. Each dimension is integrated by itself.
      *  @param : double(f)( double x0, double x1, double x2 ) : density function that depends on a three spatial dimensions.
      *  @param :  len : lenght of vector
      *  @returns double result: result of the quadrature rule (vector valued) */
-    std::vector<double> Integrate( std::vector<double>( f )( double x0, double x1, double x2 ), unsigned len );
+    virtual std::vector<double> Integrate( std::vector<double>( f )( double x0, double x1, double x2 ), unsigned len );
 
     // Quadrature Hub
     /*! @brief Creates a quadrature rule with a given name and a given order.
@@ -60,8 +65,10 @@ class QuadratureBase
     inline unsigned GetOrder() const { return _order; }       /*! @returns unsigned _order:  order of the quadrature */
     inline unsigned GetNq() const { return _nq; }             /*! @returns unsigned _nq:  number of gridpoints of the quadrature */
     inline VectorVector GetPoints() const { return _points; } /*! @returns VectorVector _points:  coordinates of gridpoints of the quadrature */
-    virtual VectorVector GetPointsSphere() const;             /*! @returns VectorVector _pointsSphere:  "---- " in spherical coordinates (my, phi)*/
-    inline Vector GetWeights() const { return _weights; }     /*! @returns Vector _weights:  weights of gridpoints of the quadrature */
+    inline VectorVector GetPointsSphere() const {
+        return _pointsSphere;
+    }                                                     /*! @returns VectorVector _pointsSphere:  "---- " in spherical coordinates (my, phi)*/
+    inline Vector GetWeights() const { return _weights; } /*! @returns Vector _weights:  weights of gridpoints of the quadrature */
     inline VectorVectorU GetConnectivity() const {
         return _connectivity;
     } /*! @returns VectorVectorU _connectivity:  connectivity of gridpoints of the quadrature */
@@ -84,10 +91,9 @@ class QuadratureBase
     unsigned _order;             /*! @brief order of the quadrature */
     unsigned _nq;                /*! @brief number of gridpoints of the quadrature */
     VectorVector _points;        /*! @brief gridpoints of the quadrature */
+    VectorVector _pointsSphere;  /*! @brief (my,phi)gridpoints of the quadrature in spherical cordinates */
     Vector _weights;             /*! @brief weights of the gridpoints of the quadrature */
     VectorVectorU _connectivity; /*! @brief connectivity of the gripoints of the quadrature */
-
-    VectorVector _pointsSphere; /*! @brief (my,phi)gridpoints of the quadrature in spherical cordinates */
 };
 
 #endif    // QUADRATURE_H

code/include/solvers/mnsolver.h

@@ -4,7 +4,7 @@
 #include "solverbase.h"
 
 class EntropyBase;
-class SphericalHarmonics;
+class SphericalBase;
 class OptimizerBase;
 
 class MNSolver : public Solver
@@ -22,11 +22,11 @@ class MNSolver : public Solver
   private:
     // --- Private member variables ---
     unsigned _nTotalEntries;    /*! @brief: Total number of equations in the system */
-    unsigned short _LMaxDegree; /*! @brief: Max Order of Moments */
+    unsigned short _LMaxDegree; /*! @brief: Max Order of Spherical Harmonics  */
 
     // Moment basis
-    SphericalHarmonics* _basis; /*! @brief: Class to compute and store current spherical harmonics basis */
-    VectorVector _moments;      /*! @brief: Moment Vector pre-computed at each quadrature point: dim= _nq x _nTotalEntries */
+    SphericalBase* _basis; /*! @brief: Class to compute and store current spherical harmonics basis */
+    VectorVector _moments; /*! @brief: Moment Vector pre-computed at each quadrature point: dim= _nq x _nTotalEntries */
 
     // Scattering
     Vector _scatterMatDiag; /*! @brief: Diagonal of the scattering matrix (its a diagonal matrix by construction) */
@@ -48,9 +48,6 @@ class MNSolver : public Solver
     // ---- Private Member functions ---
 
     // IO
-    /*! @brief Function that writes NN Training Data in a .csv file */
-    void WriteNNTrainingData( unsigned idx_pseudoTime );
-
     void PrepareVolumeOutput() override;
     void WriteVolumeOutput( unsigned idx_pseudoTime ) override;
 
@@ -74,12 +71,7 @@ class MNSolver : public Solver
     void ComputeScatterMatrix();
 
     // Helper
+    /*! @brief:  Computes the radiative flux from the solution vector of the moment system */
     void ComputeRadFlux();
-    /*! @brief : computes the global index of the moment corresponding to basis function (l,k)
-     *  @param : degree l, it must hold: 0 <= l <=_nq
-     *  @param : order k, it must hold: -l <=k <= l
-     *  @returns : global index
-     */
-    int GlobalIndex( int l, int k ) const;
 };
 #endif    // MNSOLVER_H

code/include/solvers/solverbase.h

@@ -26,12 +26,12 @@ class Solver
 
     // --------- Often used variables of member classes for faster access ----
 
-    unsigned _nEnergies;             /*! @brief number of energy/time steps, number of nodal energy values for CSD */
-    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
-    std::vector<VectorVector> _Q;    /*!  @brief  external source term */
+    unsigned _nEnergies;          /*! @brief number of energy/time steps, number of nodal energy values for CSD */
+    double _dE;                   /*! @brief energy/time step size */
+    Vector _energies;             /*! @brief energy groups used in the simulation [keV] */
+    std::vector<double> _density; /*! @brief patient density, dim(_density) = _nCells */
+    Vector _s;                    /*! @brief stopping power, dim(_s) = _nTimeSteps    */
+    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 */
@@ -40,9 +40,6 @@ class Solver
     QuadratureBase* _quadrature; /*! @brief quadrature to create members below */
     unsigned _nq;                /*! @brief number of quadrature points */
 
-    // VectorVector _quadPoints;    /*!  @brief quadrature points, dim(_quadPoints) = (_nSystem,spatialDim) */
-    // Vector _weights;             /*!  @brief quadrature weights, dim(_weights) = (_NCells) */
-
     // Mesh related members
     unsigned _nCells;                          /*! @brief number of spatial cells */
     std::vector<BOUNDARY_TYPE> _boundaryCells; /*! boundary type for all cells, dim(_boundary) = (_NCells) */
@@ -53,11 +50,11 @@ class Solver
     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 */
+    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
@@ -120,6 +117,10 @@ class Solver
     void PrepareHistoryOutput();
     /*! @brief Prints HistoryOutputFields to logger */
     void PrintHistoryOutput( unsigned iteration );
+    /*! @brief Pre Solver Screen and Logger Output */
+    void DrawPreSolverOutput();
+    /*! @brief Post Solver Screen and Logger Output */
+    void DrawPostSolverOutput();
 
   public:
     /*! @brief Solver constructor

code/include/toolboxes/datagenerator.h

+/*!
+ * \file datagenerator.h
+ * \brief Class to generate data for the neural entropy closure
+ * \author S. Schotthoefer
+ */
+
+#ifndef DATAGENERATOR_H
+#define DATAGENERATOR_H
+
+#include "common/typedef.h"
+#include <vector>
+
+class SphericalBase;
+class QuadratureBase;
+class Config;
+class NewtonOptimizer;
+class EntropyBase;
+
+class nnDataGenerator
+{
+  public:
+    /*! @brief: Class constructor. Generates training data for neural network approaches using
+     *          spherical harmonics and an entropy functional and the quadrature specified by
+     *          the options file.
+     *   @param: setSize: number of elements in training set
+     *           basisSize: length of spherical harmonics basis (maybe redundant)*/
+    nnDataGenerator( Config* settings );
+    ~nnDataGenerator();
+
+    /*! @brief: computes the training data set.
+     *          Realizable set is sampled uniformly.
+     *          Prototype: 1D, u\in[0,100] */
+    void computeTrainingData();
+
+    /*!   @brief: Writes the training data to file
+     *            Filename encryption: [TODO] */
+    void writeTrainingDataToCSV();
+
+  private:
+    Config* _settings; /*! @brief config class for global information */
+
+    VectorVector _uSol;            /*! @brief: vector with moments. Size: (setSize,basisSize)*/
+    VectorVector _alpha;           /*! @brief: vector with Lagrange multipliers. Size: (setSize,basisSize)*/
+    std::vector<double> _hEntropy; /*! @brief: vector with entropy values. Size: (setSize) */
+
+    unsigned long _setSize;
+    unsigned short _LMaxDegree; /*! @brief: Max Order of Spherical Harmonics */
+    unsigned _nTotalEntries;    /*! @brief: Total number of equations in the system */
+
+    QuadratureBase* _quadrature;    /*! @brief quadrature to create members below */
+    unsigned _nq;                   /*! @brief number of quadrature points */
+    VectorVector _quadPoints;       /*!  @brief quadrature points, dim(_quadPoints) = (_nq,spatialDim) */
+    Vector _weights;                /*!  @brief quadrature weights, dim(_weights) = (_nq) */
+    VectorVector _quadPointsSphere; /*!  @brief (my,phi), dim(_quadPoints) = (_nq,2) */
+
+    SphericalBase* _basis; /*! @brief: Class to compute and store current spherical harmonics basis */
+    VectorVector _moments; /*! @brief: Moment Vector pre-computed at each quadrature point: dim= _nq x _nTotalEntries */
+
+    NewtonOptimizer* _optimizer; /*! @brief: Class to solve minimal entropy problem */
+    EntropyBase* _entropy;       /*! @brief: Class to handle entropy functional evaluations */
+
+    // Main methods
+    void SampleSolutionU();        /*! @brief: Samples solution vectors u */
+    void ComputeEntropyH_dual();   /*! @brief: Compute the entropy functional at (u,alpha) in dual formulation */
+    void ComputeEntropyH_primal(); /*! @brief:  Compute the entropy functional at (u,alpha) in primal formulation */
+
+    // IO routines
+    void PrintTrainingData(); /*! @brief : Print computed training data to csv file and screen */
+    void PrintLoadScreen();   /*! @brief: Print screen IO*/
+    // Helper functions
+    void ComputeMoments();        /*! @brief: Pre-Compute Moments at all quadrature points. */
+    void CheckRealizability();    // Debugging helper
+
+    inline VectorVector GetuSol() { return _uSol; }                /*! @brief: Get the computed solution vector uSol */
+    inline VectorVector GetAlpha() { return _alpha; }              /*! @brief: Get the computed vector alpha */
+    inline std::vector<double> GethEntropy() { return _hEntropy; } /*! @brief: Get the computed entropy value h */
+};
+#endif    // DATAGENERATOR_H

code/include/toolboxes/errormessages.h

@@ -7,8 +7,6 @@
 #ifndef CRTSNERROR_H
 #define CRTSNERROR_H
 
-#include <iostream>
-
 #include "spdlog/spdlog.h"
 
 class ErrorMessages

code/include/toolboxes/sphericalbase.h

+/*!
+ * @file sphericalbase.h
+ * @brief Base Class to handle basis classes on the unit sphere
+ * @author S. Schotthöfer
+ *
+ */
+
+#ifndef SPHERICALBASE_H
+#define SPHERICALBASE_H
+
+#include "common/typedef.h"
+class Config;
+
+class SphericalBase
+{
+  public:
+    SphericalBase() {}
+    virtual ~SphericalBase() {}
+
+    /*! @brief: Create a set of basis functions on the unit sphere defined in settings
+     *  @param: Pointer to the config file
+     *  @returns: Pointer to the createt basis class */
+    static SphericalBase* Create( Config* settings );
+
+    /*! @brief  : Computes all N basis functions at point (my, phi)
+     *  @param  : my = cos(theta) - spherical coordinate, -1 <= x <= 1
+     *  @param  : phi - spherical coordinate, 0 <= phi <= 2*pi
+     *  @return : vector of basis functions at point (my, phi) with size N
+     */
+    virtual Vector ComputeSphericalBasis( double my, double phi ) = 0;
+
+    /*! @brief  : Computes all basis functions at point (x, y, z) on the unit sphere
+     *  @param  : x,y,z = coordinates on unit sphere
+     *  @return : vector of basis functions at point (x,y,z) with size N
+     */
+    virtual Vector ComputeSphericalBasis( double x, double y, double z ) = 0;
+
+    /*! @brief : Return size of complete Basisvector */
+    virtual unsigned GetBasisSize() = 0;
+
+    /*! @brief: Return number of basis functions with degree equals to currDegree
+     *  @param: currDegree must be smaller equals _LMaxDegree     */
+    virtual unsigned GetCurrDegreeSize( unsigned currDegree ) = 0;
+
+    /*! @brief: Computes global index of basis vector depending on order k and degree l
+     *  @param: l_degree = degree of polynomials l = 0,1,2,3,...
+     *  @param: k_order = order of element of degree l. !ATTENTION. Requirements are different for monomials and harmonics! */
+    virtual unsigned GetGlobalIndexBasis( int l_degree, int k_order ) = 0;
+
+  protected:
+    /*! @brief: maximal (polynomial) degree of the spherical basis (this is "L" in the comments)*/
+    unsigned _LMaxDegree;
+};
+
+#endif    // SPHERICALBASE_H