diff --git a/include/common/config.hpp b/include/common/config.hpp index 6957abe0..42e1992d 100644 --- a/include/common/config.hpp +++ b/include/common/config.hpp @@ -48,9 +48,11 @@ class Config // std::vector _1dIntegrationBounds; /*!< @brief Quadrature Order*/ // Mesh - unsigned _nCells; /*!< @brief Number of cells in the mesh */ - unsigned short _dim; /*!< @brief spatial dimensionality of the mesh/test case */ - bool _forcedConnectivityWrite; /*!< @brief If true, the meshconnectivity is always computed and written to .con file */ + unsigned _nCells; /*!< @brief Number of cells in the mesh */ + unsigned short _dim; /*!< @brief spatial dimensionality of the mesh/test case */ + bool _forcedConnectivityWrite; /*!< @brief If true, the meshconnectivity is always computed and written to .con file */ + bool _loadrestartSolution; /*!< @brief If true, the simulation loads a restart solution from file */ + unsigned long _saveRestartSolutionFrequency; /*!< @brief Frequency for saving restart solution to file */ // Boundary Conditions /*!< @brief List of all Pairs (marker, BOUNDARY_TYPE), e.g. (farfield,DIRICHLET). @@ -312,6 +314,8 @@ class Config unsigned GetNCells() const { return _nCells; } unsigned short GetDim() const { return _dim; } bool inline GetForcedConnectivity() const { return _forcedConnectivityWrite; } + bool inline GetLoadRestartSolution() const { return _loadrestartSolution; } + unsigned long inline GetSaveRestartSolutionFrequency() const { return _saveRestartSolutionFrequency; } // Solver Structure bool inline GetHPC() const { return _HPC; } diff --git a/include/common/io.hpp b/include/common/io.hpp index 377c1332..0c302183 100644 --- a/include/common/io.hpp +++ b/include/common/io.hpp @@ -43,6 +43,12 @@ std::string ParseArguments( int argc, char* argv[] ); void PrintLogHeader( std::string inputFile ); +void WriteRestartSolution( + const std::string& baseOutputFile, const std::vector& solution, const std::vector& scalarFlux, int rank, int idx_iter ); + +int LoadRestartSolution( + const std::string& baseInputFile, std::vector& solution, std::vector& scalarFlux, int rank, unsigned long nCells ); + // Matrix createSU2MeshFromImage( std::string imageName, std::string SU2Filename ); Deprecated #endif // IO_H diff --git a/include/solvers/snsolver_hpc.hpp b/include/solvers/snsolver_hpc.hpp index f37bb98b..7f25b882 100644 --- a/include/solvers/snsolver_hpc.hpp +++ b/include/solvers/snsolver_hpc.hpp @@ -32,7 +32,7 @@ class SNSolverHPC // Time unsigned long _nIter; /*!< @brief number of time steps, for non CSD, this equals _nEnergies, for _csd, _maxIter =_nEnergies-1*/ double _dT; /*!< @brief energy/time step size */ - + int _idx_start_iter; /*!< @brief index of first iteration */ // Mesh related members, memory optimized unsigned long _nCells; /*!< @brief number of spatial cells */ unsigned long _nSys; /*!< @brief number of equations in the transport system, i.e. num quad pts */ diff --git a/src/common/config.cpp b/src/common/config.cpp index 60ae732a..9ad9b107 100644 --- a/src/common/config.cpp +++ b/src/common/config.cpp @@ -226,6 +226,10 @@ void Config::SetConfigOptions() { /*! @brief FORCE_CONNECTIVITY_RECOMPUTE \n DESCRIPTION:If true, mesh recomputes connectivity instead of loading from file \n DEFAULT false * \ingroup Config.*/ AddBoolOption( "FORCE_CONNECTIVITY_RECOMPUTE", _forcedConnectivityWrite, false ); + /*! @brief RESTART_SOLUTION \n DESCRIPTION:If true, simulation loads a restart solution from file \n DEFAULT false + * \ingroup Config.*/ + AddBoolOption( "LOAD_RESTART_SOLUTION", _loadrestartSolution, false ); + AddUnsignedLongOption( "SAVE_RESTART_SOLUTION_FREQUENCY", _saveRestartSolutionFrequency, false ); // Quadrature relatated options /*! @brief QUAD_TYPE \n DESCRIPTION: Type of Quadrature rule \n Options: see @link QUAD_NAME \endlink \n DEFAULT: QUAD_MonteCarlo diff --git a/src/common/io.cpp b/src/common/io.cpp index 8020d396..cdc8a307 100644 --- a/src/common/io.cpp +++ b/src/common/io.cpp @@ -11,9 +11,11 @@ #include "toolboxes/errormessages.hpp" #include "toolboxes/textprocessingtoolbox.hpp" +#include #include #include #include +#include #ifdef IMPORT_MPI #include @@ -356,7 +358,7 @@ void LoadConnectivityFromFile( const std::string inputFile, for( unsigned j = 0; j < correctedNodesPerCell; ++j ) { std::getline( iss, line, ',' ); std::istringstream converter( line ); - converter >> std::fixed >> setprecision( 12 ) >> cellNeighbors[i][j]; + converter >> std::fixed >> setprecision( 15 ) >> cellNeighbors[i][j]; } // Load cellInterfaceMidPoints @@ -366,8 +368,8 @@ void LoadConnectivityFromFile( const std::string inputFile, for( unsigned k = 0; k < nDim; ++k ) { std::getline( iss, line, ',' ); std::istringstream converter( line ); - converter >> std::fixed >> setprecision( 12 ) >> cellInterfaceMidPoints[i][j][k]; // Replace with appropriate member of Vector - // std::cout << std::fixed << setprecision( 12 ) << cellInterfaceMidPoints[i][j][k] << std::endl; + converter >> std::fixed >> setprecision( 15 ) >> cellInterfaceMidPoints[i][j][k]; // Replace with appropriate member of Vector + // std::cout << std::fixed << setprecision( 15 ) << cellInterfaceMidPoints[i][j][k] << std::endl; } } // Load cellNormals @@ -377,7 +379,7 @@ void LoadConnectivityFromFile( const std::string inputFile, for( unsigned k = 0; k < nDim; ++k ) { std::getline( iss, line, ',' ); std::istringstream converter( line ); - converter >> std::fixed >> setprecision( 12 ) >> cellNormals[i][j][k]; // Replace with appropriate member of Vector + converter >> std::fixed >> setprecision( 15 ) >> cellNormals[i][j][k]; // Replace with appropriate member of Vector } } // Load cellBoundaryTypes @@ -404,7 +406,7 @@ void WriteConnecitivityToFile( const std::string outputFile, // cellBoundaryTypes (1 element), (tranlated from BOUNDARY_TYPE to unsigned) std::ofstream outFile( outputFile ); - outFile << std::fixed << setprecision( 12 ); + outFile << std::fixed << setprecision( 15 ); // const std::size_t bufferSize = 10000; // Adjust as needed // outFile.rdbuf()->pubsetbuf( 0, bufferSize ); if( outFile.is_open() ) { @@ -441,6 +443,83 @@ void WriteConnecitivityToFile( const std::string outputFile, } } +void WriteRestartSolution( + const std::string& baseOutputFile, const std::vector& solution, const std::vector& scalarFlux, int rank, int idx_iter ) { + // Generate filename with rank number + std::string outputFile = baseOutputFile + "_restart_rank_" + std::to_string( rank ); + + // Check if the file exists and delete it + if( std::filesystem::exists( outputFile ) ) { + std::filesystem::remove( outputFile ); + } + + // Open the file for writing (automatically creates a new file) + std::ofstream outFile( outputFile, std::ios::out ); + if( !outFile ) { + ErrorMessages::Error( "Error opening restart solution file.", CURRENT_FUNCTION ); + return; + } + + // Write the iteration number as the first line + outFile << idx_iter << "\n"; + outFile << std::fixed << std::setprecision( 15 ); + + // Write each element of the solution to the file on a new line + for( const double& value : solution ) { + outFile << value << "\n"; + } + + // Append each element of the scalarFlux to the file on a new line + for( const double& fluxValue : scalarFlux ) { + outFile << fluxValue << "\n"; + } + + // Close the file + outFile.close(); +} + +int LoadRestartSolution( + const std::string& baseInputFile, std::vector& solution, std::vector& scalarFlux, int rank, unsigned long nCells ) { + // Generate filename with rank number + std::string inputFile = baseInputFile + "_restart_rank_" + std::to_string( rank ); + + // Open the file for reading + std::ifstream inFile( inputFile ); + if( !inFile ) { + std::cerr << "Error opening restart solution file for reading." << std::endl; + return 0; // Signal failure + } + + // Read the iteration number from the first line + std::string line; + std::getline( inFile, line ); + int idx_iter = std::stoi( line ); // Convert the line to an integer + + // Read data into a temporary vector + std::vector tempData; + double value; + while( std::getline( inFile, line ) ) { + std::istringstream converter( line ); + converter >> std::fixed >> std::setprecision( 15 ) >> value; + tempData.push_back( value ); + } + + // Close the file + inFile.close(); + + // Verify that we have at least nCells entries to populate scalarFlux + if( tempData.size() < nCells ) { + std::cerr << "Not enough data to populate scalar flux vector." << std::endl; + return 0; // Signal failure + } + + // Allocate the last nCells entries to scalarFlux and the rest to solution + solution.assign( tempData.begin(), tempData.end() - nCells ); + scalarFlux.assign( tempData.end() - nCells, tempData.end() ); + + return idx_iter; // Return the iteration index +} + std::string ParseArguments( int argc, char* argv[] ) { std::string inputFile; std::string usage_help = "\nUsage: " + std::string( argv[0] ) + " inputfile\n"; diff --git a/src/solvers/snsolver_hpc.cpp b/src/solvers/snsolver_hpc.cpp index c361bb51..ac5759e9 100644 --- a/src/solvers/snsolver_hpc.cpp +++ b/src/solvers/snsolver_hpc.cpp @@ -21,9 +21,9 @@ SNSolverHPC::SNSolverHPC( Config* settings ) { _rank = 0; #endif - _settings = settings; - _currTime = 0.0; - + _settings = settings; + _currTime = 0.0; + _idx_start_iter = 0; _nOutputMoments = 2; // Currently only u_1 (x direction) and u_1 (y direction) are supported // Create Mesh _mesh = LoadSU2MeshFromFile( settings ); @@ -159,6 +159,9 @@ SNSolverHPC::SNSolverHPC( Config* settings ) { // _mass += _scalarFlux[idx_cell] * _areas[idx_cell]; } + if( _settings->GetLoadRestartSolution() ) + _idx_start_iter = LoadRestartSolution( _settings->GetOutputFile(), _sol, _scalarFlux, _rank, _nCells ) + 1; + #ifdef IMPORT_MPI MPI_Barrier( MPI_COMM_WORLD ); #endif @@ -265,7 +268,7 @@ void SNSolverHPC::Solve() { std::chrono::duration duration; // Loop over energies (pseudo-time of continuous slowing down approach) - for( unsigned iter = 0; iter < _nIter; iter++ ) { + for( unsigned iter = (unsigned)_idx_start_iter; iter < _nIter; iter++ ) { if( iter == _nIter - 1 ) { // last iteration _dT = _settings->GetTEnd() - iter * _dT; @@ -310,6 +313,7 @@ void SNSolverHPC::Solve() { #ifdef BUILD_MPI MPI_Barrier( MPI_COMM_WORLD ); #endif + PrintVolumeOutput( iter ); #ifdef BUILD_MPI MPI_Barrier( MPI_COMM_WORLD ); @@ -352,7 +356,7 @@ void SNSolverHPC::FluxOrder2() { for( unsigned long idx_cell = 0; idx_cell < _nCells; ++idx_cell ) { // Compute Limiter if( _cellBoundaryTypes[idx_cell] == BOUNDARY_TYPE::NONE ) { // skip ghost cells - // #pragma omp simd +#pragma omp simd for( unsigned long idx_sys = 0; idx_sys < _localNSys; idx_sys++ ) { double gaussPoint = 0; @@ -400,7 +404,7 @@ void SNSolverHPC::FluxOrder2() { } } else { - // #pragma omp simd +#pragma omp simd for( unsigned long idx_sys = 0; idx_sys < _localNSys; idx_sys++ ) { _limiter[Idx2D( idx_cell, idx_sys, _localNSys )] = 0.; // limiter should be zero at boundary _solDx[Idx3D( idx_cell, idx_sys, 0, _localNSys, _nDim )] = 0.; @@ -1282,7 +1286,12 @@ void SNSolverHPC::WriteVolumeOutput( unsigned idx_iter ) { } void SNSolverHPC::PrintVolumeOutput( int idx_iter ) { - if( _settings->GetVolumeOutputFrequency() != 0 && idx_iter % (unsigned)_settings->GetVolumeOutputFrequency() == 0 ) { + if( _settings->GetSaveRestartSolutionFrequency() != 0 && idx_iter % (int)_settings->GetSaveRestartSolutionFrequency() == 0 ) { + std::cout << "Saving restart solution at iteration " << idx_iter << std::endl; + WriteRestartSolution( _settings->GetOutputFile(), _sol, _scalarFlux, _rank, idx_iter ); + } + + if( _settings->GetVolumeOutputFrequency() != 0 && idx_iter % (int)_settings->GetVolumeOutputFrequency() == 0 ) { WriteVolumeOutput( idx_iter ); if( _rank == 0 ) { ExportVTK( _settings->GetOutputFile() + "_" + std::to_string( idx_iter ), _outputFields, _outputFieldNames, _mesh ); // slow