small tidy up in mnsolver. Started Data Generator

code/include/common/config.h

@@ -67,7 +67,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 */
 
@@ -104,6 +104,11 @@ 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 */
+
     // --- Parsing Functionality and Initializing of Options ---
     /*!
      * @brief Set default values for all options not yet set.
@@ -282,6 +287,11 @@ class Config
     std::vector<SCALAR_OUTPUT> inline GetHistoryOutput() { return _historyOutput; }
     unsigned short inline GetNHistoryOutput() { return _nHistoryOutput; }
     unsigned short inline GetHistoryOutputFrequency() { return _historyOutputFrequency; }
+
+    // Data generator
+    bool inline GetDataGeneratorMode() { return _dataGeneratorMode; }
+    unsigned long inline GetTrainingDataSetSize() { return _tainingSetSize; }
+
     // ---- Setters for option structure
 
     // Quadrature Structure
@@ -289,6 +299,7 @@ class Config
     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; }
 };

code/include/solvers/mnsolver.h

@@ -22,7 +22,7 @@ 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 */

code/include/solvers/solverbase.h

@@ -25,12 +25,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 */
@@ -39,9 +39,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) */

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 SphericalHarmonics;
+class QuadratureBase;
+class Config;
+
+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 */
+
+    std::vector<std::vector<double>> _uSol;  /*! @brief: vector with moments. Size: (setSize,basisSize)*/
+    std::vector<std::vector<double>> _alpha; /*! @brief: vector with Lagrange multipliers. Size: (setSize,basisSize)*/
+    std::vector<double> _hEntropy;           /*! @brief: vector with entropy values. Size: (setSize) */
+
+    unsigned _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) */
+
+    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 */
+
+    // Helper functions
+    /*! @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;
+    void ComputeMoments(); /*! @brief : Pre-Compute Moments at all quadrature points. */
+
+    // Main methods
+    void sampleSolutionU(); /*! @brief : Samples solution vectors u */
+};
+#endif    // DATAGENERATOR_H

code/input/exampleMN.cfg

@@ -33,7 +33,8 @@ TIME_FINAL = 0.3
 % Maximal Moment degree
 MAX_MOMENT_SOLVER = 0
 %
-%% Entropy settings
+% ---- Entropy settings ----
+%
 ENTROPY_FUNCTIONAL = MAXWELL_BOLTZMANN
 ENTROPY_OPTIMIZER = NEWTON
 %
@@ -46,6 +47,7 @@ NEWTON_STEP_SIZE = 0.7
 NEWTON_LINE_SEARCH_ITER = 1000
 %
 % ---- Boundary Conditions ----
+%
 % Example: BC_DIRICLET = (dummyMarker1, dummyMarker2)
 % Dirichlet Boundary
 BC_DIRICHLET = ( void )
@@ -53,9 +55,7 @@ BC_DIRICHLET = ( void )
 % ----- Quadrature ----
 % 
 % Quadrature Rule
-%QUAD_TYPE = MONTE_CARLO
 QUAD_TYPE = GAUSS_LEGENDRE_TENSORIZED
-%
 % Quadrature Order
 QUAD_ORDER = 8
 %

code/src/common/config.cpp

@@ -291,6 +291,13 @@ void Config::SetConfigOptions() {
     AddEnumListOption( "HISTORY_OUTPUT", _nHistoryOutput, _historyOutput, ScalarOutput_Map );
     /*! @brief History output Frequency \n DESCRIPTION: Describes history output write frequency \n DEFAULT 1 \ingroup Config */
     AddUnsignedShortOption( "HISTORY_OUTPUT_FREQUENCY", _historyOutputFrequency, 1 );
+
+    // Data generator related options
+    /*! @brief Size of training data set \n DESCRIPTION: Size of training data set  \n DEFAULT 1 \ingroup Config */
+    AddUnsignedLongOption( "TRAINING_SET_SIZE", _tainingSetSize, 1 );
+    /*! @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 );
 }
 
 void Config::SetConfigParsing( string case_filename ) {

code/src/main.cpp

@@ -17,6 +17,8 @@
 #include "quadratures/qmontecarlo.h"
 #include "solvers/sphericalharmonics.h"
 
+#include "toolboxes/datagenerator.h"
+
 double testFunc( double my, double phi ) { return my * my + phi; }
 
 double testFunc2( double x, double y, double z ) { return x + y + z; }
@@ -30,17 +32,23 @@ int main( int argc, char** argv ) {
 
     // CD  Load Settings from File
     Config* config = new Config( filename );
-
     // Print input file and run info to file
     PrintLogHeader( filename );
 
-    // Build solver
-    Solver* solver = Solver::Create( config );
-
-    // Run solver and export
-    solver->Solve();
-    solver->PrintVolumeOutput();
-
+    if( config->GetDataGeneratorMode() ) {
+        // Build Data generator
+        nnDataGenerator* datagen = new nnDataGenerator( config );
+        // Generate Data and export
+        datagen->computeTrainingData();
+    }
+    else {
+        // Build solver
+        Solver* solver = Solver::Create( config );
+
+        // Run solver and export
+        solver->Solve();
+        solver->PrintVolumeOutput();
+    }
     MPI_Finalize();
     return EXIT_SUCCESS;
 }

code/src/solvers/mnsolver.cpp

@@ -16,20 +16,17 @@
 #include <mpi.h>
 
 #include <fstream>
-//#include <chrono>
 
 MNSolver::MNSolver( Config* settings ) : Solver( settings ) {
-
-    // Is this good (fast) code using a constructor list?
     _LMaxDegree    = settings->GetMaxMomentDegree();
     _nTotalEntries = GlobalIndex( _LMaxDegree, int( _LMaxDegree ) ) + 1;
 
     // build quadrature object and store quadrature points and weights
-    _quadPoints       = _quadrature->GetPoints();
-    _weights          = _quadrature->GetWeights();
-    _nq               = _quadrature->GetNq();
+    _quadPoints = _quadrature->GetPoints();
+    _weights    = _quadrature->GetWeights();
+    //_nq               = _quadrature->GetNq();
     _quadPointsSphere = _quadrature->GetPointsSphere();
-    _settings->SetNQuadPoints( _nq );
+    //_settings->SetNQuadPoints( _nq );
 
     // Initialize Scatter Matrix --
     _scatterMatDiag = Vector( _nTotalEntries, 0.0 );

code/src/solvers/solverbase.cpp

@@ -11,7 +11,9 @@
 #include "solvers/pnsolver.h"
 #include "solvers/snsolver.h"
 
-Solver::Solver( Config* settings ) : _settings( settings ) {
+Solver::Solver( Config* settings ) {
+    _settings = settings;
+
     // @TODO save parameters from settings class
 
     // build mesh and store  and store frequently used params

code/src/toolboxes/datagenerator.cpp

+/*!
+ * \file datagenerator.cpp
+ * \brief Class to generate data for the neural entropy closure
+ * \author S. Schotthoefer
+ */
+
+#include "toolboxes/datagenerator.h"
+#include "common/config.h"
+#include "quadratures/quadraturebase.h"
+#include "solvers/sphericalharmonics.h"
+
+#include "spdlog/spdlog.h"
+
+#include <iostream>
+
+nnDataGenerator::nnDataGenerator( Config* settings ) {
+    _settings = settings;
+    _setSize  = settings->GetTrainingDataSetSize();
+
+    _LMaxDegree    = settings->GetMaxMomentDegree();
+    _nTotalEntries = (unsigned)GlobalIndex( _LMaxDegree, _LMaxDegree ) + 1;
+
+    // build quadrature object and store frequently used params
+    _quadrature       = QuadratureBase::CreateQuadrature( settings );
+    _nq               = _quadrature->GetNq();
+    _quadPoints       = _quadrature->GetPoints();
+    _weights          = _quadrature->GetWeights();
+    _quadPointsSphere = _quadrature->GetPointsSphere();
+    _settings->SetNQuadPoints( _nq );
+
+    _basis   = new SphericalHarmonics( _LMaxDegree );
+    _moments = VectorVector( _nq, Vector( _nTotalEntries, 0.0 ) );
+    ComputeMoments();
+
+    // Initialize Training Data
+    _uSol     = std::vector( _setSize, std::vector( _nTotalEntries, 0.0 ) );
+    _alpha    = std::vector( _setSize, std::vector( _nTotalEntries, 0.0 ) );
+    _hEntropy = std::vector( _setSize, 0.0 );
+}
+
+void nnDataGenerator::computeTrainingData() {
+    // Prototype: Only for _LMaxDegree == 1
+    // Prototype: u is sampled from [0,100]
+
+    // --- sample u ---
+    sampleSolutionU();
+}
+
+int nnDataGenerator::GlobalIndex( int l, int k ) const {
+    int numIndicesPrevLevel  = l * l;    // number of previous indices untill level l-1
+    int prevIndicesThisLevel = k + l;    // number of previous indices in current level
+    return numIndicesPrevLevel + prevIndicesThisLevel;
+}
+
+void nnDataGenerator::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 nnDataGenerator::sampleSolutionU() {
+    double du = 100.0 / (double)_setSize;    // Prototype: u is sampled from [0,100]
+    auto log  = spdlog::get( "event" );
+
+    for( unsigned idx_set = 0; idx_set < _setSize; idx_set++ ) {
+        _uSol[idx_set][0] = du * idx_set;
+        log->info( "{}", _uSol[idx_set][0] );
+    }
+    log->flush();
+}