startet sn-solver haf quadrature points option

code/include/common/config.h

@@ -61,6 +61,7 @@ class Config
     /*!< @brief If true, very low entries (10^-10 or smaller) of the flux matrices will be set to zero,
      * to improve floating point accuracy */
     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 */
     std::string _hydrogenFile; /*!< @brief Name of hydrogen cross section file */
@@ -226,7 +227,8 @@ class Config
     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 GetUseAllGaussPts() const { return _allGaussPts; }
+    unsigned inline GetReconsOrder() { return _reconsOrder; }
     bool inline IsCSD() const { return _csd; }
     unsigned inline GetMaxMomentDegree() { return _maxMomentDegree; }
 

code/include/quadratures/qgausslegendretensorized.h

@@ -12,7 +12,7 @@ class QGaussLegendreTensorized : public QuadratureBase
     bool CheckOrder();
 
   public:
-    QGaussLegendreTensorized( unsigned order );
+    QGaussLegendreTensorized( Config* settings );
     virtual ~QGaussLegendreTensorized() {}
 
     inline void SetName() override { _name = "Tensorized Gauss-Legendre quadrature"; }

code/include/quadratures/qmontecarlo.h

@@ -6,7 +6,7 @@
 class QMonteCarlo : public QuadratureBase
 {
   public:
-    QMonteCarlo( unsigned order );
+    QMonteCarlo( Config* settings );
     inline ~QMonteCarlo() {}
 
     inline void SetName() override { _name = "Monte Carlo Quadrature."; }

code/include/quadratures/quadraturebase.h

@@ -4,10 +4,12 @@
 #include "common/globalconstants.h"
 #include "common/typedef.h"
 
+class Config;
+
 class QuadratureBase
 {
   public:
-    QuadratureBase( unsigned order );
+    QuadratureBase( Config* settings );
     virtual ~QuadratureBase() {}
 
     // Aux functions
@@ -35,7 +37,7 @@ class QuadratureBase
      *  @param: std::string name: Name of the quadrature rule
      *  @param: unsigned order: Order of the quadrature rule
      *  @returns Quadrature* quadrature: returns pointer to instance of the given derived quadrature class */
-    static QuadratureBase* CreateQuadrature( QUAD_NAME name, unsigned order );
+    static QuadratureBase* CreateQuadrature( Config* settings );
 
     // Getter
     inline std::string GetName() const { return _name; }      /*! @returns std::string _name:  name of the quadrature */
@@ -61,7 +63,7 @@ class QuadratureBase
     virtual void SetPointsAndWeights() = 0;
 
     // Member variables
-    // TODO Config* _settings;           /*! @brief pointer to settings class that manages the solver */
+    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 */

code/input/exampleMN.cfg

@@ -17,6 +17,7 @@ LOG_DIR = ../result/logs
 MESH_FILE = linesource.su2
 %MESH_FILE = linesource_debug.su2
 %
+PROBLEM2 = Lines
 PROBLEM = LINESOURCE
 %
 % ---- Solver specifications ----

code/src/common/config.cpp

@@ -214,7 +214,7 @@ void Config::SetConfigOptions() {
     AddEnumOption( "PROBLEM", _problemName, Problem_Map, PROBLEM_ElectronRT );
     /*! @brief Solver \n DESCRIPTION: Solver used for problem \n DEFAULT SN_SOLVER @ingroup Config. */
     AddEnumOption( "SOLVER", _solverName, Solver_Map, SN_SOLVER );
-    /*!\brief RECONS_ORDER \n DESCRIPTION: Reconstruction order for solver \n DEFAULT 1 \ingroup Config.*/
+    /*! @brief RECONS_ORDER \n DESCRIPTION: Reconstruction order for solver \n DEFAULT 1 \ingroup Config.*/
     AddUnsignedShortOption( "RECONS_ORDER", _reconsOrder, 1 );
     /*! @brief CleanFluxMatrices \n DESCRIPTION:  If true, very low entries (10^-10 or smaller) of the flux matrices will be set to zero,
      * to improve floating point accuracy \n DEFAULT false \ingroup Config */
@@ -228,6 +228,8 @@ void Config::SetConfigOptions() {
     /*! @brief OxygenFile \n DESCRIPTION: If the continuous slowing down approximation is used, this referes to the cross section file for oxygen.
      * . \n DEFAULT "o.dat" \ingroup Config */
     AddStringOption( "OXYGEN_FILE", _oxygenFile, string( "ENDL_O.txt" ) );
+    /*! @brief SN_ALL_GAUSS_PTS \n DESCRIPTION: If true, the SN Solver uses all Gauss Quadrature Points for 2d. \n DEFAULT false \ingroup Config */
+    AddBoolOption( "SN_ALL_GAUSS_PTS", _allGaussPts, false );
 
     // Entropy related options
     /*! @brief Entropy Functional \n DESCRIPTION: Entropy functional used for the MN_Solver \n DEFAULT QUADRTATIC @ingroup Config. */

code/src/quadratures/qgausslegendretensorized.cpp

 #include "quadratures/qgausslegendretensorized.h"
 #include "toolboxes/errormessages.h"
 
-QGaussLegendreTensorized::QGaussLegendreTensorized( unsigned order ) : QuadratureBase( order ) {
+QGaussLegendreTensorized::QGaussLegendreTensorized( Config* settings ) : QuadratureBase( settings ) {
     SetName();
     CheckOrder();
     SetNq();
@@ -46,8 +46,13 @@ void QGaussLegendreTensorized::SetPointsAndWeights() {
         phi[i] = ( i + 0.5 ) * M_PI / _order;
     }
 
-    // unsigned range = std::floor( _order / 2.0 );    // comment (steffen): why do we only need half of the points:
-    //=> In 2D we would count everything twice. (not wrong with scaling
+    unsigned range             = _order;    // By default, use all quad points
+    double normalizationFactor = 1.0;
+    if( _settings->GetSolverName() == SN_SOLVER && _settings->GetUseAllGaussPts() == false ) {
+        range = std::floor( _order / 2.0 );    // comment (steffen): why do we only need half of the points:
+        //=> In 2D we would count everything twice. (not wrong with scaling)
+        normalizationFactor = 2.0;
+    }
 
     // resize points and weights
     _points.resize( _nq );
@@ -62,7 +67,7 @@ void QGaussLegendreTensorized::SetPointsAndWeights() {
     _weights.resize( _nq );
 
     // transform tensorized (x,y,z)-grid to spherical grid points
-    for( unsigned j = 0; j < _order; ++j ) {
+    for( unsigned j = 0; j < range; ++j ) {
         for( unsigned i = 0; i < 2 * _order; ++i ) {
             _points[j * ( 2 * _order ) + i][0] = sqrt( 1 - nodes1D[j] * nodes1D[j] ) * std::cos( phi[i] );
             _points[j * ( 2 * _order ) + i][1] = sqrt( 1 - nodes1D[j] * nodes1D[j] ) * std::sin( phi[i] );
@@ -71,7 +76,7 @@ void QGaussLegendreTensorized::SetPointsAndWeights() {
             _pointsSphere[j * ( 2 * _order ) + i][0] = nodes1D[j];    // my
             _pointsSphere[j * ( 2 * _order ) + i][1] = phi[i];        // phi
 
-            _weights[j * ( 2 * _order ) + i] = M_PI / _order * weights1D[j];
+            _weights[j * ( 2 * _order ) + i] = normalizationFactor * M_PI / _order * weights1D[j];
         }
     }
 }

code/src/quadratures/qmontecarlo.cpp

 #include "quadratures/qmontecarlo.h"
 
-QMonteCarlo::QMonteCarlo( unsigned order ) : QuadratureBase( order ) {
+QMonteCarlo::QMonteCarlo( Config* settings ) : QuadratureBase( settings ) {
     SetName();
     SetNq();
     SetPointsAndWeights();

code/src/quadratures/quadraturebase.cpp

 #include "quadratures/quadraturebase.h"
+#include "common/config.h"
 #include "quadratures/qgausslegendre1D.h"
 #include "quadratures/qgausslegendretensorized.h"
 #include "quadratures/qldfesa.h"
@@ -7,18 +8,22 @@
 #include "quadratures/qmontecarlo.h"
 #include "toolboxes/errormessages.h"
 
-QuadratureBase::QuadratureBase( unsigned order ) : _order( order ) {}
+QuadratureBase::QuadratureBase( Config* settings ) {
+    _settings = settings;
+    _order    = settings->GetQuadOrder();
+}
 
-QuadratureBase* QuadratureBase::CreateQuadrature( QUAD_NAME name, unsigned order ) {
+QuadratureBase* QuadratureBase::CreateQuadrature( Config* settings ) {
+    QUAD_NAME name = settings->GetQuadName();
 
     switch( name ) {
-        case QUAD_MonteCarlo: return new QMonteCarlo( order );
-        case QUAD_GaussLegendreTensorized: return new QGaussLegendreTensorized( order );
-        case QUAD_GaussLegendre1D: return new QGaussLegendre1D( order );
-        case QUAD_LevelSymmetric: return new QLevelSymmetric( order );
-        case QUAD_LDFESA: return new QLDFESA( order );
-        case QUAD_Lebedev: return new QLebedev( order );
-        default: return new QMonteCarlo( order );    // Use MonteCarlo as dummy
+        case QUAD_MonteCarlo: return new QMonteCarlo( settings );
+        case QUAD_GaussLegendreTensorized: return new QGaussLegendreTensorized( settings );
+        case QUAD_GaussLegendre1D: return new QGaussLegendre1D( settings );
+        case QUAD_LevelSymmetric: return new QLevelSymmetric( settings );
+        case QUAD_LDFESA: return new QLDFESA( settings );
+        case QUAD_Lebedev: return new QLebedev( settings );
+        default: return new QMonteCarlo( settings );    // Use MonteCarlo as dummy
     }
 }