added restart files for SN hpc solver

include/common/config.hpp

@@ -48,9 +48,11 @@ class Config
     // std::vector<double> _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; }

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<double>& solution, const std::vector<double>& scalarFlux, int rank, int idx_iter );
+
+int LoadRestartSolution(
+    const std::string& baseInputFile, std::vector<double>& solution, std::vector<double>& scalarFlux, int rank, unsigned long nCells );
+
 // Matrix createSU2MeshFromImage( std::string imageName, std::string SU2Filename ); Deprecated
 
 #endif    // IO_H

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 */

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

src/common/io.cpp

@@ -11,9 +11,11 @@
 #include "toolboxes/errormessages.hpp"
 #include "toolboxes/textprocessingtoolbox.hpp"
 
+#include <filesystem>
 #include <fstream>
 #include <iomanip>
 #include <iostream>
+#include <sstream>
 
 #ifdef IMPORT_MPI
 #include <mpi.h>
@@ -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<double>& solution, const std::vector<double>& 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<double>& solution, std::vector<double>& 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<double> 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";

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<double> 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