Commit d1914238 authored by steffen.schotthoefer's avatar steffen.schotthoefer
Browse files

added Option for Optimizer, added Optimizer, MNSolver compiles - not tested -...

added Option for Optimizer, added Optimizer, MNSolver compiles - not tested - base framework designed, added Entropy option and Entropy classes, added flux reconstruction step
parent c30bd657
Pipeline #94723 failed with stages
in 6 minutes and 4 seconds
#ifndef ENTROPYBASE_H
#define ENTROPYBASE_H
#include "settings/config.h"
class EntropyBase
{
public:
inline EntropyBase(){};
static EntropyBase* Create( Config* settings );
/*! @brief: computes the entropy functional
* @param: z = point where the functional should be evaluated.
* z must be in domain of the functional
* @returns: value of entropy functional at z */
virtual double Entropy( double z ) = 0;
/*! @brief: computes the derivative of the entropy functional
* @param: z = point where the derivative should be evaluated.
* z must be in domain of the functional
* @returns: value of derivative at z */
virtual double EntropyPrime( double z ) = 0;
/*! @brief: computes the dual of the entropy functional
* @param: z = point where the dual of the functional should be evaluated.
* z must be in domain of the duaö functional
* @returns: value of entropy functional at z */
virtual double EntropyDual( double y ) = 0;
/*! @brief: computes the dual of the derivative of the entropy functional
* @param: z = point where the dual of the derivative should be evaluated.
* z must be in domain of the dual functional
* @returns: value of dual of the derivative at z */
virtual double EntropyPrimeDual( double y ) = 0;
/*! @brief: checks, if value is in domain of entropy */
virtual bool CheckDomain( double z ) = 0;
};
#endif // ENTROPYBASE_H
#ifndef MAXWELLBOLTZMANNENTROPY_H
#define MAXWELLBOLTZMANNENTROPY_H
#include "entropybase.h"
#include <cmath>
class MaxwellBoltzmannEntropy : public EntropyBase
{
public:
inline MaxwellBoltzmannEntropy() : EntropyBase(){};
inline double Entropy( double z ) override { return z * log( z ) - z; };
inline double EntropyPrime( double z ) override { return log( z ); };
inline double EntropyDual( double y ) override { return exp( y ); };
inline double EntropyPrimeDual( double y ) override { return exp( y ); };
inline bool CheckDomain( double z ) override { return ( z >= 0 && std::isnan( z ) ); };
};
#endif // MAXWELLBOLTZMANNENTROPY_H
#ifndef QUADRATICENTROPY_H
#define QUADRATICENTROPY_H
#include "entropybase.h"
#include <cmath>
class QuadraticEntropy : public EntropyBase
{
public:
inline QuadraticEntropy() : EntropyBase(){};
inline double Entropy( double z ) override { return 0.5 * z * z; };
inline double EntropyPrime( double z ) override { return z; };
inline double EntropyDual( double y ) override { return 0.5 * y * y; };
inline double EntropyPrimeDual( double y ) override { return y; };
inline bool CheckDomain( double z ) override { return std::isnan( z ); };
};
#endif // QUADRATICENTROPY_H
#ifndef NEWTONOPTIMIZER_H
#define NEWTONOPTIMIZER_H
#include "optimizerbase.h"
class NewtonOptimizer : public OptimizerBase
{
public:
inline NewtonOptimizer( Config* settings ) : OptimizerBase( settings ){};
virtual Vector Solve( Vector u ) override;
};
#endif // NEWTONOPTIMIZER_H
#ifndef OPTIMIZERBASE_H
#define OPTIMIZERBASE_H
#include "entropies/entropybase.h"
#include "settings/config.h"
#include "settings/typedef.h"
class OptimizerBase
{
public:
OptimizerBase( Config* settings );
static OptimizerBase* Create( Config* settings );
/*! @brief : Computes the optimal Lagrange multilpiers for the dual entropy minimization problem
* @param : Vector u = pointer to vector of given moments. // Maybe use pointer for performance?
* @return : Vector alpha = optimal lagrange multipliers. Has the same length as Vector u. */
virtual Vector Solve( Vector u ) = 0;
protected:
EntropyBase* _entropy; /*! @brief: Class to handle entropy functional evaluations */
Config* _settings;
};
#endif // OPTIMIZERBASE_H
......@@ -55,6 +55,7 @@ class Config
double _tEnd; /*!< @brief Final Time for Simulation */
PROBLEM_NAME _problemName; /*!< @brief Name of predefined Problem */
SOLVER_NAME _solverName; /*!< @brief Name of the used Solver */
ENTROPY_NAME _entropyName; /*!< @brief Name of the used Entropy Functional */
unsigned short _maxMomentDegree; /*!< @brief Maximal Order of Moments for PN and MN Solver */
/*!< @brief If true, very low entries (10^-10 or smaller) of the flux matrices will be set to zero,
......@@ -73,6 +74,9 @@ class Config
// Scattering Kernel
KERNEL_NAME _kernelName; /*!< @brief Scattering Kernel Name*/
// Optimizer
OPTIMIZER_NAME _entropyOptimizerName; /*!< @brief Choice of optimizer */
// --- Parsing Functionality and Initializing of Options ---
/*!
* @brief Set default values for all options not yet set.
......@@ -208,7 +212,12 @@ class Config
double inline GetTEnd() const { return _tEnd; }
PROBLEM_NAME inline GetProblemName() const { return _problemName; }
SOLVER_NAME inline GetSolverName() const { return _solverName; }
ENTROPY_NAME inline GetEntropyName() const { return _entropyName; }
bool inline GetCleanFluxMat() const { return _cleanFluxMat; }
// Optimizer
OPTIMIZER_NAME inline GetOptimizerName() const { return _entropyOptimizerName; }
// Boundary Conditions
BOUNDARY_TYPE GetBoundaryType( std::string nameMarker ) const; /*! @brief Get Boundary Type of given marker */
......
......@@ -53,4 +53,14 @@ enum SOLVER_NAME { SN_SOLVER, PN_SOLVER, MN_SOLVER };
inline std::map<std::string, SOLVER_NAME> Solver_Map{ { "SN_SOLVER", SN_SOLVER }, { "PN_SOLVER", PN_SOLVER }, { "MN_SOLVER", MN_SOLVER } };
// Entropy functional
enum ENTROPY_NAME { QUADRATIC, MAXWELL_BOLZMANN, BOSE_EINSTEIN, FERMI_DIRAC };
inline std::map<std::string, ENTROPY_NAME> Entropy_Map{
{ "QUADRATIC", QUADRATIC }, { "MAXWELL_BOLZMANN", MAXWELL_BOLZMANN }, { "BOSE_EINSTEIN", BOSE_EINSTEIN }, { "FERMI_DIRAC", FERMI_DIRAC } };
// Otpimizer
enum OPTIMIZER_NAME { NEWTON, ML };
inline std::map<std::string, OPTIMIZER_NAME> Optimizer_Map{ { "NEWTON", NEWTON }, { "ML", ML } };
#endif // GLOBAL_CONSTANTS_H
#ifndef MNSOLVER_H
#define MNSOLVER_H
#include "entropies/entropybase.h"
#include "optimizers/optimizerbase.h"
#include "solverbase.h"
#include "sphericalharmonics.h"
......@@ -34,6 +36,11 @@ class MNSolver : public Solver
layout: Nx2 (in 2d), N = size of system */
QuadratureBase* _quadrature; /*! @brief: Quadrature rule to compute flux Jacobian */
EntropyBase* _entropy; /*! @brief: Class to handle entropy functionals */
VectorVector _alpha; /*! @brief: Lagrange Multipliers for Minimal Entropy problem for each gridCell
Layout: _nCells x _nTotalEntries*/
OptimizerBase* _optimizer; /*! @brief: Class to solve minimal entropy problem */
/*! @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
......@@ -44,5 +51,11 @@ class MNSolver : public Solver
/*! @brief : function for computing and setting up the System Matrices.
The System Matrices are diagonal, so there is no need to diagonalize*/
void ComputeSystemMatrices();
/*! @brief : Construct flux by computing the Moment of the FVM discretization at the interface of cell
* @param : idx_cell = current cell id
* @param : idx_neighbor = neighbor cell id
* @returns : flux for all moments at interface of idx_cell,idx_neighbor */
Vector ConstructFlux( unsigned idx_cell, unsigned idx_neighbor );
};
#endif // MNSOLVER_H
......@@ -11,7 +11,9 @@
#ifndef SPHERICALHARMONICS_H
#define SPHERICALHARMONICS_H
#include "settings/typedef.h"
#include <vector>
class SphericalHarmonics
{
public:
......@@ -28,13 +30,13 @@ class SphericalHarmonics
* @param : phi - spherical coordinate, 0 <= phi <= 2*pi
* @return : vector of basis functions at point (my, phi) with size N = L² +2L
*/
std::vector<double> ComputeSphericalBasis( double my, double phi );
Vector ComputeSphericalBasis( double my, double phi );
/*! @brief : Computes all N = L² +2L basis functions at point (x, y, z) on the unit sphere
* @param : x,y,z = coordinates on unit sphere
* @return : vector of basis functions at point (x,y,z) with size N = L² +2L
*/
std::vector<double> ComputeSphericalBasis( double x, double y, double z );
Vector ComputeSphericalBasis( double x, double y, double z );
private:
/*! @brief: maximal degree of the spherical harmonics basis (this is "L" in the comments)*/
......@@ -56,7 +58,7 @@ class SphericalHarmonics
* degree 0 <= l <= L and order -l <= k <= l
* length : N = L² +2L
*/
std::vector<double> _YBasis;
Vector _YBasis;
/*! @brief : helper function to get the global index for given k and l of
* the associated legendre polynomial P_k^l.
......
#include "entropies/entropybase.h"
#include "entropies/maxwellboltzmannentropy.h"
#include "entropies/quadraticentropy.h"
EntropyBase* EntropyBase::Create( Config* settings ) {
switch( settings->GetEntropyName() ) {
case QUADRATIC: return new QuadraticEntropy();
case MAXWELL_BOLZMANN: return new MaxwellBoltzmannEntropy();
default: return new QuadraticEntropy();
}
}
#include "optimizers/newtonoptimizer.h"
Vector NewtonOptimizer::Solve( Vector u ) {
// if we have quadratic entropy, then alpha = u;
if( _settings->GetEntropyName() == QUADRATIC ) return u;
return Vector( 1, 0.0 ); // dummy
}
#include "optimizers/optimizerbase.h"
#include "optimizers/newtonoptimizer.h"
OptimizerBase::OptimizerBase( Config* settings ) {
_entropy = EntropyBase::Create( settings );
_settings = settings;
}
OptimizerBase* OptimizerBase::Create( Config* settings ) {
switch( settings->GetOptimizerName() ) {
case NEWTON: return new NewtonOptimizer( settings );
// extend to other optimizers
default: return new NewtonOptimizer( settings );
}
}
......@@ -208,6 +208,10 @@ void Config::SetConfigOptions() {
/*! @brief CleanFluxMatrices \n DESCRIPTION: If true, very low entries (10^-10 or smaller) of the flux matrices will be set to zero,
* to improve floating point accuracy \n DEFAULT false \ingroup Config */
AddBoolOption( "CLEAN_FLUX_MATRICES", _cleanFluxMat, false );
/*! @brief Entropy Functional \n DESCRIPTION: Entropy functional used for the MN_Solver \n DEFAULT QUADRTATIC @ingroup Config. */
AddEnumOption( "ENTROPY_FUNCTIONAL", _entropyName, Entropy_Map, QUADRATIC );
/*! @brief Optimizer Name \n DESCRIPTION: Optimizer used to determine the minimal Entropy reconstruction \n DEFAULT NEWTON \ingroup Config */
AddEnumOption( "ENTROPY_OPTIMIZER", _entropyOptimizerName, Optimizer_Map, NEWTON );
// Mesh related options
// Boundary Markers
......
......@@ -24,9 +24,12 @@ MNSolver::MNSolver( Config* settings ) : Solver( settings ), _nMaxMomentsOrder(
// Initialize System Matrices
_A = VectorVector( _nTotalEntries );
// Fill System Matrices
ComputeSystemMatrices();
TextProcessingToolbox::PrintVectorVector( _A );
_entropy = EntropyBase::Create( settings );
// Initialize lagrange Multipliers
_alpha = VectorVector( _nCells );
for( unsigned idx_cells = 0; idx_cells < _nTotalEntries; idx_cells++ ) {
_alpha[idx_cells] = Vector( _nTotalEntries, 0.0 );
}
}
MNSolver::~MNSolver() { delete _quadrature; }
......@@ -46,7 +49,7 @@ void MNSolver::ComputeSystemMatrices() {
}
// Compute Integrals _A[idx_system][i]=<Omega_i*m^(l,k)>
std::vector<double> funcEval( _nTotalEntries, 0.0 );
Vector funcEval( _nTotalEntries, 0.0 );
for( unsigned i = 0; i < _nq; i++ ) {
double x = _quadrature->GetPoints()[i][0];
double y = _quadrature->GetPoints()[i][1];
......@@ -60,6 +63,34 @@ void MNSolver::ComputeSystemMatrices() {
}
}
Vector MNSolver::ConstructFlux( unsigned idx_cell, unsigned idx_neigh ) {
// ---- Integration of Moment of flux ----
double x, y, z, w, entropy, velocity;
Vector moment( _nTotalEntries, 0.0 );
Vector flux( _nTotalEntries, 0.0 );
for( unsigned idx_quad = 0; idx_quad < _nq; idx_quad++ ) {
x = _quadrature->GetPoints()[idx_quad][0];
y = _quadrature->GetPoints()[idx_quad][1];
z = _quadrature->GetPoints()[idx_quad][2];
w = _quadrature->GetWeights()[idx_quad];
moment = _basis.ComputeSphericalBasis( x, y, z );
// --- get upwinding direction (2D)
// DOUBLE CHECK IF WE NEED DIMENSION SPLITTIN!!!
velocity = _normals[idx_cell][idx_neigh][0] * x + _normals[idx_cell][idx_neigh][1] * y;
// Perform upwinding
( velocity > 0 ) ? ( entropy = _entropy->EntropyPrimeDual( blaze::dot( _alpha[idx_cell], moment ) ) )
: ( entropy = _entropy->EntropyPrimeDual( blaze::dot( _alpha[idx_neigh], moment ) ) );
// integrate
flux += moment * ( w * velocity * entropy );
}
return flux;
}
void MNSolver::Solve() {
int rank;
......@@ -86,62 +117,43 @@ void MNSolver::Solve() {
for( unsigned idx_energy = 0; idx_energy < _nEnergies; ++idx_energy ) {
// ------- Reconstruction Step -------
// Todo
// ------- Advection Step ------------
// Loop over the grid cells
for( unsigned idx_cell = 0; idx_cell < _nCells; idx_cell++ ) {
// Loop over all spatial cells
for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
if( _boundaryCells[idx_cell] == BOUNDARY_TYPE::DIRICHLET ) continue; // Dirichlet cells stay at IC, farfield assumption
// Loop over all equations of the system.
for( int idx_lDegree = 0; idx_lDegree <= int( _nMaxMomentsOrder ); idx_lDegree++ ) {
for( int idx_kOrder = -idx_lDegree; idx_kOrder <= idx_lDegree; idx_kOrder++ ) {
idx_system = unsigned( GlobalIndex( idx_lDegree, idx_kOrder ) );
psiNew[idx_cell][idx_system] = 0.0;
// Loop over all neighbor cells (edges) of cell j and compute numerical fluxes
for( unsigned l = 0; l < _neighbors[idx_cell].size(); ++l ) {
// store flux contribution on psiNew_sigmaS to save memory
if( _boundaryCells[idx_cell] == BOUNDARY_TYPE::NEUMANN && _neighbors[idx_cell][l] == _nCells )
psiNew[idx_cell][idx_system] +=
_g->Flux( _A[idx_system], _psi[idx_cell][idx_system], _psi[idx_cell][idx_system], _normals[idx_cell][l] );
else
psiNew[idx_cell][idx_system] += _g->Flux(
_A[idx_system], _psi[idx_cell][idx_system], _psi[_neighbors[idx_cell][l]][idx_system], _normals[idx_cell][l] );
}
// time update angular flux with numerical flux and total scattering cross section
psiNew[idx_cell][idx_system] = _psi[idx_cell][idx_system] -
( _dE / _areas[idx_cell] ) * psiNew[idx_cell][idx_system] /* cell averaged flux */
- _dE * _psi[idx_cell][idx_system] *
( _sigmaA[idx_energy][idx_cell] /* absorbtion influence */
+ _sigmaS[idx_energy][idx_cell] * _scatterMatDiag[idx_system] ); /* scattering influence */
}
// ------- Reconstruction Step -------
_alpha[idx_cell] = _optimizer->Solve( _psi[idx_cell] );
// ------- Compute Fluxes ---------
// Loop over neighbors
for( unsigned idx_neigh = 0; idx_neigh < _neighbors[idx_cell].size(); idx_neigh++ ) {
// Store fluxes in psiNew, to save memory
psiNew[idx_cell] += ConstructFlux( idx_cell, idx_neigh );
}
// TODO: figure out a more elegant way
// add external source contribution
// if( _Q.size() == 1u ) { // constant source for all energies
// if( _Q[0][idx_cell].size() == 1u ) // isotropic source
// psiNew[idx_cell] += _dE * _Q[0][idx_cell][0];
// else
// psiNew[idx_cell] += _dE * _Q[0][idx_cell];
//}
// else {
// if( _Q[0][idx_cell].size() == 1u ) // isotropic source
// psiNew[idx_cell] += _dE * _Q[idx_energy][idx_cell][0];
// else
// psiNew[idx_cell] += _dE * _Q[idx_energy][idx_cell];
//}
// ------ FVM Step ------
// NEED TO VECTORIZE
for( unsigned idx_system = 0; idx_system < _nTotalEntries; idx_system++ ) {
psiNew[idx_cell][idx_system] = _psi[idx_cell][idx_system] -
( _dE / _areas[idx_cell] ) * psiNew[idx_cell][idx_system] /* cell averaged flux */
- _dE * _psi[idx_cell][idx_system] *
( _sigmaA[idx_energy][idx_cell] /* absorbtion influence */
+ _sigmaS[idx_energy][idx_cell] * _scatterMatDiag[idx_system] ); /* scattering influence */
}
}
_psi = psiNew;
_psi = psiNew;
// pseudo time iteration output
double mass = 0.0;
for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
fluxNew[idx_cell] = _psi[idx_cell][0]; // zeroth moment is raditation densitiy we are interested in
_solverOutput[idx_cell] = _psi[idx_cell][0];
mass += _psi[idx_cell][0] * _areas[idx_cell];
}
// Save( n );
Save( idx_energy );
dFlux = blaze::l2Norm( fluxNew - fluxOld );
fluxOld = fluxNew;
if( rank == 0 ) log->info( "{:03.8f} {:01.5e}", _energies[idx_energy], dFlux );
......
......@@ -13,16 +13,16 @@ SphericalHarmonics::SphericalHarmonics( unsigned L_degree ) {
ComputeCoefficients();
unsigned basisSize = GlobalIdxBasis( _LMaxDegree, _LMaxDegree ) + 1;
_YBasis = std::vector<double>( basisSize, 0.0 );
_YBasis = Vector( basisSize, 0.0 );
}
std::vector<double> SphericalHarmonics::ComputeSphericalBasis( double my, double phi ) {
Vector SphericalHarmonics::ComputeSphericalBasis( double my, double phi ) {
ComputeAssLegendrePoly( my );
ComputeYBasis( phi );
return _YBasis;
}
std::vector<double> SphericalHarmonics::ComputeSphericalBasis( double x, double y, double z ) {
Vector SphericalHarmonics::ComputeSphericalBasis( double x, double y, double z ) {
// transform (x,y,z) into (my,phi)
double my = z;
......
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