Commit 732856ab authored by Steffen Schotthöfer's avatar Steffen Schotthöfer
Browse files

added new bounadary epsilon for datagenerator. parallelized parts of mn solver...

added new bounadary epsilon for datagenerator. parallelized parts of mn solver (openMP). Removed compiler warnings. Removed bug in datagenerator.
parent 765aa5a4
Pipeline #135419 passed with stage
in 17 minutes and 44 seconds
......@@ -119,9 +119,10 @@ 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 */
double _boundaryDistanceRealizableSet; /*! @brief Distance of the sampled moments to the boundary of the realizable set */
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 _RealizableSetEpsilonU0; /*! @brief Distance to 0 of the sampled moments to the boundary of the realizable set */
double _RealizableSetEpsilonU1; /*!< @brief: norm(u_1)/u_0 !< _RealizableSetEpsilonU1 */
// --- Parsing Functionality and Initializing of Options ---
/*!
......@@ -315,7 +316,8 @@ class Config
bool inline GetDataGeneratorMode() { return _dataGeneratorMode; }
unsigned long inline GetTrainingDataSetSize() { return _tainingSetSize; }
unsigned long inline GetMaxValFirstMoment() { return _maxValFirstMoment; }
double GetBoundaryDistanceRealizableSet() { return _boundaryDistanceRealizableSet; }
double GetRealizableSetEpsilonU0() { return _RealizableSetEpsilonU0; }
double GetRealizableSetEpsilonU1() { return _RealizableSetEpsilonU1; }
// ---- Setters for option structure
// This section is dangerous
......
......@@ -13,3 +13,12 @@ PROBLEM = CHECKERBOARD
BC_NEUMANN = ( void )
QUAD_TYPE = GAUSS_LEGENDRE_TENSORIZED
QUAD_ORDER = 10
% ----- Output ----
%
VOLUME_OUTPUT = (MINIMAL)
VOLUME_OUTPUT_FREQUENCY = 1
SCREEN_OUTPUT = (ITER, MASS, RMS_FLUX, VTK_OUTPUT, CSV_OUTPUT)
SCREEN_OUTPUT_FREQUENCY = 1
HISTORY_OUTPUT = (ITER, MASS, RMS_FLUX, VTK_OUTPUT, CSV_OUTPUT)
HISTORY_OUTPUT_FREQUENCY = 1
......@@ -7,9 +7,10 @@
%
% ---- Global settings ----
DATA_GENERATOR_MODE = YES
TRAINING_SET_SIZE = 10
MAX_VALUE_FIRST_MOMENT = 500
BOUNDARY_DISTANCE_REALIZABLE_SET = 0.01
TRAINING_SET_SIZE = 1000
MAX_VALUE_FIRST_MOMENT = 20
REALIZABLE_SET_EPSILON_U0 = 0.05
REALIZABLE_SET_EPSILON_U1 = 0.97
%
% ---- File specifications ----
%
......
......@@ -14,8 +14,8 @@ OUTPUT_FILE = exampleMN_mono
% Log directory
LOG_DIR = ../result/logs
% Mesh File
%MESH_FILE = linesource.su2
MESH_FILE = linesource_debug.su2
MESH_FILE = linesource.su2
%MESH_FILE = linesource_debug.su2
%
% ---- Problem description ---
%
......@@ -33,7 +33,7 @@ SOLVER = MN_SOLVER
% CFL number
CFL_NUMBER = 0.7
% Final time for simulation
TIME_FINAL = 0.3
TIME_FINAL = 0.5
% Reconstruction order (spatial flux)
RECONS_ORDER = 1
......@@ -41,7 +41,7 @@ RECONS_ORDER = 1
% ---- Entropy settings ----
%
ENTROPY_FUNCTIONAL = MAXWELL_BOLTZMANN
ENTROPY_OPTIMIZER = ML
ENTROPY_OPTIMIZER = NEWTON
NEURAL_MODEL = 4
%
% ----- Newton Solver Specifications ----
......
......@@ -320,9 +320,12 @@ 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 );
/*! @brief Distance to 0 of the sampled moments to the boundary of the realizable set \n DESCRIPTION: Distance to the boundary of the
* realizable set \n DEFAULT 0.1 \ingroup Config */
AddDoubleOption( "REALIZABLE_SET_EPSILON_U0", _RealizableSetEpsilonU0, 0.1 );
/*! @brief norm(u_1)/u_0 is enforced to be smaller than _RealizableSetEpsilonU1 \n DESCRIPTION: Distance to the boundary of the realizable set \n
* DEFAULT 0.1 \ingroup Config */
AddDoubleOption( "REALIZABLE_SET_EPSILON_U1", _RealizableSetEpsilonU1, 0.9 );
}
void Config::SetConfigParsing( string case_filename ) {
......
......@@ -34,4 +34,6 @@ VectorVector IsotropicSource2D::SetupIC() {
return psi;
}
std::vector<double> IsotropicSource2D::GetDensity( const VectorVector& cellMidPoints ) { return std::vector<double>( _settings->GetNCells(), 1.0 ); }
std::vector<double> IsotropicSource2D::GetDensity( const VectorVector& /*cellMidPoints*/ ) {
return std::vector<double>( _settings->GetNCells(), 1.0 );
}
......@@ -171,7 +171,7 @@ void CSDSolverTrafoFP::FluxUpdate() {
}
}
void CSDSolverTrafoFP::FVMUpdate( unsigned idx_energy ) {
void CSDSolverTrafoFP::FVMUpdate( unsigned /*idx_energy*/ ) {
#pragma omp parallel for
for( unsigned j = 0; j < _nCells; ++j ) {
if( _boundaryCells[j] == BOUNDARY_TYPE::DIRICHLET ) continue;
......
......@@ -275,7 +275,7 @@ void CSDSolverTrafoFP2D::WriteVolumeOutput( unsigned idx_pseudoTime ) {
}
// Solver
void CSDSolverTrafoFP2D::FVMUpdate( unsigned idx_energy ) {
void CSDSolverTrafoFP2D::FVMUpdate( unsigned /*idx_energy*/ ) {
// loop over all spatial cells
#pragma omp parallel for
for( unsigned j = 0; j < _nCells; ++j ) {
......@@ -394,15 +394,6 @@ void CSDSolverTrafoFP2D::SolverPreprocessing() {
for( unsigned k = 0; k < _nq; ++k ) _identity( k, k ) = 1.0;
// angular flux at next time step
VectorVector psiNew( _nCells, Vector( _nq, 0.0 ) );
double dFlux = 1e10;
Vector fluxNew( _nCells, 0.0 );
Vector fluxOld( _nCells, 0.0 );
// for( unsigned j = 0; j < _nCells; ++j ) {
// fluxOld[j] = dot( _sol[j], _weights );
//}
int rank;
MPI_Comm_rank( MPI_COMM_WORLD, &rank );
// if( rank == 0 ) log->info( "{:10} {:10}", "E", "dFlux" );
......
......@@ -13,8 +13,8 @@ CSDSolverTrafoFPSH2D::CSDSolverTrafoFPSH2D( Config* settings ) : SNSolver( setti
GenerateEnergyGrid( false );
// create quadrature
unsigned order = _quadrature->GetOrder();
unsigned nq = _settings->GetNQuadPoints();
unsigned order = _quadrature->GetOrder();
// unsigned nq = _settings->GetNQuadPoints();
_quadPoints = _quadrature->GetPoints();
_weights = _quadrature->GetWeights();
_quadPointsSphere = _quadrature->GetPointsSphere();
......@@ -107,7 +107,7 @@ CSDSolverTrafoFPSH2D::CSDSolverTrafoFPSH2D( Config* settings ) : SNSolver( setti
_S = Matrix( nSph, nSph, 0.0 );
unsigned counter = 0;
for( int l = 0; l <= orderSph; ++l ) {
for( int l = 0; l <= (int)orderSph; ++l ) {
for( int m = -l; m <= l; ++m ) {
_S( counter, counter ) = double( -l * ( l + 1 ) );
counter++;
......@@ -385,7 +385,7 @@ void CSDSolverTrafoFPSH2D::WriteVolumeOutput( unsigned idx_pseudoTime ) {
}
// Solver
void CSDSolverTrafoFPSH2D::FVMUpdate( unsigned idx_energy ) {
void CSDSolverTrafoFPSH2D::FVMUpdate( unsigned /*idx_energy*/ ) {
// loop over all spatial cells
#pragma omp parallel for
for( unsigned j = 0; j < _nCells; ++j ) {
......@@ -502,15 +502,6 @@ void CSDSolverTrafoFPSH2D::SolverPreprocessing() {
_identity = Matrix( _nq, _nq, 0.0 );
for( unsigned k = 0; k < _nq; ++k ) _identity( k, k ) = 1.0;
// angular flux at next time step
VectorVector psiNew( _nCells, Vector( _nq, 0.0 ) );
double dFlux = 1e10;
Vector fluxNew( _nCells, 0.0 );
Vector fluxOld( _nCells, 0.0 );
// for( unsigned j = 0; j < _nCells; ++j ) {
// fluxOld[j] = dot( _sol[j], _weights );
//}
int rank;
MPI_Comm_rank( MPI_COMM_WORLD, &rank );
if( rank == 0 ) log->info( "{:10} {:10}", "E", "dFlux" );
......@@ -546,5 +537,4 @@ void CSDSolverTrafoFPSH2D::SolverPreprocessing() {
// std::cout << "nEnergies = " << _nEnergies << std::endl;
// loop over energies (pseudo-time)
}
......@@ -15,7 +15,7 @@
#include "spdlog/spdlog.h"
#include <mpi.h>
#include <fstream>
#include <iostream>
MNSolver::MNSolver( Config* settings ) : Solver( settings ) {
......@@ -126,7 +126,7 @@ Vector MNSolver::ConstructFlux( unsigned idx_cell ) {
}
// Solution flux
flux += _moments[idx_quad] * ( _weights[idx_quad] * entropyFlux );
}
}
return flux;
}
......@@ -140,7 +140,7 @@ void MNSolver::ComputeRealizableSolution( unsigned idx_cell ) {
}
}
void MNSolver::IterPreprocessing( unsigned idx_pseudotime ) {
void MNSolver::IterPreprocessing( unsigned /*idx_pseudotime*/ ) {
// ------- Reconstruction Step ----------------
......@@ -148,7 +148,7 @@ void MNSolver::IterPreprocessing( unsigned idx_pseudotime ) {
// ------- Relizablity Preservation Step ----
for( unsigned idx_cell = 0; idx_cell < _nCells; idx_cell++ ) {
ComputeRealizableSolution( idx_cell );
ComputeRealizableSolution( idx_cell ); // already parallelized
}
}
......@@ -162,6 +162,7 @@ void MNSolver::IterPostprocessing() {
void MNSolver::ComputeRadFlux() {
double firstMomentScaleFactor = sqrt( 4 * M_PI );
#pragma omp parallel for
for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
_fluxNew[idx_cell] = _sol[idx_cell][0] * firstMomentScaleFactor;
}
......@@ -172,29 +173,32 @@ void MNSolver::FluxUpdate() {
_mesh->ReconstructSlopesU( _nTotalEntries, _solDx, _solDy, _alpha ); // unstructured reconstruction
}
// Loop over the grid cells
// Loop over the grid cells
#pragma omp parallel for
for( unsigned idx_cell = 0; idx_cell < _nCells; idx_cell++ ) {
// Dirichlet Boundaries stay
// Dirichlet Boundaries stayd
if( _boundaryCells[idx_cell] == BOUNDARY_TYPE::DIRICHLET ) continue;
_solNew[idx_cell] = ConstructFlux( idx_cell );
// if( norm( _solNew[idx_cell] ) > 0.001 ) {
// std::cout << _solNew[idx_cell] << "\n";
//}
}
}
void MNSolver::FVMUpdate( unsigned idx_energy ) {
// Loop over the grid cells
//#pragma omp parallel for
for( unsigned idx_cell = 0; idx_cell < _nCells; idx_cell++ ) {
// Dirichlet Boundaries stay
if( _boundaryCells[idx_cell] == BOUNDARY_TYPE::DIRICHLET ) continue;
for( unsigned idx_system = 0; idx_system < _nTotalEntries; idx_system++ ) {
_solNew[idx_cell][idx_system] = _sol[idx_cell][idx_system] -
( _dE / _areas[idx_cell] ) * _solNew[idx_cell][idx_system] /* cell averaged flux */
- _dE * _sol[idx_cell][idx_system] *
( _sigmaT[idx_energy][idx_cell] /* absorbtion influence */
+ _sigmaS[idx_energy][idx_cell] * _scatterMatDiag[idx_system] ); /* scattering influence */
// Flux update
for( unsigned idx_sys = 0; idx_sys < _nTotalEntries; idx_sys++ ) {
_solNew[idx_cell][idx_sys] = _sol[idx_cell][idx_sys] - ( _dE / _areas[idx_cell] ) * _solNew[idx_cell][idx_sys] /* cell averaged flux */
- _dE * _sol[idx_cell][idx_sys] *
( _sigmaT[idx_energy][idx_cell] /* absorbtion influence */
+ _sigmaS[idx_energy][idx_cell] * _scatterMatDiag[idx_sys] ); /* scattering influence */
}
// Source Term
_solNew[idx_cell][0] += _dE * _Q[0][idx_cell][0];
}
}
......
......@@ -53,7 +53,7 @@ PNSolver::PNSolver( Config* settings ) : Solver( settings ) {
// TODO
}
void PNSolver::IterPreprocessing( unsigned idx_pseudotime ) {
void PNSolver::IterPreprocessing( unsigned /*idx_pseudotime*/ ) {
// Nothing to preprocess for PNSolver
}
......@@ -150,6 +150,10 @@ void PNSolver::FluxUpdate() {
}
}
}
// std::cout << norm( _solNew[idx_cell] ) << "\n";
// if( norm( _solNew[idx_cell] ) > 0.001 ) {
// std::cout << _solNew[idx_cell] << "\n";
//}
}
}
......
......@@ -23,7 +23,7 @@ SNSolver::SNSolver( Config* settings ) : Solver( settings ) {
delete k;
}
void SNSolver::IterPreprocessing( unsigned idx_pseudotime ) {
void SNSolver::IterPreprocessing( unsigned /*idx_pseudotime*/ ) {
// Nothing to do for SNSolver
}
......
......@@ -111,6 +111,8 @@ Solver* Solver::Create( Config* settings ) {
case CSD_SN_FOKKERPLANCK_TRAFO_SH_SOLVER_2D: return new CSDSolverTrafoFPSH2D( settings );
default: ErrorMessages::Error( "Creator for the chosen solver does not yet exist. This is is the fault of the coder!", CURRENT_FUNCTION );
}
ErrorMessages::Error( "Creator for the chosen solver does not yet exist. This is is the fault of the coder!", CURRENT_FUNCTION );
return nullptr; // This code is never reached. Just to disable compiler warnings.
}
void Solver::Solve() {
......
......@@ -66,12 +66,12 @@ nnDataGenerator::nnDataGenerator( Config* settings ) {
else if( _LMaxDegree == 2 && _settings->GetSphericalBasisName() == SPHERICAL_MONOMIALS && _settings->GetDim() == 1 ) {
// Carefull: This computes only normalized moments, i.e. sampling for u_0 = 1
unsigned c = 1;
double N1 = -1.0 + _settings->GetBoundaryDistanceRealizableSet();
double N1 = -1.0 + _settings->GetRealizableSetEpsilonU0();
double N2;
double dN = 2.0 / (double)_gridSize;
while( N1 < 1.0 - _settings->GetBoundaryDistanceRealizableSet() ) {
N2 = N1 * N1 + _settings->GetBoundaryDistanceRealizableSet();
while( N2 < 1.0 - _settings->GetBoundaryDistanceRealizableSet() ) {
while( N1 < 1.0 - _settings->GetRealizableSetEpsilonU0() ) {
N2 = N1 * N1 + _settings->GetRealizableSetEpsilonU0();
while( N2 < 1.0 - _settings->GetRealizableSetEpsilonU0() ) {
c++;
N2 += dN;
}
......@@ -153,23 +153,36 @@ void nnDataGenerator::SampleSolutionU() {
// --- sample u in order 0 ---
// u_0 = <1*psi>
#pragma omp parallel for schedule( guided )
//#pragma omp parallel for schedule( guided )
for( unsigned long idx_set = 0; idx_set < _gridSize; idx_set++ ) {
Vector u1 = { 0, 0, 0 };
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 ) * ( 1 + 2 * _settings->GetBoundaryDistanceRealizableSet() );
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
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++ ) {
......@@ -187,12 +200,12 @@ void nnDataGenerator::SampleSolutionU() {
else if( _LMaxDegree == 2 && _settings->GetSphericalBasisName() == SPHERICAL_MONOMIALS && _settings->GetDim() == 1 ) {
// Carefull: This computes only normalized moments, i.e. sampling for u_0 = 1
unsigned c = 0;
double N1 = -1.0 + _settings->GetBoundaryDistanceRealizableSet();
double N1 = -1.0 + _settings->GetRealizableSetEpsilonU0();
double N2;
double dN = 2.0 / (double)_gridSize;
while( N1 < 1.0 - _settings->GetBoundaryDistanceRealizableSet() ) {
N2 = N1 * N1 + _settings->GetBoundaryDistanceRealizableSet();
while( N2 < 1.0 - _settings->GetBoundaryDistanceRealizableSet() ) {
while( N1 < 1.0 - _settings->GetRealizableSetEpsilonU0() ) {
N2 = N1 * N1 + _settings->GetRealizableSetEpsilonU0();
while( N2 < 1.0 - _settings->GetRealizableSetEpsilonU0() ) {
_uSol[c][0] = 1; // u0 (normalized i.e. N0) by Monreals notation
_uSol[c][1] = N1; // u1 (normalized i.e. N1) by Monreals notation
_uSol[c][2] = N2; // u2 (normalized i.e. N2) by Monreals notation
......@@ -203,7 +216,7 @@ void nnDataGenerator::SampleSolutionU() {
}
}
else {
ErrorMessages::Error( "Sampling for order higher than 1 is not yet supported", CURRENT_FUNCTION );
ErrorMessages::Error( "Sampling for this configuration is not yet supported", CURRENT_FUNCTION );
}
}
......@@ -253,7 +266,7 @@ void nnDataGenerator::PrintTrainingData() {
}
void nnDataGenerator::CheckRealizability() {
double epsilon = _settings->GetBoundaryDistanceRealizableSet();
double epsilon = _settings->GetRealizableSetEpsilonU0();
if( _LMaxDegree == 1 ) {
double normU1 = 0.0;
Vector u1( 3, 0.0 );
......@@ -269,9 +282,9 @@ void nnDataGenerator::CheckRealizability() {
else {
u1 = { _uSol[idx_set][1], _uSol[idx_set][2], _uSol[idx_set][3] };
normU1 = norm( u1 );
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";
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 );
......
......@@ -9,8 +9,8 @@
DATA_GENERATOR_MODE = YES
TRAINING_SET_SIZE = 2
MAX_VALUE_FIRST_MOMENT = 20
BOUNDARY_DISTANCE_REALIZABLE_SET = 0.01
REALIZABLE_SET_EPSILON_U0 = 0.01
REALIZABLE_SET_EPSILON_U1 = 0.97
%
% ---- File specifications ----
%
......
......@@ -189,6 +189,7 @@ TEST_CASE( "MN_SOLVER", "[validation_tests]" ) {
}
}
}
/*
TEST_CASE( "CSD_SN_FP_SOLVER", "[validation_tests]" ) {
std::string csd_sn_fileDir = "input/validation_tests/CSD_SN_FP_solver/";
......
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