From c5e88e96a85e7795c561d9618aa17038b10ddab5 Mon Sep 17 00:00:00 2001 From: Steffen Schotthoefer Date: Thu, 12 Sep 2024 21:40:57 -0400 Subject: [PATCH 01/21] start debugging --- CMakeLists.txt | 3 +++ src/common/io.cpp | 13 ++++++++++--- 2 files changed, 13 insertions(+), 3 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 9cae01d2..12689398 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -35,6 +35,7 @@ find_package( OpenMP REQUIRED ) if( BUILD_MPI ) find_package( MPI REQUIRED ) include_directories( ${MPI_INCLUDE_PATH} ) + message( STATUS "MPI: enabled" ) endif() find_package( LAPACK REQUIRED ) @@ -63,6 +64,8 @@ include_directories( ${CMAKE_SOURCE_DIR}/ext/spdlog/include ) if( BUILD_MPI ) set( CORE_LIBRARIES ${LAPACK_LIBRARIES} ${BLAS_LIBRARIES} ${MPI_LIBRARIES} ${VTK_LIBRARIES} OpenMP::OpenMP_CXX -lstdc++fs ) + message( STATUS "MPI: libraries loaded" ) + else() set( CORE_LIBRARIES ${LAPACK_LIBRARIES} ${BLAS_LIBRARIES} ${VTK_LIBRARIES} OpenMP::OpenMP_CXX -lstdc++fs ) endif() diff --git a/src/common/io.cpp b/src/common/io.cpp index b19ed8e4..47cb089d 100644 --- a/src/common/io.cpp +++ b/src/common/io.cpp @@ -50,8 +50,13 @@ void ExportVTK( const std::string fileName, const std::vector>& outputFieldNames, const Mesh* mesh ) { int rank = 0; - // MPI_Comm_rank( MPI_COMM_WORLD, &rank ); - if( rank == 0 ) { + int nprocs = 1; +#ifdef BUILD_MPI + // Initialize MPI + MPI_Comm_size( MPI_COMM_WORLD, &nprocs ); + MPI_Comm_rank( MPI_COMM_WORLD, &rank ); +#endif +if( rank == 0 ) { unsigned dim = mesh->GetDim(); unsigned numCells = mesh->GetNumCells(); unsigned numNodes = mesh->GetNumNodes(); @@ -140,7 +145,9 @@ void ExportVTK( const std::string fileName, // auto log = spdlog::get( "event" ); // log->info( "Result successfully exported to '{0}'!", fileNameWithExt ); } - // MPI_Barrier( MPI_COMM_WORLD ); +#ifdef BUILD_MPI + MPI_Barrier( MPI_COMM_WORLD ); +#endif } Mesh* LoadSU2MeshFromFile( const Config* settings ) { From 3fcc2ea995eacf0b3317f9b5a9df579db788404c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Steffen=20Schotth=C3=B6fer?= Date: Fri, 13 Sep 2024 10:07:33 -0400 Subject: [PATCH 02/21] mpi bug fixed on local machine- testing on cluster# --- CMakeLists.txt | 10 ++++++---- src/main.cpp | 3 +-- src/solvers/snsolver_hpc.cpp | 2 ++ 3 files changed, 9 insertions(+), 6 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 12689398..1680de14 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -9,7 +9,7 @@ option( BUILD_UNITY "enables unity build for faster compile times" ON ) option( BUILD_CODE_COV "enables compiler option required for code coverage analysis" OFF ) option( BUILD_ML "enables build with tensorflow backend access" OFF ) option( BUILD_MPI "enables build with MPI access" OFF ) - +add_definitions(-DBUILD_MPI) ################################################# @@ -26,16 +26,18 @@ if( BUILD_UNITY AND NOT BUILD_CODE_COV ) endif() ################################################# - ### LIBRARIES ################################### set(CMAKE_MODULE_PATH ${PROJECT_SOURCE_DIR}/cmake ${CMAKE_MODULE_PATH}) find_package( OpenMP REQUIRED ) +message(STATUS "MPI build flag: ${BUILD_MPI}") if( BUILD_MPI ) find_package( MPI REQUIRED ) include_directories( ${MPI_INCLUDE_PATH} ) - message( STATUS "MPI: enabled" ) + message(STATUS "C++ compiler: ${CMAKE_CXX_COMPILER}") + message(STATUS "MPI C++ compiler: ${MPI_CXX_COMPILER}") + endif() find_package( LAPACK REQUIRED ) @@ -64,7 +66,7 @@ include_directories( ${CMAKE_SOURCE_DIR}/ext/spdlog/include ) if( BUILD_MPI ) set( CORE_LIBRARIES ${LAPACK_LIBRARIES} ${BLAS_LIBRARIES} ${MPI_LIBRARIES} ${VTK_LIBRARIES} OpenMP::OpenMP_CXX -lstdc++fs ) - message( STATUS "MPI: libraries loaded" ) + message( STATUS "MPI: Libraries loaded" ) else() set( CORE_LIBRARIES ${LAPACK_LIBRARIES} ${BLAS_LIBRARIES} ${VTK_LIBRARIES} OpenMP::OpenMP_CXX -lstdc++fs ) diff --git a/src/main.cpp b/src/main.cpp index 07cfcc1f..55c169dc 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -17,9 +17,8 @@ #include "solvers/solverbase.hpp" #ifdef BUILD_GUI -#include - #include "mainwindow.h" +#include #endif int main( int argc, char** argv ) { diff --git a/src/solvers/snsolver_hpc.cpp b/src/solvers/snsolver_hpc.cpp index d5816bb7..5b7ef930 100644 --- a/src/solvers/snsolver_hpc.cpp +++ b/src/solvers/snsolver_hpc.cpp @@ -820,6 +820,8 @@ double SNSolverHPC::ComputeTimeStep( double cfl ) const { } // 2D case double charSize = __DBL_MAX__; // minimum char size of all mesh cells in the mesh + +#pragma omp parallel for reduction( min : charSize ) for( unsigned j = 0; j < _nCells; j++ ) { double currCharSize = sqrt( _areas[j] ); if( currCharSize < charSize ) { From 3156758c5a03909e33e9a4fa73bca333e973f11f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Steffen=20Schotth=C3=B6fer?= Date: Fri, 13 Sep 2024 10:32:52 -0400 Subject: [PATCH 03/21] santiy checks --- src/common/mesh.cpp | 4 ++-- src/main.cpp | 5 +++++ src/solvers/snsolver_hpc.cpp | 1 + 3 files changed, 8 insertions(+), 2 deletions(-) diff --git a/src/common/mesh.cpp b/src/common/mesh.cpp index 8ff2f2d7..1e4adf01 100644 --- a/src/common/mesh.cpp +++ b/src/common/mesh.cpp @@ -582,10 +582,10 @@ unsigned Mesh::GetCellOfKoordinate( const double x, const double y ) const { unsigned koordinate_cell_id = std::numeric_limits::max(); bool found = false; - // #pragma omp parallel for shared( found ) +#pragma omp parallel for shared( found ) for( unsigned idx_cell = 0; idx_cell < _numCells; idx_cell++ ) { if( IsPointInsideCell( idx_cell, x, y ) ) { - // #pragma omp critical + #pragma omp critical { if( !found ) { koordinate_cell_id = idx_cell; diff --git a/src/main.cpp b/src/main.cpp index 55c169dc..45b60870 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -32,7 +32,12 @@ int main( int argc, char** argv ) { // Py_SetProgramName( program ); #ifdef BUILD_MPI MPI_Init( &argc, &argv ); + printf( "MPI initialized\n" ); #endif +#ifndef BUILD_MPI + printf( "MPI not initialized\n" ); +#endif + std::string filename = ParseArguments( argc, argv ); // CD Load Settings from File diff --git a/src/solvers/snsolver_hpc.cpp b/src/solvers/snsolver_hpc.cpp index 5b7ef930..aac283f4 100644 --- a/src/solvers/snsolver_hpc.cpp +++ b/src/solvers/snsolver_hpc.cpp @@ -87,6 +87,7 @@ SNSolverHPC::SNSolverHPC( Config* settings ) { auto nodes = _mesh->GetNodes(); auto cells = _mesh->GetCells(); +#pragma omp parallel for for( unsigned idx_cell = 0; idx_cell < _nCells; idx_cell++ ) { _areas[idx_cell] = areas[idx_cell]; } From 1295161b2e1254774a7829aaa8aabf529cd55b98 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Steffen=20Schotth=C3=B6fer?= Date: Fri, 13 Sep 2024 10:43:53 -0400 Subject: [PATCH 04/21] improve parallelization --- src/common/mesh.cpp | 5 +++-- src/main.cpp | 4 ++-- src/solvers/snsolver_hpc.cpp | 2 +- 3 files changed, 6 insertions(+), 5 deletions(-) diff --git a/src/common/mesh.cpp b/src/common/mesh.cpp index 1e4adf01..4740685a 100644 --- a/src/common/mesh.cpp +++ b/src/common/mesh.cpp @@ -577,15 +577,16 @@ double Mesh::GetDistanceToOrigin( unsigned idx_cell ) const { return sqrt( distance ); } + unsigned Mesh::GetCellOfKoordinate( const double x, const double y ) const { // Experimental parallel implementation unsigned koordinate_cell_id = std::numeric_limits::max(); bool found = false; -#pragma omp parallel for shared( found ) +//#pragma omp parallel for shared( found ) for( unsigned idx_cell = 0; idx_cell < _numCells; idx_cell++ ) { if( IsPointInsideCell( idx_cell, x, y ) ) { - #pragma omp critical + //#pragma omp critical { if( !found ) { koordinate_cell_id = idx_cell; diff --git a/src/main.cpp b/src/main.cpp index 45b60870..98b8980d 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -32,10 +32,10 @@ int main( int argc, char** argv ) { // Py_SetProgramName( program ); #ifdef BUILD_MPI MPI_Init( &argc, &argv ); - printf( "MPI initialized\n" ); + printf( "| KiT-RT compiled with MPI and OpenMP parallelization\n" ); #endif #ifndef BUILD_MPI - printf( "MPI not initialized\n" ); + printf( "| KiT-RT compiled with OpenMP, but without MPI parallelization\n" ); #endif std::string filename = ParseArguments( argc, argv ); diff --git a/src/solvers/snsolver_hpc.cpp b/src/solvers/snsolver_hpc.cpp index aac283f4..eca7f257 100644 --- a/src/solvers/snsolver_hpc.cpp +++ b/src/solvers/snsolver_hpc.cpp @@ -1456,7 +1456,7 @@ std::vector SNSolverHPC::linspace2D( const std::vector& start, result.resize( num_points ); double stepX = ( end[0] - start[0] ) / ( num_points - 1 ); double stepY = ( end[1] - start[1] ) / ( num_points - 1 ); - +#pragma omp parallel for for( unsigned i = 0; i < num_points; ++i ) { double x = start[0] + i * stepX; double y = start[1] + i * stepY; From c06f7274499ad925809c151ff26d02205bcc9898 Mon Sep 17 00:00:00 2001 From: ScSteffen Date: Tue, 17 Sep 2024 13:05:53 -0400 Subject: [PATCH 05/21] fixed mpi import, fixed csv output --- CMakeLists.txt | 2 +- src/common/config.cpp | 6 +++--- src/common/io.cpp | 18 +++++++++--------- src/common/mesh.cpp | 11 +++++------ src/main.cpp | 10 +++++----- src/solvers/snsolver_hpc.cpp | 22 +++++++++++----------- 6 files changed, 34 insertions(+), 35 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 1680de14..58470146 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -9,7 +9,6 @@ option( BUILD_UNITY "enables unity build for faster compile times" ON ) option( BUILD_CODE_COV "enables compiler option required for code coverage analysis" OFF ) option( BUILD_ML "enables build with tensorflow backend access" OFF ) option( BUILD_MPI "enables build with MPI access" OFF ) -add_definitions(-DBUILD_MPI) ################################################# @@ -33,6 +32,7 @@ find_package( OpenMP REQUIRED ) message(STATUS "MPI build flag: ${BUILD_MPI}") if( BUILD_MPI ) + add_definitions(-DIMPORT_MPI) find_package( MPI REQUIRED ) include_directories( ${MPI_INCLUDE_PATH} ) message(STATUS "C++ compiler: ${CMAKE_CXX_COMPILER}") diff --git a/src/common/config.cpp b/src/common/config.cpp index 1483ac77..60ae732a 100644 --- a/src/common/config.cpp +++ b/src/common/config.cpp @@ -5,7 +5,7 @@ * * Disclaimer: This class structure was copied and modifed with open source permission from SU2 v7.0.3 https://su2code.github.io/ */ -#ifdef BUILD_MPI +#ifdef IMPORT_MPI #include #endif #include "common/config.hpp" @@ -1207,7 +1207,7 @@ bool Config::TokenizeString( string& str, string& option_name, vector& o void Config::InitLogger() { int rank = 0; -#ifdef BUILD_MPI +#ifdef IMPORT_MPI MPI_Comm_rank( MPI_COMM_WORLD, &rank ); // Initialize MPI #endif @@ -1366,7 +1366,7 @@ void Config::InitLogger() { spdlog::flush_every( std::chrono::seconds( 5 ) ); } } -#ifdef BUILD_MPI +#ifdef IMPORT_MPI MPI_Barrier( MPI_COMM_WORLD ); #endif } diff --git a/src/common/io.cpp b/src/common/io.cpp index 47cb089d..38804953 100644 --- a/src/common/io.cpp +++ b/src/common/io.cpp @@ -15,7 +15,7 @@ #include #include -#ifdef BUILD_MPI +#ifdef IMPORT_MPI #include #endif @@ -49,14 +49,14 @@ void ExportVTK( const std::string fileName, const std::vector>>& outputFields, const std::vector>& outputFieldNames, const Mesh* mesh ) { - int rank = 0; + int rank = 0; int nprocs = 1; -#ifdef BUILD_MPI +#ifdef IMPORT_MPI // Initialize MPI MPI_Comm_size( MPI_COMM_WORLD, &nprocs ); MPI_Comm_rank( MPI_COMM_WORLD, &rank ); -#endif -if( rank == 0 ) { +#endif + if( rank == 0 ) { unsigned dim = mesh->GetDim(); unsigned numCells = mesh->GetNumCells(); unsigned numNodes = mesh->GetNumNodes(); @@ -145,8 +145,8 @@ if( rank == 0 ) { // auto log = spdlog::get( "event" ); // log->info( "Result successfully exported to '{0}'!", fileNameWithExt ); } -#ifdef BUILD_MPI - MPI_Barrier( MPI_COMM_WORLD ); +#ifdef IMPORT_MPI + MPI_Barrier( MPI_COMM_WORLD ); #endif } @@ -472,7 +472,7 @@ std::string ParseArguments( int argc, char* argv[] ) { void PrintLogHeader( std::string inputFile ) { int nprocs = 1; int rank = 0; -#ifdef BUILD_MPI +#ifdef IMPORT_MPI // Initialize MPI MPI_Comm_size( MPI_COMM_WORLD, &nprocs ); MPI_Comm_rank( MPI_COMM_WORLD, &rank ); @@ -512,7 +512,7 @@ void PrintLogHeader( std::string inputFile ) { } // log->info( "------------------------------------------------------------------------" ); } -#ifdef BUILD_MPI +#ifdef IMPORT_MPI MPI_Barrier( MPI_COMM_WORLD ); #endif } diff --git a/src/common/mesh.cpp b/src/common/mesh.cpp index 4740685a..5fd6ce6e 100644 --- a/src/common/mesh.cpp +++ b/src/common/mesh.cpp @@ -6,7 +6,7 @@ #include #include #include -#ifdef BUILD_MPI +#ifdef IMPORT_MPI #include #endif #include @@ -28,7 +28,7 @@ Mesh::Mesh( const Config* settings, } int nprocs = 1; int rank = 0; -#ifdef BUILD_MPI +#ifdef IMPORT_MPI MPI_Comm_size( MPI_COMM_WORLD, &nprocs ); MPI_Comm_rank( MPI_COMM_WORLD, &rank ); #endif @@ -95,7 +95,7 @@ Mesh::~Mesh() {} void Mesh::ComputeConnectivity() { int nprocs = 1; int rank = 0; -#ifdef BUILD_MPI +#ifdef IMPORT_MPI MPI_Comm_size( MPI_COMM_WORLD, &nprocs ); MPI_Comm_rank( MPI_COMM_WORLD, &rank ); #endif @@ -577,16 +577,15 @@ double Mesh::GetDistanceToOrigin( unsigned idx_cell ) const { return sqrt( distance ); } - unsigned Mesh::GetCellOfKoordinate( const double x, const double y ) const { // Experimental parallel implementation unsigned koordinate_cell_id = std::numeric_limits::max(); bool found = false; -//#pragma omp parallel for shared( found ) + // #pragma omp parallel for shared( found ) for( unsigned idx_cell = 0; idx_cell < _numCells; idx_cell++ ) { if( IsPointInsideCell( idx_cell, x, y ) ) { - //#pragma omp critical + // #pragma omp critical { if( !found ) { koordinate_cell_id = idx_cell; diff --git a/src/main.cpp b/src/main.cpp index 98b8980d..6782e9b8 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -4,7 +4,7 @@ * @version: 0.1 */ -#ifdef BUILD_MPI +#ifdef IMPORT_MPI #include #endif #include @@ -30,12 +30,12 @@ int main( int argc, char** argv ) { #else // wchar_t* program = Py_DecodeLocale( argv[0], NULL ); // Py_SetProgramName( program ); -#ifdef BUILD_MPI +#ifdef IMPORT_MPI MPI_Init( &argc, &argv ); printf( "| KiT-RT compiled with MPI and OpenMP parallelization\n" ); #endif -#ifndef BUILD_MPI - printf( "| KiT-RT compiled with OpenMP, but without MPI parallelization\n" ); +#ifndef IMPORT_MPI + printf( "| KiT-RT compiled with OpenMP, but without MPI parallelization\n" ); #endif std::string filename = ParseArguments( argc, argv ); @@ -70,7 +70,7 @@ int main( int argc, char** argv ) { } delete config; -#ifdef BUILD_MPI +#ifdef IMPORT_MPI MPI_Finalize(); #endif diff --git a/src/solvers/snsolver_hpc.cpp b/src/solvers/snsolver_hpc.cpp index eca7f257..8af4b407 100644 --- a/src/solvers/snsolver_hpc.cpp +++ b/src/solvers/snsolver_hpc.cpp @@ -1,4 +1,4 @@ -#ifdef BUILD_MPI +#ifdef IMPORT_MPI #include #endif #include "common/config.hpp" @@ -11,12 +11,12 @@ #include "toolboxes/textprocessingtoolbox.hpp" SNSolverHPC::SNSolverHPC( Config* settings ) { -#ifdef BUILD_MPI +#ifdef IMPORT_MPI // Initialize MPI MPI_Comm_size( MPI_COMM_WORLD, &_numProcs ); MPI_Comm_rank( MPI_COMM_WORLD, &_rank ); #endif -#ifndef BUILD_MPI +#ifndef IMPORT_MPI _numProcs = 1; // default values _rank = 0; #endif @@ -148,7 +148,7 @@ SNSolverHPC::SNSolverHPC( Config* settings ) { } // _mass += _scalarFlux[idx_cell] * _areas[idx_cell]; } -#ifdef BUILD_MPI +#ifdef IMPORT_MPI MPI_Barrier( MPI_COMM_WORLD ); #endif SetGhostCells(); @@ -156,7 +156,7 @@ SNSolverHPC::SNSolverHPC( Config* settings ) { PrepareScreenOutput(); // Screen Output PrepareHistoryOutput(); // History Output } -#ifdef BUILD_MPI +#ifdef IMPORT_MPI MPI_Barrier( MPI_COMM_WORLD ); #endif delete quad; @@ -530,12 +530,12 @@ void SNSolverHPC::FVMUpdate() { temp_scalarFlux[idx_cell] = localScalarFlux; // set flux } // MPI Allreduce: _scalarFlux -#ifdef BUILD_MPI +#ifdef IMPORT_MPI MPI_Barrier( MPI_COMM_WORLD ); MPI_Allreduce( temp_scalarFlux.data(), _scalarFlux.data(), _nCells, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD ); MPI_Barrier( MPI_COMM_WORLD ); #endif -#ifndef BUILD_MPI +#ifndef IMPORT_MPI _scalarFlux = temp_scalarFlux; #endif } @@ -679,7 +679,7 @@ void SNSolverHPC::IterPostprocessing() { } } // MPI Allreduce -#ifdef BUILD_MPI +#ifdef IMPORT_MPI double tmp_curScalarOutflow = 0.0; double tmp_curScalarOutflowPeri1 = 0.0; double tmp_curScalarOutflowPeri2 = 0.0; @@ -755,12 +755,12 @@ void SNSolverHPC::IterPostprocessing() { // probe values green ComputeQOIsGreenProbingLine(); -#ifdef BUILD_MPI +#ifdef IMPORT_MPI MPI_Barrier( MPI_COMM_WORLD ); MPI_Allreduce( temp_probingMoments.data(), _probingMoments.data(), 3 * n_probes, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD ); MPI_Barrier( MPI_COMM_WORLD ); #endif -#ifndef BUILD_MPI +#ifndef IMPORT_MPI for( unsigned idx_probe = 0; idx_probe < n_probes; idx_probe++ ) { // Loop over probing cells _probingMoments[Idx2D( idx_probe, 0, 3 )] = temp_probingMoments[Idx2D( idx_probe, 0, 3 )]; _probingMoments[Idx2D( idx_probe, 1, 3 )] = temp_probingMoments[Idx2D( idx_probe, 1, 3 )]; @@ -1147,7 +1147,7 @@ void SNSolverHPC::PrintHistoryOutput( unsigned idx_iter ) { } lineToPrint += tmp + ","; } - tmp = TextProcessingToolbox::DoubleToScientificNotation( _historyOutputFields[_settings->GetNScreenOutput() - 1] ); + tmp = TextProcessingToolbox::DoubleToScientificNotation( _historyOutputFields[_settings->GetNHistoryOutput() - 1] ); lineToPrint += tmp; // Last element without comma if( _settings->GetHistoryOutputFrequency() != 0 && idx_iter % (unsigned)_settings->GetHistoryOutputFrequency() == 0 ) { From 339390810d88f14af80c63cc1cb46c9e6ffaf3a7 Mon Sep 17 00:00:00 2001 From: Steffen Schotthoefer Date: Tue, 17 Sep 2024 13:15:13 -0400 Subject: [PATCH 06/21] fixed mpi bug --- src/common/io.cpp | 9 +++------ src/solvers/snsolver_hpc.cpp | 4 +++- 2 files changed, 6 insertions(+), 7 deletions(-) diff --git a/src/common/io.cpp b/src/common/io.cpp index 47cb089d..79d0d831 100644 --- a/src/common/io.cpp +++ b/src/common/io.cpp @@ -49,14 +49,14 @@ void ExportVTK( const std::string fileName, const std::vector>>& outputFields, const std::vector>& outputFieldNames, const Mesh* mesh ) { - int rank = 0; + int rank = 0; int nprocs = 1; #ifdef BUILD_MPI // Initialize MPI MPI_Comm_size( MPI_COMM_WORLD, &nprocs ); MPI_Comm_rank( MPI_COMM_WORLD, &rank ); -#endif -if( rank == 0 ) { +#endif + if( rank == 0 ) { unsigned dim = mesh->GetDim(); unsigned numCells = mesh->GetNumCells(); unsigned numNodes = mesh->GetNumNodes(); @@ -145,9 +145,6 @@ if( rank == 0 ) { // auto log = spdlog::get( "event" ); // log->info( "Result successfully exported to '{0}'!", fileNameWithExt ); } -#ifdef BUILD_MPI - MPI_Barrier( MPI_COMM_WORLD ); -#endif } Mesh* LoadSU2MeshFromFile( const Config* settings ) { diff --git a/src/solvers/snsolver_hpc.cpp b/src/solvers/snsolver_hpc.cpp index eca7f257..f73be6e9 100644 --- a/src/solvers/snsolver_hpc.cpp +++ b/src/solvers/snsolver_hpc.cpp @@ -295,8 +295,10 @@ void SNSolverHPC::Solve() { PrintHistoryOutput( iter ); PrintVolumeOutput( iter ); } +#ifdef BUILD_MPI + MPI_Barrier( MPI_COMM_WORLD ); +#endif } - // --- Postprocessing --- if( _rank == 0 ) { DrawPostSolverOutput(); From 48aff8f25792d795c960bca02ac730c1e556d957 Mon Sep 17 00:00:00 2001 From: ScSteffen Date: Mon, 23 Sep 2024 14:05:39 -0400 Subject: [PATCH 07/21] update moment output --- src/main.cpp | 6 ++++- src/solvers/snsolver_hpc.cpp | 43 ++++++++++++++++++++++-------------- 2 files changed, 31 insertions(+), 18 deletions(-) diff --git a/src/main.cpp b/src/main.cpp index 6782e9b8..7a0820cb 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -32,7 +32,11 @@ int main( int argc, char** argv ) { // Py_SetProgramName( program ); #ifdef IMPORT_MPI MPI_Init( &argc, &argv ); - printf( "| KiT-RT compiled with MPI and OpenMP parallelization\n" ); + int rank; + MPI_Comm_rank( MPI_COMM_WORLD, &rank ); + if ( rank == 0 ) { + printf( "| KiT-RT compiled with MPI and OpenMP parallelization\n" ); + } #endif #ifndef IMPORT_MPI printf( "| KiT-RT compiled with OpenMP, but without MPI parallelization\n" ); diff --git a/src/solvers/snsolver_hpc.cpp b/src/solvers/snsolver_hpc.cpp index e4436a0d..97bfafd7 100644 --- a/src/solvers/snsolver_hpc.cpp +++ b/src/solvers/snsolver_hpc.cpp @@ -293,10 +293,13 @@ void SNSolverHPC::Solve() { // --- Print Output --- PrintScreenOutput( iter ); PrintHistoryOutput( iter ); - PrintVolumeOutput( iter ); } #ifdef BUILD_MPI MPI_Barrier( MPI_COMM_WORLD ); +#endif + PrintVolumeOutput( iter ); +#ifdef BUILD_MPI + MPI_Barrier( MPI_COMM_WORLD ); #endif } // --- Postprocessing --- @@ -308,11 +311,15 @@ void SNSolverHPC::Solve() { void SNSolverHPC::PrintVolumeOutput( int idx_iter ) { if( _settings->GetVolumeOutputFrequency() != 0 && idx_iter % (unsigned)_settings->GetVolumeOutputFrequency() == 0 ) { WriteVolumeOutput( idx_iter ); - ExportVTK( _settings->GetOutputFile() + "_" + std::to_string( idx_iter ), _outputFields, _outputFieldNames, _mesh ); // slow + if(_rank==0){ + ExportVTK( _settings->GetOutputFile() + "_" + std::to_string( idx_iter ), _outputFields, _outputFieldNames, _mesh ); // slow + } } if( idx_iter == (int)_nIter - 1 ) { // Last iteration write without suffix. WriteVolumeOutput( idx_iter ); - ExportVTK( _settings->GetOutputFile(), _outputFields, _outputFieldNames, _mesh ); + if(_rank==0){ + ExportVTK( _settings->GetOutputFile(), _outputFields, _outputFieldNames, _mesh ); + } } } @@ -1249,30 +1256,32 @@ void SNSolverHPC::WriteVolumeOutput( unsigned idx_iter ) { for( unsigned idx_group = 0; idx_group < nGroups; idx_group++ ) { switch( _settings->GetVolumeOutput()[idx_group] ) { case MINIMAL: - // for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) { - _outputFields[idx_group][0] = _scalarFlux; //[idx_cell]; - //} + if(_rank==0){ + _outputFields[idx_group][0] = _scalarFlux; + } break; case MOMENTS: +#ifdef BUILD_MPI + MPI_Barrier( MPI_COMM_WORLD ); +#endif #pragma omp parallel for for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) { _outputFields[idx_group][0][idx_cell] = 0.0; _outputFields[idx_group][1][idx_cell] = 0.0; - for( unsigned idx_moments = 0; idx_moments < _nOutputMoments; idx_moments++ ) { - _outputFields[idx_group][0][idx_cell] = 0.0; - _outputFields[idx_group][1][idx_cell] = 0.0; - _outputFields[idx_group][2][idx_cell] = 0.0; - // for( unsigned idx_sys = _startSysIdx; idx_sys < _endSysIdx; idx_sys++ ) { // TODO - // _outputFields[idx_group][0][idx_cell] += - // _quadPts[Idx2D( idx_sys, 0, _nDim )] * _sol[Idx2D( idx_cell, idx_sys, _localNSys )] * _quadWeights[idx_sys]; - // _outputFields[idx_group][1][idx_cell] += - // _quadPts[Idx2D( idx_sys, 1, _nDim )] * _sol[Idx2D( idx_cell, idx_sys, _localNSys )] * _quadWeights[idx_sys]; - // } + for( unsigned idx_sys = _startSysIdx; idx_sys < _endSysIdx; idx_sys++ ) { // TODO + _outputFields[idx_group][0][idx_cell] += + _quadPts[Idx2D( idx_sys, 0, _nDim )] * _sol[Idx2D( idx_cell, idx_sys, _localNSys )] * _quadWeights[idx_sys]; + _outputFields[idx_group][1][idx_cell] += + _quadPts[Idx2D( idx_sys, 1, _nDim )] * _sol[Idx2D( idx_cell, idx_sys, _localNSys )] * _quadWeights[idx_sys]; } } +#ifdef BUILD_MPI + MPI_Allreduce( _outputFields[idx_group][0].data(), _outputFields[idx_group][0].data(), _nCells, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD ); + MPI_Allreduce( _outputFields[idx_group][1].data(), _outputFields[idx_group][1].data(), _nCells, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD ); + MPI_Barrier( MPI_COMM_WORLD ); +#endif break; - default: ErrorMessages::Error( "Volume Output Group not defined for HPC SN Solver!", CURRENT_FUNCTION ); break; } } From f7ac225272834d4c695c8d496716ee50698c517a Mon Sep 17 00:00:00 2001 From: Steffen Date: Mon, 23 Sep 2024 16:38:20 -0400 Subject: [PATCH 08/21] another mpi fix --- src/solvers/snsolver_hpc.cpp | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/src/solvers/snsolver_hpc.cpp b/src/solvers/snsolver_hpc.cpp index 97bfafd7..bc629221 100644 --- a/src/solvers/snsolver_hpc.cpp +++ b/src/solvers/snsolver_hpc.cpp @@ -1262,14 +1262,11 @@ void SNSolverHPC::WriteVolumeOutput( unsigned idx_iter ) { break; case MOMENTS: -#ifdef BUILD_MPI - MPI_Barrier( MPI_COMM_WORLD ); -#endif #pragma omp parallel for for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) { _outputFields[idx_group][0][idx_cell] = 0.0; _outputFields[idx_group][1][idx_cell] = 0.0; - for( unsigned idx_sys = _startSysIdx; idx_sys < _endSysIdx; idx_sys++ ) { // TODO + for( unsigned idx_sys = 0; idx_sys < _localNSys; idx_sys++ ) { _outputFields[idx_group][0][idx_cell] += _quadPts[Idx2D( idx_sys, 0, _nDim )] * _sol[Idx2D( idx_cell, idx_sys, _localNSys )] * _quadWeights[idx_sys]; _outputFields[idx_group][1][idx_cell] += @@ -1277,6 +1274,7 @@ void SNSolverHPC::WriteVolumeOutput( unsigned idx_iter ) { } } #ifdef BUILD_MPI + MPI_Barrier( MPI_COMM_WORLD ); MPI_Allreduce( _outputFields[idx_group][0].data(), _outputFields[idx_group][0].data(), _nCells, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD ); MPI_Allreduce( _outputFields[idx_group][1].data(), _outputFields[idx_group][1].data(), _nCells, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD ); MPI_Barrier( MPI_COMM_WORLD ); From 3f5858fde0ed8012af8277e97b765ed2996f0276 Mon Sep 17 00:00:00 2001 From: Steffen Date: Tue, 24 Sep 2024 15:15:12 -0400 Subject: [PATCH 09/21] changed indices to unsigned long to prevent overflow for large meshes --- include/solvers/snsolver_hpc.hpp | 32 ++--- src/common/mesh.cpp | 6 +- src/problems/symmetrichohlraum.cpp | 21 +-- src/solvers/snsolver_hpc.cpp | 201 +++++++++++++++-------------- 4 files changed, 135 insertions(+), 125 deletions(-) diff --git a/include/solvers/snsolver_hpc.hpp b/include/solvers/snsolver_hpc.hpp index 0046621e..f37bb98b 100644 --- a/include/solvers/snsolver_hpc.hpp +++ b/include/solvers/snsolver_hpc.hpp @@ -20,9 +20,9 @@ class SNSolverHPC private: int _rank; int _numProcs; - unsigned _localNSys; - unsigned _startSysIdx; - unsigned _endSysIdx; + unsigned long _localNSys; + unsigned long _startSysIdx; + unsigned long _endSysIdx; double _currTime; /*!< @brief wall-time after current iteration */ Config* _settings; /*!< @brief config class for global information */ @@ -30,16 +30,16 @@ class SNSolverHPC ProblemBase* _problem; // Time - unsigned _nIter; /*!< @brief number of time steps, for non CSD, this equals _nEnergies, for _csd, _maxIter =_nEnergies-1*/ - double _dT; /*!< @brief energy/time step size */ + unsigned long _nIter; /*!< @brief number of time steps, for non CSD, this equals _nEnergies, for _csd, _maxIter =_nEnergies-1*/ + double _dT; /*!< @brief energy/time step size */ // Mesh related members, memory optimized - unsigned _nCells; /*!< @brief number of spatial cells */ - unsigned _nSys; /*!< @brief number of equations in the transport system, i.e. num quad pts */ - unsigned _nq; /*!< @brief number of quadrature points */ - unsigned _nDim; - unsigned _nNbr; - unsigned _nNodes; + unsigned long _nCells; /*!< @brief number of spatial cells */ + unsigned long _nSys; /*!< @brief number of equations in the transport system, i.e. num quad pts */ + unsigned long _nq; /*!< @brief number of quadrature points */ + unsigned long _nDim; + unsigned long _nNbr; + unsigned long _nNodes; std::vector _areas; /*!< @brief surface area of all spatial cells, dim(_areas) = _NCells */ @@ -111,9 +111,9 @@ class SNSolverHPC double _curAbsorptionHohlraumVertical; double _curAbsorptionHohlraumHorizontal; double _varAbsorptionHohlraumGreen; - std::vector< std::vector> _probingCellsHohlraum; /*!< @brief Indices of cells that contain a probing sensor */ - std::vector _probingMoments; /*!< @brief Solution Momnets at the probing cells that contain a probing sensor */ - unsigned _probingMomentsTimeIntervals; /*!< @brief Solution Momnets at the probing cells that contain a probing sensor */ + std::vector> _probingCellsHohlraum; /*!< @brief Indices of cells that contain a probing sensor */ + std::vector _probingMoments; /*!< @brief Solution Momnets at the probing cells that contain a probing sensor */ + unsigned _probingMomentsTimeIntervals; /*!< @brief Solution Momnets at the probing cells that contain a probing sensor */ unsigned _nProbingCellsLineGreen; /*!< @brief Number of sampling cells that contain a probing sensor for the sliding window */ std::vector _probingCellsLineGreen; /*!< @brief Indices of cells that contain a probing sensor for the sliding window */ @@ -198,8 +198,8 @@ class SNSolverHPC void DrawPostSolverOutput(); // Helper - unsigned Idx2D( unsigned idx1, unsigned idx2, unsigned len2 ); - unsigned Idx3D( unsigned idx1, unsigned idx2, unsigned idx3, unsigned len2, unsigned len3 ); + unsigned long Idx2D( unsigned long idx1, unsigned long idx2, unsigned long len2 ); + unsigned long Idx3D( unsigned long idx1, unsigned long idx2, unsigned long idx3, unsigned long len2, unsigned long len3 ); bool IsAbsorptionLattice( double x, double y ) const; void ComputeCellsPerimeterLattice(); diff --git a/src/common/mesh.cpp b/src/common/mesh.cpp index 5fd6ce6e..02935194 100644 --- a/src/common/mesh.cpp +++ b/src/common/mesh.cpp @@ -604,8 +604,8 @@ unsigned Mesh::GetCellOfKoordinate( const double x, const double y ) const { std::vector Mesh::GetCellsofBall( const double x, const double y, const double r ) const { std::vector cells_in_ball; - // Experimental parallel implementation - // #pragma omp parallel for +// Experimental parallel implementation +#pragma omp parallel for for( unsigned idx_cell = 0; idx_cell < _numCells; idx_cell++ ) { // Assume GetCellCenter returns the center coordinates of the cell double cell_x = _cellMidPoints[idx_cell][0]; @@ -615,7 +615,7 @@ std::vector Mesh::GetCellsofBall( const double x, const double y, cons double distance = std::sqrt( ( cell_x - x ) * ( cell_x - x ) + ( cell_y - y ) * ( cell_y - y ) ); if( distance <= r ) { - // #pragma omp critical +#pragma omp critical { cells_in_ball.push_back( idx_cell ); } } } diff --git a/src/problems/symmetrichohlraum.cpp b/src/problems/symmetrichohlraum.cpp index 73a65bc0..e73c104b 100644 --- a/src/problems/symmetrichohlraum.cpp +++ b/src/problems/symmetrichohlraum.cpp @@ -130,19 +130,23 @@ void SymmetricHohlraum::SetGhostCells() { if( vq[idx_q][0] < 0.0 ) right_inflow[idx_q] = 1.0; } +#pragma omp parallel for for( unsigned idx_cell = 0; idx_cell < _mesh->GetNumCells(); idx_cell++ ) { double x = _mesh->GetCellMidPoints()[idx_cell][0]; double y = _mesh->GetCellMidPoints()[idx_cell][1]; if( cellBoundaries[idx_cell] == BOUNDARY_TYPE::NEUMANN || cellBoundaries[idx_cell] == BOUNDARY_TYPE::DIRICHLET ) { - if( y < -0.6 ) - ghostCellMap.insert( { idx_cell, vertical_flow } ); - else if( y > 0.6 ) - ghostCellMap.insert( { idx_cell, vertical_flow } ); - else if( x < -0.6 ) - ghostCellMap.insert( { idx_cell, left_inflow } ); - else if( x > 0.6 ) - ghostCellMap.insert( { idx_cell, right_inflow } ); +#pragma omp critical + { + if( y < -0.6 ) + ghostCellMap.insert( { idx_cell, vertical_flow } ); + else if( y > 0.6 ) + ghostCellMap.insert( { idx_cell, vertical_flow } ); + else if( x < -0.6 ) + ghostCellMap.insert( { idx_cell, left_inflow } ); + else if( x > 0.6 ) + ghostCellMap.insert( { idx_cell, right_inflow } ); + } } } _ghostCells = ghostCellMap; @@ -318,6 +322,7 @@ std::vector SymmetricHohlraum::linspace2D( const std::vector& double stepX = ( end[0] - start[0] ) / ( num_points - 1 ); double stepY = ( end[1] - start[1] ) / ( num_points - 1 ); +#pragma omp parallel for for( unsigned i = 0; i < num_points; ++i ) { double x = start[0] + i * stepX; double y = start[1] + i * stepY; diff --git a/src/solvers/snsolver_hpc.cpp b/src/solvers/snsolver_hpc.cpp index bc629221..6430b268 100644 --- a/src/solvers/snsolver_hpc.cpp +++ b/src/solvers/snsolver_hpc.cpp @@ -31,13 +31,18 @@ SNSolverHPC::SNSolverHPC( Config* settings ) { auto quad = QuadratureBase::Create( settings ); _settings->SetNQuadPoints( quad->GetNq() ); - _nCells = _mesh->GetNumCells(); - _nNbr = _mesh->GetNumNodesPerCell(); - _nDim = _mesh->GetDim(); - _nNodes = _mesh->GetNumNodes(); + _nCells = static_cast( _mesh->GetNumCells() ); + _nNbr = static_cast( _mesh->GetNumNodesPerCell() ); + _nDim = static_cast( _mesh->GetDim() ); + _nNodes = static_cast( _mesh->GetNumNodes() ); + + _nq = static_cast( quad->GetNq() ); + _nSys = _nq; + + if( _numProcs > _nq ) { + ErrorMessages::Error( "The number of processors must be less than or equal to the number of quadrature points.", CURRENT_FUNCTION ); + } - _nq = quad->GetNq(); - _nSys = _nq; _localNSys = _nq / _numProcs; _startSysIdx = _rank * _localNSys; _endSysIdx = _rank * _localNSys + _localNSys; @@ -59,8 +64,8 @@ SNSolverHPC::SNSolverHPC( Config* settings ) { _relativeCellVertices = std::vector( _nCells * _nNbr * _nDim ); // Slope - _solDx = std::vector( _nCells * _localNSys * _nDim ); - _limiter = std::vector( _nCells * _localNSys ); + _solDx = std::vector( _nCells * _localNSys * _nDim, 0.0 ); + _limiter = std::vector( _nCells * _localNSys, 0.0 ); // Physics _sigmaS = std::vector( _nCells ); @@ -88,7 +93,7 @@ SNSolverHPC::SNSolverHPC( Config* settings ) { auto cells = _mesh->GetCells(); #pragma omp parallel for - for( unsigned idx_cell = 0; idx_cell < _nCells; idx_cell++ ) { + for( unsigned long idx_cell = 0; idx_cell < _nCells; idx_cell++ ) { _areas[idx_cell] = areas[idx_cell]; } @@ -106,8 +111,8 @@ SNSolverHPC::SNSolverHPC( Config* settings ) { // Copy to everything to solver _mass = 0; #pragma omp parallel for - for( unsigned idx_sys = 0; idx_sys < _localNSys; idx_sys++ ) { - for( unsigned idx_dim = 0; idx_dim < _nDim; idx_dim++ ) { + for( unsigned long idx_sys = 0; idx_sys < _localNSys; idx_sys++ ) { + for( unsigned long idx_dim = 0; idx_dim < _nDim; idx_dim++ ) { _quadPts[Idx2D( idx_sys, idx_dim, _nDim )] = quadPoints[idx_sys + _startSysIdx][idx_dim]; } _quadWeights[idx_sys] = @@ -115,18 +120,18 @@ SNSolverHPC::SNSolverHPC( Config* settings ) { } #pragma omp parallel for - for( unsigned idx_cell = 0; idx_cell < _nCells; idx_cell++ ) { + for( unsigned long idx_cell = 0; idx_cell < _nCells; idx_cell++ ) { _cellBoundaryTypes[idx_cell] = boundaryTypes[idx_cell]; - for( unsigned idx_dim = 0; idx_dim < _nDim; idx_dim++ ) { + for( unsigned long idx_dim = 0; idx_dim < _nDim; idx_dim++ ) { _cellMidPoints[Idx2D( idx_cell, idx_dim, _nDim )] = cellMidPts[idx_cell][idx_dim]; } - for( unsigned idx_nbr = 0; idx_nbr < _nNbr; idx_nbr++ ) { + for( unsigned long idx_nbr = 0; idx_nbr < _nNbr; idx_nbr++ ) { _neighbors[Idx2D( idx_cell, idx_nbr, _nNbr )] = neighbors[idx_cell][idx_nbr]; - for( unsigned idx_dim = 0; idx_dim < _nDim; idx_dim++ ) { + for( unsigned long idx_dim = 0; idx_dim < _nDim; idx_dim++ ) { _normals[Idx3D( idx_cell, idx_nbr, idx_dim, _nNbr, _nDim )] = normals[idx_cell][idx_nbr][idx_dim]; _interfaceMidPoints[Idx3D( idx_cell, idx_nbr, idx_dim, _nNbr, _nDim )] = interfaceMidPts[idx_cell][idx_nbr][idx_dim]; _relativeInterfaceMidPt[Idx3D( idx_cell, idx_nbr, idx_dim, _nNbr, _nDim )] = @@ -140,7 +145,7 @@ SNSolverHPC::SNSolverHPC( Config* settings ) { _sigmaT[idx_cell] = sigmaT[0][idx_cell]; _scalarFlux[idx_cell] = 0; - for( unsigned idx_sys = 0; idx_sys < _localNSys; idx_sys++ ) { + for( unsigned long idx_sys = 0; idx_sys < _localNSys; idx_sys++ ) { _source[Idx2D( idx_cell, idx_sys, _localNSys )] = source[0][idx_cell][0]; // CAREFUL HERE hardcoded to isotropic source _sol[Idx2D( idx_cell, idx_sys, _localNSys )] = 0.0; // initial condition is zero @@ -148,6 +153,7 @@ SNSolverHPC::SNSolverHPC( Config* settings ) { } // _mass += _scalarFlux[idx_cell] * _areas[idx_cell]; } + #ifdef IMPORT_MPI MPI_Barrier( MPI_COMM_WORLD ); #endif @@ -255,6 +261,7 @@ void SNSolverHPC::Solve() { std::chrono::duration duration; // Loop over energies (pseudo-time of continuous slowing down approach) for( unsigned iter = 0; iter < _nIter; iter++ ) { + if( iter == _nIter - 1 ) { // last iteration _dT = _settings->GetTEnd() - iter * _dT; } @@ -265,9 +272,9 @@ void SNSolverHPC::Solve() { ( _spatialOrder == 2 ) ? FluxOrder2() : FluxOrder1(); FVMUpdate(); #pragma omp parallel for - for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) { + for( unsigned long idx_cell = 0; idx_cell < _nCells; ++idx_cell ) { #pragma omp simd - for( unsigned idx_sys = 0; idx_sys < _localNSys; idx_sys++ ) { + for( unsigned long idx_sys = 0; idx_sys < _localNSys; idx_sys++ ) { _sol[Idx2D( idx_cell, idx_sys, _localNSys )] = 0.5 * ( solRK0[Idx2D( idx_cell, idx_sys, _localNSys )] + _sol[Idx2D( idx_cell, idx_sys, _localNSys )] ); // Solution averaging with HEUN @@ -275,6 +282,7 @@ void SNSolverHPC::Solve() { } } else { + ( _spatialOrder == 2 ) ? FluxOrder2() : FluxOrder1(); FVMUpdate(); } @@ -311,13 +319,13 @@ void SNSolverHPC::Solve() { void SNSolverHPC::PrintVolumeOutput( int idx_iter ) { if( _settings->GetVolumeOutputFrequency() != 0 && idx_iter % (unsigned)_settings->GetVolumeOutputFrequency() == 0 ) { WriteVolumeOutput( idx_iter ); - if(_rank==0){ + if( _rank == 0 ) { ExportVTK( _settings->GetOutputFile() + "_" + std::to_string( idx_iter ), _outputFields, _outputFieldNames, _mesh ); // slow } } if( idx_iter == (int)_nIter - 1 ) { // Last iteration write without suffix. WriteVolumeOutput( idx_iter ); - if(_rank==0){ + if( _rank == 0 ) { ExportVTK( _settings->GetOutputFile(), _outputFields, _outputFieldNames, _mesh ); } } @@ -328,53 +336,48 @@ void SNSolverHPC::FluxOrder2() { double const eps = 1e-10; #pragma omp parallel for - for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) { // Compute Slopes - if( _cellBoundaryTypes[idx_cell] == BOUNDARY_TYPE::NONE ) { // skip ghost cells -#pragma omp simd - for( unsigned idx_sys = 0; idx_sys < _localNSys; idx_sys++ ) { - - _limiter[Idx2D( idx_cell, idx_sys, _localNSys )] = 1.; // limiter should be zero at boundary - + for( unsigned long idx_cell = 0; idx_cell < _nCells; ++idx_cell ) { // Compute Slopes + if( _cellBoundaryTypes[idx_cell] == BOUNDARY_TYPE::NONE ) { // skip ghost cells + // #pragma omp simd + for( unsigned long idx_sys = 0; idx_sys < _localNSys; idx_sys++ ) { // + _limiter[Idx2D( idx_cell, idx_sys, _localNSys )] = 1.; // limiter should be zero at boundary// _solDx[Idx3D( idx_cell, idx_sys, 0, _localNSys, _nDim )] = 0.; - _solDx[Idx3D( idx_cell, idx_sys, 1, _localNSys, _nDim )] = 0.; - - double solInterfaceAvg = 0.0; - for( unsigned idx_nbr = 0; idx_nbr < _nNbr; ++idx_nbr ) { // Compute slopes and mininum and maximum - unsigned idx_nbr_glob = _neighbors[Idx2D( idx_cell, idx_nbr, _nNbr )]; - + _solDx[Idx3D( idx_cell, idx_sys, 1, _localNSys, _nDim )] = 0.; // + double solInterfaceAvg = 0.0; + for( unsigned long idx_nbr = 0; idx_nbr < _nNbr; ++idx_nbr ) { // Compute slopes and mininum and maximum + unsigned nbr_glob = _neighbors[Idx2D( idx_cell, idx_nbr, _nNbr )]; // // Slopes - solInterfaceAvg = 0.5 * ( _sol[Idx2D( idx_cell, idx_sys, _localNSys )] + _sol[Idx2D( idx_nbr_glob, idx_sys, _localNSys )] ); + solInterfaceAvg = 0.5 * ( _sol[Idx2D( idx_cell, idx_sys, _localNSys )] + _sol[Idx2D( nbr_glob, idx_sys, _localNSys )] ); _solDx[Idx3D( idx_cell, idx_sys, 0, _localNSys, _nDim )] += solInterfaceAvg * _normals[Idx3D( idx_cell, idx_nbr, 0, _nNbr, _nDim )]; _solDx[Idx3D( idx_cell, idx_sys, 1, _localNSys, _nDim )] += solInterfaceAvg * _normals[Idx3D( idx_cell, idx_nbr, 1, _nNbr, _nDim )]; - } - + } // _solDx[Idx3D( idx_cell, idx_sys, 0, _localNSys, _nDim )] /= _areas[idx_cell]; _solDx[Idx3D( idx_cell, idx_sys, 1, _localNSys, _nDim )] /= _areas[idx_cell]; } } } #pragma omp parallel for - for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) { // Compute Limiter - if( _cellBoundaryTypes[idx_cell] == BOUNDARY_TYPE::NONE ) { // skip ghost cells + for( unsigned long idx_cell = 0; idx_cell < _nCells; ++idx_cell ) { // Compute Limiter + if( _cellBoundaryTypes[idx_cell] == BOUNDARY_TYPE::NONE ) { // skip ghost cells -#pragma omp simd - for( unsigned idx_sys = 0; idx_sys < _localNSys; idx_sys++ ) { + // #pragma omp simd + for( unsigned long idx_sys = 0; idx_sys < _localNSys; idx_sys++ ) { double gaussPoint = 0; double r = 0; double minSol = _sol[Idx2D( idx_cell, idx_sys, _localNSys )]; double maxSol = _sol[Idx2D( idx_cell, idx_sys, _localNSys )]; - for( unsigned idx_nbr = 0; idx_nbr < _nNbr; ++idx_nbr ) { // Compute slopes and mininum and maximum - unsigned idx_nbr_glob = _neighbors[Idx2D( idx_cell, idx_nbr, _nNbr )]; + for( unsigned long idx_nbr = 0; idx_nbr < _nNbr; ++idx_nbr ) { // Compute slopes and mininum and maximum + unsigned nbr_glob = _neighbors[Idx2D( idx_cell, idx_nbr, _nNbr )]; // Compute ptswise max and minimum solultion values of current and neighbor cells - maxSol = std::max( _sol[Idx2D( idx_nbr_glob, idx_sys, _localNSys )], maxSol ); - minSol = std::min( _sol[Idx2D( idx_nbr_glob, idx_sys, _localNSys )], minSol ); + maxSol = std::max( _sol[Idx2D( nbr_glob, idx_sys, _localNSys )], maxSol ); + minSol = std::min( _sol[Idx2D( nbr_glob, idx_sys, _localNSys )], minSol ); } - for( unsigned idx_nbr = 0; idx_nbr < _nNbr; idx_nbr++ ) { // Compute limiter, see https://arxiv.org/pdf/1710.07187.pdf + for( unsigned long idx_nbr = 0; idx_nbr < _nNbr; idx_nbr++ ) { // Compute limiter, see https://arxiv.org/pdf/1710.07187.pdf // Compute test value at cell vertex, called gaussPt gaussPoint = @@ -407,27 +410,26 @@ void SNSolverHPC::FluxOrder2() { } } else { -#pragma omp simd - for( unsigned idx_sys = 0; idx_sys < _localNSys; idx_sys++ ) { + // #pragma omp simd + for( unsigned long idx_sys = 0; idx_sys < _localNSys; idx_sys++ ) { _limiter[Idx2D( idx_cell, idx_sys, _localNSys )] = 0.; // limiter should be zero at boundary _solDx[Idx3D( idx_cell, idx_sys, 0, _localNSys, _nDim )] = 0.; _solDx[Idx3D( idx_cell, idx_sys, 1, _localNSys, _nDim )] = 0.; } } } - #pragma omp parallel for - for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) { // Compute Flux + for( unsigned long idx_cell = 0; idx_cell < _nCells; ++idx_cell ) { // Compute Flux #pragma omp simd - for( unsigned idx_sys = 0; idx_sys < _localNSys; idx_sys++ ) { + for( unsigned long idx_sys = 0; idx_sys < _localNSys; idx_sys++ ) { _flux[Idx2D( idx_cell, idx_sys, _localNSys )] = 0.; } // Fluxes - for( unsigned idx_nbr = 0; idx_nbr < _nNbr; ++idx_nbr ) { + for( unsigned long idx_nbr = 0; idx_nbr < _nNbr; ++idx_nbr ) { if( _cellBoundaryTypes[idx_cell] == BOUNDARY_TYPE::NEUMANN && _neighbors[Idx2D( idx_cell, idx_nbr, _nNbr )] == _nCells ) { -#pragma omp simd - for( unsigned idx_sys = 0; idx_sys < _localNSys; idx_sys++ ) { + // #pragma omp simd + for( unsigned long idx_sys = 0; idx_sys < _localNSys; idx_sys++ ) { double localInner = _quadPts[Idx2D( idx_sys, 0, _nDim )] * _normals[Idx3D( idx_cell, idx_nbr, 0, _nNbr, _nDim )] + _quadPts[Idx2D( idx_sys, 1, _nDim )] * _normals[Idx3D( idx_cell, idx_nbr, 1, _nNbr, _nDim )]; if( localInner > 0 ) { @@ -440,12 +442,13 @@ void SNSolverHPC::FluxOrder2() { } } else { -// Second order + + unsigned long nbr_glob = _neighbors[Idx2D( idx_cell, idx_nbr, _nNbr )]; // global idx of neighbor cell + // Second order #pragma omp simd - for( unsigned idx_sys = 0; idx_sys < _localNSys; idx_sys++ ) { + for( unsigned long idx_sys = 0; idx_sys < _localNSys; idx_sys++ ) { // store flux contribution on psiNew_sigmaS to save memory - unsigned idx_nbr_glob = _neighbors[Idx2D( idx_cell, idx_nbr, _nNbr )]; // global idx of neighbor cell - double localInner = _quadPts[Idx2D( idx_sys, 0, _nDim )] * _normals[Idx3D( idx_cell, idx_nbr, 0, _nNbr, _nDim )] + + double localInner = _quadPts[Idx2D( idx_sys, 0, _nDim )] * _normals[Idx3D( idx_cell, idx_nbr, 0, _nNbr, _nDim )] + _quadPts[Idx2D( idx_sys, 1, _nDim )] * _normals[Idx3D( idx_cell, idx_nbr, 1, _nNbr, _nDim )]; _flux[Idx2D( idx_cell, idx_sys, _localNSys )] += @@ -455,14 +458,14 @@ void SNSolverHPC::FluxOrder2() { _relativeInterfaceMidPt[Idx3D( idx_cell, idx_nbr, 0, _nNbr, _nDim )] + _solDx[Idx3D( idx_cell, idx_sys, 1, _localNSys, _nDim )] * _relativeInterfaceMidPt[Idx3D( idx_cell, idx_nbr, 1, _nNbr, _nDim )] ) ) - : localInner * ( _sol[Idx2D( idx_nbr_glob, idx_sys, _localNSys )] + - _limiter[Idx2D( idx_nbr_glob, idx_sys, _localNSys )] * - ( _solDx[Idx3D( idx_nbr_glob, idx_sys, 0, _localNSys, _nDim )] * + : localInner * ( _sol[Idx2D( nbr_glob, idx_sys, _localNSys )] + + _limiter[Idx2D( nbr_glob, idx_sys, _localNSys )] * + ( _solDx[Idx3D( nbr_glob, idx_sys, 0, _localNSys, _nDim )] * ( _interfaceMidPoints[Idx3D( idx_cell, idx_nbr, 0, _nNbr, _nDim )] - - _cellMidPoints[Idx2D( idx_nbr_glob, 0, _nDim )] ) + - _solDx[Idx3D( idx_nbr_glob, idx_sys, 1, _localNSys, _nDim )] * + _cellMidPoints[Idx2D( nbr_glob, 0, _nDim )] ) + + _solDx[Idx3D( nbr_glob, idx_sys, 1, _localNSys, _nDim )] * ( _interfaceMidPoints[Idx3D( idx_cell, idx_nbr, 1, _nNbr, _nDim )] - - _cellMidPoints[Idx2D( idx_nbr_glob, 1, _nDim )] ) ) ); + _cellMidPoints[Idx2D( nbr_glob, 1, _nDim )] ) ) ); } } } @@ -472,18 +475,18 @@ void SNSolverHPC::FluxOrder2() { void SNSolverHPC::FluxOrder1() { #pragma omp parallel for - for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) { + for( unsigned long idx_cell = 0; idx_cell < _nCells; ++idx_cell ) { #pragma omp simd - for( unsigned idx_sys = 0; idx_sys < _localNSys; idx_sys++ ) { + for( unsigned long idx_sys = 0; idx_sys < _localNSys; idx_sys++ ) { _flux[Idx2D( idx_cell, idx_sys, _localNSys )] = 0.0; // Reset temporary variable } // Fluxes - for( unsigned idx_nbr = 0; idx_nbr < _nNbr; ++idx_nbr ) { + for( unsigned long idx_nbr = 0; idx_nbr < _nNbr; ++idx_nbr ) { if( _cellBoundaryTypes[idx_cell] == BOUNDARY_TYPE::NEUMANN && _neighbors[Idx2D( idx_cell, idx_nbr, _nNbr )] == _nCells ) { #pragma omp simd - for( unsigned idx_sys = 0; idx_sys < _localNSys; idx_sys++ ) { + for( unsigned long idx_sys = 0; idx_sys < _localNSys; idx_sys++ ) { double localInner = _quadPts[Idx2D( idx_sys, 0, _nDim )] * _normals[Idx3D( idx_cell, idx_nbr, 0, _nNbr, _nDim )] + _quadPts[Idx2D( idx_sys, 1, _nDim )] * _normals[Idx3D( idx_cell, idx_nbr, 1, _nNbr, _nDim )]; if( localInner > 0 ) { @@ -496,9 +499,9 @@ void SNSolverHPC::FluxOrder1() { } } else { - unsigned nbr_glob = _neighbors[Idx2D( idx_cell, idx_nbr, _nNbr )]; // global idx of neighbor cell + unsigned long nbr_glob = _neighbors[Idx2D( idx_cell, idx_nbr, _nNbr )]; // global idx of neighbor cell #pragma omp simd - for( unsigned idx_sys = 0; idx_sys < _localNSys; idx_sys++ ) { + for( unsigned long idx_sys = 0; idx_sys < _localNSys; idx_sys++ ) { double localInner = _quadPts[Idx2D( idx_sys, 0, _nDim )] * _normals[Idx3D( idx_cell, idx_nbr, 0, _nNbr, _nDim )] + _quadPts[Idx2D( idx_sys, 1, _nDim )] * _normals[Idx3D( idx_cell, idx_nbr, 1, _nNbr, _nDim )]; @@ -517,10 +520,10 @@ void SNSolverHPC::FVMUpdate() { std::vector temp_scalarFlux( _nCells ); // for MPI allreduce #pragma omp parallel for reduction( + : _mass, _rmsFlux ) - for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) { + for( unsigned long idx_cell = 0; idx_cell < _nCells; ++idx_cell ) { #pragma omp simd - for( unsigned idx_sys = 0; idx_sys < _localNSys; idx_sys++ ) { + for( unsigned long idx_sys = 0; idx_sys < _localNSys; idx_sys++ ) { // Update _sol[Idx2D( idx_cell, idx_sys, _localNSys )] = ( 1 - _dT * _sigmaT[idx_cell] ) * _sol[Idx2D( idx_cell, idx_sys, _localNSys )] - @@ -530,7 +533,7 @@ void SNSolverHPC::FVMUpdate() { double localScalarFlux = 0; #pragma omp simd reduction( + : localScalarFlux ) - for( unsigned idx_sys = 0; idx_sys < _localNSys; idx_sys++ ) { + for( unsigned long idx_sys = 0; idx_sys < _localNSys; idx_sys++ ) { _sol[Idx2D( idx_cell, idx_sys, _localNSys )] = std::max( _sol[Idx2D( idx_cell, idx_sys, _localNSys )], 0.0 ); localScalarFlux += _sol[Idx2D( idx_cell, idx_sys, _localNSys )] * _quadWeights[idx_sys]; } @@ -570,7 +573,7 @@ void SNSolverHPC::IterPostprocessing() { _curAbsorptionHohlraumVertical, \ _curAbsorptionHohlraumHorizontal, \ a_g ) reduction( max : _curMaxOrdinateOutflow, _curMaxAbsorptionLattice ) - for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) { + for( unsigned long idx_cell = 0; idx_cell < _nCells; ++idx_cell ) { if( _settings->GetProblemName() == PROBLEM_Lattice || _settings->GetProblemName() == PROBLEM_HalfLattice ) { if( IsAbsorptionLattice( _cellMidPoints[Idx2D( idx_cell, 0, _nDim )], _cellMidPoints[Idx2D( idx_cell, 1, _nDim )] ) ) { @@ -620,9 +623,9 @@ void SNSolverHPC::IterPostprocessing() { if( _settings->GetProblemName() == PROBLEM_Lattice ) { // Outflow out of inner and middle perimeter if( _isPerimeterLatticeCell1[idx_cell] ) { // inner - for( unsigned idx_nbr_helper = 0; idx_nbr_helper < _cellsLatticePerimeter1[idx_cell].size(); ++idx_nbr_helper ) { + for( unsigned long idx_nbr_helper = 0; idx_nbr_helper < _cellsLatticePerimeter1[idx_cell].size(); ++idx_nbr_helper ) { #pragma omp simd reduction( + : _curScalarOutflowPeri1 ) - for( unsigned idx_sys = 0; idx_sys < _localNSys; idx_sys++ ) { + for( unsigned long idx_sys = 0; idx_sys < _localNSys; idx_sys++ ) { double localInner = _quadPts[Idx2D( idx_sys, 0, _nDim )] * _normals[Idx3D( idx_cell, _cellsLatticePerimeter1[idx_cell][idx_nbr_helper], 0, _nNbr, _nDim )] + _quadPts[Idx2D( idx_sys, 1, _nDim )] * @@ -637,9 +640,9 @@ void SNSolverHPC::IterPostprocessing() { } } if( _isPerimeterLatticeCell2[idx_cell] ) { // middle - for( unsigned idx_nbr_helper = 0; idx_nbr_helper < _cellsLatticePerimeter2[idx_cell].size(); ++idx_nbr_helper ) { + for( unsigned long idx_nbr_helper = 0; idx_nbr_helper < _cellsLatticePerimeter2[idx_cell].size(); ++idx_nbr_helper ) { #pragma omp simd reduction( + : _curScalarOutflowPeri2 ) - for( unsigned idx_sys = 0; idx_sys < _localNSys; idx_sys++ ) { + for( unsigned long idx_sys = 0; idx_sys < _localNSys; idx_sys++ ) { double localInner = _quadPts[Idx2D( idx_sys, 0, _nDim )] * _normals[Idx3D( idx_cell, _cellsLatticePerimeter2[idx_cell][idx_nbr_helper], 0, _nNbr, _nDim )] + _quadPts[Idx2D( idx_sys, 1, _nDim )] * @@ -659,11 +662,11 @@ void SNSolverHPC::IterPostprocessing() { // Iterate over face cell faces double currOrdinatewiseOutflow = 0.0; - for( unsigned idx_nbr = 0; idx_nbr < _nNbr; ++idx_nbr ) { + for( unsigned long idx_nbr = 0; idx_nbr < _nNbr; ++idx_nbr ) { // Find face that points outward if( _neighbors[Idx2D( idx_cell, idx_nbr, _nNbr )] == _nCells ) { #pragma omp simd reduction( + : _curScalarOutflow ) reduction( max : _curMaxOrdinateOutflow ) - for( unsigned idx_sys = 0; idx_sys < _localNSys; idx_sys++ ) { + for( unsigned long idx_sys = 0; idx_sys < _localNSys; idx_sys++ ) { double localInner = _quadPts[Idx2D( idx_sys, 0, _nDim )] * _normals[Idx3D( idx_cell, idx_nbr, 0, _nNbr, _nDim )] + _quadPts[Idx2D( idx_sys, 1, _nDim )] * _normals[Idx3D( idx_cell, idx_nbr, 1, _nNbr, _nDim )]; @@ -717,7 +720,7 @@ void SNSolverHPC::IterPostprocessing() { std::vector temp_probingMoments( 3 * n_probes ); // for MPI allreduce #pragma omp parallel for reduction( + : _varAbsorptionHohlraumGreen ) - for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) { + for( unsigned long idx_cell = 0; idx_cell < _nCells; ++idx_cell ) { double x = _cellMidPoints[Idx2D( idx_cell, 0, _nDim )]; double y = _cellMidPoints[Idx2D( idx_cell, 1, _nDim )]; @@ -740,15 +743,15 @@ void SNSolverHPC::IterPostprocessing() { } // Probes value moments // #pragma omp parallel for - for( unsigned idx_probe = 0; idx_probe < n_probes; idx_probe++ ) { // Loop over probing cells + for( unsigned long idx_probe = 0; idx_probe < n_probes; idx_probe++ ) { // Loop over probing cells temp_probingMoments[Idx2D( idx_probe, 0, 3 )] = 0.0; temp_probingMoments[Idx2D( idx_probe, 1, 3 )] = 0.0; temp_probingMoments[Idx2D( idx_probe, 2, 3 )] = 0.0; - // for( unsigned idx_ball = 0; idx_ball < _probingCellsHohlraum[idx_probe].size(); idx_ball++ ) { + // for( unsigned long idx_ball = 0; idx_ball < _probingCellsHohlraum[idx_probe].size(); idx_ball++ ) { // std::cout << idx_ball << _areas[_probingCellsHohlraum[idx_probe][idx_ball]] / ( 0.01 * 0.01 * M_PI ) << std::endl; //} - for( unsigned idx_sys = 0; idx_sys < _localNSys; idx_sys++ ) { - for( unsigned idx_ball = 0; idx_ball < _probingCellsHohlraum[idx_probe].size(); idx_ball++ ) { + for( unsigned long idx_sys = 0; idx_sys < _localNSys; idx_sys++ ) { + for( unsigned long idx_ball = 0; idx_ball < _probingCellsHohlraum[idx_probe].size(); idx_ball++ ) { temp_probingMoments[Idx2D( idx_probe, 0, 3 )] += _sol[Idx2D( _probingCellsHohlraum[idx_probe][idx_ball], idx_sys, _localNSys )] * _quadWeights[idx_sys] * _areas[_probingCellsHohlraum[idx_probe][idx_ball]] / ( 0.01 * 0.01 * M_PI ); @@ -770,7 +773,7 @@ void SNSolverHPC::IterPostprocessing() { MPI_Barrier( MPI_COMM_WORLD ); #endif #ifndef IMPORT_MPI - for( unsigned idx_probe = 0; idx_probe < n_probes; idx_probe++ ) { // Loop over probing cells + for( unsigned long idx_probe = 0; idx_probe < n_probes; idx_probe++ ) { // Loop over probing cells _probingMoments[Idx2D( idx_probe, 0, 3 )] = temp_probingMoments[Idx2D( idx_probe, 0, 3 )]; _probingMoments[Idx2D( idx_probe, 1, 3 )] = temp_probingMoments[Idx2D( idx_probe, 1, 3 )]; _probingMoments[Idx2D( idx_probe, 2, 3 )] = temp_probingMoments[Idx2D( idx_probe, 2, 3 )]; @@ -832,7 +835,7 @@ double SNSolverHPC::ComputeTimeStep( double cfl ) const { double charSize = __DBL_MAX__; // minimum char size of all mesh cells in the mesh #pragma omp parallel for reduction( min : charSize ) - for( unsigned j = 0; j < _nCells; j++ ) { + for( unsigned long j = 0; j < _nCells; j++ ) { double currCharSize = sqrt( _areas[j] ); if( currCharSize < charSize ) { charSize = currCharSize; @@ -1243,9 +1246,9 @@ void SNSolverHPC::DrawPostSolverOutput() { log->info( "--------------------------- Solver Finished ----------------------------" ); } -unsigned SNSolverHPC::Idx2D( unsigned idx1, unsigned idx2, unsigned len2 ) { return idx1 * len2 + idx2; } +unsigned long SNSolverHPC::Idx2D( unsigned long idx1, unsigned long idx2, unsigned long len2 ) { return idx1 * len2 + idx2; } -unsigned SNSolverHPC::Idx3D( unsigned idx1, unsigned idx2, unsigned idx3, unsigned len2, unsigned len3 ) { +unsigned long SNSolverHPC::Idx3D( unsigned long idx1, unsigned long idx2, unsigned long idx3, unsigned long len2, unsigned long len3 ) { return ( idx1 * len2 + idx2 ) * len3 + idx3; } @@ -1256,8 +1259,8 @@ void SNSolverHPC::WriteVolumeOutput( unsigned idx_iter ) { for( unsigned idx_group = 0; idx_group < nGroups; idx_group++ ) { switch( _settings->GetVolumeOutput()[idx_group] ) { case MINIMAL: - if(_rank==0){ - _outputFields[idx_group][0] = _scalarFlux; + if( _rank == 0 ) { + _outputFields[idx_group][0] = _scalarFlux; } break; @@ -1267,17 +1270,19 @@ void SNSolverHPC::WriteVolumeOutput( unsigned idx_iter ) { _outputFields[idx_group][0][idx_cell] = 0.0; _outputFields[idx_group][1][idx_cell] = 0.0; for( unsigned idx_sys = 0; idx_sys < _localNSys; idx_sys++ ) { - _outputFields[idx_group][0][idx_cell] += - _quadPts[Idx2D( idx_sys, 0, _nDim )] * _sol[Idx2D( idx_cell, idx_sys, _localNSys )] * _quadWeights[idx_sys]; - _outputFields[idx_group][1][idx_cell] += - _quadPts[Idx2D( idx_sys, 1, _nDim )] * _sol[Idx2D( idx_cell, idx_sys, _localNSys )] * _quadWeights[idx_sys]; + _outputFields[idx_group][0][idx_cell] += + _quadPts[Idx2D( idx_sys, 0, _nDim )] * _sol[Idx2D( idx_cell, idx_sys, _localNSys )] * _quadWeights[idx_sys]; + _outputFields[idx_group][1][idx_cell] += + _quadPts[Idx2D( idx_sys, 1, _nDim )] * _sol[Idx2D( idx_cell, idx_sys, _localNSys )] * _quadWeights[idx_sys]; } } #ifdef BUILD_MPI - MPI_Barrier( MPI_COMM_WORLD ); - MPI_Allreduce( _outputFields[idx_group][0].data(), _outputFields[idx_group][0].data(), _nCells, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD ); - MPI_Allreduce( _outputFields[idx_group][1].data(), _outputFields[idx_group][1].data(), _nCells, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD ); - MPI_Barrier( MPI_COMM_WORLD ); + MPI_Barrier( MPI_COMM_WORLD ); + MPI_Allreduce( + _outputFields[idx_group][0].data(), _outputFields[idx_group][0].data(), _nCells, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD ); + MPI_Allreduce( + _outputFields[idx_group][1].data(), _outputFields[idx_group][1].data(), _nCells, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD ); + MPI_Barrier( MPI_COMM_WORLD ); #endif break; default: ErrorMessages::Error( "Volume Output Group not defined for HPC SN Solver!", CURRENT_FUNCTION ); break; From f079854fb4881ad7df8272a9a091d4dbefb3390a Mon Sep 17 00:00:00 2001 From: ScSteffen Date: Mon, 14 Oct 2024 19:32:49 -0400 Subject: [PATCH 10/21] added tessalation quadrature --- include/common/globalconstants.hpp | 2 + include/quadratures/qsphericaltessalation.hpp | 38 ++++ src/quadratures/qsphericaltessalation.cpp | 205 ++++++++++++++++++ src/quadratures/quadraturebase.cpp | 2 + src/solvers/snsolver_hpc.cpp | 9 +- tests/test_quadrature.cpp | 18 +- 6 files changed, 268 insertions(+), 6 deletions(-) create mode 100644 include/quadratures/qsphericaltessalation.hpp create mode 100644 src/quadratures/qsphericaltessalation.cpp diff --git a/include/common/globalconstants.hpp b/include/common/globalconstants.hpp index 826bd338..debd3830 100644 --- a/include/common/globalconstants.hpp +++ b/include/common/globalconstants.hpp @@ -62,6 +62,7 @@ enum QUAD_NAME { QUAD_LevelSymmetric, QUAD_Lebedev, QUAD_LDFESA, + QUAD_SphericalTessalation2D, QUAD_Product, QUAD_Rectangular1D, QUAD_Rectangular2D, @@ -82,6 +83,7 @@ inline std::map Quadrature_Map{ { "MONTE_CARLO", QUAD_Mo { "LEVEL_SYMMETRIC", QUAD_LevelSymmetric }, { "LEBEDEV", QUAD_Lebedev }, { "LDFESA", QUAD_LDFESA }, + { "TESSALATION", QUAD_SphericalTessalation2D }, { "MIDPOINT_1D", QUAD_Midpoint1D }, { "MIDPOINT_2D", QUAD_Midpoint2D }, { "MIDPOINT_3D", QUAD_Midpoint3D }, diff --git a/include/quadratures/qsphericaltessalation.hpp b/include/quadratures/qsphericaltessalation.hpp new file mode 100644 index 00000000..c7721d3e --- /dev/null +++ b/include/quadratures/qsphericaltessalation.hpp @@ -0,0 +1,38 @@ +#ifndef QSPHERICALTRIANGLE_H +#define QSPHERICALTRIANGLE_H + +#include "quadraturebase.hpp" + +class QSphericalTessalation : public QuadratureBase +{ + // Implementation is done accordingly to Kendall Atkinson 1981, Australian Matematical Society. + + public: + QSphericalTessalation( Config* settings ); + + virtual ~QSphericalTessalation() {} + + protected: + virtual bool CheckOrder(); + virtual inline void SetName() override { _name = "Tesselated Spherical Triangle Quadrature"; } + void SetNq() override; + void SetPointsAndWeights() override; + void SetConnectivity() override; + + private: + std::vector, 3>> subdivide_triangle( const std::array, 3>& triangle ); + std::vector, 3>> recursive_subdivide( const std::vector, 3>>& triangles, + int order ); + std::array compute_centroid( const std::array, 3>& triangle ); + std::array map_to_unit_sphere( const std::array& point ); + double dot_product( const std::array& v1, const std::array& v2 ); + double angle_between_vectors( const std::array& a, const std::array& b, const std::array& c ); + double spherical_triangle_area( const std::array& a, const std::array& b, const std::array& c ); + std::array midpoint( const std::array& p1, const std::array& p2 ); + void reflect_and_permute( const std::vector>& points, + const std::vector& weights, + std::vector>& full_points, + std::vector& full_weights ); +}; + +#endif // QSPHERICALTRIANGLE_H diff --git a/src/quadratures/qsphericaltessalation.cpp b/src/quadratures/qsphericaltessalation.cpp new file mode 100644 index 00000000..678c7438 --- /dev/null +++ b/src/quadratures/qsphericaltessalation.cpp @@ -0,0 +1,205 @@ +#include "quadratures/qsphericaltessalation.hpp" + +#include +#include +#include +#include + +QSphericalTessalation::QSphericalTessalation( Config* settings ) : QuadratureBase( settings ) { + SetName(); + CheckOrder(); + SetNq(); + SetPointsAndWeights(); + SetConnectivity(); + _supportedDimensions = { 2 }; // Extension to 3d is straightforward +} + +void QSphericalTessalation::SetPointsAndWeights() { + // Define basal triangle vertices + std::array A = { 1.0, 0.0, 0.0 }; + std::array B = { 0.0, 1.0, 0.0 }; + std::array C = { 0.0, 0.0, 1.0 }; + std::array, 3> basal_triangle = { A, B, C }; + + // Perform recursive subdivision + std::vector, 3>> final_triangles = recursive_subdivide( { basal_triangle }, _order ); + + // Compute centroids and map to unit sphere + std::vector> centroids; + std::vector> mapped_points; + std::vector weights; + + for( const auto& tri : final_triangles ) { + auto centroid = compute_centroid( tri ); + auto mapped_centroid = map_to_unit_sphere( centroid ); + centroids.push_back( centroid ); + mapped_points.push_back( mapped_centroid ); + + // Map triangle vertices to unit sphere + std::array a = map_to_unit_sphere( tri[0] ); + std::array b = map_to_unit_sphere( tri[1] ); + std::array c = map_to_unit_sphere( tri[2] ); + + // Calculate area using spherical excess + double area = spherical_triangle_area( a, b, c ); + weights.push_back( area ); + } + + // Vectors to hold the full set of points and weights for the upper hemisphere + std::vector> full_points; + std::vector full_weights; + + // Perform reflection and permutation + reflect_and_permute( mapped_points, weights, full_points, full_weights ); + + _nq = full_points.size(); + _pointsKarth.resize( _nq ); + _pointsSphere.resize( _nq ); + _weights.resize( _nq ); + for( size_t i = 0; i < _nq; ++i ) { + _pointsKarth[i] = { full_points[i][0], full_points[i][1], full_points[i][2] }; + _weights[i] = full_weights[i]; + } + double w = 0.0; + // Transform _points to _pointsSphere ==>transform (x,y,z) into (my,phi) + for( unsigned idx = 0; idx < _nq; idx++ ) { + _pointsSphere[idx].resize( 3 ); // (my,phi) + _pointsSphere[idx][0] = _pointsKarth[idx][2]; // my = z + _pointsSphere[idx][1] = atan2( _pointsKarth[idx][1], _pointsKarth[idx][0] ); // phi in [-pi,pi] + _pointsSphere[idx][2] = 1.0; // radius r + + // adapt intervall s.t. phi in [0,2pi] + if( _pointsSphere[idx][1] < 0 ) { + _pointsSphere[idx][1] = 2 * M_PI + _pointsSphere[idx][1]; + } + + // std::cout << _pointsKarth[idx][0] << " " << _pointsKarth[idx][1] << " " << _pointsKarth[idx][2] << std::endl; + // std::cout << _pointsSphere[idx][0] << " " << _pointsSphere[idx][1] << " " << _pointsSphere[idx][2] << std::endl; + // std::cout << _weights[idx] << std::endl; + w += _weights[idx]; + } + // std::cout << w << std::endl; + // exit( 1 ); +} + +void QSphericalTessalation::SetNq() { _nq = pow( GetOrder(), 4 ); } + +void QSphericalTessalation::SetConnectivity() { // TODO + // Not initialized for this quadrature. + VectorVectorU connectivity; + _connectivity = connectivity; +} + +bool QSphericalTessalation::CheckOrder() { return true; } + +// Helper function to compute midpoint of two vectors +std::array QSphericalTessalation::midpoint( const std::array& p1, const std::array& p2 ) { + return { ( p1[0] + p2[0] ) / 2, ( p1[1] + p2[1] ) / 2, ( p1[2] + p2[2] ) / 2 }; +} + +// Subdivide a triangle into four smaller triangles +std::vector, 3>> QSphericalTessalation::subdivide_triangle( const std::array, 3>& triangle ) { + std::array A = triangle[0]; + std::array B = triangle[1]; + std::array C = triangle[2]; + + // Compute midpoints of edges + std::array AB_mid = midpoint( A, B ); + std::array BC_mid = midpoint( B, C ); + std::array CA_mid = midpoint( C, A ); + + // Define four new triangles + std::vector, 3>> new_triangles = { + { A, AB_mid, CA_mid }, { AB_mid, B, BC_mid }, { CA_mid, BC_mid, C }, { AB_mid, BC_mid, CA_mid } }; + + return new_triangles; +} + +// Recursively subdivide triangles to the desired order +std::vector, 3>> +QSphericalTessalation::recursive_subdivide( const std::vector, 3>>& triangles, int order ) { + if( order == 1 ) { + return triangles; + } + else { + std::vector, 3>> new_triangles; + for( const auto& triangle : triangles ) { + std::vector, 3>> subdivided = subdivide_triangle( triangle ); + new_triangles.insert( new_triangles.end(), subdivided.begin(), subdivided.end() ); + } + return recursive_subdivide( new_triangles, order - 1 ); + } +} + +// Compute the centroid of a triangle +std::array QSphericalTessalation::compute_centroid( const std::array, 3>& triangle ) { + std::array A = triangle[0]; + std::array B = triangle[1]; + std::array C = triangle[2]; + + return { ( A[0] + B[0] + C[0] ) / 3, ( A[1] + B[1] + C[1] ) / 3, ( A[2] + B[2] + C[2] ) / 3 }; +} + +// Normalize a vector to map it to the unit sphere +std::array QSphericalTessalation::map_to_unit_sphere( const std::array& point ) { + double norm = std::sqrt( point[0] * point[0] + point[1] * point[1] + point[2] * point[2] ); + return { point[0] / norm, point[1] / norm, point[2] / norm }; +} + +// Helper function to calculate dot product between two vectors +double QSphericalTessalation::dot_product( const std::array& v1, const std::array& v2 ) { + return v1[0] * v2[0] + v1[1] * v2[1] + v1[2] * v2[2]; // Correct dot product. +} + +// Calculate angle at vertex A of spherical triangle ABC +double +QSphericalTessalation::angle_between_vectors( const std::array& a, const std::array& b, const std::array& c ) { + double ab = dot_product( b, a ); + double ac = dot_product( c, a ); + double bc = dot_product( b, c ); + return std::acos( ( ab - ac * bc ) / ( std::sqrt( 1 - ac * ac ) * std::sqrt( 1 - bc * bc ) ) ); +} + +// Calculate area of spherical triangle using spherical excess formula +double +QSphericalTessalation::spherical_triangle_area( const std::array& a, const std::array& b, const std::array& c ) { + double angle_A = angle_between_vectors( b, a, c ); + double angle_B = angle_between_vectors( c, b, a ); + double angle_C = angle_between_vectors( a, c, b ); + + // Spherical excess + return ( angle_A + angle_B + angle_C ) - M_PI; +} + +// Function to reflect and permute points from one octant to generate points and weights for the upper half-sphere +void QSphericalTessalation::reflect_and_permute( const std::vector>& points, + const std::vector& weights, + std::vector>& full_points, + std::vector& full_weights ) { + // Loop through the points in the first octant + for( size_t i = 0; i < points.size(); ++i ) { + const auto& point = points[i]; + double weight = weights[i]; + + double x = point[0]; + double y = point[1]; + double z = point[2]; + + // Generate the four permutations for the upper hemisphere + std::array, 4> permutations = { std::array{ x, y, z }, + std::array{ -x, y, z }, + std::array{ x, -y, z }, + std::array{ -x, -y, z } }; + + // Append each reflected point and its corresponding weight + for( const auto& perm : permutations ) { + full_points.push_back( perm ); + full_weights.push_back( weight ); + } + } + + // Optional: Project the points onto the z=0 plane (for further use, if necessary) + // for( auto& p : full_points ) { + // p[2] = 0.0; // Set z to 0 for all points + //} +} \ No newline at end of file diff --git a/src/quadratures/quadraturebase.cpp b/src/quadratures/quadraturebase.cpp index 7c5a94e7..dd312ca7 100644 --- a/src/quadratures/quadraturebase.cpp +++ b/src/quadratures/quadraturebase.cpp @@ -10,6 +10,7 @@ #include "quadratures/qmontecarlo.hpp" #include "quadratures/qproduct.hpp" #include "quadratures/qrectangular.hpp" +#include "quadratures/qsphericaltessalation.hpp" #include "toolboxes/errormessages.hpp" #include "toolboxes/textprocessingtoolbox.hpp" @@ -31,6 +32,7 @@ QuadratureBase* QuadratureBase::Create( Config* settings ) { case QUAD_GaussChebyshev1D: return new QGaussChebyshev1D( settings ); case QUAD_LevelSymmetric: return new QLevelSymmetric( settings ); case QUAD_LDFESA: return new QLDFESA( settings ); + case QUAD_SphericalTessalation2D: return new QSphericalTessalation( settings ); case QUAD_Lebedev: return new QLebedev( settings ); case QUAD_Product: return new QProduct( settings ); // case QUAD_Midpoint1D: return new QMidpoint1D( settings ); diff --git a/src/solvers/snsolver_hpc.cpp b/src/solvers/snsolver_hpc.cpp index 6430b268..807a5195 100644 --- a/src/solvers/snsolver_hpc.cpp +++ b/src/solvers/snsolver_hpc.cpp @@ -115,8 +115,13 @@ SNSolverHPC::SNSolverHPC( Config* settings ) { for( unsigned long idx_dim = 0; idx_dim < _nDim; idx_dim++ ) { _quadPts[Idx2D( idx_sys, idx_dim, _nDim )] = quadPoints[idx_sys + _startSysIdx][idx_dim]; } - _quadWeights[idx_sys] = - 2.0 * quadWeights[idx_sys + _startSysIdx]; // Rescaling of quadweights TODO: Check if this needs general refactoring + if( _settings->GetQuadName() == QUAD_GaussLegendreTensorized2D ) { + _quadWeights[idx_sys] = + 2.0 * quadWeights[idx_sys + _startSysIdx]; // Rescaling of quadweights TODO: Check if this needs general refactoring + } + else { + _quadWeights[idx_sys] = quadWeights[idx_sys + _startSysIdx]; // Rescaling of quadweights TODO: Check if this needs general refactoring} + } } #pragma omp parallel for diff --git a/tests/test_quadrature.cpp b/tests/test_quadrature.cpp index 89994ca0..330e2cb9 100644 --- a/tests/test_quadrature.cpp +++ b/tests/test_quadrature.cpp @@ -12,6 +12,7 @@ std::vector quadraturenames = { QUAD_MonteCarlo, QUAD_LevelSymmetric, QUAD_Lebedev, QUAD_LDFESA, + QUAD_SphericalTessalation2D, QUAD_Product, QUAD_Rectangular1D, QUAD_Rectangular2D, @@ -26,6 +27,7 @@ std::vector> quadratureorders = { { 3, 5, 7, 9, 11, 13, 15, 17, 19, 21, 23, 25, 27, 29, 31, 35, 41, 47, 53, 59, 65, 71, 77, 83, 89, 95, 101, 107, 113, 119, 125, 131 }, // Available orders for Lebedev { 1, 2, 3 }, // Available Orders for LDFESA + { 1, 2, 3, 4, 5 }, // Available Orders for SphericalTessalation { 4, 6, 8, 10 }, // Available Orders for Product { 60, 80, 100 }, // Rectangular 1D { 60, 80, 100 }, // Rectangular 2D @@ -119,6 +121,12 @@ TEST_CASE( "Quadrature Tests", "[quadrature]" ) { PrintErrorMsg( config, std::abs( Q->SumUpWeights() - 8.0 ), Q->SumUpWeights(), lowAccuracyTesting ); } } + else if( quadraturename == QUAD_SphericalTessalation2D ) { + if( !approxequal( Q->SumUpWeights(), 2 * M_PI, lowAccuracyTesting ) ) { + testPassed = false; + PrintErrorMsg( config, std::abs( Q->SumUpWeights() - 2 * M_PI ), Q->SumUpWeights(), lowAccuracyTesting ); + } + } else { if( !approxequal( Q->SumUpWeights(), 4 * M_PI, lowAccuracyTesting ) ) { testPassed = false; @@ -330,7 +338,8 @@ TEST_CASE( "Quadrature Tests", "[quadrature]" ) { PrintErrorMsg( config, std::abs( result - 0 ), result, lowAccuracyTesting ); } } - else if( quadraturename == QUAD_GaussLegendreTensorized2D || quadraturename == QUAD_Midpoint2D ) { + else if( quadraturename == QUAD_SphericalTessalation2D || quadraturename == QUAD_GaussLegendreTensorized2D || + quadraturename == QUAD_Midpoint2D ) { result = Q->Integrate( sin ); if( !approxequal( result, 0, lowAccuracyTesting ) ) { testPassed = false; @@ -459,9 +468,10 @@ TEST_CASE( "Quadrature Tests", "[quadrature]" ) { QuadratureBase* Q = QuadratureBase::Create( config ); // Note: Leaving out Quad_GaussLegendreTensorized with half weights... (to be added) - if( quadraturename != QUAD_GaussLegendreTensorized2D && quadraturename != QUAD_GaussLegendre1D && quadraturename != QUAD_MonteCarlo && - quadraturename != QUAD_Midpoint2D && quadraturename != QUAD_Midpoint1D && quadraturename != QUAD_Rectangular1D && - quadraturename != QUAD_Rectangular2D && quadraturename != QUAD_Rectangular3D ) // MonteCarlo is too low order... + if( quadraturename != QUAD_SphericalTessalation2D && quadraturename != QUAD_GaussLegendreTensorized2D && + quadraturename != QUAD_GaussLegendre1D && quadraturename != QUAD_MonteCarlo && quadraturename != QUAD_Midpoint2D && + quadraturename != QUAD_Midpoint1D && quadraturename != QUAD_Rectangular1D && quadraturename != QUAD_Rectangular2D && + quadraturename != QUAD_Rectangular3D ) // MonteCarlo is too low order... { if( quadraturename == QUAD_LevelSymmetric && quadratureorder == 20 ) continue; // Order 20 is somehow errorous result = Q->IntegrateSpherical( Omega_0 ); From 197e5309cfe6aa5bca80adc96cf7e1f985da87b4 Mon Sep 17 00:00:00 2001 From: ScSteffen Date: Tue, 15 Oct 2024 09:56:28 -0400 Subject: [PATCH 11/21] tidy code --- src/solvers/snsolver_hpc.cpp | 29 ++++++++++++++--------------- 1 file changed, 14 insertions(+), 15 deletions(-) diff --git a/src/solvers/snsolver_hpc.cpp b/src/solvers/snsolver_hpc.cpp index 807a5195..89b69347 100644 --- a/src/solvers/snsolver_hpc.cpp +++ b/src/solvers/snsolver_hpc.cpp @@ -321,21 +321,6 @@ void SNSolverHPC::Solve() { } } -void SNSolverHPC::PrintVolumeOutput( int idx_iter ) { - if( _settings->GetVolumeOutputFrequency() != 0 && idx_iter % (unsigned)_settings->GetVolumeOutputFrequency() == 0 ) { - WriteVolumeOutput( idx_iter ); - if( _rank == 0 ) { - ExportVTK( _settings->GetOutputFile() + "_" + std::to_string( idx_iter ), _outputFields, _outputFieldNames, _mesh ); // slow - } - } - if( idx_iter == (int)_nIter - 1 ) { // Last iteration write without suffix. - WriteVolumeOutput( idx_iter ); - if( _rank == 0 ) { - ExportVTK( _settings->GetOutputFile(), _outputFields, _outputFieldNames, _mesh ); - } - } -} - void SNSolverHPC::FluxOrder2() { double const eps = 1e-10; @@ -1296,6 +1281,20 @@ void SNSolverHPC::WriteVolumeOutput( unsigned idx_iter ) { } } +void SNSolverHPC::PrintVolumeOutput( int idx_iter ) { + if( _settings->GetVolumeOutputFrequency() != 0 && idx_iter % (unsigned)_settings->GetVolumeOutputFrequency() == 0 ) { + WriteVolumeOutput( idx_iter ); + if( _rank == 0 ) { + ExportVTK( _settings->GetOutputFile() + "_" + std::to_string( idx_iter ), _outputFields, _outputFieldNames, _mesh ); // slow + } + } + if( idx_iter == (int)_nIter - 1 ) { // Last iteration write without suffix. + WriteVolumeOutput( idx_iter ); + if( _rank == 0 ) { + ExportVTK( _settings->GetOutputFile(), _outputFields, _outputFieldNames, _mesh ); + } + } +} void SNSolverHPC::PrepareVolumeOutput() { unsigned nGroups = (unsigned)_settings->GetNVolumeOutput(); From 77bf64f917af295f221074842ba6bc194d301fe5 Mon Sep 17 00:00:00 2001 From: ScSteffen Date: Tue, 15 Oct 2024 14:29:30 -0400 Subject: [PATCH 12/21] changed tesslation implementation from recursive to direct --- include/quadratures/qsphericaltessalation.hpp | 5 +- src/quadratures/qsphericaltessalation.cpp | 132 +++++++++--------- src/solvers/snsolver_hpc.cpp | 1 + 3 files changed, 67 insertions(+), 71 deletions(-) diff --git a/include/quadratures/qsphericaltessalation.hpp b/include/quadratures/qsphericaltessalation.hpp index c7721d3e..02f4aaaa 100644 --- a/include/quadratures/qsphericaltessalation.hpp +++ b/include/quadratures/qsphericaltessalation.hpp @@ -20,9 +20,6 @@ class QSphericalTessalation : public QuadratureBase void SetConnectivity() override; private: - std::vector, 3>> subdivide_triangle( const std::array, 3>& triangle ); - std::vector, 3>> recursive_subdivide( const std::vector, 3>>& triangles, - int order ); std::array compute_centroid( const std::array, 3>& triangle ); std::array map_to_unit_sphere( const std::array& point ); double dot_product( const std::array& v1, const std::array& v2 ); @@ -33,6 +30,8 @@ class QSphericalTessalation : public QuadratureBase const std::vector& weights, std::vector>& full_points, std::vector& full_weights ); + + std::vector, 3>> generate_tessellation( const std::array, 3>& triangle, int order ); }; #endif // QSPHERICALTRIANGLE_H diff --git a/src/quadratures/qsphericaltessalation.cpp b/src/quadratures/qsphericaltessalation.cpp index 678c7438..a454086b 100644 --- a/src/quadratures/qsphericaltessalation.cpp +++ b/src/quadratures/qsphericaltessalation.cpp @@ -14,33 +14,72 @@ QSphericalTessalation::QSphericalTessalation( Config* settings ) : QuadratureBas _supportedDimensions = { 2 }; // Extension to 3d is straightforward } +std::vector, 3>> QSphericalTessalation::generate_tessellation( const std::array, 3>& triangle, + int order ) { + std::vector, 3>> triangles; + + // Vertices of the triangle + std::array A = triangle[0]; + std::array B = triangle[1]; + std::array C = triangle[2]; + + // Generate grid points along edges and inside the triangle + std::vector> grid_points; + for( int i = 0; i <= order; ++i ) { + for( int j = 0; j <= order - i; ++j ) { + double alpha = static_cast( i ) / order; + double beta = static_cast( j ) / order; + double gamma = 1.0 - alpha - beta; + std::array point = { + alpha * A[0] + beta * B[0] + gamma * C[0], alpha * A[1] + beta * B[1] + gamma * C[1], alpha * A[2] + beta * B[2] + gamma * C[2] }; + grid_points.push_back( point ); + } + } + + // Create triangles from the grid points + auto idx = [&]( int i, int j ) -> int { return i * ( order + 1 ) - ( i * ( i - 1 ) ) / 2 + j; }; + + for( int i = 0; i < order; ++i ) { + for( int j = 0; j < order - i; ++j ) { + std::array p1 = grid_points[idx( i, j )]; + std::array p2 = grid_points[idx( i + 1, j )]; + std::array p3 = grid_points[idx( i, j + 1 )]; + triangles.push_back( { p1, p2, p3 } ); + + if( j < order - i - 1 ) { + std::array p4 = grid_points[idx( i + 1, j + 1 )]; + triangles.push_back( { p2, p4, p3 } ); + } + } + } + + return triangles; +} + void QSphericalTessalation::SetPointsAndWeights() { // Define basal triangle vertices - std::array A = { 1.0, 0.0, 0.0 }; - std::array B = { 0.0, 1.0, 0.0 }; - std::array C = { 0.0, 0.0, 1.0 }; - std::array, 3> basal_triangle = { A, B, C }; + std::array A = { 1.0, 0.0, 0.0 }; + std::array B = { 0.0, 1.0, 0.0 }; + std::array C = { 0.0, 0.0, 1.0 }; - // Perform recursive subdivision - std::vector, 3>> final_triangles = recursive_subdivide( { basal_triangle }, _order ); + // Generate the tessellation for the given order + std::vector, 3>> final_triangles = generate_tessellation( { A, B, C }, _order ); - // Compute centroids and map to unit sphere + // Compute centroids of triangles and map them to the unit sphere std::vector> centroids; std::vector> mapped_points; std::vector weights; for( const auto& tri : final_triangles ) { + // Compute the centroid auto centroid = compute_centroid( tri ); auto mapped_centroid = map_to_unit_sphere( centroid ); - centroids.push_back( centroid ); - mapped_points.push_back( mapped_centroid ); + centroids.push_back( mapped_centroid ); - // Map triangle vertices to unit sphere - std::array a = map_to_unit_sphere( tri[0] ); - std::array b = map_to_unit_sphere( tri[1] ); - std::array c = map_to_unit_sphere( tri[2] ); - - // Calculate area using spherical excess + // Compute the area of the spherical triangle + auto a = map_to_unit_sphere( tri[0] ); + auto b = map_to_unit_sphere( tri[1] ); + auto c = map_to_unit_sphere( tri[2] ); double area = spherical_triangle_area( a, b, c ); weights.push_back( area ); } @@ -49,40 +88,36 @@ void QSphericalTessalation::SetPointsAndWeights() { std::vector> full_points; std::vector full_weights; - // Perform reflection and permutation - reflect_and_permute( mapped_points, weights, full_points, full_weights ); + // Perform reflection and permutation (as done in the original code) + reflect_and_permute( centroids, weights, full_points, full_weights ); + // Set the points and weights for the quadrature _nq = full_points.size(); _pointsKarth.resize( _nq ); _pointsSphere.resize( _nq ); _weights.resize( _nq ); + for( size_t i = 0; i < _nq; ++i ) { _pointsKarth[i] = { full_points[i][0], full_points[i][1], full_points[i][2] }; _weights[i] = full_weights[i]; } + double w = 0.0; - // Transform _points to _pointsSphere ==>transform (x,y,z) into (my,phi) for( unsigned idx = 0; idx < _nq; idx++ ) { - _pointsSphere[idx].resize( 3 ); // (my,phi) - _pointsSphere[idx][0] = _pointsKarth[idx][2]; // my = z - _pointsSphere[idx][1] = atan2( _pointsKarth[idx][1], _pointsKarth[idx][0] ); // phi in [-pi,pi] - _pointsSphere[idx][2] = 1.0; // radius r + _pointsSphere[idx].resize( 3 ); + _pointsSphere[idx][0] = _pointsKarth[idx][2]; + _pointsSphere[idx][1] = atan2( _pointsKarth[idx][1], _pointsKarth[idx][0] ); + _pointsSphere[idx][2] = 1.0; - // adapt intervall s.t. phi in [0,2pi] if( _pointsSphere[idx][1] < 0 ) { _pointsSphere[idx][1] = 2 * M_PI + _pointsSphere[idx][1]; } - // std::cout << _pointsKarth[idx][0] << " " << _pointsKarth[idx][1] << " " << _pointsKarth[idx][2] << std::endl; - // std::cout << _pointsSphere[idx][0] << " " << _pointsSphere[idx][1] << " " << _pointsSphere[idx][2] << std::endl; - // std::cout << _weights[idx] << std::endl; w += _weights[idx]; } - // std::cout << w << std::endl; - // exit( 1 ); } -void QSphericalTessalation::SetNq() { _nq = pow( GetOrder(), 4 ); } +void QSphericalTessalation::SetNq() { _nq = 4 * pow( GetOrder(), 2 ); } void QSphericalTessalation::SetConnectivity() { // TODO // Not initialized for this quadrature. @@ -92,45 +127,6 @@ void QSphericalTessalation::SetConnectivity() { // TODO bool QSphericalTessalation::CheckOrder() { return true; } -// Helper function to compute midpoint of two vectors -std::array QSphericalTessalation::midpoint( const std::array& p1, const std::array& p2 ) { - return { ( p1[0] + p2[0] ) / 2, ( p1[1] + p2[1] ) / 2, ( p1[2] + p2[2] ) / 2 }; -} - -// Subdivide a triangle into four smaller triangles -std::vector, 3>> QSphericalTessalation::subdivide_triangle( const std::array, 3>& triangle ) { - std::array A = triangle[0]; - std::array B = triangle[1]; - std::array C = triangle[2]; - - // Compute midpoints of edges - std::array AB_mid = midpoint( A, B ); - std::array BC_mid = midpoint( B, C ); - std::array CA_mid = midpoint( C, A ); - - // Define four new triangles - std::vector, 3>> new_triangles = { - { A, AB_mid, CA_mid }, { AB_mid, B, BC_mid }, { CA_mid, BC_mid, C }, { AB_mid, BC_mid, CA_mid } }; - - return new_triangles; -} - -// Recursively subdivide triangles to the desired order -std::vector, 3>> -QSphericalTessalation::recursive_subdivide( const std::vector, 3>>& triangles, int order ) { - if( order == 1 ) { - return triangles; - } - else { - std::vector, 3>> new_triangles; - for( const auto& triangle : triangles ) { - std::vector, 3>> subdivided = subdivide_triangle( triangle ); - new_triangles.insert( new_triangles.end(), subdivided.begin(), subdivided.end() ); - } - return recursive_subdivide( new_triangles, order - 1 ); - } -} - // Compute the centroid of a triangle std::array QSphericalTessalation::compute_centroid( const std::array, 3>& triangle ) { std::array A = triangle[0]; diff --git a/src/solvers/snsolver_hpc.cpp b/src/solvers/snsolver_hpc.cpp index 89b69347..c361bb51 100644 --- a/src/solvers/snsolver_hpc.cpp +++ b/src/solvers/snsolver_hpc.cpp @@ -1295,6 +1295,7 @@ void SNSolverHPC::PrintVolumeOutput( int idx_iter ) { } } } + void SNSolverHPC::PrepareVolumeOutput() { unsigned nGroups = (unsigned)_settings->GetNVolumeOutput(); From 430e4c54c92f2d559aaa9d9facc72ac89a0df4f1 Mon Sep 17 00:00:00 2001 From: Steffen Date: Tue, 15 Oct 2024 15:40:51 -0400 Subject: [PATCH 13/21] removed some auto calls --- src/quadratures/qsphericaltessalation.cpp | 51 ++++++++++++----------- 1 file changed, 26 insertions(+), 25 deletions(-) diff --git a/src/quadratures/qsphericaltessalation.cpp b/src/quadratures/qsphericaltessalation.cpp index a454086b..5f9e7b89 100644 --- a/src/quadratures/qsphericaltessalation.cpp +++ b/src/quadratures/qsphericaltessalation.cpp @@ -37,8 +37,12 @@ std::vector, 3>> QSphericalTessalation::generat } // Create triangles from the grid points - auto idx = [&]( int i, int j ) -> int { return i * ( order + 1 ) - ( i * ( i - 1 ) ) / 2 + j; }; - + auto idx = [&](int i, int j) -> int { + if (i < 0 || j < 0 || i > order || j > order - i) { + throw std::out_of_range("Index out of range in idx calculation"); + } + return i * (order + 1) - (i * (i - 1)) / 2 + j; + }; for( int i = 0; i < order; ++i ) { for( int j = 0; j < order - i; ++j ) { std::array p1 = grid_points[idx( i, j )]; @@ -52,7 +56,6 @@ std::vector, 3>> QSphericalTessalation::generat } } } - return triangles; } @@ -69,52 +72,50 @@ void QSphericalTessalation::SetPointsAndWeights() { std::vector> centroids; std::vector> mapped_points; std::vector weights; - for( const auto& tri : final_triangles ) { - // Compute the centroid auto centroid = compute_centroid( tri ); auto mapped_centroid = map_to_unit_sphere( centroid ); - centroids.push_back( mapped_centroid ); - - // Compute the area of the spherical triangle - auto a = map_to_unit_sphere( tri[0] ); - auto b = map_to_unit_sphere( tri[1] ); - auto c = map_to_unit_sphere( tri[2] ); + centroids.push_back( centroid ); + mapped_points.push_back( mapped_centroid ); + // Map triangle vertices to unit sphere + std::array a = map_to_unit_sphere( tri[0] ); + std::array b = map_to_unit_sphere( tri[1] ); + std::array c = map_to_unit_sphere( tri[2] ); + // Calculate area using spherical excess double area = spherical_triangle_area( a, b, c ); weights.push_back( area ); } - // Vectors to hold the full set of points and weights for the upper hemisphere std::vector> full_points; std::vector full_weights; - - // Perform reflection and permutation (as done in the original code) - reflect_and_permute( centroids, weights, full_points, full_weights ); - - // Set the points and weights for the quadrature + // Perform reflection and permutation + reflect_and_permute( mapped_points, weights, full_points, full_weights ); _nq = full_points.size(); _pointsKarth.resize( _nq ); _pointsSphere.resize( _nq ); _weights.resize( _nq ); - for( size_t i = 0; i < _nq; ++i ) { _pointsKarth[i] = { full_points[i][0], full_points[i][1], full_points[i][2] }; _weights[i] = full_weights[i]; } - double w = 0.0; + // Transform _points to _pointsSphere ==>transform (x,y,z) into (my,phi) for( unsigned idx = 0; idx < _nq; idx++ ) { - _pointsSphere[idx].resize( 3 ); - _pointsSphere[idx][0] = _pointsKarth[idx][2]; - _pointsSphere[idx][1] = atan2( _pointsKarth[idx][1], _pointsKarth[idx][0] ); - _pointsSphere[idx][2] = 1.0; - + _pointsSphere[idx].resize( 3 ); // (my,phi) + _pointsSphere[idx][0] = _pointsKarth[idx][2]; // my = z + _pointsSphere[idx][1] = atan2( _pointsKarth[idx][1], _pointsKarth[idx][0] ); // phi in [-pi,pi] + _pointsSphere[idx][2] = 1.0; // radius r + // adapt intervall s.t. phi in [0,2pi] if( _pointsSphere[idx][1] < 0 ) { _pointsSphere[idx][1] = 2 * M_PI + _pointsSphere[idx][1]; } - + // std::cout << _pointsKarth[idx][0] << " " << _pointsKarth[idx][1] << " " << _pointsKarth[idx][2] << std::endl; + // std::cout << _pointsSphere[idx][0] << " " << _pointsSphere[idx][1] << " " << _pointsSphere[idx][2] << std::endl; + // std::cout << _weights[idx] << std::endl; w += _weights[idx]; } + // std::cout << w << std::endl; + // exit( 1 ); } void QSphericalTessalation::SetNq() { _nq = 4 * pow( GetOrder(), 2 ); } From ed277b470b909f069755e250ad5c9c53f1fbf078 Mon Sep 17 00:00:00 2001 From: Steffen Date: Thu, 17 Oct 2024 10:27:54 -0400 Subject: [PATCH 14/21] added restart files for SN hpc solver --- include/common/config.hpp | 10 ++-- include/common/io.hpp | 6 +++ include/solvers/snsolver_hpc.hpp | 2 +- src/common/config.cpp | 4 ++ src/common/io.cpp | 89 ++++++++++++++++++++++++++++++-- src/solvers/snsolver_hpc.cpp | 23 ++++++--- 6 files changed, 118 insertions(+), 16 deletions(-) diff --git a/include/common/config.hpp b/include/common/config.hpp index 6957abe0..42e1992d 100644 --- a/include/common/config.hpp +++ b/include/common/config.hpp @@ -48,9 +48,11 @@ class Config // std::vector _1dIntegrationBounds; /*!< @brief Quadrature Order*/ // Mesh - unsigned _nCells; /*!< @brief Number of cells in the mesh */ - unsigned short _dim; /*!< @brief spatial dimensionality of the mesh/test case */ - bool _forcedConnectivityWrite; /*!< @brief If true, the meshconnectivity is always computed and written to .con file */ + unsigned _nCells; /*!< @brief Number of cells in the mesh */ + unsigned short _dim; /*!< @brief spatial dimensionality of the mesh/test case */ + bool _forcedConnectivityWrite; /*!< @brief If true, the meshconnectivity is always computed and written to .con file */ + bool _loadrestartSolution; /*!< @brief If true, the simulation loads a restart solution from file */ + unsigned long _saveRestartSolutionFrequency; /*!< @brief Frequency for saving restart solution to file */ // Boundary Conditions /*!< @brief List of all Pairs (marker, BOUNDARY_TYPE), e.g. (farfield,DIRICHLET). @@ -312,6 +314,8 @@ class Config unsigned GetNCells() const { return _nCells; } unsigned short GetDim() const { return _dim; } bool inline GetForcedConnectivity() const { return _forcedConnectivityWrite; } + bool inline GetLoadRestartSolution() const { return _loadrestartSolution; } + unsigned long inline GetSaveRestartSolutionFrequency() const { return _saveRestartSolutionFrequency; } // Solver Structure bool inline GetHPC() const { return _HPC; } diff --git a/include/common/io.hpp b/include/common/io.hpp index 377c1332..0c302183 100644 --- a/include/common/io.hpp +++ b/include/common/io.hpp @@ -43,6 +43,12 @@ std::string ParseArguments( int argc, char* argv[] ); void PrintLogHeader( std::string inputFile ); +void WriteRestartSolution( + const std::string& baseOutputFile, const std::vector& solution, const std::vector& scalarFlux, int rank, int idx_iter ); + +int LoadRestartSolution( + const std::string& baseInputFile, std::vector& solution, std::vector& scalarFlux, int rank, unsigned long nCells ); + // Matrix createSU2MeshFromImage( std::string imageName, std::string SU2Filename ); Deprecated #endif // IO_H diff --git a/include/solvers/snsolver_hpc.hpp b/include/solvers/snsolver_hpc.hpp index f37bb98b..7f25b882 100644 --- a/include/solvers/snsolver_hpc.hpp +++ b/include/solvers/snsolver_hpc.hpp @@ -32,7 +32,7 @@ class SNSolverHPC // Time unsigned long _nIter; /*!< @brief number of time steps, for non CSD, this equals _nEnergies, for _csd, _maxIter =_nEnergies-1*/ double _dT; /*!< @brief energy/time step size */ - + int _idx_start_iter; /*!< @brief index of first iteration */ // Mesh related members, memory optimized unsigned long _nCells; /*!< @brief number of spatial cells */ unsigned long _nSys; /*!< @brief number of equations in the transport system, i.e. num quad pts */ diff --git a/src/common/config.cpp b/src/common/config.cpp index 60ae732a..9ad9b107 100644 --- a/src/common/config.cpp +++ b/src/common/config.cpp @@ -226,6 +226,10 @@ void Config::SetConfigOptions() { /*! @brief FORCE_CONNECTIVITY_RECOMPUTE \n DESCRIPTION:If true, mesh recomputes connectivity instead of loading from file \n DEFAULT false * \ingroup Config.*/ AddBoolOption( "FORCE_CONNECTIVITY_RECOMPUTE", _forcedConnectivityWrite, false ); + /*! @brief RESTART_SOLUTION \n DESCRIPTION:If true, simulation loads a restart solution from file \n DEFAULT false + * \ingroup Config.*/ + AddBoolOption( "LOAD_RESTART_SOLUTION", _loadrestartSolution, false ); + AddUnsignedLongOption( "SAVE_RESTART_SOLUTION_FREQUENCY", _saveRestartSolutionFrequency, false ); // Quadrature relatated options /*! @brief QUAD_TYPE \n DESCRIPTION: Type of Quadrature rule \n Options: see @link QUAD_NAME \endlink \n DEFAULT: QUAD_MonteCarlo diff --git a/src/common/io.cpp b/src/common/io.cpp index 8020d396..cdc8a307 100644 --- a/src/common/io.cpp +++ b/src/common/io.cpp @@ -11,9 +11,11 @@ #include "toolboxes/errormessages.hpp" #include "toolboxes/textprocessingtoolbox.hpp" +#include #include #include #include +#include #ifdef IMPORT_MPI #include @@ -356,7 +358,7 @@ void LoadConnectivityFromFile( const std::string inputFile, for( unsigned j = 0; j < correctedNodesPerCell; ++j ) { std::getline( iss, line, ',' ); std::istringstream converter( line ); - converter >> std::fixed >> setprecision( 12 ) >> cellNeighbors[i][j]; + converter >> std::fixed >> setprecision( 15 ) >> cellNeighbors[i][j]; } // Load cellInterfaceMidPoints @@ -366,8 +368,8 @@ void LoadConnectivityFromFile( const std::string inputFile, for( unsigned k = 0; k < nDim; ++k ) { std::getline( iss, line, ',' ); std::istringstream converter( line ); - converter >> std::fixed >> setprecision( 12 ) >> cellInterfaceMidPoints[i][j][k]; // Replace with appropriate member of Vector - // std::cout << std::fixed << setprecision( 12 ) << cellInterfaceMidPoints[i][j][k] << std::endl; + converter >> std::fixed >> setprecision( 15 ) >> cellInterfaceMidPoints[i][j][k]; // Replace with appropriate member of Vector + // std::cout << std::fixed << setprecision( 15 ) << cellInterfaceMidPoints[i][j][k] << std::endl; } } // Load cellNormals @@ -377,7 +379,7 @@ void LoadConnectivityFromFile( const std::string inputFile, for( unsigned k = 0; k < nDim; ++k ) { std::getline( iss, line, ',' ); std::istringstream converter( line ); - converter >> std::fixed >> setprecision( 12 ) >> cellNormals[i][j][k]; // Replace with appropriate member of Vector + converter >> std::fixed >> setprecision( 15 ) >> cellNormals[i][j][k]; // Replace with appropriate member of Vector } } // Load cellBoundaryTypes @@ -404,7 +406,7 @@ void WriteConnecitivityToFile( const std::string outputFile, // cellBoundaryTypes (1 element), (tranlated from BOUNDARY_TYPE to unsigned) std::ofstream outFile( outputFile ); - outFile << std::fixed << setprecision( 12 ); + outFile << std::fixed << setprecision( 15 ); // const std::size_t bufferSize = 10000; // Adjust as needed // outFile.rdbuf()->pubsetbuf( 0, bufferSize ); if( outFile.is_open() ) { @@ -441,6 +443,83 @@ void WriteConnecitivityToFile( const std::string outputFile, } } +void WriteRestartSolution( + const std::string& baseOutputFile, const std::vector& solution, const std::vector& scalarFlux, int rank, int idx_iter ) { + // Generate filename with rank number + std::string outputFile = baseOutputFile + "_restart_rank_" + std::to_string( rank ); + + // Check if the file exists and delete it + if( std::filesystem::exists( outputFile ) ) { + std::filesystem::remove( outputFile ); + } + + // Open the file for writing (automatically creates a new file) + std::ofstream outFile( outputFile, std::ios::out ); + if( !outFile ) { + ErrorMessages::Error( "Error opening restart solution file.", CURRENT_FUNCTION ); + return; + } + + // Write the iteration number as the first line + outFile << idx_iter << "\n"; + outFile << std::fixed << std::setprecision( 15 ); + + // Write each element of the solution to the file on a new line + for( const double& value : solution ) { + outFile << value << "\n"; + } + + // Append each element of the scalarFlux to the file on a new line + for( const double& fluxValue : scalarFlux ) { + outFile << fluxValue << "\n"; + } + + // Close the file + outFile.close(); +} + +int LoadRestartSolution( + const std::string& baseInputFile, std::vector& solution, std::vector& scalarFlux, int rank, unsigned long nCells ) { + // Generate filename with rank number + std::string inputFile = baseInputFile + "_restart_rank_" + std::to_string( rank ); + + // Open the file for reading + std::ifstream inFile( inputFile ); + if( !inFile ) { + std::cerr << "Error opening restart solution file for reading." << std::endl; + return 0; // Signal failure + } + + // Read the iteration number from the first line + std::string line; + std::getline( inFile, line ); + int idx_iter = std::stoi( line ); // Convert the line to an integer + + // Read data into a temporary vector + std::vector tempData; + double value; + while( std::getline( inFile, line ) ) { + std::istringstream converter( line ); + converter >> std::fixed >> std::setprecision( 15 ) >> value; + tempData.push_back( value ); + } + + // Close the file + inFile.close(); + + // Verify that we have at least nCells entries to populate scalarFlux + if( tempData.size() < nCells ) { + std::cerr << "Not enough data to populate scalar flux vector." << std::endl; + return 0; // Signal failure + } + + // Allocate the last nCells entries to scalarFlux and the rest to solution + solution.assign( tempData.begin(), tempData.end() - nCells ); + scalarFlux.assign( tempData.end() - nCells, tempData.end() ); + + return idx_iter; // Return the iteration index +} + std::string ParseArguments( int argc, char* argv[] ) { std::string inputFile; std::string usage_help = "\nUsage: " + std::string( argv[0] ) + " inputfile\n"; diff --git a/src/solvers/snsolver_hpc.cpp b/src/solvers/snsolver_hpc.cpp index c361bb51..ac5759e9 100644 --- a/src/solvers/snsolver_hpc.cpp +++ b/src/solvers/snsolver_hpc.cpp @@ -21,9 +21,9 @@ SNSolverHPC::SNSolverHPC( Config* settings ) { _rank = 0; #endif - _settings = settings; - _currTime = 0.0; - + _settings = settings; + _currTime = 0.0; + _idx_start_iter = 0; _nOutputMoments = 2; // Currently only u_1 (x direction) and u_1 (y direction) are supported // Create Mesh _mesh = LoadSU2MeshFromFile( settings ); @@ -159,6 +159,9 @@ SNSolverHPC::SNSolverHPC( Config* settings ) { // _mass += _scalarFlux[idx_cell] * _areas[idx_cell]; } + if( _settings->GetLoadRestartSolution() ) + _idx_start_iter = LoadRestartSolution( _settings->GetOutputFile(), _sol, _scalarFlux, _rank, _nCells ) + 1; + #ifdef IMPORT_MPI MPI_Barrier( MPI_COMM_WORLD ); #endif @@ -265,7 +268,7 @@ void SNSolverHPC::Solve() { std::chrono::duration duration; // Loop over energies (pseudo-time of continuous slowing down approach) - for( unsigned iter = 0; iter < _nIter; iter++ ) { + for( unsigned iter = (unsigned)_idx_start_iter; iter < _nIter; iter++ ) { if( iter == _nIter - 1 ) { // last iteration _dT = _settings->GetTEnd() - iter * _dT; @@ -310,6 +313,7 @@ void SNSolverHPC::Solve() { #ifdef BUILD_MPI MPI_Barrier( MPI_COMM_WORLD ); #endif + PrintVolumeOutput( iter ); #ifdef BUILD_MPI MPI_Barrier( MPI_COMM_WORLD ); @@ -352,7 +356,7 @@ void SNSolverHPC::FluxOrder2() { for( unsigned long idx_cell = 0; idx_cell < _nCells; ++idx_cell ) { // Compute Limiter if( _cellBoundaryTypes[idx_cell] == BOUNDARY_TYPE::NONE ) { // skip ghost cells - // #pragma omp simd +#pragma omp simd for( unsigned long idx_sys = 0; idx_sys < _localNSys; idx_sys++ ) { double gaussPoint = 0; @@ -400,7 +404,7 @@ void SNSolverHPC::FluxOrder2() { } } else { - // #pragma omp simd +#pragma omp simd for( unsigned long idx_sys = 0; idx_sys < _localNSys; idx_sys++ ) { _limiter[Idx2D( idx_cell, idx_sys, _localNSys )] = 0.; // limiter should be zero at boundary _solDx[Idx3D( idx_cell, idx_sys, 0, _localNSys, _nDim )] = 0.; @@ -1282,7 +1286,12 @@ void SNSolverHPC::WriteVolumeOutput( unsigned idx_iter ) { } void SNSolverHPC::PrintVolumeOutput( int idx_iter ) { - if( _settings->GetVolumeOutputFrequency() != 0 && idx_iter % (unsigned)_settings->GetVolumeOutputFrequency() == 0 ) { + if( _settings->GetSaveRestartSolutionFrequency() != 0 && idx_iter % (int)_settings->GetSaveRestartSolutionFrequency() == 0 ) { + std::cout << "Saving restart solution at iteration " << idx_iter << std::endl; + WriteRestartSolution( _settings->GetOutputFile(), _sol, _scalarFlux, _rank, idx_iter ); + } + + if( _settings->GetVolumeOutputFrequency() != 0 && idx_iter % (int)_settings->GetVolumeOutputFrequency() == 0 ) { WriteVolumeOutput( idx_iter ); if( _rank == 0 ) { ExportVTK( _settings->GetOutputFile() + "_" + std::to_string( idx_iter ), _outputFields, _outputFieldNames, _mesh ); // slow From 3e23828f19f076c59364db0e7af0ae1fb181fcd9 Mon Sep 17 00:00:00 2001 From: ScSteffen Date: Fri, 18 Oct 2024 17:47:36 -0400 Subject: [PATCH 15/21] changed green segment absorption output --- include/common/globalconstants.hpp | 6 +- include/common/mesh.hpp | 6 +- include/solvers/snsolver_hpc.hpp | 12 ++- src/common/config.cpp | 62 +++++++------ src/common/mesh.cpp | 39 ++++++++ src/solvers/snsolver_hpc.cpp | 141 ++++++++++++++++++++++++----- src/solvers/solverbase.cpp | 7 +- 7 files changed, 216 insertions(+), 57 deletions(-) diff --git a/include/common/globalconstants.hpp b/include/common/globalconstants.hpp index debd3830..b0ca20e3 100644 --- a/include/common/globalconstants.hpp +++ b/include/common/globalconstants.hpp @@ -197,7 +197,8 @@ enum SCALAR_OUTPUT { TOTAL_PARTICLE_ABSORPTION_HORIZONTAL, PROBE_MOMENT_TIME_TRACE, VAR_ABSORPTION_GREEN, - VAR_ABSORPTION_GREEN_LINE + ABSORPTION_GREEN_BLOCK, + ABSORPTION_GREEN_LINE }; inline std::map ScalarOutput_Map{ { "ITER", ITER }, @@ -221,7 +222,8 @@ inline std::map ScalarOutput_Map{ { "ITER", ITER }, { "TOTAL_PARTICLE_ABSORPTION_HORIZONTAL", TOTAL_PARTICLE_ABSORPTION_HORIZONTAL }, { "PROBE_MOMENT_TIME_TRACE", PROBE_MOMENT_TIME_TRACE }, { "VAR_ABSORPTION_GREEN", VAR_ABSORPTION_GREEN }, - { "VAR_ABSORPTION_GREEN_LINE", VAR_ABSORPTION_GREEN_LINE } }; + { "ABSORPTION_GREEN_BLOCK", ABSORPTION_GREEN_BLOCK }, + { "ABSORPTION_GREEN_LINE", ABSORPTION_GREEN_LINE } }; // Spherical Basis Name enum SPHERICAL_BASIS_NAME { SPHERICAL_HARMONICS, SPHERICAL_MONOMIALS, SPHERICAL_MONOMIALS_ROTATED }; diff --git a/include/common/mesh.hpp b/include/common/mesh.hpp index 04c5ff7e..92f42f41 100644 --- a/include/common/mesh.hpp +++ b/include/common/mesh.hpp @@ -129,10 +129,14 @@ class Mesh * @return cell_idx: unsigned */ unsigned GetCellOfKoordinate( const double x, const double y ) const; - /*! @brief Returns index of cells contained in the ball around the coordinate (x,y) with radius r + /*! @brief Returns index of cells contained in the ball around the coordinate (x,y) with radius r * @return cell_idxs: std::vector */ std::vector GetCellsofBall( const double x, const double y, const double r ) const; + /*! @brief Returns index of cells contained in the rectangle with specified corner coordinates/* + * @return cell_idxs: std::vector */ + std::vector GetCellsofRectangle( const std::vector>& cornercoordinates ) const; + /*! @brief ComputeSlopes calculates the slope in every cell into x and y direction using the divergence theorem. * @param nq is number of quadrature points * @param psiDerX is slope in x direction (gets computed. Slope is stored here) diff --git a/include/solvers/snsolver_hpc.hpp b/include/solvers/snsolver_hpc.hpp index 7f25b882..4bdb32ed 100644 --- a/include/solvers/snsolver_hpc.hpp +++ b/include/solvers/snsolver_hpc.hpp @@ -111,14 +111,18 @@ class SNSolverHPC double _curAbsorptionHohlraumVertical; double _curAbsorptionHohlraumHorizontal; double _varAbsorptionHohlraumGreen; + std::vector> _probingCellsHohlraum; /*!< @brief Indices of cells that contain a probing sensor */ std::vector _probingMoments; /*!< @brief Solution Momnets at the probing cells that contain a probing sensor */ unsigned _probingMomentsTimeIntervals; /*!< @brief Solution Momnets at the probing cells that contain a probing sensor */ - unsigned _nProbingCellsLineGreen; /*!< @brief Number of sampling cells that contain a probing sensor for the sliding window */ - std::vector _probingCellsLineGreen; /*!< @brief Indices of cells that contain a probing sensor for the sliding window */ - std::vector _absorptionValsIntegrated; /*!< @brief Avg Absorption value at the sampleing points of lineGreen */ - std::vector _varAbsorptionValsIntegrated; /*!< @brief Var in Avg Absorption value at the sampleing points of lineGreen */ + unsigned _nProbingCellsLineGreen; /*!< @brief Number of sampling cells that contain a probing sensor for the sliding window */ + std::vector _probingCellsLineGreen; /*!< @brief Indices of cells that contain a probing sensor for the sliding window */ + std::vector _absorptionValsLineSegment; /*!< @brief Avg Absorption value at the sampleing points of lineGreen */ + + unsigned _nProbingCellsBlocksGreen; + std::vector> _probingCellsBlocksGreen; /*!< @brief Indices of cells that contain a probing sensor blocks */ + std::vector _absorptionValsBlocksGreen; /*!< @brief Avg Absorption value at the sampleing blocks of lineGreen */ // Design parameters std::vector _cornerUpperLeftGreen; /*!< @brief Coord of corner of the green area (minus thickness/2 of it) relative to the green center */ diff --git a/src/common/config.cpp b/src/common/config.cpp index 9ad9b107..4e189192 100644 --- a/src/common/config.cpp +++ b/src/common/config.cpp @@ -788,21 +788,23 @@ void Config::SetPostprocessing() { break; case PROBLEM_QuarterHohlraum: case PROBLEM_SymmetricHohlraum: - legalOutputs = { ITER, - WALL_TIME, - MASS, - RMS_FLUX, - VTK_OUTPUT, - CSV_OUTPUT, - CUR_OUTFLOW, - TOTAL_OUTFLOW, - MAX_OUTFLOW, - TOTAL_PARTICLE_ABSORPTION_CENTER, - TOTAL_PARTICLE_ABSORPTION_VERTICAL, - TOTAL_PARTICLE_ABSORPTION_HORIZONTAL, - TOTAL_PARTICLE_ABSORPTION, - PROBE_MOMENT_TIME_TRACE, - VAR_ABSORPTION_GREEN }; + legalOutputs = { + ITER, + WALL_TIME, + MASS, + RMS_FLUX, + VTK_OUTPUT, + CSV_OUTPUT, + CUR_OUTFLOW, + TOTAL_OUTFLOW, + MAX_OUTFLOW, + TOTAL_PARTICLE_ABSORPTION_CENTER, + TOTAL_PARTICLE_ABSORPTION_VERTICAL, + TOTAL_PARTICLE_ABSORPTION_HORIZONTAL, + TOTAL_PARTICLE_ABSORPTION, + PROBE_MOMENT_TIME_TRACE, + VAR_ABSORPTION_GREEN, + }; it = std::find( legalOutputs.begin(), legalOutputs.end(), _screenOutput[idx_screenOutput] ); @@ -898,7 +900,7 @@ void Config::SetPostprocessing() { } // Check if the choice of history output is compatible to the problem - for( unsigned short idx_screenOutput = 0; idx_screenOutput < _nHistoryOutput; idx_screenOutput++ ) { + for( unsigned short idx_historyOutput = 0; idx_historyOutput < _nHistoryOutput; idx_historyOutput++ ) { std::vector legalOutputs; std::vector::iterator it; @@ -921,9 +923,9 @@ void Config::SetPostprocessing() { CUR_PARTICLE_ABSORPTION, TOTAL_PARTICLE_ABSORPTION, MAX_PARTICLE_ABSORPTION }; - it = std::find( legalOutputs.begin(), legalOutputs.end(), _historyOutput[idx_screenOutput] ); + it = std::find( legalOutputs.begin(), legalOutputs.end(), _historyOutput[idx_historyOutput] ); if( it == legalOutputs.end() ) { - std::string foundKey = findKey( ScalarOutput_Map, _historyOutput[idx_screenOutput] ); + std::string foundKey = findKey( ScalarOutput_Map, _historyOutput[idx_historyOutput] ); ErrorMessages::Error( "Illegal output field <" + foundKey + "> for option HISTORY_OUTPUT for this test case.\n" @@ -950,18 +952,19 @@ void Config::SetPostprocessing() { TOTAL_PARTICLE_ABSORPTION, PROBE_MOMENT_TIME_TRACE, VAR_ABSORPTION_GREEN, - VAR_ABSORPTION_GREEN_LINE }; + ABSORPTION_GREEN_BLOCK, + ABSORPTION_GREEN_LINE }; - it = std::find( legalOutputs.begin(), legalOutputs.end(), _historyOutput[idx_screenOutput] ); + it = std::find( legalOutputs.begin(), legalOutputs.end(), _historyOutput[idx_historyOutput] ); if( it == legalOutputs.end() ) { - std::string foundKey = findKey( ScalarOutput_Map, _historyOutput[idx_screenOutput] ); + std::string foundKey = findKey( ScalarOutput_Map, _historyOutput[idx_historyOutput] ); ErrorMessages::Error( "Illegal output field <" + foundKey + "> for option HISTORY_OUTPUT for this test case.\n" "Supported fields are: ITER, MASS, RMS_FLUX, VTK_OUTPUT, CSV_OUTPUT, TOTAL_PARTICLE_ABSORPTION_CENTER, \n " "TOTAL_PARTICLE_ABSORPTION_VERTICAL, TOTAL_PARTICLE_ABSORPTION_HORIZONTAL,PROBE_MOMENT_TIME_TRACE, CUR_OUTFLOW, \n" - "TOTAL_OUTFLOW, MAX_OUTFLOW , VAR_ABSORPTION_GREEN, VAR_ABSORPTION_GREEN_LINE \n" + "TOTAL_OUTFLOW, MAX_OUTFLOW , VAR_ABSORPTION_GREEN, ABSORPTION_GREEN_BLOCK, ABSORPTION_GREEN_LINE \n" "Please check your .cfg file.", CURRENT_FUNCTION ); } @@ -969,10 +972,10 @@ void Config::SetPostprocessing() { default: legalOutputs = { ITER, WALL_TIME, MASS, RMS_FLUX, VTK_OUTPUT, CSV_OUTPUT, CUR_OUTFLOW, TOTAL_OUTFLOW, MAX_OUTFLOW }; - it = std::find( legalOutputs.begin(), legalOutputs.end(), _historyOutput[idx_screenOutput] ); + it = std::find( legalOutputs.begin(), legalOutputs.end(), _historyOutput[idx_historyOutput] ); if( it == legalOutputs.end() ) { - std::string foundKey = findKey( ScalarOutput_Map, _historyOutput[idx_screenOutput] ); + std::string foundKey = findKey( ScalarOutput_Map, _historyOutput[idx_historyOutput] ); ErrorMessages::Error( "Illegal output field <" + foundKey + "> for option SCREEN_OUTPUT for this test case.\n" @@ -1023,11 +1026,18 @@ void Config::SetPostprocessing() { if( _problemName == PROBLEM_SymmetricHohlraum || _problemName == PROBLEM_QuarterHohlraum ) { std::vector::iterator it; - it = find( _historyOutput.begin(), _historyOutput.end(), VAR_ABSORPTION_GREEN_LINE ); + it = find( _historyOutput.begin(), _historyOutput.end(), ABSORPTION_GREEN_LINE ); if( it != _historyOutput.end() ) { _historyOutput.erase( it ); _nHistoryOutput += _nProbingCellsLineGreenHohlraum - 1; // extend the screen output by the number of probing points - for( unsigned i = 0; i < _nProbingCellsLineGreenHohlraum; i++ ) _historyOutput.push_back( VAR_ABSORPTION_GREEN_LINE ); + for( unsigned i = 0; i < _nProbingCellsLineGreenHohlraum; i++ ) _historyOutput.push_back( ABSORPTION_GREEN_LINE ); + } + + it = find( _historyOutput.begin(), _historyOutput.end(), ABSORPTION_GREEN_BLOCK ); + if( it != _historyOutput.end() ) { + _historyOutput.erase( it ); + _nHistoryOutput += 44 - 1; // extend the screen output by the number of probing points + for( unsigned i = 0; i < 44; i++ ) _historyOutput.push_back( ABSORPTION_GREEN_BLOCK ); } } } diff --git a/src/common/mesh.cpp b/src/common/mesh.cpp index 02935194..ed7ef3c3 100644 --- a/src/common/mesh.cpp +++ b/src/common/mesh.cpp @@ -337,6 +337,7 @@ _cellMidPoints[i] ); interfaceMidFlatPart[pos] = ComputeCellInterfaceMidpoints( } } */ + // assign boundary types to all cells _cellBoundaryTypes.resize( _numCells, BOUNDARY_TYPE::NONE ); #pragma omp parallel for @@ -632,6 +633,44 @@ std::vector Mesh::GetCellsofBall( const double x, const double y, cons return cells_in_ball; } +std::vector Mesh::GetCellsofRectangle( const std::vector>& cornercoordinates ) const { + + std::vector cells_in_rectangle; + + // Compute x_min, x_max, y_min, y_max from the corner coordinates + double x_min = cornercoordinates[0][0]; + double x_max = cornercoordinates[0][0]; + double y_min = cornercoordinates[0][1]; + double y_max = cornercoordinates[0][1]; + + // Loop over the cornercoordinates to find the minimum and maximum x and y values + for( const auto& corner : cornercoordinates ) { + x_min = std::min( x_min, corner[0] ); + x_max = std::max( x_max, corner[0] ); + y_min = std::min( y_min, corner[1] ); + y_max = std::max( y_max, corner[1] ); + } + // Loop over all cells and check if their midpoints fall within the rectangle +#pragma omp parallel for + for( unsigned idx_cell = 0; idx_cell < _numCells; ++idx_cell ) { + double cell_x = _cellMidPoints[idx_cell][0]; + double cell_y = _cellMidPoints[idx_cell][1]; + + // Check if the cell's midpoint is inside the rectangle + if( cell_x >= x_min && cell_x <= x_max && cell_y >= y_min && cell_y <= y_max ) { +#pragma omp critical + { cells_in_rectangle.push_back( idx_cell ); } + } + } + if( cells_in_rectangle.empty() ) { // take the only cell that contains the point + std::cout << "No cells found within rectangle centered at defined by corner coordinates. " << std::endl; + std::cout << "Taking the only cell that contains the point." << std::endl; + cells_in_rectangle.push_back( GetCellOfKoordinate( x_min + x_max / 2.0, y_min + y_max / 2.0 ) ); + } + + return cells_in_rectangle; +} + bool Mesh::IsPointInsideCell( unsigned idx_cell, double x, double y ) const { bool inside = false; if( _numNodesPerCell == 3 ) { diff --git a/src/solvers/snsolver_hpc.cpp b/src/solvers/snsolver_hpc.cpp index ac5759e9..6ebcff27 100644 --- a/src/solvers/snsolver_hpc.cpp +++ b/src/solvers/snsolver_hpc.cpp @@ -206,7 +206,6 @@ SNSolverHPC::SNSolverHPC( Config* settings ) { unsigned n_probes = 4; if( _settings->GetProblemName() == PROBLEM_SymmetricHohlraum ) n_probes = 4; if( _settings->GetProblemName() == PROBLEM_QuarterHohlraum ) n_probes = 2; - //_probingMomentsTimeIntervals = 10; _probingMoments = std::vector( n_probes * 3, 0. ); // 10 time horizons if( _settings->GetProblemName() == PROBLEM_SymmetricHohlraum ) { @@ -244,10 +243,11 @@ SNSolverHPC::SNSolverHPC( Config* settings ) { _nProbingCellsLineGreen = _settings->GetNumProbingCellsLineHohlraum(); - SetProbingCellsLineGreen(); // ONLY FOR HOHLRAUM + _nProbingCellsBlocksGreen = 44; + _absorptionValsBlocksGreen = std::vector( _nProbingCellsBlocksGreen, 0. ); + _absorptionValsLineSegment = std::vector( _nProbingCellsLineGreen, 0.0 ); - _absorptionValsIntegrated = std::vector( _nProbingCellsLineGreen, 0.0 ); - _varAbsorptionValsIntegrated = std::vector( _nProbingCellsLineGreen, 0.0 ); + SetProbingCellsLineGreen(); // ONLY FOR HOHLRAUM } } @@ -610,7 +610,8 @@ void SNSolverHPC::IterPostprocessing() { bool green4 = x > 0.15 + _settings->GetPosXCenterGreenHohlraum() && x < 0.2 + _settings->GetPosXCenterGreenHohlraum() && y > -0.35 + _settings->GetPosYCenterGreenHohlraum() && y < 0.35 + _settings->GetPosYCenterGreenHohlraum(); if( green1 || green2 || green3 || green4 ) { - a_g += ( _sigmaT[idx_cell] - _sigmaS[idx_cell] ) * _scalarFlux[idx_cell] * _areas[idx_cell]; + a_g += ( _sigmaT[idx_cell] - _sigmaS[idx_cell] ) * _scalarFlux[idx_cell] * _areas[idx_cell] / + ( 44 * 0.05 * 0.05 ); // divisor is area of the green } } @@ -759,8 +760,6 @@ void SNSolverHPC::IterPostprocessing() { } } - // probe values green - ComputeQOIsGreenProbingLine(); #ifdef IMPORT_MPI MPI_Barrier( MPI_COMM_WORLD ); MPI_Allreduce( temp_probingMoments.data(), _probingMoments.data(), 3 * n_probes, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD ); @@ -774,7 +773,8 @@ void SNSolverHPC::IterPostprocessing() { } #endif } - + // probe values green + ComputeQOIsGreenProbingLine(); // Update time integral values on rank 0 if( _rank == 0 ) { _totalScalarOutflow += _curScalarOutflow * _dT; @@ -889,6 +889,7 @@ void SNSolverHPC::PrepareScreenOutput() { } break; case VAR_ABSORPTION_GREEN: _screenOutputFieldNames[idx_field] = "Var. absorption green"; break; + default: ErrorMessages::Error( "Screen output field not defined!", CURRENT_FUNCTION ); break; } } @@ -944,7 +945,7 @@ void SNSolverHPC::WriteScalarOutput( unsigned idx_iter ) { } idx_field--; break; - case VAR_ABSORPTION_GREEN: _screenOutputFields[idx_field] = _varAbsorptionHohlraumGreen; break; + case VAR_ABSORPTION_GREEN: _screenOutputFields[idx_field] = _absorptionValsBlocksGreen[0]; break; default: ErrorMessages::Error( "Screen output group not defined!", CURRENT_FUNCTION ); break; } } @@ -1002,9 +1003,17 @@ void SNSolverHPC::WriteScalarOutput( unsigned idx_iter ) { idx_field--; break; case VAR_ABSORPTION_GREEN: _historyOutputFields[idx_field] = _varAbsorptionHohlraumGreen; break; - case VAR_ABSORPTION_GREEN_LINE: + case ABSORPTION_GREEN_LINE: for( unsigned i = 0; i < _settings->GetNumProbingCellsLineHohlraum(); i++ ) { - _historyOutputFieldNames[idx_field] = _varAbsorptionValsIntegrated[i]; + _historyOutputFields[idx_field] = _absorptionValsLineSegment[i]; + idx_field++; + } + idx_field--; + break; + case ABSORPTION_GREEN_BLOCK: + for( unsigned i = 0; i < 44; i++ ) { + _historyOutputFields[idx_field] = _absorptionValsBlocksGreen[i]; + // std::cout << _absorptionValsBlocksGreen[i] << "/" << _historyOutputFields[idx_field] << std::endl; idx_field++; } idx_field--; @@ -1046,7 +1055,8 @@ void SNSolverHPC::PrintScreenOutput( unsigned idx_iter ) { TOTAL_PARTICLE_ABSORPTION_HORIZONTAL, PROBE_MOMENT_TIME_TRACE, VAR_ABSORPTION_GREEN, - VAR_ABSORPTION_GREEN_LINE }; + ABSORPTION_GREEN_BLOCK, + ABSORPTION_GREEN_LINE }; std::vector booleanFields = { VTK_OUTPUT, CSV_OUTPUT }; if( !( integerFields.end() == std::find( integerFields.begin(), integerFields.end(), _settings->GetScreenOutput()[idx_field] ) ) ) { @@ -1125,7 +1135,15 @@ void SNSolverHPC::PrepareHistoryOutput() { idx_field--; break; case VAR_ABSORPTION_GREEN: _historyOutputFieldNames[idx_field] = "Var. absorption green"; break; - case VAR_ABSORPTION_GREEN_LINE: + case ABSORPTION_GREEN_BLOCK: + for( unsigned i = 0; i < 44; i++ ) { + _historyOutputFieldNames[idx_field] = "Probe Green Block " + std::to_string( i ); + idx_field++; + } + idx_field--; + break; + + case ABSORPTION_GREEN_LINE: for( unsigned i = 0; i < _settings->GetNumProbingCellsLineHohlraum(); i++ ) { _historyOutputFieldNames[idx_field] = "Probe Green Line " + std::to_string( i ); idx_field++; @@ -1155,7 +1173,7 @@ void SNSolverHPC::PrintHistoryOutput( unsigned idx_iter ) { } tmp = TextProcessingToolbox::DoubleToScientificNotation( _historyOutputFields[_settings->GetNHistoryOutput() - 1] ); lineToPrint += tmp; // Last element without comma - + // std::cout << lineToPrint << std::endl; if( _settings->GetHistoryOutputFrequency() != 0 && idx_iter % (unsigned)_settings->GetHistoryOutputFrequency() == 0 ) { log->info( lineToPrint ); } @@ -1287,7 +1305,7 @@ void SNSolverHPC::WriteVolumeOutput( unsigned idx_iter ) { void SNSolverHPC::PrintVolumeOutput( int idx_iter ) { if( _settings->GetSaveRestartSolutionFrequency() != 0 && idx_iter % (int)_settings->GetSaveRestartSolutionFrequency() == 0 ) { - std::cout << "Saving restart solution at iteration " << idx_iter << std::endl; + // std::cout << "Saving restart solution at iteration " << idx_iter << std::endl; WriteRestartSolution( _settings->GetOutputFile(), _sol, _scalarFlux, _rank, idx_iter ); } @@ -1447,6 +1465,79 @@ void SNSolverHPC::SetProbingCellsLineGreen() { _probingCellsLineGreen.insert( _probingCellsLineGreen.end(), side2.begin(), side2.end() ); _probingCellsLineGreen.insert( _probingCellsLineGreen.end(), side3.begin(), side3.end() ); _probingCellsLineGreen.insert( _probingCellsLineGreen.end(), side4.begin(), side4.end() ); + + // Block-wise approach + // Initialize the vector to store the corner points of each block + std::vector>> block_corners; + + double block_size = 0.05; + + // Loop to fill the blocks + for( int i = 0; i < 8; ++i ) { // 8 blocks in the x-direction (horizontal) + + // Top row + double x1 = -0.2 + i * block_size; + double y1 = 0.4; + double x2 = x1 + block_size; + double y2 = y1 - block_size; + + std::vector> corners = { + { x1, y1 }, // top-left + { x2, y1 }, // top-right + { x2, y2 }, // bottom-right + { x1, y2 } // bottom-left + }; + block_corners.push_back( corners ); + + // bottom row + y1 = -0.35; + y2 = y1 - block_size; + + corners = { + { x1, y1 }, // top-left + { x2, y1 }, // top-right + { x2, y2 }, // bottom-right + { x1, y2 } // bottom-left + }; + block_corners.push_back( corners ); + } + for( int j = 0; j < 14; ++j ) { // 14 blocks in the y-direction (vertical) + + // left column + double x1 = -0.2; + double y1 = 0.35 - j * block_size; + double x2 = x1 + block_size; + double y2 = y1 - block_size; + + // Store the four corner points for this block + std::vector> corners = { + { x1, y1 }, // top-left + { x2, y1 }, // top-right + { x2, y2 }, // bottom-right + { x1, y2 } // bottom-left + }; + + block_corners.push_back( corners ); + + // right column double x1 = 0.15; + x1 = 0.15; + x2 = x1 + block_size; + + // Store the four corner points for this block + corners = { + { x1, y1 }, // top-left + { x2, y1 }, // top-right + { x2, y2 }, // bottom-right + { x1, y2 } // bottom-left + }; + + block_corners.push_back( corners ); + } + + // Compute the probing cells for each block + for( int i = 0; i < _nProbingCellsBlocksGreen; i++ ) { + _probingCellsBlocksGreen.push_back( _mesh->GetCellsofRectangle( block_corners[i] ) ); + } } } @@ -1456,17 +1547,25 @@ void SNSolverHPC::ComputeQOIsGreenProbingLine() { double dl = 2 * ( horizontalLineWidth + verticalLineWidth ) / ( (double)_nProbingCellsLineGreen ); double area = dl * _thicknessGreen; - double a_g = 0; double l_max = _nProbingCellsLineGreen * dl; +#pragma omp parallel for for( unsigned i = 0; i < _nProbingCellsLineGreen; i++ ) { // Loop over probing cells - _absorptionValsIntegrated[i] = - ( _sigmaT[_probingCellsLineGreen[i]] - _sigmaS[_probingCellsLineGreen[i]] ) * _scalarFlux[_probingCellsLineGreen[i]] * area; - a_g += _absorptionValsIntegrated[i] / (double)_nProbingCellsLineGreen; + _absorptionValsLineSegment[i] = ( _sigmaT[_probingCellsLineGreen[i]] - _sigmaS[_probingCellsLineGreen[i]] ) * + _scalarFlux[_probingCellsLineGreen[i]] * _areas[_probingCellsLineGreen[i]] / dl; } - for( unsigned i = 0; i < _nProbingCellsLineGreen; i++ ) { // Loop over probing cells - _varAbsorptionValsIntegrated[i] = dl / l_max * ( a_g - _absorptionValsIntegrated[i] ) * ( a_g - _absorptionValsIntegrated[i] ); + + // Block-wise approach + // #pragma omp parallel for + for( unsigned i = 0; i < _nProbingCellsBlocksGreen; i++ ) { + _absorptionValsBlocksGreen[i] = 0.0; + for( unsigned j = 0; j < _probingCellsBlocksGreen[i].size(); j++ ) { + _absorptionValsBlocksGreen[i] += ( _sigmaT[_probingCellsBlocksGreen[i][j]] - _sigmaS[_probingCellsBlocksGreen[i][j]] ) * + _scalarFlux[_probingCellsBlocksGreen[i][j]] * _areas[_probingCellsBlocksGreen[i][j]]; + } } + // std::cout << _absorptionValsBlocksGreen[1] << std::endl; + // std::cout << _absorptionValsLineSegment[1] << std::endl; } std::vector SNSolverHPC::linspace2D( const std::vector& start, const std::vector& end, unsigned num_points ) { diff --git a/src/solvers/solverbase.cpp b/src/solvers/solverbase.cpp index fbae12f6..18e8c2fd 100644 --- a/src/solvers/solverbase.cpp +++ b/src/solvers/solverbase.cpp @@ -380,7 +380,7 @@ void SolverBase::WriteScalarOutput( unsigned idx_iter ) { idx_field--; break; case VAR_ABSORPTION_GREEN: _historyOutputFields[idx_field] = _problem->GetVarAbsorptionHohlraumGreen(); break; - case VAR_ABSORPTION_GREEN_LINE: + case ABSORPTION_GREEN_LINE: for( unsigned i = 0; i < _settings->GetNumProbingCellsLineHohlraum(); i++ ) { _historyOutputFieldNames[idx_field] = _problem->GetCurrentVarProbeValuesGreenLine()[i]; idx_field++; @@ -420,7 +420,8 @@ void SolverBase::PrintScreenOutput( unsigned idx_iter ) { TOTAL_PARTICLE_ABSORPTION_HORIZONTAL, PROBE_MOMENT_TIME_TRACE, VAR_ABSORPTION_GREEN, - VAR_ABSORPTION_GREEN_LINE }; + ABSORPTION_GREEN_BLOCK, + ABSORPTION_GREEN_LINE }; std::vector booleanFields = { VTK_OUTPUT, CSV_OUTPUT }; if( !( integerFields.end() == std::find( integerFields.begin(), integerFields.end(), _settings->GetScreenOutput()[idx_field] ) ) ) { @@ -495,7 +496,7 @@ void SolverBase::PrepareHistoryOutput() { idx_field--; break; case VAR_ABSORPTION_GREEN: _historyOutputFieldNames[idx_field] = "Var. absorption green"; break; - case VAR_ABSORPTION_GREEN_LINE: + case ABSORPTION_GREEN_LINE: for( unsigned i = 0; i < _settings->GetNumProbingCellsLineHohlraum(); i++ ) { _historyOutputFieldNames[idx_field] = "Probe Green Line " + std::to_string( i ); idx_field++; From 682c3bbca318017e3c49b0a7f490f0351e8373b7 Mon Sep 17 00:00:00 2001 From: ScSteffen Date: Tue, 22 Oct 2024 10:11:53 -0400 Subject: [PATCH 16/21] change restart file to binary --- src/common/io.cpp | 60 ++++++++++++++++++++++++----------------------- 1 file changed, 31 insertions(+), 29 deletions(-) diff --git a/src/common/io.cpp b/src/common/io.cpp index cdc8a307..f89c88cb 100644 --- a/src/common/io.cpp +++ b/src/common/io.cpp @@ -453,26 +453,27 @@ void WriteRestartSolution( std::filesystem::remove( outputFile ); } - // Open the file for writing (automatically creates a new file) - std::ofstream outFile( outputFile, std::ios::out ); + // Open the file for binary writing + std::ofstream outFile( outputFile, std::ios::out | std::ios::binary ); if( !outFile ) { - ErrorMessages::Error( "Error opening restart solution file.", CURRENT_FUNCTION ); + std::cerr << "Error opening restart solution file." << std::endl; return; } - // Write the iteration number as the first line - outFile << idx_iter << "\n"; - outFile << std::fixed << std::setprecision( 15 ); + // Write the iteration number as the first item (in binary) + outFile.write( reinterpret_cast( &idx_iter ), sizeof( idx_iter ) ); - // Write each element of the solution to the file on a new line - for( const double& value : solution ) { - outFile << value << "\n"; - } + // Write the size of the solution and scalarFlux vectors (optional but useful for reading) + size_t solution_size = solution.size(); + size_t scalarFlux_size = scalarFlux.size(); + outFile.write( reinterpret_cast( &solution_size ), sizeof( solution_size ) ); + outFile.write( reinterpret_cast( &scalarFlux_size ), sizeof( scalarFlux_size ) ); - // Append each element of the scalarFlux to the file on a new line - for( const double& fluxValue : scalarFlux ) { - outFile << fluxValue << "\n"; - } + // Write each element of the solution vector in binary + outFile.write( reinterpret_cast( solution.data() ), solution_size * sizeof( double ) ); + + // Write each element of the scalarFlux vector in binary + outFile.write( reinterpret_cast( scalarFlux.data() ), scalarFlux_size * sizeof( double ) ); // Close the file outFile.close(); @@ -483,31 +484,32 @@ int LoadRestartSolution( // Generate filename with rank number std::string inputFile = baseInputFile + "_restart_rank_" + std::to_string( rank ); - // Open the file for reading - std::ifstream inFile( inputFile ); + // Open the file for binary reading + std::ifstream inFile( inputFile, std::ios::in | std::ios::binary ); if( !inFile ) { std::cerr << "Error opening restart solution file for reading." << std::endl; return 0; // Signal failure } - // Read the iteration number from the first line - std::string line; - std::getline( inFile, line ); - int idx_iter = std::stoi( line ); // Convert the line to an integer + // Read the iteration number + int idx_iter = 0; + inFile.read( reinterpret_cast( &idx_iter ), sizeof( idx_iter ) ); - // Read data into a temporary vector - std::vector tempData; - double value; - while( std::getline( inFile, line ) ) { - std::istringstream converter( line ); - converter >> std::fixed >> std::setprecision( 15 ) >> value; - tempData.push_back( value ); - } + // Read the size of the solution and scalarFlux vectors + size_t solution_size, scalarFlux_size; + inFile.read( reinterpret_cast( &solution_size ), sizeof( solution_size ) ); + inFile.read( reinterpret_cast( &scalarFlux_size ), sizeof( scalarFlux_size ) ); + + // Temporary vector to hold the full data + std::vector tempData( solution_size + scalarFlux_size ); + + // Read the entire data block in binary form + inFile.read( reinterpret_cast( tempData.data() ), ( solution_size + scalarFlux_size ) * sizeof( double ) ); // Close the file inFile.close(); - // Verify that we have at least nCells entries to populate scalarFlux + // Verify that we have enough entries to populate scalarFlux if( tempData.size() < nCells ) { std::cerr << "Not enough data to populate scalar flux vector." << std::endl; return 0; // Signal failure From bf3915882084da38aad41badf41f483d86592e1b Mon Sep 17 00:00:00 2001 From: Steffen Date: Tue, 5 Nov 2024 13:00:34 -0500 Subject: [PATCH 17/21] fixed probe green line --- src/solvers/snsolver_hpc.cpp | 134 +++++++++++++++++++---------------- 1 file changed, 74 insertions(+), 60 deletions(-) diff --git a/src/solvers/snsolver_hpc.cpp b/src/solvers/snsolver_hpc.cpp index 6ebcff27..df701061 100644 --- a/src/solvers/snsolver_hpc.cpp +++ b/src/solvers/snsolver_hpc.cpp @@ -20,7 +20,6 @@ SNSolverHPC::SNSolverHPC( Config* settings ) { _numProcs = 1; // default values _rank = 0; #endif - _settings = settings; _currTime = 0.0; _idx_start_iter = 0; @@ -43,13 +42,15 @@ SNSolverHPC::SNSolverHPC( Config* settings ) { ErrorMessages::Error( "The number of processors must be less than or equal to the number of quadrature points.", CURRENT_FUNCTION ); } - _localNSys = _nq / _numProcs; + _localNSys = _nSys / _numProcs;//(_numProcs-1); _startSysIdx = _rank * _localNSys; _endSysIdx = _rank * _localNSys + _localNSys; if( _rank == _numProcs - 1 ) { - _localNSys += _nq % _numProcs; - _endSysIdx = _nSys; + _localNSys = _nSys - _startSysIdx; + _endSysIdx = _nSys; } + + //std::cout << "Rank: " << _rank << " startSysIdx: " << _startSysIdx << " endSysIdx: " << _endSysIdx << " localNSys: " << _localNSys << std::endl; _spatialOrder = _settings->GetSpatialOrder(); _temporalOrder = _settings->GetTemporalOrder(); @@ -205,7 +206,8 @@ SNSolverHPC::SNSolverHPC( Config* settings ) { unsigned n_probes = 4; if( _settings->GetProblemName() == PROBLEM_SymmetricHohlraum ) n_probes = 4; - if( _settings->GetProblemName() == PROBLEM_QuarterHohlraum ) n_probes = 2; + //if( _settings->GetProblemName() == PROBLEM_QuarterHohlraum ) n_probes = 2; + _probingMoments = std::vector( n_probes * 3, 0. ); // 10 time horizons if( _settings->GetProblemName() == PROBLEM_SymmetricHohlraum ) { @@ -216,12 +218,12 @@ SNSolverHPC::SNSolverHPC( Config* settings ) { _mesh->GetCellsofBall( 0., 0.5, 0.01 ), }; } - else if( _settings->GetProblemName() == PROBLEM_QuarterHohlraum ) { - _probingCellsHohlraum = { - _mesh->GetCellsofBall( 0.4, 0., 0.01 ), - _mesh->GetCellsofBall( 0., 0.5, 0.01 ), - }; - } + //else if( _settings->GetProblemName() == PROBLEM_QuarterHohlraum ) { + // _probingCellsHohlraum = { + // _mesh->GetCellsofBall( 0.4, 0., 0.01 ), + // _mesh->GetCellsofBall( 0., 0.5, 0.01 ), + // }; + //} // Green _thicknessGreen = 0.05; @@ -233,13 +235,13 @@ SNSolverHPC::SNSolverHPC( Config* settings ) { _cornerUpperRightGreen = { 0.2 + _centerGreen[0], 0.4 + _centerGreen[1] }; _cornerLowerRightGreen = { 0.2 + _centerGreen[0], -0.4 + _centerGreen[1] }; } - else { - _centerGreen = { 0.0, 0.0 }; - _cornerUpperLeftGreen = { 0., 0.4 - _thicknessGreen / 2.0 }; - _cornerLowerLeftGreen = { 0., +_thicknessGreen / 2.0 }; - _cornerUpperRightGreen = { 0.2 - _thicknessGreen / 2.0, 0.4 - _thicknessGreen / 2.0 }; - _cornerLowerRightGreen = { 0.2 - _thicknessGreen / 2.0, 0. + _thicknessGreen / 2.0 }; - } + //else { + // _centerGreen = { 0.0, 0.0 }; + // _cornerUpperLeftGreen = { 0., 0.4 - _thicknessGreen / 2.0 }; + // _cornerLowerLeftGreen = { 0., +_thicknessGreen / 2.0 }; + // _cornerUpperRightGreen = { 0.2 - _thicknessGreen / 2.0, 0.4 - _thicknessGreen / 2.0 }; + // _cornerLowerRightGreen = { 0.2 - _thicknessGreen / 2.0, 0. + _thicknessGreen / 2.0 }; + //} _nProbingCellsLineGreen = _settings->GetNumProbingCellsLineHohlraum(); @@ -577,7 +579,7 @@ void SNSolverHPC::IterPostprocessing() { } } - if( _settings->GetProblemName() == PROBLEM_SymmetricHohlraum || _settings->GetProblemName() == PROBLEM_QuarterHohlraum ) { + if( _settings->GetProblemName() == PROBLEM_SymmetricHohlraum){//} || _settings->GetProblemName() == PROBLEM_QuarterHohlraum ) { double x = _cellMidPoints[Idx2D( idx_cell, 0, _nDim )]; double y = _cellMidPoints[Idx2D( idx_cell, 1, _nDim )]; @@ -938,7 +940,7 @@ void SNSolverHPC::WriteScalarOutput( unsigned idx_iter ) { case TOTAL_PARTICLE_ABSORPTION_HORIZONTAL: _screenOutputFields[idx_field] = _totalAbsorptionHohlraumHorizontal; break; case PROBE_MOMENT_TIME_TRACE: if( _settings->GetProblemName() == PROBLEM_SymmetricHohlraum ) n_probes = 4; - if( _settings->GetProblemName() == PROBLEM_QuarterHohlraum ) n_probes = 2; + //if( _settings->GetProblemName() == PROBLEM_QuarterHohlraum ) n_probes = 2; for( unsigned i = 0; i < n_probes; i++ ) { _screenOutputFields[idx_field] = _probingMoments[Idx2D( i, 0, 3 )]; idx_field++; @@ -993,7 +995,7 @@ void SNSolverHPC::WriteScalarOutput( unsigned idx_iter ) { case TOTAL_PARTICLE_ABSORPTION_HORIZONTAL: _historyOutputFields[idx_field] = _totalAbsorptionHohlraumHorizontal; break; case PROBE_MOMENT_TIME_TRACE: if( _settings->GetProblemName() == PROBLEM_SymmetricHohlraum ) n_probes = 4; - if( _settings->GetProblemName() == PROBLEM_QuarterHohlraum ) n_probes = 2; + //if( _settings->GetProblemName() == PROBLEM_QuarterHohlraum ) n_probes = 2; for( unsigned i = 0; i < n_probes; i++ ) { for( unsigned j = 0; j < 3; j++ ) { _historyOutputFields[idx_field] = _probingMoments[Idx2D( i, j, 3 )]; @@ -1125,7 +1127,7 @@ void SNSolverHPC::PrepareHistoryOutput() { case TOTAL_PARTICLE_ABSORPTION_HORIZONTAL: _historyOutputFieldNames[idx_field] = "Cumulated_absorption_horizontal_wall"; break; case PROBE_MOMENT_TIME_TRACE: if( _settings->GetProblemName() == PROBLEM_SymmetricHohlraum ) n_probes = 4; - if( _settings->GetProblemName() == PROBLEM_QuarterHohlraum ) n_probes = 2; + //if( _settings->GetProblemName() == PROBLEM_QuarterHohlraum ) n_probes = 2; for( unsigned i = 0; i < n_probes; i++ ) { for( unsigned j = 0; j < 3; j++ ) { _historyOutputFieldNames[idx_field] = "Probe " + std::to_string( i ) + " u_" + std::to_string( j ); @@ -1417,55 +1419,67 @@ void SNSolverHPC::SetGhostCells() { void SNSolverHPC::SetProbingCellsLineGreen() { - if( _settings->GetProblemName() == PROBLEM_QuarterHohlraum ) { - double verticalLineWidth = std::abs( _cornerUpperLeftGreen[1] - _cornerLowerLeftGreen[1] ); - double horizontalLineWidth = std::abs( _cornerUpperLeftGreen[0] - _cornerUpperRightGreen[0] ); - - // double dx = 2 * ( horizontalLineWidth + verticalLineWidth ) / ( (double)_nProbingCellsLineGreen ); - - unsigned nHorizontalProbingCells = - (unsigned)std::ceil( _nProbingCellsLineGreen * ( horizontalLineWidth / ( horizontalLineWidth + verticalLineWidth ) ) ); - unsigned nVerticalProbingCells = _nProbingCellsLineGreen - nHorizontalProbingCells; - - _probingCellsLineGreen = std::vector( _nProbingCellsLineGreen ); - - // Sample points on each side of the rectangle - std::vector side3 = linspace2D( _cornerLowerRightGreen, _cornerUpperRightGreen, nVerticalProbingCells ); - std::vector side4 = linspace2D( _cornerUpperRightGreen, _cornerUpperLeftGreen, nHorizontalProbingCells ); - - // Combine the points from each side - _probingCellsLineGreen.insert( _probingCellsLineGreen.end(), side3.begin(), side3.end() ); - _probingCellsLineGreen.insert( _probingCellsLineGreen.end(), side4.begin(), side4.end() ); - } - else if( _settings->GetProblemName() == PROBLEM_SymmetricHohlraum ) { - double verticalLineWidth = std::abs( _cornerUpperLeftGreen[1] - _cornerLowerLeftGreen[1] ); - double horizontalLineWidth = std::abs( _cornerUpperLeftGreen[0] - _cornerUpperRightGreen[0] ); - - // double dx = 2 * ( horizontalLineWidth + verticalLineWidth ) / ( (double)_nProbingCellsLineGreen ); - - unsigned nHorizontalProbingCells = - (unsigned)std::ceil( _nProbingCellsLineGreen / 2 * ( horizontalLineWidth / ( horizontalLineWidth + verticalLineWidth ) ) ); - unsigned nVerticalProbingCells = _nProbingCellsLineGreen - nHorizontalProbingCells; - - _probingCellsLineGreen = std::vector( _nProbingCellsLineGreen ); - + //if( _settings->GetProblemName() == PROBLEM_QuarterHohlraum ) { + // double verticalLineWidth = std::abs( _cornerUpperLeftGreen[1] - _cornerLowerLeftGreen[1] ); + // double horizontalLineWidth = std::abs( _cornerUpperLeftGreen[0] - _cornerUpperRightGreen[0] ); + // + // // double dx = 2 * ( horizontalLineWidth + verticalLineWidth ) / ( (double)_nProbingCellsLineGreen ); + // + // unsigned nHorizontalProbingCells = + // (unsigned)std::ceil( _nProbingCellsLineGreen * ( horizontalLineWidth / ( horizontalLineWidth + verticalLineWidth ) ) ); + // unsigned nVerticalProbingCells = _nProbingCellsLineGreen - nHorizontalProbingCells; + // + // _probingCellsLineGreen = std::vector( _nProbingCellsLineGreen ); + // + // // Sample points on each side of the rectangle + // std::vector side3 = linspace2D( _cornerLowerRightGreen, _cornerUpperRightGreen, nVerticalProbingCells ); + // std::vector side4 = linspace2D( _cornerUpperRightGreen, _cornerUpperLeftGreen, nHorizontalProbingCells ); + // + // // Combine the points from each side + // _probingCellsLineGreen.insert( _probingCellsLineGreen.end(), side3.begin(), side3.end() ); + // _probingCellsLineGreen.insert( _probingCellsLineGreen.end(), side4.begin(), side4.end() ); + //} + //else + if( _settings->GetProblemName() == PROBLEM_SymmetricHohlraum ) { + + std::vector p1 = { _cornerUpperLeftGreen[0] + _thicknessGreen / 2.0, _cornerUpperLeftGreen[1] - _thicknessGreen / 2.0 }; std::vector p2 = { _cornerLowerLeftGreen[0] + _thicknessGreen / 2.0, _cornerLowerLeftGreen[1] + _thicknessGreen / 2.0 }; std::vector p3 = { _cornerUpperRightGreen[0] - _thicknessGreen / 2.0, _cornerUpperRightGreen[1] - _thicknessGreen / 2.0 }; std::vector p4 = { _cornerLowerRightGreen[0] - _thicknessGreen / 2.0, _cornerLowerRightGreen[1] + _thicknessGreen / 2.0 }; + double verticalLineWidth = std::abs( p1[1] - p2[1]); + double horizontalLineWidth = std::abs( p1[0] - p3[0] ); + + double pt_ratio_h = horizontalLineWidth / ( horizontalLineWidth + verticalLineWidth ); + double pt_ratio_v = verticalLineWidth / ( horizontalLineWidth + verticalLineWidth ); + + unsigned nHorizontalProbingCells = (unsigned)std::ceil( _nProbingCellsLineGreen / 2 * pt_ratio_h ); + unsigned nVerticalProbingCells = _nProbingCellsLineGreen / 2 - nHorizontalProbingCells; + + _probingCellsLineGreen = std::vector( 2*(nVerticalProbingCells + nHorizontalProbingCells) ); + // Sample points on each side of the rectangle std::vector side1 = linspace2D( p1, p2, nVerticalProbingCells ); std::vector side2 = linspace2D( p2, p3, nHorizontalProbingCells ); std::vector side3 = linspace2D( p3, p4, nVerticalProbingCells ); std::vector side4 = linspace2D( p4, p1, nHorizontalProbingCells ); - // Combine the points from each side - _probingCellsLineGreen.insert( _probingCellsLineGreen.end(), side1.begin(), side1.end() ); - _probingCellsLineGreen.insert( _probingCellsLineGreen.end(), side2.begin(), side2.end() ); - _probingCellsLineGreen.insert( _probingCellsLineGreen.end(), side3.begin(), side3.end() ); - _probingCellsLineGreen.insert( _probingCellsLineGreen.end(), side4.begin(), side4.end() ); + for ( unsigned i = 0; i < nVerticalProbingCells; ++i ) { + _probingCellsLineGreen[i] = side1[i]; + } + for ( unsigned i = 0; i < nHorizontalProbingCells; ++i ) { + _probingCellsLineGreen[i+nVerticalProbingCells] = side2[i]; + } + for ( unsigned i = 0; i < nVerticalProbingCells; ++i ) { + _probingCellsLineGreen[i+nVerticalProbingCells+nHorizontalProbingCells] = side3[i]; + } + for ( unsigned i = 0; i < nHorizontalProbingCells; ++i ) { + _probingCellsLineGreen[i+2*nVerticalProbingCells+nHorizontalProbingCells] = side4[i]; + } + + // Block-wise approach // Initialize the vector to store the corner points of each block std::vector>> block_corners; @@ -1552,7 +1566,7 @@ void SNSolverHPC::ComputeQOIsGreenProbingLine() { #pragma omp parallel for for( unsigned i = 0; i < _nProbingCellsLineGreen; i++ ) { // Loop over probing cells _absorptionValsLineSegment[i] = ( _sigmaT[_probingCellsLineGreen[i]] - _sigmaS[_probingCellsLineGreen[i]] ) * - _scalarFlux[_probingCellsLineGreen[i]] * _areas[_probingCellsLineGreen[i]] / dl; + _scalarFlux[_probingCellsLineGreen[i]]; } // Block-wise approach From 8af674b2c1b34a8f52640fd0845bd29557768f14 Mon Sep 17 00:00:00 2001 From: Steffen Date: Tue, 5 Nov 2024 14:02:51 -0500 Subject: [PATCH 18/21] update mpi parallal partition --- src/solvers/snsolver_hpc.cpp | 126 ++++++++++++++++++----------------- 1 file changed, 66 insertions(+), 60 deletions(-) diff --git a/src/solvers/snsolver_hpc.cpp b/src/solvers/snsolver_hpc.cpp index df701061..97fd04dc 100644 --- a/src/solvers/snsolver_hpc.cpp +++ b/src/solvers/snsolver_hpc.cpp @@ -42,15 +42,24 @@ SNSolverHPC::SNSolverHPC( Config* settings ) { ErrorMessages::Error( "The number of processors must be less than or equal to the number of quadrature points.", CURRENT_FUNCTION ); } - _localNSys = _nSys / _numProcs;//(_numProcs-1); - _startSysIdx = _rank * _localNSys; - _endSysIdx = _rank * _localNSys + _localNSys; - if( _rank == _numProcs - 1 ) { - _localNSys = _nSys - _startSysIdx; - _endSysIdx = _nSys; + if( _numProcs == 1 ) { + _localNSys = _nSys; + _startSysIdx = 0; + _endSysIdx = _nSys; } - - //std::cout << "Rank: " << _rank << " startSysIdx: " << _startSysIdx << " endSysIdx: " << _endSysIdx << " localNSys: " << _localNSys << std::endl; + else { + _localNSys = _nSys / ( _numProcs - 1 ); + _startSysIdx = _rank * _localNSys; + _endSysIdx = _rank * _localNSys + _localNSys; + + if( _rank == _numProcs - 1 ) { + _localNSys = _nSys - _startSysIdx; + _endSysIdx = _nSys; + } + } + + // std::cout << "Rank: " << _rank << " startSysIdx: " << _startSysIdx << " endSysIdx: " << _endSysIdx << " localNSys: " << _localNSys << + // std::endl; _spatialOrder = _settings->GetSpatialOrder(); _temporalOrder = _settings->GetTemporalOrder(); @@ -206,8 +215,8 @@ SNSolverHPC::SNSolverHPC( Config* settings ) { unsigned n_probes = 4; if( _settings->GetProblemName() == PROBLEM_SymmetricHohlraum ) n_probes = 4; - //if( _settings->GetProblemName() == PROBLEM_QuarterHohlraum ) n_probes = 2; - + // if( _settings->GetProblemName() == PROBLEM_QuarterHohlraum ) n_probes = 2; + _probingMoments = std::vector( n_probes * 3, 0. ); // 10 time horizons if( _settings->GetProblemName() == PROBLEM_SymmetricHohlraum ) { @@ -218,12 +227,12 @@ SNSolverHPC::SNSolverHPC( Config* settings ) { _mesh->GetCellsofBall( 0., 0.5, 0.01 ), }; } - //else if( _settings->GetProblemName() == PROBLEM_QuarterHohlraum ) { - // _probingCellsHohlraum = { - // _mesh->GetCellsofBall( 0.4, 0., 0.01 ), - // _mesh->GetCellsofBall( 0., 0.5, 0.01 ), - // }; - //} + // else if( _settings->GetProblemName() == PROBLEM_QuarterHohlraum ) { + // _probingCellsHohlraum = { + // _mesh->GetCellsofBall( 0.4, 0., 0.01 ), + // _mesh->GetCellsofBall( 0., 0.5, 0.01 ), + // }; + // } // Green _thicknessGreen = 0.05; @@ -235,13 +244,13 @@ SNSolverHPC::SNSolverHPC( Config* settings ) { _cornerUpperRightGreen = { 0.2 + _centerGreen[0], 0.4 + _centerGreen[1] }; _cornerLowerRightGreen = { 0.2 + _centerGreen[0], -0.4 + _centerGreen[1] }; } - //else { - // _centerGreen = { 0.0, 0.0 }; - // _cornerUpperLeftGreen = { 0., 0.4 - _thicknessGreen / 2.0 }; - // _cornerLowerLeftGreen = { 0., +_thicknessGreen / 2.0 }; - // _cornerUpperRightGreen = { 0.2 - _thicknessGreen / 2.0, 0.4 - _thicknessGreen / 2.0 }; - // _cornerLowerRightGreen = { 0.2 - _thicknessGreen / 2.0, 0. + _thicknessGreen / 2.0 }; - //} + // else { + // _centerGreen = { 0.0, 0.0 }; + // _cornerUpperLeftGreen = { 0., 0.4 - _thicknessGreen / 2.0 }; + // _cornerLowerLeftGreen = { 0., +_thicknessGreen / 2.0 }; + // _cornerUpperRightGreen = { 0.2 - _thicknessGreen / 2.0, 0.4 - _thicknessGreen / 2.0 }; + // _cornerLowerRightGreen = { 0.2 - _thicknessGreen / 2.0, 0. + _thicknessGreen / 2.0 }; + // } _nProbingCellsLineGreen = _settings->GetNumProbingCellsLineHohlraum(); @@ -579,7 +588,7 @@ void SNSolverHPC::IterPostprocessing() { } } - if( _settings->GetProblemName() == PROBLEM_SymmetricHohlraum){//} || _settings->GetProblemName() == PROBLEM_QuarterHohlraum ) { + if( _settings->GetProblemName() == PROBLEM_SymmetricHohlraum ) { //} || _settings->GetProblemName() == PROBLEM_QuarterHohlraum ) { double x = _cellMidPoints[Idx2D( idx_cell, 0, _nDim )]; double y = _cellMidPoints[Idx2D( idx_cell, 1, _nDim )]; @@ -940,7 +949,7 @@ void SNSolverHPC::WriteScalarOutput( unsigned idx_iter ) { case TOTAL_PARTICLE_ABSORPTION_HORIZONTAL: _screenOutputFields[idx_field] = _totalAbsorptionHohlraumHorizontal; break; case PROBE_MOMENT_TIME_TRACE: if( _settings->GetProblemName() == PROBLEM_SymmetricHohlraum ) n_probes = 4; - //if( _settings->GetProblemName() == PROBLEM_QuarterHohlraum ) n_probes = 2; + // if( _settings->GetProblemName() == PROBLEM_QuarterHohlraum ) n_probes = 2; for( unsigned i = 0; i < n_probes; i++ ) { _screenOutputFields[idx_field] = _probingMoments[Idx2D( i, 0, 3 )]; idx_field++; @@ -995,7 +1004,7 @@ void SNSolverHPC::WriteScalarOutput( unsigned idx_iter ) { case TOTAL_PARTICLE_ABSORPTION_HORIZONTAL: _historyOutputFields[idx_field] = _totalAbsorptionHohlraumHorizontal; break; case PROBE_MOMENT_TIME_TRACE: if( _settings->GetProblemName() == PROBLEM_SymmetricHohlraum ) n_probes = 4; - //if( _settings->GetProblemName() == PROBLEM_QuarterHohlraum ) n_probes = 2; + // if( _settings->GetProblemName() == PROBLEM_QuarterHohlraum ) n_probes = 2; for( unsigned i = 0; i < n_probes; i++ ) { for( unsigned j = 0; j < 3; j++ ) { _historyOutputFields[idx_field] = _probingMoments[Idx2D( i, j, 3 )]; @@ -1127,7 +1136,7 @@ void SNSolverHPC::PrepareHistoryOutput() { case TOTAL_PARTICLE_ABSORPTION_HORIZONTAL: _historyOutputFieldNames[idx_field] = "Cumulated_absorption_horizontal_wall"; break; case PROBE_MOMENT_TIME_TRACE: if( _settings->GetProblemName() == PROBLEM_SymmetricHohlraum ) n_probes = 4; - //if( _settings->GetProblemName() == PROBLEM_QuarterHohlraum ) n_probes = 2; + // if( _settings->GetProblemName() == PROBLEM_QuarterHohlraum ) n_probes = 2; for( unsigned i = 0; i < n_probes; i++ ) { for( unsigned j = 0; j < 3; j++ ) { _historyOutputFieldNames[idx_field] = "Probe " + std::to_string( i ) + " u_" + std::to_string( j ); @@ -1419,45 +1428,44 @@ void SNSolverHPC::SetGhostCells() { void SNSolverHPC::SetProbingCellsLineGreen() { - //if( _settings->GetProblemName() == PROBLEM_QuarterHohlraum ) { - // double verticalLineWidth = std::abs( _cornerUpperLeftGreen[1] - _cornerLowerLeftGreen[1] ); - // double horizontalLineWidth = std::abs( _cornerUpperLeftGreen[0] - _cornerUpperRightGreen[0] ); + // if( _settings->GetProblemName() == PROBLEM_QuarterHohlraum ) { + // double verticalLineWidth = std::abs( _cornerUpperLeftGreen[1] - _cornerLowerLeftGreen[1] ); + // double horizontalLineWidth = std::abs( _cornerUpperLeftGreen[0] - _cornerUpperRightGreen[0] ); // - // // double dx = 2 * ( horizontalLineWidth + verticalLineWidth ) / ( (double)_nProbingCellsLineGreen ); + // // double dx = 2 * ( horizontalLineWidth + verticalLineWidth ) / ( (double)_nProbingCellsLineGreen ); // - // unsigned nHorizontalProbingCells = - // (unsigned)std::ceil( _nProbingCellsLineGreen * ( horizontalLineWidth / ( horizontalLineWidth + verticalLineWidth ) ) ); - // unsigned nVerticalProbingCells = _nProbingCellsLineGreen - nHorizontalProbingCells; + // unsigned nHorizontalProbingCells = + // (unsigned)std::ceil( _nProbingCellsLineGreen * ( horizontalLineWidth / ( horizontalLineWidth + verticalLineWidth ) ) ); + // unsigned nVerticalProbingCells = _nProbingCellsLineGreen - nHorizontalProbingCells; // - // _probingCellsLineGreen = std::vector( _nProbingCellsLineGreen ); + // _probingCellsLineGreen = std::vector( _nProbingCellsLineGreen ); // - // // Sample points on each side of the rectangle - // std::vector side3 = linspace2D( _cornerLowerRightGreen, _cornerUpperRightGreen, nVerticalProbingCells ); - // std::vector side4 = linspace2D( _cornerUpperRightGreen, _cornerUpperLeftGreen, nHorizontalProbingCells ); + // // Sample points on each side of the rectangle + // std::vector side3 = linspace2D( _cornerLowerRightGreen, _cornerUpperRightGreen, nVerticalProbingCells ); + // std::vector side4 = linspace2D( _cornerUpperRightGreen, _cornerUpperLeftGreen, nHorizontalProbingCells ); // - // // Combine the points from each side - // _probingCellsLineGreen.insert( _probingCellsLineGreen.end(), side3.begin(), side3.end() ); - // _probingCellsLineGreen.insert( _probingCellsLineGreen.end(), side4.begin(), side4.end() ); - //} - //else + // // Combine the points from each side + // _probingCellsLineGreen.insert( _probingCellsLineGreen.end(), side3.begin(), side3.end() ); + // _probingCellsLineGreen.insert( _probingCellsLineGreen.end(), side4.begin(), side4.end() ); + // } + // else if( _settings->GetProblemName() == PROBLEM_SymmetricHohlraum ) { - - + std::vector p1 = { _cornerUpperLeftGreen[0] + _thicknessGreen / 2.0, _cornerUpperLeftGreen[1] - _thicknessGreen / 2.0 }; std::vector p2 = { _cornerLowerLeftGreen[0] + _thicknessGreen / 2.0, _cornerLowerLeftGreen[1] + _thicknessGreen / 2.0 }; std::vector p3 = { _cornerUpperRightGreen[0] - _thicknessGreen / 2.0, _cornerUpperRightGreen[1] - _thicknessGreen / 2.0 }; std::vector p4 = { _cornerLowerRightGreen[0] - _thicknessGreen / 2.0, _cornerLowerRightGreen[1] + _thicknessGreen / 2.0 }; - double verticalLineWidth = std::abs( p1[1] - p2[1]); + double verticalLineWidth = std::abs( p1[1] - p2[1] ); double horizontalLineWidth = std::abs( p1[0] - p3[0] ); double pt_ratio_h = horizontalLineWidth / ( horizontalLineWidth + verticalLineWidth ); double pt_ratio_v = verticalLineWidth / ( horizontalLineWidth + verticalLineWidth ); unsigned nHorizontalProbingCells = (unsigned)std::ceil( _nProbingCellsLineGreen / 2 * pt_ratio_h ); - unsigned nVerticalProbingCells = _nProbingCellsLineGreen / 2 - nHorizontalProbingCells; + unsigned nVerticalProbingCells = _nProbingCellsLineGreen / 2 - nHorizontalProbingCells; - _probingCellsLineGreen = std::vector( 2*(nVerticalProbingCells + nHorizontalProbingCells) ); + _probingCellsLineGreen = std::vector( 2 * ( nVerticalProbingCells + nHorizontalProbingCells ) ); // Sample points on each side of the rectangle std::vector side1 = linspace2D( p1, p2, nVerticalProbingCells ); @@ -1465,21 +1473,19 @@ void SNSolverHPC::SetProbingCellsLineGreen() { std::vector side3 = linspace2D( p3, p4, nVerticalProbingCells ); std::vector side4 = linspace2D( p4, p1, nHorizontalProbingCells ); - for ( unsigned i = 0; i < nVerticalProbingCells; ++i ) { - _probingCellsLineGreen[i] = side1[i]; + for( unsigned i = 0; i < nVerticalProbingCells; ++i ) { + _probingCellsLineGreen[i] = side1[i]; } - for ( unsigned i = 0; i < nHorizontalProbingCells; ++i ) { - _probingCellsLineGreen[i+nVerticalProbingCells] = side2[i]; + for( unsigned i = 0; i < nHorizontalProbingCells; ++i ) { + _probingCellsLineGreen[i + nVerticalProbingCells] = side2[i]; } - for ( unsigned i = 0; i < nVerticalProbingCells; ++i ) { - _probingCellsLineGreen[i+nVerticalProbingCells+nHorizontalProbingCells] = side3[i]; + for( unsigned i = 0; i < nVerticalProbingCells; ++i ) { + _probingCellsLineGreen[i + nVerticalProbingCells + nHorizontalProbingCells] = side3[i]; } - for ( unsigned i = 0; i < nHorizontalProbingCells; ++i ) { - _probingCellsLineGreen[i+2*nVerticalProbingCells+nHorizontalProbingCells] = side4[i]; + for( unsigned i = 0; i < nHorizontalProbingCells; ++i ) { + _probingCellsLineGreen[i + 2 * nVerticalProbingCells + nHorizontalProbingCells] = side4[i]; } - - // Block-wise approach // Initialize the vector to store the corner points of each block std::vector>> block_corners; @@ -1565,8 +1571,8 @@ void SNSolverHPC::ComputeQOIsGreenProbingLine() { #pragma omp parallel for for( unsigned i = 0; i < _nProbingCellsLineGreen; i++ ) { // Loop over probing cells - _absorptionValsLineSegment[i] = ( _sigmaT[_probingCellsLineGreen[i]] - _sigmaS[_probingCellsLineGreen[i]] ) * - _scalarFlux[_probingCellsLineGreen[i]]; + _absorptionValsLineSegment[i] = + ( _sigmaT[_probingCellsLineGreen[i]] - _sigmaS[_probingCellsLineGreen[i]] ) * _scalarFlux[_probingCellsLineGreen[i]]; } // Block-wise approach From 59e24a10d73e57a7872aa128539c1d31a5087485 Mon Sep 17 00:00:00 2001 From: Steffen Date: Wed, 6 Nov 2024 15:22:07 -0500 Subject: [PATCH 19/21] update block and line wois --- include/common/io.hpp | 22 ++++-- src/common/io.cpp | 30 ++++++-- src/solvers/snsolver_hpc.cpp | 131 +++++++++++++++++++++-------------- 3 files changed, 124 insertions(+), 59 deletions(-) diff --git a/include/common/io.hpp b/include/common/io.hpp index 0c302183..a42edead 100644 --- a/include/common/io.hpp +++ b/include/common/io.hpp @@ -43,11 +43,25 @@ std::string ParseArguments( int argc, char* argv[] ); void PrintLogHeader( std::string inputFile ); -void WriteRestartSolution( - const std::string& baseOutputFile, const std::vector& solution, const std::vector& scalarFlux, int rank, int idx_iter ); +void WriteRestartSolution( const std::string& baseOutputFile, + const std::vector& solution, + const std::vector& scalarFlux, + int rank, + int idx_iter, + double totalAbsorptionHohlraumCenter, + double totalAbsorptionHohlraumVertical, + double totalAbsorptionHohlraumHorizontal, + double totalAbsorption ); -int LoadRestartSolution( - const std::string& baseInputFile, std::vector& solution, std::vector& scalarFlux, int rank, unsigned long nCells ); +int LoadRestartSolution( const std::string& baseInputFile, + std::vector& solution, + std::vector& scalarFlux, + int rank, + unsigned long nCells, + double& totalAbsorptionHohlraumCenter, + double& totalAbsorptionHohlraumVertical, + double& totalAbsorptionHohlraumHorizontal, + double& totalAbsorption ); // Matrix createSU2MeshFromImage( std::string imageName, std::string SU2Filename ); Deprecated diff --git a/src/common/io.cpp b/src/common/io.cpp index f89c88cb..aa10c89b 100644 --- a/src/common/io.cpp +++ b/src/common/io.cpp @@ -443,8 +443,15 @@ void WriteConnecitivityToFile( const std::string outputFile, } } -void WriteRestartSolution( - const std::string& baseOutputFile, const std::vector& solution, const std::vector& scalarFlux, int rank, int idx_iter ) { +void WriteRestartSolution( const std::string& baseOutputFile, + const std::vector& solution, + const std::vector& scalarFlux, + int rank, + int idx_iter, + double totalAbsorptionHohlraumCenter, + double totalAbsorptionHohlraumVertical, + double totalAbsorptionHohlraumHorizontal, + double totalAbsorption ) { // Generate filename with rank number std::string outputFile = baseOutputFile + "_restart_rank_" + std::to_string( rank ); @@ -462,6 +469,10 @@ void WriteRestartSolution( // Write the iteration number as the first item (in binary) outFile.write( reinterpret_cast( &idx_iter ), sizeof( idx_iter ) ); + outFile.write( reinterpret_cast( &totalAbsorptionHohlraumCenter ), sizeof( totalAbsorptionHohlraumCenter ) ); + outFile.write( reinterpret_cast( &totalAbsorptionHohlraumVertical ), sizeof( totalAbsorptionHohlraumVertical ) ); + outFile.write( reinterpret_cast( &totalAbsorptionHohlraumHorizontal ), sizeof( totalAbsorptionHohlraumHorizontal ) ); + outFile.write( reinterpret_cast( &totalAbsorption ), sizeof( totalAbsorption ) ); // Write the size of the solution and scalarFlux vectors (optional but useful for reading) size_t solution_size = solution.size(); @@ -479,8 +490,15 @@ void WriteRestartSolution( outFile.close(); } -int LoadRestartSolution( - const std::string& baseInputFile, std::vector& solution, std::vector& scalarFlux, int rank, unsigned long nCells ) { +int LoadRestartSolution( const std::string& baseInputFile, + std::vector& solution, + std::vector& scalarFlux, + int rank, + unsigned long nCells, + double& totalAbsorptionHohlraumCenter, + double& totalAbsorptionHohlraumVertical, + double& totalAbsorptionHohlraumHorizontal, + double& totalAbsorption ) { // Generate filename with rank number std::string inputFile = baseInputFile + "_restart_rank_" + std::to_string( rank ); @@ -494,6 +512,10 @@ int LoadRestartSolution( // Read the iteration number int idx_iter = 0; inFile.read( reinterpret_cast( &idx_iter ), sizeof( idx_iter ) ); + inFile.read( reinterpret_cast( &totalAbsorptionHohlraumCenter ), sizeof( totalAbsorptionHohlraumCenter ) ); + inFile.read( reinterpret_cast( &totalAbsorptionHohlraumVertical ), sizeof( totalAbsorptionHohlraumVertical ) ); + inFile.read( reinterpret_cast( &totalAbsorptionHohlraumHorizontal ), sizeof( totalAbsorptionHohlraumHorizontal ) ); + inFile.read( reinterpret_cast( &totalAbsorption ), sizeof( totalAbsorption ) ); // Read the size of the solution and scalarFlux vectors size_t solution_size, scalarFlux_size; diff --git a/src/solvers/snsolver_hpc.cpp b/src/solvers/snsolver_hpc.cpp index 97fd04dc..8d1c9fef 100644 --- a/src/solvers/snsolver_hpc.cpp +++ b/src/solvers/snsolver_hpc.cpp @@ -169,26 +169,6 @@ SNSolverHPC::SNSolverHPC( Config* settings ) { // _mass += _scalarFlux[idx_cell] * _areas[idx_cell]; } - if( _settings->GetLoadRestartSolution() ) - _idx_start_iter = LoadRestartSolution( _settings->GetOutputFile(), _sol, _scalarFlux, _rank, _nCells ) + 1; - -#ifdef IMPORT_MPI - MPI_Barrier( MPI_COMM_WORLD ); -#endif - SetGhostCells(); - if( _rank == 0 ) { - PrepareScreenOutput(); // Screen Output - PrepareHistoryOutput(); // History Output - } -#ifdef IMPORT_MPI - MPI_Barrier( MPI_COMM_WORLD ); -#endif - delete quad; - - // Initialiye QOIS - _mass = 0; - _rmsFlux = 0; - // Lattice { _curAbsorptionLattice = 0; @@ -212,6 +192,38 @@ SNSolverHPC::SNSolverHPC( Config* settings ) { _curAbsorptionHohlraumVertical = 0; _curAbsorptionHohlraumHorizontal = 0; _varAbsorptionHohlraumGreen = 0; + } + + if( _settings->GetLoadRestartSolution() ) + _idx_start_iter = LoadRestartSolution( _settings->GetOutputFile(), + _sol, + _scalarFlux, + _rank, + _nCells, + _totalAbsorptionHohlraumCenter, + _totalAbsorptionHohlraumVertical, + _totalAbsorptionHohlraumHorizontal, + _totalAbsorptionLattice ) + + 1; + +#ifdef IMPORT_MPI + MPI_Barrier( MPI_COMM_WORLD ); +#endif + SetGhostCells(); + if( _rank == 0 ) { + PrepareScreenOutput(); // Screen Output + PrepareHistoryOutput(); // History Output + } +#ifdef IMPORT_MPI + MPI_Barrier( MPI_COMM_WORLD ); +#endif + delete quad; + + // Initialiye QOIS + _mass = 0; + _rmsFlux = 0; + + { // Hohlraum unsigned n_probes = 4; if( _settings->GetProblemName() == PROBLEM_SymmetricHohlraum ) n_probes = 4; @@ -1317,7 +1329,15 @@ void SNSolverHPC::WriteVolumeOutput( unsigned idx_iter ) { void SNSolverHPC::PrintVolumeOutput( int idx_iter ) { if( _settings->GetSaveRestartSolutionFrequency() != 0 && idx_iter % (int)_settings->GetSaveRestartSolutionFrequency() == 0 ) { // std::cout << "Saving restart solution at iteration " << idx_iter << std::endl; - WriteRestartSolution( _settings->GetOutputFile(), _sol, _scalarFlux, _rank, idx_iter ); + WriteRestartSolution( _settings->GetOutputFile(), + _sol, + _scalarFlux, + _rank, + idx_iter, + _totalAbsorptionHohlraumCenter, + _totalAbsorptionHohlraumVertical, + _totalAbsorptionHohlraumHorizontal, + _totalAbsorptionLattice ); } if( _settings->GetVolumeOutputFrequency() != 0 && idx_iter % (int)_settings->GetVolumeOutputFrequency() == 0 ) { @@ -1452,12 +1472,12 @@ void SNSolverHPC::SetProbingCellsLineGreen() { if( _settings->GetProblemName() == PROBLEM_SymmetricHohlraum ) { std::vector p1 = { _cornerUpperLeftGreen[0] + _thicknessGreen / 2.0, _cornerUpperLeftGreen[1] - _thicknessGreen / 2.0 }; - std::vector p2 = { _cornerLowerLeftGreen[0] + _thicknessGreen / 2.0, _cornerLowerLeftGreen[1] + _thicknessGreen / 2.0 }; - std::vector p3 = { _cornerUpperRightGreen[0] - _thicknessGreen / 2.0, _cornerUpperRightGreen[1] - _thicknessGreen / 2.0 }; - std::vector p4 = { _cornerLowerRightGreen[0] - _thicknessGreen / 2.0, _cornerLowerRightGreen[1] + _thicknessGreen / 2.0 }; + std::vector p2 = { _cornerUpperRightGreen[0] - _thicknessGreen / 2.0, _cornerUpperRightGreen[1] - _thicknessGreen / 2.0 }; + std::vector p3 = { _cornerLowerRightGreen[0] - _thicknessGreen / 2.0, _cornerLowerRightGreen[1] + _thicknessGreen / 2.0 }; + std::vector p4 = { _cornerLowerLeftGreen[0] + _thicknessGreen / 2.0, _cornerLowerLeftGreen[1] + _thicknessGreen / 2.0 }; - double verticalLineWidth = std::abs( p1[1] - p2[1] ); - double horizontalLineWidth = std::abs( p1[0] - p3[0] ); + double verticalLineWidth = std::abs( p1[1] - p4[1] ); + double horizontalLineWidth = std::abs( p1[0] - p2[0] ); double pt_ratio_h = horizontalLineWidth / ( horizontalLineWidth + verticalLineWidth ); double pt_ratio_v = verticalLineWidth / ( horizontalLineWidth + verticalLineWidth ); @@ -1468,22 +1488,22 @@ void SNSolverHPC::SetProbingCellsLineGreen() { _probingCellsLineGreen = std::vector( 2 * ( nVerticalProbingCells + nHorizontalProbingCells ) ); // Sample points on each side of the rectangle - std::vector side1 = linspace2D( p1, p2, nVerticalProbingCells ); - std::vector side2 = linspace2D( p2, p3, nHorizontalProbingCells ); - std::vector side3 = linspace2D( p3, p4, nVerticalProbingCells ); - std::vector side4 = linspace2D( p4, p1, nHorizontalProbingCells ); + std::vector side1 = linspace2D( p1, p2, nHorizontalProbingCells ); // upper horizontal + std::vector side2 = linspace2D( p2, p3, nVerticalProbingCells ); // right vertical + std::vector side3 = linspace2D( p3, p4, nHorizontalProbingCells ); // lower horizontal + std::vector side4 = linspace2D( p4, p1, nVerticalProbingCells ); // left vertical - for( unsigned i = 0; i < nVerticalProbingCells; ++i ) { - _probingCellsLineGreen[i] = side1[i]; - } for( unsigned i = 0; i < nHorizontalProbingCells; ++i ) { - _probingCellsLineGreen[i + nVerticalProbingCells] = side2[i]; + _probingCellsLineGreen[i] = side1[i]; } for( unsigned i = 0; i < nVerticalProbingCells; ++i ) { - _probingCellsLineGreen[i + nVerticalProbingCells + nHorizontalProbingCells] = side3[i]; + _probingCellsLineGreen[i + nHorizontalProbingCells] = side2[i]; } for( unsigned i = 0; i < nHorizontalProbingCells; ++i ) { - _probingCellsLineGreen[i + 2 * nVerticalProbingCells + nHorizontalProbingCells] = side4[i]; + _probingCellsLineGreen[i + nVerticalProbingCells + nHorizontalProbingCells] = side3[i]; + } + for( unsigned i = 0; i < nVerticalProbingCells; ++i ) { + _probingCellsLineGreen[i + nVerticalProbingCells + 2 * nHorizontalProbingCells] = side4[i]; } // Block-wise approach @@ -1493,7 +1513,7 @@ void SNSolverHPC::SetProbingCellsLineGreen() { double block_size = 0.05; // Loop to fill the blocks - for( int i = 0; i < 8; ++i ) { // 8 blocks in the x-direction (horizontal) + for( int i = 0; i < 8; ++i ) { // 8 blocks in the x-direction (horizontal) (upper) (left to right) // Top row double x1 = -0.2 + i * block_size; @@ -1508,43 +1528,52 @@ void SNSolverHPC::SetProbingCellsLineGreen() { { x1, y2 } // bottom-left }; block_corners.push_back( corners ); + } - // bottom row - y1 = -0.35; - y2 = y1 - block_size; + for( int j = 0; j < 14; ++j ) { // 14 blocks in the y-direction (vertical) + // right column double x1 = 0.15; + double x1 = 0.15; + double y1 = 0.35 - j * block_size; + double x2 = x1 + block_size; + double y2 = y1 - block_size; - corners = { + // Store the four corner points for this block + std::vector> corners = { { x1, y1 }, // top-left { x2, y1 }, // top-right { x2, y2 }, // bottom-right { x1, y2 } // bottom-left }; + block_corners.push_back( corners ); } - for( int j = 0; j < 14; ++j ) { // 14 blocks in the y-direction (vertical) - // left column - double x1 = -0.2; - double y1 = 0.35 - j * block_size; + for( int i = 0; i < 8; ++i ) { // 8 blocks in the x-direction (horizontal) (lower) (right to left) + // bottom row + double x1 = 0.15 - i * block_size; + double y1 = -0.35; double x2 = x1 + block_size; double y2 = y1 - block_size; - // Store the four corner points for this block std::vector> corners = { { x1, y1 }, // top-left { x2, y1 }, // top-right { x2, y2 }, // bottom-right { x1, y2 } // bottom-left }; - block_corners.push_back( corners ); + } - // right column double x1 = 0.15; - x1 = 0.15; - x2 = x1 + block_size; + for( int j = 0; j < 14; ++j ) { // 14 blocks in the y-direction (vertical) (down to up) + + // left column + double x1 = -0.2; + double y1 = -0.35 + j * block_size; + double x2 = x1 + block_size; + double y2 = y1 - block_size; // Store the four corner points for this block - corners = { + std::vector> corners = { { x1, y1 }, // top-left { x2, y1 }, // top-right { x2, y2 }, // bottom-right From d462947576b966a744ae00d60ba791623ba62930 Mon Sep 17 00:00:00 2001 From: Steffen Date: Thu, 7 Nov 2024 19:43:30 -0500 Subject: [PATCH 20/21] added postprocessing fix --- src/solvers/snsolver_hpc.cpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/solvers/snsolver_hpc.cpp b/src/solvers/snsolver_hpc.cpp index 8d1c9fef..8f93f516 100644 --- a/src/solvers/snsolver_hpc.cpp +++ b/src/solvers/snsolver_hpc.cpp @@ -797,7 +797,9 @@ void SNSolverHPC::IterPostprocessing() { #endif } // probe values green - ComputeQOIsGreenProbingLine(); + if( _settings->GetProblemName() == PROBLEM_SymmetricHohlraum ) { + ComputeQOIsGreenProbingLine(); + } // Update time integral values on rank 0 if( _rank == 0 ) { _totalScalarOutflow += _curScalarOutflow * _dT; From 14f75b1de1c90073b7693af8335fb0bb62bc6c40 Mon Sep 17 00:00:00 2001 From: Steffen Date: Mon, 13 Jan 2025 16:33:13 -0500 Subject: [PATCH 21/21] fixes bug in qoi3 --- src/solvers/snsolver_hpc.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/solvers/snsolver_hpc.cpp b/src/solvers/snsolver_hpc.cpp index 8f93f516..6bec5280 100644 --- a/src/solvers/snsolver_hpc.cpp +++ b/src/solvers/snsolver_hpc.cpp @@ -1570,7 +1570,7 @@ void SNSolverHPC::SetProbingCellsLineGreen() { // left column double x1 = -0.2; - double y1 = -0.35 + j * block_size; + double y1 = -0.3 + j * block_size; double x2 = x1 + block_size; double y2 = y1 - block_size;