diff --git a/CMakeLists.txt b/CMakeLists.txt index 9cae01d2..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 ) - ################################################# @@ -26,15 +25,19 @@ 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 ) + add_definitions(-DIMPORT_MPI) find_package( MPI REQUIRED ) include_directories( ${MPI_INCLUDE_PATH} ) + message(STATUS "C++ compiler: ${CMAKE_CXX_COMPILER}") + message(STATUS "MPI C++ compiler: ${MPI_CXX_COMPILER}") + endif() find_package( LAPACK REQUIRED ) @@ -63,6 +66,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/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/globalconstants.hpp b/include/common/globalconstants.hpp index 826bd338..b0ca20e3 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 }, @@ -195,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 }, @@ -219,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/io.hpp b/include/common/io.hpp index 377c1332..a42edead 100644 --- a/include/common/io.hpp +++ b/include/common/io.hpp @@ -43,6 +43,26 @@ 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, + 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, + double& totalAbsorptionHohlraumCenter, + double& totalAbsorptionHohlraumVertical, + double& totalAbsorptionHohlraumHorizontal, + double& totalAbsorption ); + // Matrix createSU2MeshFromImage( std::string imageName, std::string SU2Filename ); Deprecated #endif // IO_H 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/quadratures/qsphericaltessalation.hpp b/include/quadratures/qsphericaltessalation.hpp new file mode 100644 index 00000000..02f4aaaa --- /dev/null +++ b/include/quadratures/qsphericaltessalation.hpp @@ -0,0 +1,37 @@ +#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::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 ); + + std::vector, 3>> generate_tessellation( const std::array, 3>& triangle, int order ); +}; + +#endif // QSPHERICALTRIANGLE_H diff --git a/include/solvers/snsolver_hpc.hpp b/include/solvers/snsolver_hpc.hpp index 0046621e..4bdb32ed 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 */ + int _idx_start_iter; /*!< @brief index of first iteration */ // 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,14 +111,18 @@ 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 */ - 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 */ + 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 _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 */ @@ -198,8 +202,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/config.cpp b/src/common/config.cpp index 1483ac77..4e189192 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" @@ -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 @@ -784,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] ); @@ -894,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; @@ -917,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" @@ -946,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 ); } @@ -965,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" @@ -1019,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 ); } } } @@ -1207,7 +1221,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 +1380,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 b19ed8e4..aa10c89b 100644 --- a/src/common/io.cpp +++ b/src/common/io.cpp @@ -11,11 +11,13 @@ #include "toolboxes/errormessages.hpp" #include "toolboxes/textprocessingtoolbox.hpp" +#include #include #include #include +#include -#ifdef BUILD_MPI +#ifdef IMPORT_MPI #include #endif @@ -49,8 +51,13 @@ void ExportVTK( const std::string fileName, const std::vector>>& outputFields, const std::vector>& outputFieldNames, const Mesh* mesh ) { - int rank = 0; - // MPI_Comm_rank( MPI_COMM_WORLD, &rank ); + int rank = 0; + int nprocs = 1; +#ifdef IMPORT_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(); @@ -140,7 +147,6 @@ void ExportVTK( const std::string fileName, // auto log = spdlog::get( "event" ); // log->info( "Result successfully exported to '{0}'!", fileNameWithExt ); } - // MPI_Barrier( MPI_COMM_WORLD ); } Mesh* LoadSU2MeshFromFile( const Config* settings ) { @@ -352,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 @@ -362,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 @@ -373,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 @@ -400,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() ) { @@ -437,6 +443,107 @@ 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, + double totalAbsorptionHohlraumCenter, + double totalAbsorptionHohlraumVertical, + double totalAbsorptionHohlraumHorizontal, + double totalAbsorption ) { + // 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 binary writing + std::ofstream outFile( outputFile, std::ios::out | std::ios::binary ); + if( !outFile ) { + std::cerr << "Error opening restart solution file." << std::endl; + return; + } + + // 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(); + size_t scalarFlux_size = scalarFlux.size(); + outFile.write( reinterpret_cast( &solution_size ), sizeof( solution_size ) ); + outFile.write( reinterpret_cast( &scalarFlux_size ), sizeof( scalarFlux_size ) ); + + // 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(); +} + +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 ); + + // 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 + 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; + 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 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 + } + + // 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"; @@ -465,7 +572,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 ); @@ -505,7 +612,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 8ff2f2d7..ed7ef3c3 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 @@ -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 @@ -604,8 +605,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 +616,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 ); } } } @@ -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/main.cpp b/src/main.cpp index 07cfcc1f..7a0820cb 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 @@ -17,9 +17,8 @@ #include "solvers/solverbase.hpp" #ifdef BUILD_GUI -#include - #include "mainwindow.h" +#include #endif int main( int argc, char** argv ) { @@ -31,9 +30,18 @@ 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 ); + 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" ); #endif + std::string filename = ParseArguments( argc, argv ); // CD Load Settings from File @@ -66,7 +74,7 @@ int main( int argc, char** argv ) { } delete config; -#ifdef BUILD_MPI +#ifdef IMPORT_MPI MPI_Finalize(); #endif 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/quadratures/qsphericaltessalation.cpp b/src/quadratures/qsphericaltessalation.cpp new file mode 100644 index 00000000..5f9e7b89 --- /dev/null +++ b/src/quadratures/qsphericaltessalation.cpp @@ -0,0 +1,202 @@ +#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 +} + +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 { + 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 )]; + 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 }; + + // Generate the tessellation for the given order + std::vector, 3>> final_triangles = generate_tessellation( { A, B, C }, _order ); + + // 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 ) { + 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 = 4 * pow( GetOrder(), 2 ); } + +void QSphericalTessalation::SetConnectivity() { // TODO + // Not initialized for this quadrature. + VectorVectorU connectivity; + _connectivity = connectivity; +} + +bool QSphericalTessalation::CheckOrder() { return true; } + +// 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 d5816bb7..6bec5280 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,19 +11,18 @@ #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 - - _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 ); @@ -31,21 +30,37 @@ 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(); - - _nq = quad->GetNq(); - _nSys = _nq; - _localNSys = _nq / _numProcs; - _startSysIdx = _rank * _localNSys; - _endSysIdx = _rank * _localNSys + _localNSys; - if( _rank == _numProcs - 1 ) { - _localNSys += _nq % _numProcs; - _endSysIdx = _nSys; + _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 ); } + if( _numProcs == 1 ) { + _localNSys = _nSys; + _startSysIdx = 0; + _endSysIdx = _nSys; + } + 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(); @@ -59,8 +74,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 ); @@ -87,7 +102,8 @@ SNSolverHPC::SNSolverHPC( Config* settings ) { auto nodes = _mesh->GetNodes(); auto cells = _mesh->GetCells(); - for( unsigned idx_cell = 0; idx_cell < _nCells; idx_cell++ ) { +#pragma omp parallel for + for( unsigned long idx_cell = 0; idx_cell < _nCells; idx_cell++ ) { _areas[idx_cell] = areas[idx_cell]; } @@ -105,27 +121,32 @@ 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] = - 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 - 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 )] = @@ -139,7 +160,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 @@ -147,22 +168,6 @@ SNSolverHPC::SNSolverHPC( Config* settings ) { } // _mass += _scalarFlux[idx_cell] * _areas[idx_cell]; } -#ifdef BUILD_MPI - MPI_Barrier( MPI_COMM_WORLD ); -#endif - SetGhostCells(); - if( _rank == 0 ) { - PrepareScreenOutput(); // Screen Output - PrepareHistoryOutput(); // History Output - } -#ifdef BUILD_MPI - MPI_Barrier( MPI_COMM_WORLD ); -#endif - delete quad; - - // Initialiye QOIS - _mass = 0; - _rmsFlux = 0; // Lattice { @@ -187,11 +192,43 @@ 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; - if( _settings->GetProblemName() == PROBLEM_QuarterHohlraum ) n_probes = 2; - //_probingMomentsTimeIntervals = 10; + // if( _settings->GetProblemName() == PROBLEM_QuarterHohlraum ) n_probes = 2; + _probingMoments = std::vector( n_probes * 3, 0. ); // 10 time horizons if( _settings->GetProblemName() == PROBLEM_SymmetricHohlraum ) { @@ -202,12 +239,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; @@ -219,20 +256,21 @@ 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(); - 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 } } @@ -253,7 +291,8 @@ 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; } @@ -264,9 +303,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 @@ -274,6 +313,7 @@ void SNSolverHPC::Solve() { } } else { + ( _spatialOrder == 2 ) ? FluxOrder2() : FluxOrder1(); FVMUpdate(); } @@ -292,79 +332,69 @@ 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 --- if( _rank == 0 ) { DrawPostSolverOutput(); } } -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( idx_iter == (int)_nIter - 1 ) { // Last iteration write without suffix. - WriteVolumeOutput( idx_iter ); - ExportVTK( _settings->GetOutputFile(), _outputFields, _outputFieldNames, _mesh ); - } -} - 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++ ) { + 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 = @@ -398,26 +428,25 @@ void SNSolverHPC::FluxOrder2() { } else { #pragma omp simd - for( unsigned idx_sys = 0; idx_sys < _localNSys; idx_sys++ ) { + 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 ) { @@ -430,12 +459,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 )] += @@ -445,14 +475,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 )] ) ) ); } } } @@ -462,18 +492,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 ) { @@ -486,9 +516,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 )]; @@ -507,10 +537,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 )] - @@ -520,7 +550,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]; } @@ -529,12 +559,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 } @@ -560,7 +590,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 )] ) ) { @@ -570,7 +600,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 )]; @@ -603,16 +633,17 @@ 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 } } 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 )] * @@ -627,9 +658,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 )] * @@ -649,11 +680,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 )]; @@ -678,7 +709,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; @@ -707,7 +738,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 )]; @@ -730,15 +761,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 ); @@ -752,22 +783,23 @@ 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 - for( unsigned idx_probe = 0; idx_probe < n_probes; idx_probe++ ) { // Loop over probing cells +#ifndef IMPORT_MPI + 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 )]; } #endif } - + // probe values green + if( _settings->GetProblemName() == PROBLEM_SymmetricHohlraum ) { + ComputeQOIsGreenProbingLine(); + } // Update time integral values on rank 0 if( _rank == 0 ) { _totalScalarOutflow += _curScalarOutflow * _dT; @@ -820,7 +852,9 @@ double SNSolverHPC::ComputeTimeStep( double cfl ) const { } // 2D case double charSize = __DBL_MAX__; // minimum char size of all mesh cells in the mesh - for( unsigned j = 0; j < _nCells; j++ ) { + +#pragma omp parallel for reduction( min : charSize ) + for( unsigned long j = 0; j < _nCells; j++ ) { double currCharSize = sqrt( _areas[j] ); if( currCharSize < charSize ) { charSize = currCharSize; @@ -880,6 +914,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; } } @@ -928,14 +963,14 @@ 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++; } 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; } } @@ -983,7 +1018,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 )]; @@ -993,9 +1028,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--; @@ -1037,7 +1080,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] ) ) ) { @@ -1106,7 +1150,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 ); @@ -1116,7 +1160,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++; @@ -1144,9 +1196,9 @@ 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 - + // std::cout << lineToPrint << std::endl; if( _settings->GetHistoryOutputFrequency() != 0 && idx_iter % (unsigned)_settings->GetHistoryOutputFrequency() == 0 ) { log->info( lineToPrint ); } @@ -1231,9 +1283,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; } @@ -1244,9 +1296,9 @@ 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: @@ -1254,26 +1306,56 @@ void SNSolverHPC::WriteVolumeOutput( unsigned idx_iter ) { 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 = 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]; } } +#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 ); +#endif break; - default: ErrorMessages::Error( "Volume Output Group not defined for HPC SN Solver!", CURRENT_FUNCTION ); break; } } } } +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, + _totalAbsorptionHohlraumCenter, + _totalAbsorptionHohlraumVertical, + _totalAbsorptionHohlraumHorizontal, + _totalAbsorptionLattice ); + } + + 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 + } + } + 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(); @@ -1368,54 +1450,145 @@ 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 ); + // + // 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 = { _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 dx = 2 * ( horizontalLineWidth + verticalLineWidth ) / ( (double)_nProbingCellsLineGreen ); + double verticalLineWidth = std::abs( p1[1] - p4[1] ); + double horizontalLineWidth = std::abs( p1[0] - p2[0] ); - unsigned nHorizontalProbingCells = - (unsigned)std::ceil( _nProbingCellsLineGreen * ( horizontalLineWidth / ( horizontalLineWidth + verticalLineWidth ) ) ); - unsigned nVerticalProbingCells = _nProbingCellsLineGreen - nHorizontalProbingCells; + double pt_ratio_h = horizontalLineWidth / ( horizontalLineWidth + verticalLineWidth ); + double pt_ratio_v = verticalLineWidth / ( horizontalLineWidth + verticalLineWidth ); - _probingCellsLineGreen = std::vector( _nProbingCellsLineGreen ); + 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 side3 = linspace2D( _cornerLowerRightGreen, _cornerUpperRightGreen, nVerticalProbingCells ); - std::vector side4 = linspace2D( _cornerUpperRightGreen, _cornerUpperLeftGreen, 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 - // 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] ); + for( unsigned i = 0; i < nHorizontalProbingCells; ++i ) { + _probingCellsLineGreen[i] = side1[i]; + } + for( unsigned i = 0; i < nVerticalProbingCells; ++i ) { + _probingCellsLineGreen[i + nHorizontalProbingCells] = side2[i]; + } + for( unsigned i = 0; i < nHorizontalProbingCells; ++i ) { + _probingCellsLineGreen[i + nVerticalProbingCells + nHorizontalProbingCells] = side3[i]; + } + for( unsigned i = 0; i < nVerticalProbingCells; ++i ) { + _probingCellsLineGreen[i + nVerticalProbingCells + 2 * nHorizontalProbingCells] = side4[i]; + } - // double dx = 2 * ( horizontalLineWidth + verticalLineWidth ) / ( (double)_nProbingCellsLineGreen ); + // Block-wise approach + // Initialize the vector to store the corner points of each block + std::vector>> block_corners; - unsigned nHorizontalProbingCells = - (unsigned)std::ceil( _nProbingCellsLineGreen / 2 * ( horizontalLineWidth / ( horizontalLineWidth + verticalLineWidth ) ) ); - unsigned nVerticalProbingCells = _nProbingCellsLineGreen - nHorizontalProbingCells; + double block_size = 0.05; - _probingCellsLineGreen = std::vector( _nProbingCellsLineGreen ); + // Loop to fill the blocks + for( int i = 0; i < 8; ++i ) { // 8 blocks in the x-direction (horizontal) (upper) (left to right) - 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 }; + // Top row + double x1 = -0.2 + i * block_size; + double y1 = 0.4; + double x2 = x1 + block_size; + double y2 = y1 - block_size; - // 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() ); + 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) + // 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; + + // 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 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; + + 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) (down to up) + + // left column + double x1 = -0.2; + double y1 = -0.3 + 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 ); + } + + // Compute the probing cells for each block + for( int i = 0; i < _nProbingCellsBlocksGreen; i++ ) { + _probingCellsBlocksGreen.push_back( _mesh->GetCellsofRectangle( block_corners[i] ) ); + } } } @@ -1425,17 +1598,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]]; } - 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 ) { @@ -1453,7 +1634,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; 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++; 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 );