Merge branch 'MN_solver' into 'master'

Mn solver

Closes #58 and #57

See merge request !19

code/include/entropies/entropybase.h

+#ifndef ENTROPYBASE_H
+#define ENTROPYBASE_H
+
+// Foward declaration
+class Config;
+
+class EntropyBase
+{
+  public:
+    inline EntropyBase() {}
+
+    virtual inline ~EntropyBase() {}
+
+    static EntropyBase* Create( Config* settings );
+
+    /*! @brief: computes the entropy functional
+     *  @param: z = point where the functional should be evaluated.
+     *          z must be in domain of the functional
+     *  @returns: value of entropy functional at z */
+    virtual double Entropy( double z ) = 0;
+
+    /*! @brief: computes eta'(z).z must be in domain of the functional.
+     *  @param: z = point where the derivative should be evaluated. */
+    virtual double EntropyPrime( double z ) = 0;
+
+    /*! @brief: computes the dual of the  entropy functional
+     *  @param: z = point where the dual of the functional should be evaluated.
+     *          z must be in domain of the duaö functional
+     *  @returns: value of entropy functional at z */
+    virtual double EntropyDual( double y ) = 0;
+
+    // /*! @brief: computes m*eta_*'(alpha*m). alpha*m must be in domain of the functional
+    //  *          alpha, m and grad must be of same size.
+    //  *  @param: alpha = point where the derivative should be evaluated.
+    //  *  @param: m = moment vector
+    //  *  @param: grad: vector in which the resulting gradient is stored */
+    // virtual void EntropyPrimeDual( Vector& alpha, Vector& m, Vector& grad ) = 0;
+
+    /*! @brief: computes eta_*'(y).
+     *  @param: y = point where the derivative should be evaluated.
+     *  @returns: value of the derivative of the  entropy functional at y */
+    virtual double EntropyPrimeDual( double z ) = 0;
+
+    /*! @brief: computes the hessian of the dual entropy functional
+     *  @param: y = point where the hessian should be evaluated;
+     *  @returns: value of the hessian at y */
+    virtual double EntropyHessianDual( double y ) = 0;
+
+    /*! @brief: checks, if value is in domain of entropy */
+    virtual bool CheckDomain( double z ) = 0;
+};
+
+#endif    // ENTROPYBASE_H

code/include/entropies/maxwellboltzmannentropy.h

+#ifndef MAXWELLBOLTZMANNENTROPY_H
+#define MAXWELLBOLTZMANNENTROPY_H
+
+#include "entropybase.h"
+#include <cmath>
+
+class MaxwellBoltzmannEntropy : public EntropyBase
+{
+  public:
+    inline MaxwellBoltzmannEntropy() : EntropyBase() {}
+
+    inline ~MaxwellBoltzmannEntropy() {}
+
+    inline double Entropy( double z ) override { return z * log( z ) - z; }
+
+    inline double EntropyPrime( double z ) override { return log( z ); }
+
+    // inline void EntropyPrime( Vector& alpha, Vector& m, Vector& grad ) override { grad = m * log( dot( alpha, m ) ); }
+
+    inline double EntropyDual( double y ) override { return exp( y ); }
+
+    // inline void EntropyPrimeDual( Vector& alpha, Vector& m, Vector& grad ) override { grad = m * exp( dot( alpha, m ) ); }
+
+    inline double EntropyPrimeDual( double y ) override { return exp( y ); }
+
+    inline double EntropyHessianDual( double y ) override { return exp( y ); }
+
+    inline bool CheckDomain( double z ) override { return ( z >= 0 && std::isnan( z ) ); }
+};
+
+#endif    // MAXWELLBOLTZMANNENTROPY_H

code/include/entropies/quadraticentropy.h

+#ifndef QUADRATICENTROPY_H
+#define QUADRATICENTROPY_H
+
+#include "entropybase.h"
+#include <cmath>
+
+class QuadraticEntropy : public EntropyBase
+{
+  public:
+    inline QuadraticEntropy() : EntropyBase() {}
+
+    inline ~QuadraticEntropy() {}
+
+    inline double Entropy( double z ) override { return 0.5 * z * z; }
+
+    inline double EntropyPrime( double z ) override { return z; }
+
+    inline double EntropyDual( double y ) override { return 0.5 * y * y; }
+
+    // inline void EntropyPrimeDual( Vector& alpha, Vector& m, Vector& grad ) override { grad = m * dot( alpha, m ); }
+
+    inline double EntropyPrimeDual( double y ) override { return y; }
+
+    inline double EntropyHessianDual( double y ) override { return 1.0; }
+
+    inline bool CheckDomain( double z ) override { return std::isnan( z ); }
+};
+
+#endif    // QUADRATICENTROPY_H

code/include/fluxes/laxfriedrichsflux.h

@@ -12,7 +12,7 @@ class LaxFriedrichsFlux : public NumericalFlux
      * @brief LaxFriedrichsFlux
      * @param settings
      */
-    LaxFriedrichsFlux( Config* settings );
+    LaxFriedrichsFlux();
 
     /**
      * @brief Flux computes flux on edge for fixed ordinate at a given edge

code/include/fluxes/numericalflux.h

 #ifndef NUMERICALFLUX_H
 #define NUMERICALFLUX_H
 
-#include "settings/config.h"
 #include "settings/typedef.h"
 
 class NumericalFlux
 {
   public:
-    NumericalFlux( Config* settings );
-    static NumericalFlux* Create( Config* settings );
+    NumericalFlux();
+    static NumericalFlux* Create();
+
     /**
      * @brief Flux computes flux on edge for fixed ordinate at a given edge
      * @param Omega fixed ordinate for flux computation

code/include/fluxes/upwindflux.h

@@ -10,7 +10,7 @@ class UpwindFlux : public NumericalFlux
      * @brief UpwindFlux
      * @param settings
      */
-    UpwindFlux( Config* settings );
+    UpwindFlux();
 
     /**
      * @brief Flux computes flux on edge for fixed ordinate at a given edge

code/include/io.h

 #ifndef IO_H
 #define IO_H
 
-#include <chrono>
-#include <filesystem>
-#include <iostream>
+#include "settings/typedef.h"
 #include <string>
+#include <vector>
 
-#include <mpi.h>
-#include <omp.h>
-#include <vtkCellArray.h>
-#include <vtkCellData.h>
-#include <vtkCellDataToPointData.h>
-#include <vtkDoubleArray.h>
-#include <vtkPointData.h>
-#include <vtkPointDataToCellData.h>
-#include <vtkQuad.h>
-#include <vtkSmartPointer.h>
-#include <vtkTriangle.h>
-#include <vtkUnstructuredGrid.h>
-#include <vtkUnstructuredGridWriter.h>
-
-#include <Python.h>
-#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION
-#include <numpy/arrayobject.h>
-
-#include "mesh.h"
-#include "settings/config.h"
-
-using vtkPointsSP                 = vtkSmartPointer<vtkPoints>;
-using vtkUnstructuredGridSP       = vtkSmartPointer<vtkUnstructuredGrid>;
-using vtkTriangleSP               = vtkSmartPointer<vtkTriangle>;
-using vtkCellArraySP              = vtkSmartPointer<vtkCellArray>;
-using vtkDoubleArraySP            = vtkSmartPointer<vtkDoubleArray>;
-using vtkUnstructuredGridWriterSP = vtkSmartPointer<vtkUnstructuredGridWriter>;
-using vtkCellDataToPointDataSP    = vtkSmartPointer<vtkCellDataToPointData>;
-using vtkPointDataToCellDataSP    = vtkSmartPointer<vtkPointDataToCellData>;
+// Forward Declarations
+class Config;
+class Mesh;
+// class Matrix;
 
 void ExportVTK( const std::string fileName,
                 const std::vector<std::vector<std::vector<double>>>& results,

code/include/kernels/scatteringkernelbase.h

 #ifndef SCATTERINGKERNELBASE_CPP
 #define SCATTERINGKERNELBASE_CPP
 
-#include "quadratures/quadraturebase.h"
 #include "settings/globalconstants.h"
 #include "settings/typedef.h"
 
+// Forward Declaration
+class QuadratureBase;
+
 class ScatteringKernel
 {
   private:
@@ -15,7 +17,7 @@ class ScatteringKernel
 
   public:
     ScatteringKernel( QuadratureBase* quad );
-    ~ScatteringKernel();
+    virtual ~ScatteringKernel();
 
     virtual Matrix GetScatteringKernel() = 0;
 

code/include/mesh.h

 #ifndef MESH_H
 #define MESH_H
 
-#include <algorithm>
-#include <mpi.h>
-#include <omp.h>
-#include <vector>
-
 #include "blaze/math/CompressedMatrix.h"
-#include "metis.h"
-#include "parmetis.h"
-
 #include "settings/globalconstants.h"
 #include "settings/typedef.h"
-#include "toolboxes/errormessages.h"
 
-#include "reconstructor.h"
+#include <vector>
 
 class Mesh
 {
@@ -132,7 +123,6 @@ class Mesh
 
     void ReconstructSlopesS( unsigned nq, VectorVector& psiDerX, VectorVector& psiDerY, const VectorVector& psi ) const;
     void ReconstructSlopesU( unsigned nq, VectorVector& psiDerX, VectorVector& psiDerY, const VectorVector& psi ) const;
-
 };
 
 #endif    // MESH_H

code/include/optimizers/newtonoptimizer.h

+#ifndef NEWTONOPTIMIZER_H
+#define NEWTONOPTIMIZER_H
+
+#include "optimizerbase.h"
+
+class QuadratureBase;
+
+class NewtonOptimizer : public OptimizerBase
+{
+  public:
+    NewtonOptimizer( Config* settings );
+
+    inline ~NewtonOptimizer() {}
+
+    void Solve( Vector& lambda, Vector& u, VectorVector& moments ) override;
+
+  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 );
+
+    /*! @brief: Computes hessian of objective function and stores it in hessian
+                grad = <mXm*eta*'(alpha*m)> */
+    void ComputeHessian( Vector& alpha, VectorVector& moments, Matrix& hessian );
+
+    double ComputeObjFunc( Vector& alpha, Vector& sol, VectorVector& moments );
+
+    QuadratureBase* _quadrature; /*! @brief: used quadrature */    // THis is memory doubling! Try to use a pointer.
+    unsigned _nq;                                                  /*! @brief: number of quadrature points */
+    Vector _weights;                                               /*!  @brief quadrature weights, dim(_weights) = (_nq) */
+    VectorVector _quadPointsSphere;                                /*!  @brief (my,phi), dim(_quadPoints) = (_nq,2) */
+
+    double _epsilon;                 /*! @brief: Termination criterion for newton optimizer */
+    unsigned short _maxIterations;   /*! @brief: Max iterations of the newton solver */
+    double _alpha;                   /*! @brief: Newton Step Size */
+    unsigned short _maxLineSearches; /*! @brief: Max amount of line searches for Newton Algo */
+};
+
+#endif    // NEWTONOPTIMIZER_H

code/include/optimizers/optimizerbase.h

+#ifndef OPTIMIZERBASE_H
+#define OPTIMIZERBASE_H
+
+#include "entropies/entropybase.h"
+#include "settings/typedef.h"
+
+// Foward declaration
+class Config;
+
+class OptimizerBase
+{
+  public:
+    OptimizerBase( Config* settings );
+
+    virtual inline ~OptimizerBase() { delete _entropy; }
+
+    static OptimizerBase* Create( Config* settings );
+
+    /*! @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 ) = 0;
+
+  protected:
+    EntropyBase* _entropy; /*! @brief: Class to handle entropy functional evaluations */
+    Config* _settings;
+};
+
+#endif    // OPTIMIZERBASE_H

code/include/physics.h

 #ifndef PHYSICS_H
 #define PHYSICS_H
 
-// include Matrix, Vector definitions
-#include "math.h"
-#include "settings/config.h"
 #include "settings/typedef.h"
-#include "spline.h"
-
-#include <fstream>
-#include <list>
-#include <tuple>
 
 class Physics
 {

code/include/problems/problembase.h

 #ifndef PROBLEMBASE_H
 #define PROBLEMBASE_H
 
-#include "mesh.h"
-#include "physics.h"
-#include "settings/config.h"
 #include "settings/typedef.h"
 
+// Forward Declaration
+class Config;
+class Physics;
+class Mesh;
+
 class ProblemBase
 {
 

code/include/quadratures/qgausslegendretensorized.h

@@ -5,6 +5,7 @@
 
 class QGaussLegendreTensorized : 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 );
@@ -14,10 +15,12 @@ class QGaussLegendreTensorized : public QuadratureBase
     QGaussLegendreTensorized( unsigned order );
     virtual ~QGaussLegendreTensorized() {}
 
-    inline void SetName() override { _name = "Tensorized Gauss-Legendre quadrature."; }
-    inline void SetNq() override { _nq = pow( GetOrder(), 2 ); }
+    inline void SetName() override { _name = "Tensorized Gauss-Legendre quadrature"; }
+    inline void SetNq() override { _nq = 2 * pow( GetOrder(), 2 ); }
     void SetPointsAndWeights() override;
     void SetConnectivity() override;
+
+    inline VectorVector GetPointsSphere() const override { return _pointsSphere; }
 };
 
 #endif    // QGAUSSLEGENDRETENSORIZED_H

code/include/quadratures/qmontecarlo.h

@@ -7,12 +7,13 @@ class QMonteCarlo : public QuadratureBase
 {
   public:
     QMonteCarlo( unsigned order );
-    virtual ~QMonteCarlo() {}
+    inline ~QMonteCarlo() {}
 
     inline void SetName() override { _name = "Monte Carlo Quadrature."; }
     inline void SetNq() override { _nq = GetOrder(); }
     void SetPointsAndWeights() override;
     void SetConnectivity() override;
+    VectorVector GetPointsSphere() const override;
 };
 
 #endif    // QMONTECARLO_H

code/include/quadratures/quadraturebase.h

@@ -3,7 +3,6 @@
 
 #include "settings/globalconstants.h"
 #include "settings/typedef.h"
-#include "toolboxes/errormessages.h"
 
 class QuadratureBase
 {
@@ -43,6 +42,7 @@ 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 VectorVectorU GetConnectivity() const {
         return _connectivity;
@@ -61,12 +61,15 @@ class QuadratureBase
     virtual void SetPointsAndWeights() = 0;
 
     // Member variables
+    // TODO Config* _settings;           /*! @brief pointer to settings class that manages the solver */
     std::string _name;           /*! @brief name of the quadrature */
     unsigned _order;             /*! @brief order of the quadrature */
     unsigned _nq;                /*! @brief number of gridpoints of the quadrature */
     VectorVector _points;        /*! @brief gridpoints of the quadrature */
     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/settings/config.h

@@ -18,7 +18,9 @@
 #include "spdlog/spdlog.h"
 
 #include "globalconstants.h"
-#include "optionstructure.h"
+
+// Forward declaration
+class OptionBase;
 
 /*!
  * @class Config
@@ -56,6 +58,7 @@ class Config
     double _tEnd;                    /*!< @brief Final Time for Simulation */
     PROBLEM_NAME _problemName;       /*!< @brief Name of predefined Problem   */
     SOLVER_NAME _solverName;         /*!< @brief Name of the used Solver */
+    ENTROPY_NAME _entropyName;       /*!< @brief Name of the used Entropy Functional */
     unsigned short _maxMomentDegree; /*!< @brief Maximal Order of Moments for PN and MN Solver */
     unsigned short _reconsOrder;     /*!< @brief Spatial Order of Accuracy for Solver */
 
@@ -78,6 +81,14 @@ class Config
     // Scattering Kernel
     KERNEL_NAME _kernelName; /*!< @brief Scattering Kernel Name*/
 
+    // 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 */
+    double _newtonStepSize;               /*!< @brief Stepsize factor for newton optimizer */
+    unsigned short _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 */
+
     // --- Parsing Functionality and Initializing of Options ---
     /*!
      * @brief Set default values for all options not yet set.
@@ -209,15 +220,25 @@ class Config
     unsigned GetNCells() { return _nCells; }
 
     // Solver Structure
+    unsigned short inline GetMaxMomentDegree() const { return _maxMomentDegree; }
     double inline GetCFL() const { return _CFL; }
     double inline GetTEnd() const { return _tEnd; }
     PROBLEM_NAME inline GetProblemName() const { return _problemName; }
     SOLVER_NAME inline GetSolverName() const { return _solverName; }
+    ENTROPY_NAME inline GetEntropyName() const { return _entropyName; }
     bool inline GetCleanFluxMat() const { return _cleanFluxMat; }
     unsigned GetReconsOrder() { return _reconsOrder; }
     bool inline IsCSD() const { return _csd; }
     unsigned inline GetMaxMomentDegree() { return _maxMomentDegree; }
 
+    //  Optimizer
+    OPTIMIZER_NAME inline GetOptimizerName() const { return _entropyOptimizerName; }
+    double inline GetNewtonOptimizerEpsilon() const { return _optimizerEpsilon; }
+    unsigned inline GetNewtonIter() const { return _newtonIter; }
+    double inline GetNewtonStepSize() const { return _newtonStepSize; }
+    unsigned inline GetMaxLineSearches() const { return _newtonLineSearchIter; }
+    bool inline GetNewtonFastMode() const { return _newtonFastMode; }
+
     // Boundary Conditions
     BOUNDARY_TYPE GetBoundaryType( std::string nameMarker ) const; /*! @brief Get Boundary Type of given marker */
 

code/include/settings/globalconstants.h

@@ -53,8 +53,18 @@ enum KERNEL_NAME { KERNEL_Isotropic, KERNEL_Isotropic1D };
 inline std::map<std::string, KERNEL_NAME> Kernel_Map{ { "ISOTROPIC", KERNEL_Isotropic }, { "ISOTROPIC_1D", KERNEL_Isotropic1D } };
 
 // Solver name
-enum SOLVER_NAME { SN_SOLVER, PN_SOLVER };
+enum SOLVER_NAME { SN_SOLVER, PN_SOLVER, MN_SOLVER };
 
-inline std::map<std::string, SOLVER_NAME> Solver_Map{ { "SN_SOLVER", SN_SOLVER }, { "PN_SOLVER", PN_SOLVER } };
+inline std::map<std::string, SOLVER_NAME> Solver_Map{ { "SN_SOLVER", SN_SOLVER }, { "PN_SOLVER", PN_SOLVER }, { "MN_SOLVER", MN_SOLVER } };
 
+// Entropy functional
+enum ENTROPY_NAME { QUADRATIC, MAXWELL_BOLZMANN, BOSE_EINSTEIN, FERMI_DIRAC };
+
+inline std::map<std::string, ENTROPY_NAME> Entropy_Map{
+    { "QUADRATIC", QUADRATIC }, { "MAXWELL_BOLZMANN", MAXWELL_BOLZMANN }, { "BOSE_EINSTEIN", BOSE_EINSTEIN }, { "FERMI_DIRAC", FERMI_DIRAC } };
+
+// Otpimizer
+enum OPTIMIZER_NAME { NEWTON, ML };
+
+inline std::map<std::string, OPTIMIZER_NAME> Optimizer_Map{ { "NEWTON", NEWTON }, { "ML", ML } };
 #endif    // GLOBAL_CONSTANTS_H

code/include/settings/optionstructure.h

@@ -25,35 +25,17 @@ class OptionBase
     OptionBase(){};
     virtual ~OptionBase() = 0;
 
-    virtual std::string SetValue( std::vector<std::string> value ) {
-        this->_value = value;
-        return "";
-    }
-    std::vector<std::string> GetValue() { return _value; }
+    virtual std::string SetValue( std::vector<std::string> value );
+
+    std::vector<std::string> GetValue();
+
     virtual void SetDefault() = 0;
 
-    std::string optionCheckMultipleValues( std::vector<std::string>& option_value, std::string type_id, std::string option_name ) {
-        if( option_value.size() != 1 ) {
-            std::string newString;
-            newString.append( option_name );
-            newString.append( ": multiple values for type " );
-            newString.append( type_id );
-            return newString;
-        }
-        return "";
-    }
+    std::string optionCheckMultipleValues( std::vector<std::string>& option_value, std::string type_id, std::string option_name );
 
-    std::string badValue( std::vector<std::string>& option_value, std::string type_id, std::string option_name ) {
-        std::string newString;
-        newString.append( option_name );
-        newString.append( ": improper option value for type " );
-        newString.append( type_id );
-        return newString;
-    }
+    std::string badValue( std::vector<std::string>& option_value, std::string type_id, std::string option_name );
 };
 
-inline OptionBase::~OptionBase() {}
-
 template <class Tenum> class OptionEnum : public OptionBase
 {
 
@@ -71,6 +53,7 @@ template <class Tenum> class OptionEnum : public OptionBase
     }
 
     ~OptionEnum() override{};
+
     std::string SetValue( std::vector<std::string> option_value ) override {
         OptionBase::SetValue( option_value );
         // Check if there is more than one string
@@ -78,7 +61,6 @@ template <class Tenum> class OptionEnum : public OptionBase
         if( out.compare( "" ) != 0 ) {
             return out;
         }
-
         // Check to see if the enum value is in the map
         if( this->_map.find( option_value[0] ) == _map.end() ) {
             std::string str;
@@ -104,28 +86,13 @@ class OptionDouble : public OptionBase
     std::string _name;    // identifier for the option
 
   public:
-    OptionDouble( std::string option_field_name, double& option_field, double default_value ) : _field( option_field ) {
-        this->_def  = default_value;
-        this->_name = option_field_name;
-    }
+    OptionDouble( std::string option_field_name, double& option_field, double default_value );
 
     ~OptionDouble() override{};
-    std::string SetValue( std::vector<std::string> option_value ) override {
-        OptionBase::SetValue( option_value );
-        // check if there is more than one value
-        std::string out = optionCheckMultipleValues( option_value, "double", this->_name );
-        if( out.compare( "" ) != 0 ) {
-            return out;
-        }
-        std::istringstream is( option_value[0] );
-        double val;
-        if( is >> val ) {
-            this->_field = val;
-            return "";
-        }
-        return badValue( option_value, "double", this->_name );
-    }
-    void SetDefault() override { this->_field = this->_def; }
+
+    std::string SetValue( std::vector<std::string> option_value ) override;
+
+    void SetDefault() override;
 };
 
 class OptionString : public OptionBase
@@ -135,24 +102,13 @@ class OptionString : public OptionBase
     std::string _name;      // identifier for the option
 
   public: