Skip to content

GitLab

Commit 7abca39a authored by Steffen Schotthöfer's avatar Steffen Schotthöfer
Browse files

started monomial basis data generator and mn solver 


Former-commit-id: 64319962
parent 16e9da8f
Hide whitespace changes
Inline Side-by-side
code/include/solvers/mnsolver.h
View file @ 7abca39a
......@@ -4,7 +4,7 @@
#include "solverbase.h"
class EntropyBase;
class SphericalHarmonics;
class SphericalBase;
class OptimizerBase;
class MNSolver : public Solver
......@@ -25,8 +25,8 @@ class MNSolver : public Solver
unsigned short _LMaxDegree; /*! @brief: Max Order of Spherical Harmonics */
// Moment basis
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 */
SphericalBase* _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 */
// Scattering
Vector _scatterMatDiag; /*! @brief: Diagonal of the scattering matrix (its a diagonal matrix by construction) */
......
code/include/toolboxes/datagenerator.h
View file @ 7abca39a
......@@ -43,7 +43,7 @@ class nnDataGenerator
VectorVector _alpha; /*! @brief: vector with Lagrange multipliers. Size: (setSize,basisSize)*/
std::vector<double> _hEntropy; /*! @brief: vector with entropy values. Size: (setSize) */
unsigned _setSize;
unsigned long _setSize;
unsigned short _LMaxDegree; /*! @brief: Max Order of Spherical Harmonics */
unsigned _nTotalEntries; /*! @brief: Total number of equations in the system */
......@@ -59,20 +59,15 @@ class nnDataGenerator
NewtonOptimizer* _optimizer; /*! @brief: Class to solve minimal entropy problem */
EntropyBase* _entropy; /*! @brief: Class to handle entropy functional evaluations */
// 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 */
void computeEntropyH_dual(); /*! @brief : Compute the entropy functional at (u,alpha) in dual formulation */
void computeEntropyH_primal(); /*! @brief: Compute the entropy functional at (u,alpha) in primal formulation */
void SampleSolutionU(); /*! @brief: Samples solution vectors u */
void ComputeEntropyH_dual(); /*! @brief: Compute the entropy functional at (u,alpha) in dual formulation */
void ComputeEntropyH_primal(); /*! @brief: Compute the entropy functional at (u,alpha) in primal formulation */
// IO routines
void printTrainingData(); /*! @brief : Print computed training data to csv file and screen */
void PrintTrainingData(); /*! @brief : Print computed training data to csv file and screen */
// Helper functions
void ComputeMoments(); /*! @brief: Pre-Compute Moments at all quadrature points. */
};
#endif // DATAGENERATOR_H
code/input/exampleMN.cfg
View file @ 7abca39a
......@@ -24,6 +24,9 @@ SCATTER_COEFF = 1
%
% ---- Solver specifications ----
%
% Choice of basis functions
SPHERICAL_BASIS = SPHERICAL_HARMONICS
%
% Solver type
SOLVER = MN_SOLVER
% CFL number
......@@ -31,7 +34,7 @@ CFL_NUMBER = 0.7
% Final time for simulation
TIME_FINAL = 0.3
% Maximal Moment degree
MAX_MOMENT_SOLVER = 0
MAX_MOMENT_SOLVER = 1
%
% ---- Entropy settings ----
%
......
code/src/solvers/mnsolver.cpp
View file @ 7abca39a
......@@ -8,7 +8,7 @@
#include "problems/problembase.h"
#include "quadratures/quadraturebase.h"
#include "toolboxes/errormessages.h"
#include "toolboxes/sphericalharmonics.h"
#include "toolboxes/sphericalbase.h"
#include "toolboxes/textprocessingtoolbox.h"
// externals
......@@ -19,7 +19,8 @@
MNSolver::MNSolver( Config* settings ) : Solver( settings ) {
_LMaxDegree = settings->GetMaxMomentDegree();
_nTotalEntries = GlobalIndex( _LMaxDegree, int( _LMaxDegree ) ) + 1;
_basis = SphericalBase::Create( _settings );
_nTotalEntries = _basis->GetBasisSize();
// build quadrature object and store quadrature points and weights
_quadPoints = _quadrature->GetPoints();
......@@ -46,8 +47,6 @@ MNSolver::MNSolver( Config* settings ) : Solver( settings ) {
_alpha = VectorVector( _nCells, Vector( _nTotalEntries, 0.0 ) );
// Initialize and Pre-Compute Moments at quadrature points
_basis = new SphericalHarmonics( _LMaxDegree );
_moments = VectorVector( _nq, Vector( _nTotalEntries, 0.0 ) );
ComputeMoments();
......@@ -103,10 +102,9 @@ Vector MNSolver::ConstructFlux( unsigned idx_cell ) {
for( unsigned idx_neigh = 0; idx_neigh < _neighbors[idx_cell].size(); idx_neigh++ ) {
// Left side reconstruction
if( _reconsOrder > 1 ) {
alphaL = _alpha[idx_cell] +
_solDx[idx_cell] * ( _interfaceMidPoints[idx_cell][idx_neigh][0] - _cellMidPoints[idx_cell][0] ) +
_solDy[idx_cell] * ( _interfaceMidPoints[idx_cell][idx_neigh][1] - _cellMidPoints[idx_cell][1] );
}
alphaL = _alpha[idx_cell] + _solDx[idx_cell] * ( _interfaceMidPoints[idx_cell][idx_neigh][0] - _cellMidPoints[idx_cell][0] ) +
_solDy[idx_cell] * ( _interfaceMidPoints[idx_cell][idx_neigh][1] - _cellMidPoints[idx_cell][1] );
}
else {
alphaL = _alpha[idx_cell];
}
......@@ -118,8 +116,10 @@ Vector MNSolver::ConstructFlux( unsigned idx_cell ) {
else {
if( _reconsOrder > 1 ) {
alphaR = _alpha[_neighbors[idx_cell][idx_neigh]] +
_solDx[_neighbors[idx_cell][idx_neigh]] * ( _interfaceMidPoints[idx_cell][idx_neigh][0] - _cellMidPoints[_neighbors[idx_cell][idx_neigh]][0] ) +
_solDy[_neighbors[idx_cell][idx_neigh]] * ( _interfaceMidPoints[idx_cell][idx_neigh][1] - _cellMidPoints[_neighbors[idx_cell][idx_neigh]][1] );
_solDx[_neighbors[idx_cell][idx_neigh]] *
( _interfaceMidPoints[idx_cell][idx_neigh][0] - _cellMidPoints[_neighbors[idx_cell][idx_neigh]][0] ) +
_solDy[_neighbors[idx_cell][idx_neigh]] *
( _interfaceMidPoints[idx_cell][idx_neigh][1] - _cellMidPoints[_neighbors[idx_cell][idx_neigh]][1] );
}
else {
alphaR = _alpha[_neighbors[idx_cell][idx_neigh]];
......
code/src/toolboxes/datagenerator.cpp
View file @ 7abca39a
......@@ -8,6 +8,7 @@
#include "common/config.h"
#include "entropies/entropybase.h"
#include "optimizers/newtonoptimizer.h"
#include "quadratures/qlebedev.h"
#include "quadratures/quadraturebase.h"
#include "toolboxes/errormessages.h"
#include "toolboxes/sphericalbase.h"
......@@ -59,22 +60,16 @@ void nnDataGenerator::computeTrainingData() {
// Prototype: u is sampled from [0,100]
// --- sample u ---
sampleSolutionU();
SampleSolutionU();
// --- compute alphas ---
_optimizer->SolveMultiCell( _alpha, _uSol, _moments );
// --- compute entropy functional ---
computeEntropyH_primal();
ComputeEntropyH_primal();
// --- Print everything ----
printTrainingData();
}
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;
PrintTrainingData();
}
void nnDataGenerator::ComputeMoments() {
......@@ -88,29 +83,83 @@ void nnDataGenerator::ComputeMoments() {
}
}
void nnDataGenerator::sampleSolutionU() {
void nnDataGenerator::SampleSolutionU() {
// Use necessary conditions from Monreal, Dissertation, Chapter 3.2.1, Page 26
// --- sample u in order 0 ---
// u_0 = <1*psi>
double du = 100.0 / (double)_setSize; // Prototype: u is sampled from [0,100]
// --- Determine stepsizes etc ---
double du0 = 10.0 / (double)_setSize; // Prototype: u is sampled from [0,10]
for( unsigned idx_set = 0; idx_set < _setSize; idx_set++ ) {
_uSol[idx_set][0] = du * idx_set;
}
// different processes for different
if( _LMaxDegree == 0 ) {
// --- sample u in order 0 ---
// u_0 = <1*psi>
// --- 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) */
for( unsigned idx_set = 0; idx_set < _setSize; idx_set++ ) {
_uSol[idx_set][0] = du0 * idx_set;
}
}
if( _LMaxDegree == 1 ) {
// Sample points on unit sphere. (Use Lebedev quadrature)
unsigned order = 5; // is avail.
QLebedev* quad = new QLebedev( order );
VectorVector qpoints = quad->GetPoints(); // carthesian coordinates.
unsigned long nq = (unsigned long)quad->GetNq();
// Allocate memory.
unsigned long setSizeU0 = _setSize;
unsigned long setSizeU1 = nq * setSizeU0 * ( setSizeU0 + 1 ) / 2;
_setSize = setSizeU1;
// REFACTOR THIS
_uSol = VectorVector( _setSize, Vector( _nTotalEntries, 0.0 ) );
_alpha = VectorVector( _setSize, Vector( _nTotalEntries, 0.0 ) );
_hEntropy = std::vector<double>( _setSize, 0.0 );
// --- sample u in order 0 ---
// u_0 = <1*psi>
for( unsigned long idx_set = 0; idx_set < setSizeU0; 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) */
// loop over all radii
for( unsigned long idx_subset = 0; idx_subset < idx_set; idx_subset++ ) {
double radius = du0 * ( idx_subset ); // dont use radius 0
unsigned long localIdx = outerIdx + idx_subset * nq;
for( unsigned quad_idx = 0; quad_idx < nq; quad_idx++ ) {
unsigned long radiusIdx = localIdx + quad_idx; // gives the global index
_uSol[radiusIdx][0] = radius;
// scale quadpoints with radius
_uSol[radiusIdx][1] = radius * qpoints[quad_idx][0];
_uSol[radiusIdx][2] = radius * qpoints[quad_idx][1];
_uSol[radiusIdx][3] = radius * qpoints[quad_idx][2];
std::cout << " radiusIdx: " << radiusIdx << "\n";
// std::cout << " _uSol[radiusIdx][0]: " << _uSol[radiusIdx][0] << "\n";
// std::cout << " _uSol[radiusIdx][1]: " << _uSol[radiusIdx][1] << "\n";
// std::cout << " _uSol[radiusIdx][2]: " << _uSol[radiusIdx][2] << "\n";
// std::cout << " _uSol[radiusIdx][3]: " << _uSol[radiusIdx][3] << "\n";
}
}
}
}
else {
ErrorMessages::Error( "Sampling for order higher than 1 is not yet supported", CURRENT_FUNCTION );
}
}
void nnDataGenerator::computeEntropyH_dual() {
void nnDataGenerator::ComputeEntropyH_dual() {
for( unsigned idx_set = 0; idx_set < _setSize; idx_set++ ) {
_hEntropy[idx_set] = _optimizer->ComputeObjFunc( _alpha[idx_set], _uSol[idx_set], _moments );
}
}
void nnDataGenerator::computeEntropyH_primal() {
void nnDataGenerator::ComputeEntropyH_primal() {
double result = 0.0;
for( unsigned idx_set = 0; idx_set < _setSize; idx_set++ ) {
......@@ -123,7 +172,7 @@ void nnDataGenerator::computeEntropyH_primal() {
}
}
void nnDataGenerator::printTrainingData() {
void nnDataGenerator::PrintTrainingData() {
auto log = spdlog::get( "event" );
auto logCSV = spdlog::get( "tabular" );
......
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

Imprint