Commit 4aa14039 authored by steffen.schotthoefer's avatar steffen.schotthoefer Committed by jannick.wolters
Browse files

now using spdlog everywhere ( closes #35 ); replaced hard exit calls without...

now using spdlog everywhere ( closes #35 ); replaced hard exit calls without message by ErrorMessage class ( updates #31 )
parent 4d332829
......@@ -6,14 +6,14 @@
#include <omp.h>
#include <vector>
#include "blaze/Blaze.h"
#include "blaze/math/CompressedMatrix.h"
#include "metis.h"
#include "parmetis.h"
#include "spdlog/spdlog.h"
#include "typedef.h"
#include "settings/globalconstants.h"
#include "toolboxes/errormessages.h"
#include "typedef.h"
class Mesh
{
......@@ -52,18 +52,57 @@ class Mesh
std::vector<std::pair<BOUNDARY_TYPE, std::vector<unsigned>>> boundaries );
~Mesh();
unsigned GetDim() const;
unsigned GetNumCells() const;
unsigned GetNumNodes() const;
unsigned GetNumNodesPerCell() const;
inline unsigned GetDim() const { return _dim; }
inline unsigned GetNumCells() const { return _numCells; }
inline unsigned GetNumNodes() const { return _numNodes; }
inline unsigned GetNumNodesPerCell() const { return _numNodesPerCell; }
/**
* @brief Returns all node coordinates
* @return dimension: numNodes x dim
*/
const std::vector<Vector>& GetNodes() const;
/**
* @brief Returns the mid point coordinates of each cell
* @return dimension: numCells x dim
*/
const std::vector<Vector>& GetCellMidPoints() const;
/**
* @brief Returns all node IDs that construct up each cell
* @return dimension: numCells x numNodes
*/
const std::vector<std::vector<unsigned>>& GetCells() const;
/**
* @brief Returns the cell area of each cell
* @return dimension: numCells
*/
const std::vector<double>& GetCellAreas() const;
/**
* @brief Return the color/ID of the mesh partition
* @return dimension: numCells
*/
const std::vector<unsigned>& GetPartitionIDs() const;
/**
* @brief Returns the neighbor cell IDs for every cell
* @return dimension: numCells x numNodes
*/
const std::vector<std::vector<unsigned>>& GetNeighbours() const;
/**
* @brief Returns the edge length scaled normal vectors of each cell
* @return dimension: numCells x numNodes x dim
*/
const std::vector<std::vector<Vector>>& GetNormals() const;
/**
* @brief Returns the boundary enum for each cell. BOUNDARY_TYPE::NONE is the default.
* @return dimension: numCells
*/
const std::vector<BOUNDARY_TYPE>& GetBoundaryTypes() const;
/**
......
......@@ -5,6 +5,10 @@
class QGaussLegendreTensorized : public QuadratureBase
{
private:
double Pythag( const double a, const double b );
std::pair<Vector, Matrix> ComputeEigenValTriDiagMatrix( const Matrix& mat );
public:
QGaussLegendreTensorized( unsigned order );
virtual ~QGaussLegendreTensorized() {}
......
......@@ -2,6 +2,7 @@
#define QUADRATURE_H
#include "settings/globalconstants.h"
#include "toolboxes/errormessages.h"
#include "typedef.h"
#include <iostream>
#include <string>
......
......@@ -21,90 +21,90 @@
#include "optionstructure.h"
/*!
* \class Config
* \brief Main class for defining the problem; basically this class reads the configuration file, and
* @class Config
* @brief Main class for defining the problem; basically this class reads the configuration file, and
* stores all the information.
*/
class Config
{
private:
std::string _fileName; /*!< \brief Name of the current file without extension */
std::string _fileName; /*!< @rief Name of the current file without extension */
bool _baseConfig;
int _commRank, _commSize; /*!< \brief MPI rank and size.*/ // Not yet used!!
int _commRank, _commSize; /*!< @brief MPI rank and size.*/ // Not yet used!!
// --- Options ---
// File Structure
std::string _inputDir; /*!< \brief Directory for input files*/
std::string _outputDir; /*!< \brief Directory for output files*/
std::string _outputFile; /*!< \brief Name of output file*/
std::string _logDir; /*!< \brief Directory of log file*/
std::string _meshFile; /*!< \brief Name of mesh file*/
std::string _inputDir; /*!< @brief Directory for input files*/
std::string _outputDir; /*!< @brief Directory for output files*/
std::string _outputFile; /*!< @brief Name of output file*/
std::string _logDir; /*!< @brief Directory of log file*/
std::string _meshFile; /*!< @brief Name of mesh file*/
// Quadrature
QUAD_NAME _quadName; /*!< \brief Quadrature Name*/
unsigned short _quadOrder; /*!< \brief Quadrature Order*/
QUAD_NAME _quadName; /*!< @brief Quadrature Name*/
unsigned short _quadOrder; /*!< @brief Quadrature Order*/
unsigned _nQuadPoints;
// Mesh
unsigned _nCells;
// Solver
double _CFL; /*!< \brief CFL Number for Solver*/
double _tEnd; /*!< \brief Final Time for Simulation */
double _CFL; /*!< @brief CFL Number for Solver*/
double _tEnd; /*!< @brief Final Time for Simulation */
PROBLEM_NAME _problemName;
// Boundary Conditions
/*!< \brief List of all Pairs (marker, BOUNDARY_TYPE), e.g. (farfield,DIRICHLET).
/*!< @brief List of all Pairs (marker, BOUNDARY_TYPE), e.g. (farfield,DIRICHLET).
Each Boundary Conditions must have an entry in enum BOUNDARY_TYPE*/
std::vector<std::pair<std::string, BOUNDARY_TYPE>> _boundaries;
unsigned short _nMarkerDirichlet; /*!< \brief Number of Dirichlet BC markers. Enum entry: DIRICHLET */
unsigned short _nMarkerNeumann; /*!< \brief Number of Neumann BC markers. Enum entry: Neumann */
std::vector<std::string> _MarkerDirichlet; /*!< \brief Dirichlet BC markers. */
std::vector<std::string> _MarkerNeumann; /*!< \brief Neumann BC markers. */
unsigned short _nMarkerDirichlet; /*!< @brief Number of Dirichlet BC markers. Enum entry: DIRICHLET */
unsigned short _nMarkerNeumann; /*!< @brief Number of Neumann BC markers. Enum entry: Neumann */
std::vector<std::string> _MarkerDirichlet; /*!< @brief Dirichlet BC markers. */
std::vector<std::string> _MarkerNeumann; /*!< @brief Neumann BC markers. */
// Scattering Kernel
KERNEL_NAME _kernelName; /*!< \brief Scattering Kernel Name*/
KERNEL_NAME _kernelName; /*!< @brief Scattering Kernel Name*/
// --- Parsing Functionality and Initializing of Options ---
/*!
* \brief Set default values for all options not yet set.
* @brief Set default values for all options not yet set.
*/
void SetDefault( void );
/*!
* \brief Set the config options.
* @brief Set the config options.
* ==> Set new config options here.
*/
void SetConfigOptions( void );
/*!
* \brief Set the config file parsing.
* @brief Set the config file parsing.
*/
void SetConfigParsing( char case_filename[MAX_STRING_SIZE] );
/*!
* \brief Config file screen output.
* @brief Config file screen output.
*/
void SetOutput( void );
/*!
* \brief Initializes pointers to null
* @brief Initializes pointers to null
*/
void SetPointersNull( void );
/*!
* \brief Config file postprocessing.
* @brief Config file postprocessing.
*/
void SetPostprocessing( void );
/*!
* \brief breaks an input line from the config file into a set of tokens
* \param[in] str - the input line string
* \param[out] option_name - the name of the option found at the beginning of the line
* \param[out] option_value - the tokens found after the "=" sign on the line
* \return false if the line is empty or a commment, true otherwise
* @brief breaks an input line from the config file into a set of tokens
* @param[in] str - the input line string
* @param[out] option_name - the name of the option found at the beginning of the line
* @param[out] option_value - the tokens found after the "=" sign on the line
* @return false if the line is empty or a commment, true otherwise
*/
bool TokenizeString( std::string& str, std::string& option_name, std::vector<std::string>& option_value );
......@@ -134,7 +134,7 @@ class Config
// List and Array options should also be able to be specified with the string "NONE" indicating that there
// are no elements. This allows the option to be present in a config file but left blank.
/*!< \brief addDoubleOption creates a config file parser for an option with the given name whose
/*!< @brief addDoubleOption creates a config file parser for an option with the given name whose
value can be represented by a su2double.*/
// Simple Options
......@@ -165,19 +165,19 @@ class Config
public:
/*!
* \brief Constructor of the class which reads the input file.
* @brief Constructor of the class which reads the input file.
*/
Config( char case_filename[MAX_STRING_SIZE] );
/*!
* \brief Destructor of the class.
* @brief Destructor of the class.
*/
~Config( void );
// ---- Getters for option values ----
/*!
* \brief Get Value of this option.
* @brief Get Value of this option.
* Please keep alphabetical order within each subcategory
*/
// File structure
......@@ -202,7 +202,7 @@ class Config
PROBLEM_NAME inline GetProblemName() const { return _problemName; }
// Boundary Conditions
BOUNDARY_TYPE GetBoundaryType( std::string nameMarker ) const; /*! \brief Get Boundary Type of given marker */
BOUNDARY_TYPE GetBoundaryType( std::string nameMarker ) const; /*! @brief Get Boundary Type of given marker */
// Scattering Kernel
KERNEL_NAME inline GetKernelName() const { return _kernelName; }
......
......@@ -15,8 +15,8 @@
// --- Definition for global constants goes here ---
const double PI_NUMBER = M_PI; /*!< \brief Pi number. */
const unsigned int MAX_STRING_SIZE = 200; /*!< \brief Maximum size for strings. */
const double PI_NUMBER = M_PI; /*!< @brief Pi number. */
const unsigned int MAX_STRING_SIZE = 200; /*!< @brief Maximum size for strings. */
// --- Definition of enums goes here ---
......@@ -24,12 +24,12 @@ enum BOUNDARY_TYPE { DIRICHLET, NEUMANN, NONE, INVALID };
// --- Definition of enums for EnumOptions goes here ---
/*! \brief Enum for all currently available quadratures in rtsn.
/*! @brief Enum for all currently available quadratures in rtsn.
* Option enums are written in capital letters with underscores as spaces (e.g option "time integration" has option enum "TIME_INTEGRATION")
*/
enum QUAD_NAME { QUAD_MonteCarlo, QUAD_GaussLegendreTensorized, QUAD_LevelSymmetric, QUAD_Lebedev, QUAD_LDFESA };
/*! \brief Conversion Map String to enum
/*! @brief Conversion Map String to enum
*/
inline std::map<std::string, QUAD_NAME> Quadrature_Map{ { "MONTE_CARLO", QUAD_MonteCarlo },
{ "GAUSS_LEGENDRE_TENSORIZED", QUAD_GaussLegendreTensorized },
......
#ifndef TYPEDEFS_H
#define TYPEDEFS_H
#include <blaze/math/CompressedMatrix.h>
#include <blaze/math/DynamicMatrix.h>
#include <blaze/math/DynamicVector.h>
......
......@@ -11,5 +11,5 @@ PROBLEM = CHECKERBOARD
% ---- Boundary Conditions ----
BC_NEUMANN = ( void )
QUAD_TYPE = LEBEDEV
QUAD_ORDER = 35
QUAD_TYPE = GAUSS_LEGENDRE_TENSORIZED
QUAD_ORDER = 10
......@@ -34,7 +34,7 @@ void ExportVTK( const std::string fileName,
pts->SetPoint( nodeID++, node[0], node[1], node[2] );
}
else {
exit( EXIT_FAILURE );
ErrorMessages::Error( "Unsupported dimension (d=" + std::to_string( dim ) + ")!", CURRENT_FUNCTION );
}
}
vtkCellArraySP cellArray = vtkCellArraySP::New();
......@@ -72,8 +72,8 @@ void ExportVTK( const std::string fileName,
break;
default:
auto log = spdlog::get( "event" );
log->error( "[ERROR][IO::ExportVTK] Invalid dimension" );
exit( EXIT_FAILURE );
ErrorMessages::Error( "Please implement output for results of size " + std::to_string( results[i].size() ) + "!",
CURRENT_FUNCTION );
}
grid->GetCellData()->AddArray( cellData );
}
......@@ -106,7 +106,8 @@ Mesh* LoadSU2MeshFromFile( const Config* settings ) {
std::vector<std::vector<unsigned>> cells;
std::vector<std::pair<BOUNDARY_TYPE, std::vector<unsigned>>> boundaries;
if( !std::filesystem::exists( settings->GetMeshFile() ) ) exit( EXIT_FAILURE );
if( !std::filesystem::exists( settings->GetMeshFile() ) )
ErrorMessages::Error( "Cannot find mesh file '" + settings->GetMeshFile() + "!", CURRENT_FUNCTION );
std::ifstream ifs( settings->GetMeshFile(), std::ios::in );
std::string line;
......@@ -181,7 +182,7 @@ Mesh* LoadSU2MeshFromFile( const Config* settings ) {
}
}
else {
exit( EXIT_FAILURE );
ErrorMessages::Error( "Invalid mesh file detected! Make sure boundaries are provided.'", CURRENT_FUNCTION );
}
}
std::sort( bnodes.begin(), bnodes.end() );
......@@ -215,8 +216,7 @@ Mesh* LoadSU2MeshFromFile( const Config* settings ) {
break;
}
default: {
log->error( "[LoadSU2MeshFromFile] Unsupported mesh type!" );
exit( EXIT_FAILURE );
ErrorMessages::Error( "Unsupported mesh type!'", CURRENT_FUNCTION );
}
}
}
......@@ -225,8 +225,7 @@ Mesh* LoadSU2MeshFromFile( const Config* settings ) {
}
bool mixedElementMesh = !std::equal( numNodesPerCell.begin() + 1, numNodesPerCell.end(), numNodesPerCell.begin() );
if( mixedElementMesh ) {
log->error( "[LoadSU2MeshFromFile] Mixed element meshes are currently not supported!" );
exit( EXIT_FAILURE );
ErrorMessages::Error( "Mixed element meshes are currently not supported!'", CURRENT_FUNCTION );
}
ifs.clear();
ifs.seekg( 0, std::ios::beg );
......@@ -254,8 +253,7 @@ Mesh* LoadSU2MeshFromFile( const Config* settings ) {
}
}
else {
log->error( "[LoadSU2MeshFromFile] File not found" );
exit( EXIT_FAILURE );
ErrorMessages::Error( "Cannot open mesh file '" + settings->GetMeshFile() + "!", CURRENT_FUNCTION );
}
ifs.close();
return new Mesh( nodes, cells, boundaries );
......@@ -281,8 +279,7 @@ std::string ParseArguments( int argc, char* argv[] ) {
inputFile = std::string( argv[i] );
std::ifstream f( inputFile );
if( !f.is_open() ) {
std::cerr << "[ERROR] Unable to open inputfile '" << inputFile << "'!" << std::endl;
exit( EXIT_FAILURE );
ErrorMessages::OptionNotSetError( "Unable to open inputfile '" + inputFile + "' !", CURRENT_FUNCTION );
}
}
}
......
......@@ -6,7 +6,6 @@
#include "settings/config.h"
int main( int argc, char** argv ) {
MPI_Init( &argc, &argv );
std::string filename = "default.cfg";
......
......@@ -82,10 +82,9 @@ void Mesh::ComputeConnectivity() {
if( std::any_of( neighborsFlat.begin(), neighborsFlat.end(), []( int i ) { return i == -1; } ) ) {
for( unsigned idx = 0; idx < neighborsFlat.size(); ++idx ) {
if( neighborsFlat[idx] == -1 ) _log->error( "[mesh] Detected unassigned faces at index {0}", idx );
break;
if( neighborsFlat[idx] == -1 )
ErrorMessages::Error( "Detected unassigned faces at index " + std::to_string( idx ) + " !", CURRENT_FUNCTION );
}
exit( EXIT_FAILURE );
}
_cellNeighbors.resize( _numCells );
......@@ -153,7 +152,8 @@ void Mesh::ComputeCellAreas() {
break;
}
default: {
exit( EXIT_FAILURE );
ErrorMessages::Error( "Area computation for cells with " + std::to_string( _numNodesPerCell ) + " nodes is not implemented yet!",
CURRENT_FUNCTION );
}
}
}
......@@ -320,11 +320,6 @@ void Mesh::ComputeSlopes( unsigned nq, VectorVector& psiDerX, VectorVector& psiD
}
}
unsigned Mesh::GetDim() const { return _dim; }
unsigned Mesh::GetNumCells() const { return _numCells; }
unsigned Mesh::GetNumNodes() const { return _numNodes; }
unsigned Mesh::GetNumNodesPerCell() const { return _numNodesPerCell; }
const std::vector<Vector>& Mesh::GetNodes() const { return _nodes; }
const std::vector<Vector>& Mesh::GetCellMidPoints() const { return _cellMidPoints; }
const std::vector<std::vector<unsigned>>& Mesh::GetCells() const { return _cells; }
......
......@@ -7,13 +7,60 @@ QGaussLegendreTensorized::QGaussLegendreTensorized( unsigned order ) : Quadratur
SetConnectivity();
}
void QGaussLegendreTensorized::SetPointsAndWeights() { // TODO
// Compute Points
// Nq random points on the sphere.
_points = VectorVector( GetNq() );
void QGaussLegendreTensorized::SetPointsAndWeights() {
Vector nodes( _order ), weights( _order );
// Compute Weights
_weights = Vector( GetNq(), 4.0 * M_PI / GetNq() );
// construct companion matrix
Matrix CM( _order, _order, 0.0 );
for( unsigned i = 0; i < _order - 1; ++i ) {
CM( i + 1, i ) = std::sqrt( 1 / ( 4 - 1 / std::pow( static_cast<double>( i + 1 ), 2 ) ) );
CM( i, i + 1 ) = std::sqrt( 1 / ( 4 - 1 / std::pow( static_cast<double>( i + 1 ), 2 ) ) );
}
auto evSys = ComputeEigenValTriDiagMatrix( CM );
for( unsigned i = 0; i < _order; ++i ) {
if( std::fabs( evSys.first[i] ) < 1e-15 )
nodes[i] = 0;
else
nodes[i] = evSys.first[i];
weights[i] = 2 * std::pow( evSys.second( 0, i ), 2 );
}
for( unsigned i = 0; i < _order; ++i ) {
nodes[i] = ( nodes[i] + 1.0 ) * 0.5;
}
std::vector<unsigned> p( nodes.size() );
std::iota( p.begin(), p.end(), 0 );
std::sort( p.begin(), p.end(), [&]( unsigned i, unsigned j ) { return nodes[i] < nodes[j]; } );
Vector sorted_nodes( static_cast<unsigned>( p.size() ) ), sorted_weights( static_cast<unsigned>( p.size() ) );
std::transform( p.begin(), p.end(), sorted_nodes.begin(), [&]( unsigned i ) { return nodes[i]; } );
std::transform( p.begin(), p.end(), sorted_weights.begin(), [&]( unsigned i ) { return weights[i]; } );
nodes = sorted_nodes;
weights = sorted_weights;
Vector phi( 2 * _order );
for( unsigned i = 0; i < 2 * _order; ++i ) {
phi[i] = ( i + 0.5 ) * M_PI / _order;
}
unsigned range = std::floor( _order / 2.0 );
_points.resize( _nq );
for( auto& p : _points ) {
p.resize( 3 );
}
_weights.resize( _nq );
for( unsigned j = 0; j < range; ++j ) {
for( unsigned i = 0; i < 2 * _order; ++i ) {
_points[j * ( 2 * _order ) + i][0] = sqrt( 1 - nodes[j] * nodes[j] ) * std::cos( phi[i] );
_points[j * ( 2 * _order ) + i][1] = sqrt( 1 - nodes[j] * nodes[j] ) * std::sin( phi[i] );
_points[j * ( 2 * _order ) + i][2] = nodes[j];
_weights[j * ( 2 * _order ) + i] = 2.0 * M_PI / _order * weights[j];
}
}
}
void QGaussLegendreTensorized::SetConnectivity() { // TODO
......@@ -21,3 +68,74 @@ void QGaussLegendreTensorized::SetConnectivity() { // TODO
VectorVectorU connectivity;
_connectivity = connectivity;
}
std::pair<Vector, Matrix> QGaussLegendreTensorized::ComputeEigenValTriDiagMatrix( const Matrix& mat ) {
unsigned n = mat.rows();
Vector d( n, 0.0 ), e( n, 0.0 );
Matrix z( n, n, 0.0 );
for( unsigned i = 0; i < n; ++i ) {
d[i] = mat( i, i );
z( i, i ) = 1.0;
i == 0 ? e[i] = 0.0 : e[i] = mat( i, i - 1 );
}
int m, l, iter, i, k;
m = l = iter = i = k = 0;
double s, r, p, g, f, dd, c, b;
s = r = p = g = f = dd = c = b = 0.0;
const double eps = std::numeric_limits<double>::epsilon();
for( i = 1; i < static_cast<int>( n ); i++ ) e[i - 1] = e[i];
e[n - 1] = 0.0;
for( l = 0; l < static_cast<int>( n ); l++ ) {
iter = 0;
do {
for( m = l; m < static_cast<int>( n ) - 1; m++ ) {
dd = std::fabs( d[m] ) + std::fabs( d[m + 1] );
if( std::fabs( e[m] ) <= eps * dd ) break;
}
if( m != l ) {
if( iter++ == 30 ) ErrorMessages::Error( "Solving the tridiagonal matrix took too many iterations!", CURRENT_FUNCTION );
g = ( d[l + 1] - d[l] ) / ( 2.0 * e[l] );
r = Pythag( g, 1.0 );
g = d[m] - d[l] + e[l] / ( g + std::copysign( r, g ) );
s = c = 1.0;
p = 0.0;
for( i = m - 1; i >= l; i-- ) {
f = s * e[i];
b = c * e[i];
e[i + 1] = ( r = Pythag( f, g ) );
if( r == 0.0 ) {
d[i + 1] -= p;
e[m] = 0.0;
break;
}
s = f / r;
c = g / r;
g = d[i + 1] - p;
r = ( d[i] - g ) * s + 2.0 * c * b;
d[i + 1] = g + ( p = s * r );
g = c * r - b;
for( k = 0; k < static_cast<int>( n ); k++ ) {
f = z( static_cast<unsigned>( k ), static_cast<unsigned>( i ) + 1 );
z( static_cast<unsigned>( k ), static_cast<unsigned>( i ) + 1 ) =
s * z( static_cast<unsigned>( k ), static_cast<unsigned>( i ) ) + c * f;
z( static_cast<unsigned>( k ), static_cast<unsigned>( i ) ) =
c * z( static_cast<unsigned>( k ), static_cast<unsigned>( i ) ) - s * f;
}
}
if( r == 0.0 && i >= l ) continue;
d[l] -= p;
e[l] = g;
e[m] = 0.0;
}
} while( m != l );
}
return std::make_pair( d, z );
}
double QGaussLegendreTensorized::Pythag( const double a, const double b ) {
double absa = std::fabs( a ), absb = std::fabs( b );
return ( absa > absb ? absa * std::sqrt( 1.0 + ( absb / absa ) * ( absb / absa ) )
: ( absb == 0.0 ? 0.0 : absb * std::sqrt( 1.0 + ( absa / absb ) * ( absa / absb ) ) ) );
}
......@@ -30,7 +30,7 @@ std::string QLDFESA::GetLookupTable() {
case 1: lookupTable = __1_ldfesa_txt; break;
case 2: lookupTable = __2_ldfesa_txt; break;
case 3: lookupTable = __3_ldfesa_txt; break;
default: std::cerr << "Error: Invalid order chosen" << std::endl; exit( EXIT_FAILURE );
default: ErrorMessages::Error( "Invalid quadrature order chosen!", CURRENT_FUNCTION );
}
std::string lookupTableString( reinterpret_cast<char*>( lookupTable ) );
......
......@@ -61,7 +61,7 @@ std::string QLebedev::GetLookupTable() {
case 119: lookupTable = __119_lebedev_txt; break;
case 125: lookupTable = __125_lebedev_txt; break;
case 131: lookupTable = __131_lebedev_txt; break;
default: std::cerr << "Error: Invalid order chosen" << std::endl; exit( EXIT_FAILURE );
default: ErrorMessages::Error( "Invalid quadrature order chosen!", CURRENT_FUNCTION );
}
std::string lookupTableString( reinterpret_cast<char*>( lookupTable ) );
......
......@@ -36,7 +36,7 @@ std::string QLevelSymmetric::GetLookupTable() {
case 16: lookupTable = __16_levelsym_txt; break;
case 18: lookupTable = __18_levelsym_txt; break;
case 20: lookupTable = __20_levelsym_txt; break;