Commit 03bec8b6 authored by Steffen Schotthöfer's avatar Steffen Schotthöfer
Browse files

added 2D Data Generator and example cfg file

parent 67e3aada
/*!
* \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
......@@ -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)*/
......
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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
%
......@@ -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
......
/*!
* \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 );
}
}
......@@ -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 );
}
......
......@@ -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;
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment