added 2D Data Generator and example cfg file

Former-commit-id: 03bec8b6

code/include/toolboxes/datagenerator2D.h

+/*!
+ * \file datagenerator3D.h
+ * \brief Class to generate data for the neural entropy closure in 3D spatial dimensions
+ * \author S. Schotthoefer
+ */
+
+#ifndef DATAGENERATOR2D_H
+#define DATAGENERATOR2D_H
+
+#include "toolboxes/datageneratorbase.h"
+
+class DataGenerator2D : public DataGeneratorBase
+{
+  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 settings config class with global information*/
+    DataGenerator2D( Config* settings );
+    ~DataGenerator2D();
+
+  private:
+    // Main methods
+    void SampleSolutionU() override; /*!< @brief Samples solution vectors u */
+
+    // Helper functions
+    void ComputeMoments() override;        /*!< @brief Pre-Compute Moments at all quadrature points. */
+    void ComputeSetSize() override;        /*!< @brief Computes the size of the training set, depending on the chosen settings.*/
+    void CheckRealizability() override;    // Debugging helper
+};
+#endif    // DATAGENERATOR2D_H

code/include/toolboxes/sphericalharmonics.h

@@ -49,17 +49,17 @@ class SphericalHarmonics : public SphericalBase
     std::vector<double> GetAssLegendrePoly( const double my );
 
     /*! @brief Returns length of the basis, i.e. number of elements of the basis */
-    unsigned GetBasisSize() override;
+    unsigned GetBasisSize() override final;
 
     /*! @brief Return number of basis functions with degree equals to currDegree
      *  @param currDegreeL = degree of polynomials that are counted   */
-    unsigned GetCurrDegreeSize( unsigned currDegree ) override;
+    unsigned GetCurrDegreeSize( unsigned currDegree ) override final;
 
     /*! @brief  helper function to get the global index for given k and l of
      *          the basis function Y_k^l.
      *  @param  l_degree - current degree of basis function, 0 <= l <= L
      *  @param  k_order  - current order of basis function,  -l <= k <= l      */
-    unsigned GetGlobalIndexBasis( int l_degree, int k_order ) override;
+    unsigned GetGlobalIndexBasis( int l_degree, int k_order ) override final;
 
   private:
     /*! @brief maximal degree of the spherical harmonics basis (this is "L" in the comments)*/

code/input/DataGenerator2D.cfg

+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%  Data Generator for                  % 
+%   Neural Entropy Closure             %
+%  Author <Steffen Schotthöfer>		   %
+%  Date   17.12.2020                   %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% ---- Global settings ----
+DATA_GENERATOR_MODE = YES
+TRAINING_SET_SIZE = 20
+MAX_VALUE_FIRST_MOMENT = 5
+REALIZABLE_SET_EPSILON_U0 = 0.05
+REALIZABLE_SET_EPSILON_U1 = 0.97
+SPATIAL_DIM = 2
+%
+% ---- File specifications ----
+% 
+% Output directory
+OUTPUT_DIR = ../result
+% Output file
+OUTPUT_FILE = DataGenerator2D
+% 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.0001
+NEWTON_STEP_SIZE = 0.7
+NEWTON_LINE_SEARCH_ITER = 100000
+%
+% ----- Quadrature ----
+% 
+% Quadrature Rule
+QUAD_TYPE = GAUSS_LEGENDRE_TENSORIZED
+% Quadrature Order
+QUAD_ORDER = 8
+%

code/src/solvers/csdsolvertrafofpsh2d.cpp

@@ -95,7 +95,7 @@ CSDSolverTrafoFPSH2D::CSDSolverTrafoFPSH2D( Config* settings ) : SNSolver( setti
         blaze::column( Y, q ) = sph.ComputeSphericalBasis( _quadPointsSphere[q][0], _quadPointsSphere[q][1] );
     }
 
-    _O = Y;
+    //_O = Y;
     _O = blaze::trans( Y );
     _M = _O;    // transposed to simplify weight multiplication. We will transpose _M later again
 

code/src/toolboxes/datagenerator2D.cpp

+/*!
+ * \file datagenerator3D.cpp
+ * \brief Class to generate data for the neural entropy closure
+ * \author S. Schotthoefer
+ */
+
+#include "toolboxes/datagenerator2D.h"
+#include "common/config.h"
+#include "quadratures/quadraturebase.h"
+#include "toolboxes/errormessages.h"
+#include "toolboxes/sphericalbase.h"
+
+#include <iostream>
+#include <omp.h>
+
+DataGenerator2D::DataGenerator2D( Config* settings ) : DataGeneratorBase( settings ) {
+    ComputeMoments();
+
+    // Initialize Training Data
+    ComputeSetSize();
+
+    _uSol     = VectorVector( _setSize, Vector( _nTotalEntries, 0.0 ) );
+    _alpha    = VectorVector( _setSize, Vector( _nTotalEntries, 0.0 ) );
+    _hEntropy = std::vector<double>( _setSize, 0.0 );
+}
+
+DataGenerator2D::~DataGenerator2D() {}
+
+void DataGenerator2D::ComputeMoments() {
+    double my, phi;
+
+    for( unsigned idx_quad = 0; idx_quad < _nq; idx_quad++ ) {
+        my                 = _quadPointsSphere[idx_quad][0];
+        phi                = _quadPointsSphere[idx_quad][1];
+        _moments[idx_quad] = _basis->ComputeSphericalBasis( my, phi );
+    }
+}
+
+void DataGenerator2D::SampleSolutionU() {
+    // Use necessary conditions from Monreal, Dissertation, Chapter 3.2.1, Page 26
+
+    // --- Determine stepsizes etc ---
+    double du0 = _settings->GetMaxValFirstMoment() / (double)_gridSize;
+
+    // different processes for different
+    if( _LMaxDegree == 0 ) {
+        // --- sample u in order 0 ---
+        // u_0 = <1*psi>
+
+        for( unsigned idx_set = 0; idx_set < _setSize; idx_set++ ) {
+            _uSol[idx_set][0] = du0 * idx_set;
+        }
+    }
+    else if( _LMaxDegree == 1 && _settings->GetSphericalBasisName() == SPHERICAL_MONOMIALS ) {
+        // Sample points on unit sphere.
+        QuadratureBase* quad = QuadratureBase::Create( _settings );
+        VectorVector qpoints = quad->GetPoints();    // carthesian coordinates.
+        unsigned long nq     = (unsigned long)quad->GetNq();
+
+        // --- sample u in order 0 ---
+        // u_0 = <1*psi>
+
+#pragma omp parallel for schedule( guided )
+        for( unsigned long idx_set = 0; idx_set < _gridSize; idx_set++ ) {
+            unsigned long outerIdx = ( idx_set - 1 ) * ( idx_set ) / 2;    // sum over all radii up to current
+            outerIdx *= nq;                                                // in each radius step, use all quad points
+
+            //  --- sample u in order 1 ---
+            /* order 1 has 3 elements. (omega_x, omega_y,  omega_z) = omega, let u_1 = (u_x, u_y, u_z) = <omega*psi>
+             * Condition u_0 >= norm(u_1)   */
+            double radiusU0 = du0 * ( idx_set + 1 );
+            // Boundary correction
+            if( radiusU0 < _settings->GetRealizableSetEpsilonU0() ) {
+                radiusU0 = _settings->GetRealizableSetEpsilonU0();    // Boundary close to 0
+            }
+
+            unsigned long localIdx = 0;
+            double radius          = 0.0;
+            unsigned long innerIdx = 0;
+            // loop over all radii
+            for( unsigned long idx_subset = 0; idx_subset < idx_set; idx_subset++ ) {
+                radius = du0 * ( idx_subset + 1 );    // dont use radius 0 ==> shift by one
+                // Boundary correction
+                if( radius < _settings->GetRealizableSetEpsilonU0() ) {
+                    radius = _settings->GetRealizableSetEpsilonU0();    // Boundary close to 0
+                }
+                if( radius / radiusU0 > 0.95 * _settings->GetRealizableSetEpsilonU1() ) {
+                    // small offset to take care of rounding errors in quadrature
+                    radius = radiusU0 * 0.95 * _settings->GetRealizableSetEpsilonU1();    // "outer" boundary
+                }
+
+                localIdx = outerIdx + idx_subset * nq;
+
+                for( unsigned long quad_idx = 0; quad_idx < nq; quad_idx++ ) {
+                    innerIdx = localIdx + quad_idx;    // gives the global index
+
+                    _uSol[innerIdx][0] = radiusU0;    // Prevent 1st order moments to be too close to the boundary
+                    // scale quadpoints with radius
+                    _uSol[innerIdx][1] = radius * qpoints[quad_idx][0];
+                    _uSol[innerIdx][2] = radius * qpoints[quad_idx][1];
+                    //_uSol[innerIdx][3] = radius * qpoints[quad_idx][2];
+                }
+            }
+        }
+    }
+    else {
+        ErrorMessages::Error( "Sampling for this configuration is not yet supported", CURRENT_FUNCTION );
+    }
+}
+
+void DataGenerator2D::CheckRealizability() {
+    double epsilon = _settings->GetRealizableSetEpsilonU0();
+    if( _LMaxDegree == 1 ) {
+#pragma omp parallel for schedule( guided )
+        for( unsigned idx_set = 0; idx_set < _setSize; idx_set++ ) {
+            double normU1 = 0.0;
+            Vector u1( 3, 0.0 );
+            if( _uSol[idx_set][0] < epsilon ) {
+                if( std::abs( _uSol[idx_set][1] ) > 0 || std::abs( _uSol[idx_set][2] ) > 0 ) {
+                    ErrorMessages::Error( "Moment not realizable [code 0]. Values: (" + std::to_string( _uSol[idx_set][0] ) + "|" +
+                                              std::to_string( _uSol[idx_set][1] ) + "|" + std::to_string( _uSol[idx_set][2] ) + ")",
+                                          CURRENT_FUNCTION );
+                }
+            }
+            else {
+                u1     = { _uSol[idx_set][1], _uSol[idx_set][2] };
+                normU1 = norm( u1 );
+                if( normU1 / _uSol[idx_set][0] > _settings->GetRealizableSetEpsilonU1() ) {    // Danger Hardcoded
+                    std::cout << "normU1 / _uSol[" << idx_set << "][0]: " << normU1 / _uSol[idx_set][0] << "\n";
+                    std::cout << "normU1: " << normU1 << " | _uSol[idx_set][0] " << _uSol[idx_set][0] << "\n";
+                    ErrorMessages::Error( "Moment to close to boundary of realizable set [code 1].\nBoundary ratio: " +
+                                              std::to_string( normU1 / _uSol[idx_set][0] ),
+                                          CURRENT_FUNCTION );
+                }
+                if( normU1 / _uSol[idx_set][0] <= 0 /*+ 0.5 * epsilon*/ ) {
+                    // std::cout << "_uSol" << _uSol[idx_set][1] << " | " << _uSol[idx_set][2] << " | " << _uSol[idx_set][3] << " \n";
+                    // std::cout << "normU1 / _uSol[" << idx_set << "][0]: " << normU1 / _uSol[idx_set][0] << "\n";
+                    // std::cout << "normU1: " << normU1 << " | _uSol[idx_set][0] " << _uSol[idx_set][0] << "\n";
+                    ErrorMessages::Error( "Moment to close to boundary of realizable set [code 2].\nBoundary ratio: " +
+                                              std::to_string( normU1 / _uSol[idx_set][0] ),
+                                          CURRENT_FUNCTION );
+                }
+            }
+        }
+    }
+}
+
+void DataGenerator2D::ComputeSetSize() {
+    if( _LMaxDegree == 0 ) {
+    }
+    else if( _LMaxDegree == 1 && _settings->GetSphericalBasisName() == SPHERICAL_MONOMIALS ) {
+        // Sample points on unit sphere.
+        QuadratureBase* quad = QuadratureBase::Create( _settings );
+        unsigned long nq     = (unsigned long)quad->GetNq();
+
+        // Allocate memory.
+        _setSize = nq * _gridSize * ( _gridSize - 1 ) / 2;
+
+        delete quad;
+    }
+    else {
+        ErrorMessages::Error( "Sampling of training data of degree higher than 1 is not yet implemented.", CURRENT_FUNCTION );
+    }
+}

code/src/toolboxes/datageneratorbase.cpp

@@ -11,6 +11,7 @@
 #include "quadratures/quadraturebase.h"
 #include "spdlog/spdlog.h"
 #include "toolboxes/datagenerator1D.h"
+#include "toolboxes/datagenerator2D.h"
 #include "toolboxes/datagenerator3D.h"
 #include "toolboxes/errormessages.h"
 #include "toolboxes/sphericalbase.h"
@@ -69,7 +70,7 @@ DataGeneratorBase::~DataGeneratorBase() {
 DataGeneratorBase* DataGeneratorBase::Create( Config* settings ) {
     switch( settings->GetDim() ) {
         case 1: return new DataGenerator1D( settings );
-        case 2: ErrorMessages::Error( "2D Sampling is not yet supported.", CURRENT_FUNCTION );
+        case 2: return new DataGenerator2D( settings );
         case 3: return new DataGenerator3D( settings );
         default: ErrorMessages::Error( "Sampling for more than 3 dimensions is not yet supported.", CURRENT_FUNCTION );
     }

code/src/toolboxes/sphericalmonomials.cpp

@@ -130,9 +130,9 @@ double SphericalMonomials::Omega_y( double my, double phi ) { return sqrt( 1 - m
 double SphericalMonomials::Omega_z( double my ) { return my; }
 
 double SphericalMonomials::Power( double basis, unsigned exponent ) {
-    if( exponent == 0 ) return 1.0;
-    double result = basis;
-    for( unsigned i = 1; i < exponent; i++ ) {
+    if( exponent == 0 ) return 1.0;               // basis^0
+    double result = basis;                        // basis^1
+    for( unsigned i = 1; i < exponent; i++ ) {    // exp> 1
         result = result * basis;
     }
     return result;