Commit 7d240211 authored by jannick.wolters's avatar jannick.wolters
Browse files

removed debug functions; readded logger; improved cxx20 to cxx17 port

parent d967febb
......@@ -10,7 +10,7 @@ TBD
## Build
### Required dependencies
- Compiler with C++20 support (e.g. gcc/g++ version 9)
- Compiler with C++17 support
- cmake >= v3.12.4
- LAPACK
- OpenMP
......
Subproject commit a51b4856377a71f81b6d74b9af459305c4c644f8
Subproject commit 83b9149930f392d7797b54fe97a66ab3f2120671
......@@ -30,7 +30,6 @@
#include "settings/CConfig.h"
using vtkPointsSP = vtkSmartPointer<vtkPoints>;
using vtkUnstructuredGridSP = vtkSmartPointer<vtkUnstructuredGrid>;
using vtkTriangleSP = vtkSmartPointer<vtkTriangle>;
......
......@@ -168,6 +168,7 @@ class CConfig
return _outputDir;
}
std::string inline GetOutputFile() const { return _outputFile; }
std::string inline GetLogDir() const { return _logDir; }
// Quadrature Structure
QUAD_NAME inline GetQuadName() const { return _quadName; }
......
......@@ -15,8 +15,8 @@
// --- Definition for global constants goes here ---
const double PI_NUMBER = 4.0 * std::atan( 1.0 ); /*!< \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 ---
......
#ifndef SNSOLVER_H
#define SNSOLVER_H
#include <mpi.h>
#include "solver.h"
#include "typedef.h"
......
......@@ -9,51 +9,53 @@
#include <iostream>
#include "spdlog/spdlog.h"
class CRTSNError
{
public:
CRTSNError();
inline static void Error(std::string ErrorMsg, std::string FunctionName){
//if (Rank == 0){ //For MPI implementation later
std::cout << std::endl << std::endl;
std::cout << "Error in \"" << FunctionName << "\": " << std::endl;
std::cout << "-------------------------------------------------------------------------" << std::endl;
std::cout << ErrorMsg << std::endl;
std::cout << "------------------------------ Error Exit -------------------------------" << std::endl;
std::cout << std::endl << std::endl;
//}
exit(EXIT_FAILURE);
}
inline static void OptionNotSetError(std::string OptionName, std::string FunctionName){
//if (Rank == 0){
std::cout << std::endl << std::endl;
std::cout << "Error in \"" << FunctionName << "\": " << std::endl;
std::cout << "-------------------------------------------------------------------------" << std::endl;
std::cout << "Option "<< OptionName << " not set. Please check your config file." << std::endl;
std::cout << "------------------------------ Error Exit -------------------------------" << std::endl;
std::cout << std::endl << std::endl;
//}
exit(EXIT_FAILURE);
}
public:
CRTSNError();
inline static void Error( std::string ErrorMsg, std::string FunctionName ) {
auto log = spdlog::get( "event" );
// if (Rank == 0){ //For MPI implementation later
log->error( "\n" );
log->error( "Error in \"{0}\": ", FunctionName );
log->error( "-------------------------------------------------------------------------" );
log->error( ErrorMsg );
log->error( "------------------------------ Error Exit -------------------------------" );
log->error( "\n" );
//}
exit( EXIT_FAILURE );
}
inline static void OptionNotSetError( std::string OptionName, std::string FunctionName ) {
// if (Rank == 0){
std::cout << std::endl << std::endl;
std::cout << "Error in \"" << FunctionName << "\": " << std::endl;
std::cout << "-------------------------------------------------------------------------" << std::endl;
std::cout << "Option " << OptionName << " not set. Please check your config file." << std::endl;
std::cout << "------------------------------ Error Exit -------------------------------" << std::endl;
std::cout << std::endl << std::endl;
//}
exit( EXIT_FAILURE );
}
};
#endif // CRTSNERROR_H
#endif // CRTSNERROR_H
/* Depending on the compiler, define the correct macro to get the current function name */
#if defined(__GNUC__) || (defined(__ICC) && (__ICC >= 600))
# define CURRENT_FUNCTION __PRETTY_FUNCTION__
#elif defined(__FUNCSIG__)
# define CURRENT_FUNCTION __FUNCSIG__
#elif (defined(__INTEL_COMPILER) && (__INTEL_COMPILER >= 600))
# define CURRENT_FUNCTION __FUNCTION__
#elif defined(__STDC_VERSION__) && (__STDC_VERSION__ >= 199901)
# define CURRENT_FUNCTION __func__
#elif defined(__cplusplus) && (__cplusplus >= 201103)
# define CURRENT_FUNCTION __func__
#if defined( __GNUC__ ) || ( defined( __ICC ) && ( __ICC >= 600 ) )
#define CURRENT_FUNCTION __PRETTY_FUNCTION__
#elif defined( __FUNCSIG__ )
#define CURRENT_FUNCTION __FUNCSIG__
#elif( defined( __INTEL_COMPILER ) && ( __INTEL_COMPILER >= 600 ) )
#define CURRENT_FUNCTION __FUNCTION__
#elif defined( __STDC_VERSION__ ) && ( __STDC_VERSION__ >= 199901 )
#define CURRENT_FUNCTION __func__
#elif defined( __cplusplus ) && ( __cplusplus >= 201103 )
#define CURRENT_FUNCTION __func__
#else
# define CURRENT_FUNCTION "(unknown)"
#define CURRENT_FUNCTION "(unknown)"
#endif
......@@ -16,7 +16,7 @@ MESH_FILE = checkerboard.su2
% ---- Solver specifications ----
%
% CFL number
CFL_NUMBER = 0.9
CFL_NUMBER = 0.4
% Final time for simulation
TIME_FINAL = 0.3
......@@ -29,4 +29,4 @@ BC_DIRICHLET = ( void )
% Quadrature Type
QUAD_TYPE = MONTE_CARLO
% Quadrature Order
QUAD_ORDER = 42
QUAD_ORDER = 5000
......@@ -95,7 +95,10 @@ void ExportVTK( const std::string fileName,
cellData->SetTuple3( i, results[i][0][j], results[i][1][j], results[i][2][j] );
}
break;
default: std::cout << "[ERROR][IO::ExportVTK] Invalid dimension" << std::endl;
default:
auto log = spdlog::get( "event" );
log->error( "[ERROR][IO::ExportVTK] Invalid dimension" );
exit( EXIT_FAILURE );
}
grid->GetCellData()->AddArray( cellData );
}
......@@ -411,7 +414,7 @@ Settings* ReadInputFile( std::string inputFile ) {
auto cwd = std::filesystem::current_path();
std::string tmp = std::filesystem::path( inputFile ).parent_path().string();
if( tmp.substr( tmp.size() - 1 ) != "/" ) tmp.append( "/" );
if( tmp[tmp.size() - 1] != '/' ) tmp.append( "/" );
settings->_inputDir = tmp;
// section IO
......@@ -428,7 +431,7 @@ Settings* ReadInputFile( std::string inputFile ) {
auto outputDir = io->get_as<std::string>( "outputDir" );
if( outputDir ) {
std::string tmp = *outputDir;
if( tmp.substr( tmp.size() - 1 ) != "/" ) tmp.append( "/" );
if( tmp[tmp.size() - 1] != '/' ) tmp.append( "/" );
settings->_outputDir = std::filesystem::path( tmp );
}
else {
......@@ -448,7 +451,7 @@ Settings* ReadInputFile( std::string inputFile ) {
auto logDir = io->get_as<std::string>( "logDir" );
if( logDir ) {
std::string tmp = *logDir;
if( tmp.substr( tmp.size() - 1 ) != "/" ) tmp.append( "/" );
if( tmp[tmp.size() - 1] != '/' ) tmp.append( "/" );
settings->_logDir = std::filesystem::path( tmp );
}
else {
......
......@@ -9,17 +9,6 @@ int main( int argc, char** argv ) {
MPI_Init( &argc, &argv );
/*
Settings* settings = ReadInputFile( inputFile );
InitLogger( settings->GetLogDir(), spdlog::level::info, spdlog::level::info );
PrintLogHeader( settings->GetInputFile() );
// build solver
Solver* solver = Solver::Create( settings );
solver->Solve();
solver->Save();
*/
char config_file_name[MAX_STRING_SIZE];
std::string filename = ParseArguments( argc, argv );
......@@ -31,12 +20,19 @@ int main( int argc, char** argv ) {
// Load Settings from File
CConfig* config = new CConfig( config_file_name );
// build solver
// Init stdout and file logger
InitLogger( config->GetLogDir(), spdlog::level::info, spdlog::level::info );
// Print input file and run info to file
PrintLogHeader( config_file_name );
// Build solver
Solver* solver = Solver::Create( config );
// Run solver and export
solver->Solve();
solver->Save();
// TODO: call solver
MPI_Finalize();
return EXIT_SUCCESS;
}
......@@ -18,17 +18,11 @@ void Mesh::ComputeConnectivity() {
int comm_size, comm_rank;
MPI_Comm_size( MPI_COMM_WORLD, &comm_size );
MPI_Comm_rank( MPI_COMM_WORLD, &comm_rank );
// unsigned chunkSize = std::ceil( static_cast<float>( _numCells ) / static_cast<float>( comm_size ) );
// unsigned mpiCellStart = comm_rank * chunkSize;
// unsigned mpiCellEnd = std::min( ( comm_rank + 1 ) * chunkSize, _numCells );
// std::vector<int> neighborsFlatPart( _numNodesPerCell * chunkSize, -1 );
// std::vector<Vector> normalsFlatPart( _numNodesPerCell * chunkSize, Vector( _dim, -1.0 ) );
unsigned chunkSize = _numCells;
unsigned mpiCellStart = 0u;
unsigned mpiCellEnd = _numCells;
std::vector<int> neighborsFlatPart( _numNodesPerCell * _numCells, -1 );
std::vector<Vector> normalsFlatPart( _numNodesPerCell * _numCells, Vector( _dim, -1.0 ) );
unsigned chunkSize = std::ceil( static_cast<float>( _numCells ) / static_cast<float>( comm_size ) );
unsigned mpiCellStart = comm_rank * chunkSize;
unsigned mpiCellEnd = std::min( ( comm_rank + 1 ) * chunkSize, _numCells );
std::vector<int> neighborsFlatPart( _numNodesPerCell * chunkSize, -1 );
std::vector<Vector> normalsFlatPart( _numNodesPerCell * chunkSize, Vector( _dim, -1.0 ) );
auto sortedCells( _cells );
for( unsigned i = 0; i < _numCells; ++i ) {
......@@ -114,32 +108,6 @@ void Mesh::ComputeConnectivity() {
}
}
unsigned ctr = 0;
for( unsigned i = 0; i < _numCells; ++i ) {
if( _cellNeighbors[i].size() > _numNodesPerCell ) {
std::cout << i << "\t\t";
for( unsigned j = 0; j < _cellNeighbors[i].size(); ++j ) {
std::cout << _cellNeighbors[i][j] << "\t";
}
std::cout << std::endl;
ctr++;
}
}
std::cout << ctr << std::endl;
ctr = 0;
for( unsigned i = 0; i < _numCells; ++i ) {
if( _cellNeighbors[i].size() < _numNodesPerCell ) {
std::cout << i << "\t\t";
for( unsigned j = 0; j < _cellNeighbors[i].size(); ++j ) {
std::cout << _cellNeighbors[i][j] << "\t";
}
std::cout << std::endl;
ctr++;
}
}
std::cout << ctr << std::endl << std::endl << std::endl;
_isBoundaryCell.resize( _numCells, false );
for( unsigned i = 0; i < _numCells; ++i ) {
if( std::any_of( _cellNeighbors[i].begin(), _cellNeighbors[i].end(), [this]( unsigned i ) { return i == this->_ghostCellID; } ) )
......
#include "quadratures/qlookupquadrature.h"
#include "settings/GlobalConstants.h" // for PI_NUMBER
#include <fstream>
#include <sstream>
#include "quadratures/qlookupquadrature.h"
#include "settings/GlobalConstants.h" // for PI_NUMBER
QLookupQuadrature::QLookupQuadrature( unsigned order ) : Quadrature( order ){
}
QLookupQuadrature::QLookupQuadrature( unsigned order ) : Quadrature( order ) {}
void QLookupQuadrature::printAvailOrders() const{
void QLookupQuadrature::printAvailOrders() const {
std::cout << "Available orders: (";
for(unsigned i = 0; i < _availableOrders.size()-1 ; i++){
for( unsigned i = 0; i < _availableOrders.size() - 1; i++ ) {
std::cout << _availableOrders[i] << "|";
}
std::cout << _availableOrders[_availableOrders.size()-1] << ")\n";
std::cout << _availableOrders[_availableOrders.size() - 1] << ")\n";
}
bool QLookupQuadrature::CheckOrder(){
std::vector<unsigned>::iterator it = std::find(_availableOrders.begin(), _availableOrders.end(), _order);
if (it == _availableOrders.end()){
std::cerr << "ERROR! Order "<< _order << " for " << GetName() << " not available. " << std::endl;
bool QLookupQuadrature::CheckOrder() {
std::vector<unsigned>::iterator it = std::find( _availableOrders.begin(), _availableOrders.end(), _order );
auto log = spdlog::get( "event" );
if( it == _availableOrders.end() ) {
log->error( "ERROR! Order {0} for {1} not available.", _order, GetName() );
printAvailOrders();
exit(EXIT_FAILURE);
exit( EXIT_FAILURE );
}
return true;
}
void QLookupQuadrature::SetNq() {
//find iterator of current order
std::vector<unsigned>::iterator it = std::find(_availableOrders.begin(), _availableOrders.end(), _order);
// find iterator of current order
std::vector<unsigned>::iterator it = std::find( _availableOrders.begin(), _availableOrders.end(), _order );
// Get index of element from iterator
unsigned index = std::distance(_availableOrders.begin(), it);
unsigned index = std::distance( _availableOrders.begin(), it );
// Get _nq
_nq = _nqByOrder[index];
}
void QLookupQuadrature::SetPointsAndWeights() {
_weights.resize(_nq);
_points.resize(_nq);
_weights.resize( _nq );
_points.resize( _nq );
double sumWeights = 0;
......@@ -47,34 +46,35 @@ void QLookupQuadrature::SetPointsAndWeights() {
std::stringstream in( lookupTable );
for(unsigned idx_point = 0; idx_point < _nq ; idx_point++){
for( unsigned idx_point = 0; idx_point < _nq; idx_point++ ) {
count = 0;
_points[idx_point].resize(3);
_points[idx_point].resize( 3 );
line.clear();
in >> line; //Get line of lookupTable string
in >> line; // Get line of lookupTable string
if(line.empty()){
if( line.empty() ) {
std::cout << _name << " with order " << _order << ".\n";
std::cerr << "Internal Error: Length of lookup table vector does not fit requested point vector length.\n";
exit(EXIT_FAILURE);
exit( EXIT_FAILURE );
}
std::stringstream ss(line); //give line to stringstream
std::stringstream ss( line ); // give line to stringstream
for (double double_in; ss >> double_in;) { //parse line
if(count <3) _points[idx_point][count] = double_in;
else{
for( double double_in; ss >> double_in; ) { // parse line
if( count < 3 )
_points[idx_point][count] = double_in;
else {
sumWeights += double_in;
_weights[idx_point] = double_in;
}
count ++;
if (ss.peek() == ',') ss.ignore();
count++;
if( ss.peek() == ',' ) ss.ignore();
}
}
//Correct the scaling of the weights
for (unsigned idx =0; idx < _nq; idx ++){
// Correct the scaling of the weights
for( unsigned idx = 0; idx < _nq; idx++ ) {
_weights[idx] /= sumWeights;
_weights[idx] *= 4.0 * PI_NUMBER;
}
......
......@@ -4,41 +4,17 @@
SNSolver::SNSolver( CConfig* settings ) : Solver( settings ) {}
void SNSolver::Solve() {
// TODO: set up initial condition
auto log = spdlog::get( "event" );
// angular flux at next time step (maybe store angular flux at all time steps, since time becomes energy?)
VectorVector psiNew = _psi;
unsigned ctr = 0;
unsigned ctr2 = 0;
// DEBUG: Summation to zero property
for( unsigned j = 0; j < _nCells; ++j ) {
Vector tmp( 2, 0.0 );
for( unsigned l = 0; l < _neighbors[j].size(); ++l ) {
if( norm( _normals[j][l] ) < 1e-7 ) ctr2++;
tmp = tmp + _normals[j][l];
unsigned nghPosJ;
unsigned I = _neighbors[j][l];
if( I >= _nCells ) continue;
for( unsigned m = 0; m < _neighbors[I].size(); ++m ) {
if( _neighbors[I][m] == j ) nghPosJ = m;
}
if( norm( _normals[j][l] + _normals[I][nghPosJ] ) > 1e-7 ) {
ctr++;
std::cout << _normals[j][l] << std::endl << _normals[I][nghPosJ];
std::cout << "---------" << std::endl;
}
}
// std::cout << tmp << std::endl; // needs to sum up to 0
}
std::cout << ctr << "\t" << ctr2 << std::endl;
std::cout << _nTimeSteps << std::endl;
// return;
int rank;
MPI_Comm_rank( MPI_COMM_WORLD, &rank );
if( rank == 0 ) log->info( "{:10} {:10}", "t", "#placeholder" );
// loop over energies (pseudo-time)
for( unsigned n = 0; n < _nTimeSteps; ++n ) {
std::cout << "time " << n * _dt << std::endl;
// loop over all spatial cells
for( unsigned j = 0; j < _nCells; ++j ) {
if( _boundaryCells[j] ) continue;
......@@ -49,10 +25,6 @@ void SNSolver::Solve() {
for( unsigned l = 0; l < _neighbors[j].size(); ++l ) {
// store flux contribution on psiNew_sigmaSH20 to save memory
psiNew[j][k] -= _g->Flux( _quadPoints[k], _psi[j][k], _psi[_neighbors[j][l]][k], _normals[j][l] );
if( _g->Flux( _quadPoints[k], _psi[j][k], _psi[_neighbors[j][l]][k], _normals[j][l] ) !=
-_g->Flux( _quadPoints[k], _psi[_neighbors[j][l]][k], _psi[j][k], -_normals[j][l] ) )
std::cout << "flux missmatch" << std::endl;
}
// time update angular flux with numerical flux and total scattering cross section
psiNew[j][k] = _psi[j][k] + ( _dt / _areas[j] ) * psiNew[j][k] + _dt * _sigmaTH20[n] * _psi[j][k];
......@@ -61,6 +33,7 @@ void SNSolver::Solve() {
psiNew[j] += _sigmaSH20[n] * _psi[j] * _weights; // multiply scattering matrix with psi
}
_psi = psiNew;
if( rank == 0 ) log->info( "{:03.8f} {:01.5e}", n * _dt, 0.0 );
}
}
......@@ -132,4 +105,6 @@ void SNSolver::Save() const {
std::vector<std::vector<double>> scalarField( 1, flux );
std::vector<std::vector<std::vector<double>>> results{ scalarField };
ExportVTK( _settings->GetOutputFile(), results, fieldNames, _settings, _mesh );
auto log = spdlog::get( "event" );
log->info( "Result successfully exported to '{0}'!", _settings->GetOutputFile() );
}
......@@ -7,8 +7,6 @@
Solver::Solver( CConfig* settings ) : _settings( settings ) {
// @TODO save parameters from settings class
// std::cout << "In Solver..." << std::endl;
// build quadrature object and store quadrature points and weights
Quadrature* q = Quadrature::CreateQuadrature( settings->GetQuadName(), settings->GetQuadOrder() );
_quadPoints = q->GetPoints();
......@@ -22,8 +20,6 @@ Solver::Solver( CConfig* settings ) : _settings( settings ) {
_normals = _mesh->GetNormals();
_nCells = _mesh->GetNumCells();
// std::cout << "After Mesh..." << std::endl;
// setup angular flux array (maybe directly call SetupIC() from physics class? )
_psi = std::vector( _nCells, Vector( _nq, 1e-7 ) );
......@@ -32,9 +28,7 @@ Solver::Solver( CConfig* settings ) : _settings( settings ) {
// write some IC
for( unsigned j = 0; j < nodes.size(); ++j ) {
// std::cout << norm( nodes[j] - midPoint ) << std::endl;
if( norm( nodes[j] - midPoint ) <= 0.1 ) {
// std::cout << nodes[j] << std::endl;
for( unsigned k = 0; k < _nq; ++k ) {
_psi[j][k] = 1.0;
}
......@@ -62,14 +56,8 @@ Solver::Solver( CConfig* settings ) : _settings( settings ) {
double Solver::ComputeTimeStep( double cfl ) const {
double maxEdge = -1.0;
// std::cout << _nCells << std::endl;
// std::cout << _areas.size() << std::endl;
// std::cout << _normals.size() << std::endl;
for( unsigned j = 0; j < _nCells; ++j ) {
// std::cout << j;
// std::cout << " " << _areas[j] << " " << _normals[j].size() << std::endl;
for( unsigned l = 0; l < _normals[j].size(); ++l ) {
// std::cout << _normals[j][l] << std::endl;
double currentEdge = _areas[j] / norm( _normals[j][l] );
if( currentEdge > maxEdge ) maxEdge = currentEdge;
}
......
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