fixed bug in dataGenerator. Added option to determine min distance to boundary of realizable set

Former-commit-id: 015a1b18

code/include/common/config.h

@@ -112,8 +112,9 @@ class Config
     // 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 */
+    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 ---
     /*!
@@ -298,6 +299,7 @@ class Config
     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

code/input/exampleMN.cfg

@@ -66,7 +66,7 @@ QUAD_ORDER = 8
 %
 % ----- Output ---- 
 %
-VOLUME_OUTPUT = (MINIMAL, MOMENTS)
+VOLUME_OUTPUT = (MINIMAL, MOMENTS, DUAL_MOMENTS)
 VOLUME_OUTPUT_FREQUENCY = 1
 SCREEN_OUTPUT = (ITER, MASS, RMS_FLUX, VTK_OUTPUT, CSV_OUTPUT)
 SCREEN_OUTPUT_FREQUENCY = 1

code/src/common/config.cpp

@@ -309,6 +309,9 @@ void Config::SetConfigOptions() {
     /*! @brief Data generator mode \n DESCRIPTION: Check, if data generator mode is active. If yes, no solver is called, but instead the data
      *         generator is executed \n DEFAULT false \ingroup Config */
     AddBoolOption( "DATA_GENERATOR_MODE", _dataGeneratorMode, false );
+    /*! @brief Distance to the boundary of the realizable set \n DESCRIPTION: Distance to the boundary of the realizable set  \n DEFAULT 0.1 \ingroup
+     * Config */
+    AddDoubleOption( "BOUNDARY_DISTANCE_REALIZABLE_SET", _boundaryDistanceRealizableSet, 0.1 );
 }
 
 void Config::SetConfigParsing( string case_filename ) {

code/src/optimizers/newtonoptimizer.cpp

@@ -175,8 +175,16 @@ void NewtonOptimizer::Solve( Vector& lambda, Vector& sol, const VectorVector& mo
             return;
         }
     }
-    ErrorMessages::Error( " Newton did not converge! Norm of gradient is: " + std::to_string( norm( dlambdaNew ) ) + " at cell " +
-                              std::to_string( idx_cell ) + ".\nObjective function value is " +
-                              std::to_string( ComputeObjFunc( lambda, sol, moments ) ) + " .",
+    std::string uSolString = "At moment: (" + std::to_string( sol[0] );
+    for( unsigned i = 1; i < nSize; i++ ) {
+        uSolString += " | " + std::to_string( sol[i] );
+    }
+    uSolString += ").";
+
+    Vector u1     = { sol[1], sol[2], sol[3] };
+    double normU1 = norm( u1 );
+    ErrorMessages::Error( "Newton did not converge at cell " + std::to_string( idx_cell ) + "\n" + uSolString +
+                              "\nNorm of gradient: " + std::to_string( norm( dlambdaNew ) ) + "\nObjective function value: " +
+                              std::to_string( ComputeObjFunc( lambda, sol, moments ) ) + "\nBoundary Ratio: " + std::to_string( normU1 / sol[0] ),
                           CURRENT_FUNCTION );
 }

code/src/toolboxes/datagenerator.cpp

@@ -109,9 +109,7 @@ void nnDataGenerator::SampleSolutionU() {
         }
     }
     if( _LMaxDegree == 1 ) {
-        // Sample points on unit sphere. (Use Lebedev quadrature)
-        // unsigned order       = 5;         // is avail.
-        double epsilon       = 0.5;
+        // Sample points on unit sphere.
         QuadratureBase* quad = QuadratureBase::Create( _settings );
         VectorVector qpoints = quad->GetPoints();    // carthesian coordinates.
         unsigned long nq     = (unsigned long)quad->GetNq();
@@ -138,6 +136,7 @@ void nnDataGenerator::SampleSolutionU() {
             //  --- 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 ) * ( 1 + 2 * _settings->GetBoundaryDistanceRealizableSet() );
 
             unsigned long localIdx = 0;
             double radius          = 0.0;
@@ -150,7 +149,7 @@ void nnDataGenerator::SampleSolutionU() {
                 for( unsigned long quad_idx = 0; quad_idx < nq; quad_idx++ ) {
                     innerIdx = localIdx + quad_idx;    // gives the global index
 
-                    _uSol[innerIdx][0] = radius + 2 * epsilon;
+                    _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];
@@ -210,28 +209,33 @@ void nnDataGenerator::PrintTrainingData() {
 }
 
 void nnDataGenerator::CheckRealizability() {
+    double epsilon = _settings->GetBoundaryDistanceRealizableSet();
     if( _LMaxDegree == 1 ) {
         double normU1 = 0.0;
         Vector u1( 3, 0.0 );
         for( unsigned idx_set = 0; idx_set < _setSize; idx_set++ ) {
-            if( _uSol[idx_set][0] == 0 ) {
-                if( _uSol[idx_set][1] > 0 || _uSol[idx_set][2] > 0 || _uSol[idx_set][3] > 0 ) {
+            if( _uSol[idx_set][0] < epsilon ) {
+                if( std::abs( _uSol[idx_set][1] ) > 0 || std::abs( _uSol[idx_set][2] ) > 0 || std::abs( _uSol[idx_set][3] ) > 0 ) {
                     ErrorMessages::Error( "Moment not realizable [code 0].", CURRENT_FUNCTION );
                 }
             }
             else {
                 u1     = { _uSol[idx_set][1], _uSol[idx_set][2], _uSol[idx_set][3] };
                 normU1 = norm( u1 );
-                if( normU1 / _uSol[idx_set][0] >= 1 ) {
+                if( normU1 / _uSol[idx_set][0] > 1 - 0.99 * epsilon ) {
                     // 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 not realizable [code 1].", CURRENT_FUNCTION );
+                    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 ) {
+                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 not realizable [code 2].", CURRENT_FUNCTION );
+                    ErrorMessages::Error( "Moment to close to boundary of realizable set [code 2].\nBoundary ratio: " +
+                                              std::to_string( normU1 / _uSol[idx_set][0] ),
+                                          CURRENT_FUNCTION );
                 }
             }
         }