diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp b/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp index 362e2d51c4b0..051f8d86b073 100644 --- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp +++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp @@ -362,6 +362,35 @@ class CFVMFlowSolverBase : public CSolver { void LoadRestart_impl(CGeometry **geometry, CSolver ***solver, CConfig *config, int iter, bool update_geo, su2double* RestartSolution = nullptr, unsigned short nVar_Restart = 0); + /*! + * \brief Load restart solution fields using template specialization based on flow regime. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] solver - Container vector with all of the solvers. + * \param[in] config - Definition of the particular problem. + * \param[in] update_geo - Flag for updating coords and grid velocity. + * \param[in] static_fsi - Flag for static FSI problems. + * \param[in] steady_restart - Flag for steady restart. + * \param[in] SolutionRestart - Optional solution restart buffer. + * \param[in] nVar_Restart - Number of restart variables. + * \param[in] restart_filename - Name of the restart file. + */ + void LoadRestartSolutionFields(CGeometry **geometry, CSolver ***solver, CConfig *config, + bool update_geo, bool static_fsi, bool steady_restart, + su2double* SolutionRestart, unsigned short nVar_Restart, + const string& restart_filename); + + /*! + * \brief Helper function to load grid data (coordinates and velocities) from restart file. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] iPoint_Local - Local point index. + * \param[in] baseIndex - Base index in restart data array. + * \param[in] update_geo - Flag for updating coords and grid velocity. + * \param[in] static_fsi - Flag for static FSI problems. + * \param[in] steady_restart - Flag for steady restart. + */ + void LoadRestartGridData(CGeometry **geometry, long iPoint_Local, unsigned long baseIndex, + bool update_geo, bool static_fsi, bool steady_restart); + /*! * \brief Generic implementation to compute the time step based on CFL and conv/visc eigenvalues. * \param[in] geometry - Geometrical definition of the problem. diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl index 8dc0e4b63931..e7a39a98b99d 100644 --- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl +++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl @@ -30,7 +30,10 @@ #include "../gradients/computeGradientsLeastSquares.hpp" #include "../limiters/computeLimiters.hpp" #include "../numerics_simd/CNumericsSIMD.hpp" +#include "../../Common/include/CConfig.hpp" +#include "../../Common/include/geometry/CGeometry.hpp" #include "CFVMFlowSolverBase.hpp" +#include "CRestartFieldNames.hpp" template void CFVMFlowSolverBase::AeroCoeffsArray::allocate(int size) { @@ -821,134 +824,276 @@ void CFVMFlowSolverBase::UpdateCustomBoundaryConditions( } template -void CFVMFlowSolverBase::LoadRestart_impl(CGeometry **geometry, CSolver ***solver, CConfig *config, int iter, - bool update_geo, su2double* SolutionRestart, - unsigned short nVar_Restart) { - /*--- Restart the solution from file information ---*/ +void CFVMFlowSolverBase::LoadRestartSolutionFields(CGeometry** geometry, CSolver*** solver, CConfig* config, + bool update_geo, bool static_fsi, bool steady_restart, + su2double* SolutionRestart, unsigned short nVar_Restart, + const string& restart_filename) { + /*--- Helper lambda to find field index ---*/ + auto FindFieldIndex = [&](const string& name) -> int { + auto it = std::find(this->fields.begin(), this->fields.end(), name); + if (it != this->fields.end()) return (int)std::distance(this->fields.begin(), it); + string qname = "\"" + name + "\""; + it = std::find(this->fields.begin(), this->fields.end(), qname); + if (it != this->fields.end()) return (int)std::distance(this->fields.begin(), it); + return -1; + }; + + unsigned long counter = 0; + + /*--- Incompressible Flow Logic ---*/ + if constexpr (R == ENUM_REGIME::INCOMPRESSIBLE) { + const bool energy = config->GetEnergy_Equation(); + const bool weakly_coupled_heat = config->GetWeakly_Coupled_Heat(); + const bool flamelet = (config->GetKind_FluidModel() == FLUID_FLAMELET); + + /*--- Find field indices for solution variables ---*/ + int fieldIdx_Pressure = FindFieldIndex(RestartFieldNames::PRESSURE); + vector fieldIdx_Velocity(nDim, -1); + for (auto iDim = 0u; iDim < nDim; iDim++) { + fieldIdx_Velocity[iDim] = FindFieldIndex(RestartFieldNames::GetVelocityName(iDim)); + } + int fieldIdx_Temperature = -1; + if (energy || weakly_coupled_heat || flamelet) { + fieldIdx_Temperature = FindFieldIndex(RestartFieldNames::TEMPERATURE); + if (fieldIdx_Temperature < 0 && rank == MASTER_NODE) { + cout << "\nWARNING: The restart file does not contain " << RestartFieldNames::TEMPERATURE << " field." << endl; + cout << "Temperature will be initialized from existing solution values." << endl; + } + } - string restart_filename = config->GetSolution_FileName(); - const bool static_fsi = ((config->GetTime_Marching() == TIME_MARCHING::STEADY) && config->GetFSI_Simulation()); + /*--- Load data from the restart into correct containers. ---*/ + for (auto iPoint_Global = 0ul; iPoint_Global < geometry[MESH_0]->GetGlobal_nPointDomain(); iPoint_Global++) { + const auto iPoint_Local = geometry[MESH_0]->GetGlobal_to_Local_Point(iPoint_Global); - /*--- To make this routine safe to call in parallel most of it can only be executed by one thread. ---*/ - BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS { + if (iPoint_Local > -1) { + auto baseIndex = counter * Restart_Vars[1]; - if (nVar_Restart == 0) nVar_Restart = nVar; + su2double* targetBuffer; - /*--- Skip coordinates ---*/ + /*--- If SolutionRestart is provided (non-null), use it. + Otherwise, write directly to nodes (which already contain defaults). ---*/ + if (SolutionRestart != nullptr) { + targetBuffer = SolutionRestart; // Caution: assumes SolutionRestart is per-point or advanced externally? + // But if SolutionRestart is provided, we should probably initialize it with defaults too. + const su2double* existingSolution = nodes->GetSolution(iPoint_Local); + for (auto iVar = 0u; iVar < nVar; iVar++) targetBuffer[iVar] = existingSolution[iVar]; + } else { + targetBuffer = nodes->GetSolution(iPoint_Local); + } - unsigned short skipVars = nDim; + /*--- Load Pressure (index 0) ---*/ + if (fieldIdx_Pressure >= 0) { + targetBuffer[0] = Restart_Data[baseIndex + fieldIdx_Pressure]; + } else { + if (rank == MASTER_NODE && counter == 0) { + cout << "WARNING: " << RestartFieldNames::PRESSURE + << " field not found in restart file, using existing solution value." << endl; + } + } - /*--- Read the restart data from either an ASCII or binary SU2 file. ---*/ - if (config->GetRead_Binary_Restart()) { - restart_filename = config->GetFilename(restart_filename, ".dat", iter); - Read_SU2_Restart_Binary(geometry[MESH_0], config, restart_filename); - } else { - restart_filename = config->GetFilename(restart_filename, ".csv", iter); - Read_SU2_Restart_ASCII(geometry[MESH_0], config, restart_filename); - } + /*--- Load Velocity components (indices 1..nDim) ---*/ + for (auto iDim = 0u; iDim < nDim; iDim++) { + if (fieldIdx_Velocity[iDim] >= 0) { + targetBuffer[iDim + 1] = Restart_Data[baseIndex + fieldIdx_Velocity[iDim]]; + } else { + if (rank == MASTER_NODE && counter == 0) { + cout << "WARNING: " << RestartFieldNames::GetVelocityName(iDim) + << " field not found in restart file, using existing solution value." << endl; + } + } + } - bool steady_restart = config->GetSteadyRestart(); - if (update_geo && dynamic_grid) { - auto notFound = fields.end(); - if (find(fields.begin(), notFound, string("\"Grid_Velocity_x\"")) == notFound) { - if (rank == MASTER_NODE) - cout << "\nWARNING: The restart file does not contain grid velocities, these will be set to zero.\n" << endl; - steady_restart = true; + /*--- Load Temperature (index nDim+1) if energy equation is enabled ---*/ + if (energy || weakly_coupled_heat || flamelet) { + if (fieldIdx_Temperature >= 0) { + targetBuffer[nDim + 1] = Restart_Data[baseIndex + fieldIdx_Temperature]; + } else { + if (rank == MASTER_NODE && counter == 0) { + cout << "WARNING: " << RestartFieldNames::TEMPERATURE + << " field not found in restart file, using existing solution value." << endl; + } + } + } + + /*--- If we used a separate buffer, copy it back if needed (only if SolutionRestart logic requires it) ---*/ + if (SolutionRestart != nullptr) { + nodes->SetSolution(iPoint_Local, targetBuffer); + } + + /*--- Handle grid velocities and coordinates ---*/ + LoadRestartGridData(geometry, iPoint_Local, baseIndex, update_geo, static_fsi, steady_restart); + + counter++; } } + } - /*--- Load data from the restart into correct containers. ---*/ + /*--- Compressible Flow Logic ---*/ + else if constexpr (R == ENUM_REGIME::COMPRESSIBLE) { + /*--- Find field indices for solution variables ---*/ + int fieldIdx_Density = FindFieldIndex(RestartFieldNames::DENSITY); + vector fieldIdx_Momentum(nDim, -1); + for (auto iDim = 0u; iDim < nDim; iDim++) { + fieldIdx_Momentum[iDim] = FindFieldIndex(RestartFieldNames::GetMomentumName(iDim)); + } + int fieldIdx_Energy = FindFieldIndex(RestartFieldNames::ENERGY); - unsigned long counter = 0; + /*--- Load data from the restart into correct containers. ---*/ for (auto iPoint_Global = 0ul; iPoint_Global < geometry[MESH_0]->GetGlobal_nPointDomain(); iPoint_Global++) { - - /*--- Retrieve local index. If this node from the restart file lives - on the current processor, we will load and instantiate the vars. ---*/ - const auto iPoint_Local = geometry[MESH_0]->GetGlobal_to_Local_Point(iPoint_Global); if (iPoint_Local > -1) { + auto baseIndex = counter * Restart_Vars[1]; - /*--- We need to store this point's data, so jump to the correct - offset in the buffer of data from the restart file and load it. ---*/ + su2double* targetBuffer; - auto index = counter * Restart_Vars[1] + skipVars; + if (SolutionRestart != nullptr) { + targetBuffer = SolutionRestart; + const su2double* existingSolution = nodes->GetSolution(iPoint_Local); + for (auto iVar = 0u; iVar < nVar; iVar++) targetBuffer[iVar] = existingSolution[iVar]; + } else { + targetBuffer = nodes->GetSolution(iPoint_Local); + } - if (SolutionRestart == nullptr) { - for (auto iVar = 0u; iVar < nVar_Restart; iVar++) - nodes->SetSolution(iPoint_Local, iVar, Restart_Data[index+iVar]); + /*--- Load Density (index 0) ---*/ + if (fieldIdx_Density >= 0) { + targetBuffer[0] = Restart_Data[baseIndex + fieldIdx_Density]; + } else { + if (rank == MASTER_NODE && counter == 0) { + cout << "WARNING: " << RestartFieldNames::DENSITY + << " field not found in restart file, using existing solution value." << endl; + } } - else { - /*--- Used as buffer, allows defaults for nVar > nVar_Restart. ---*/ - for (auto iVar = 0u; iVar < nVar_Restart; iVar++) - SolutionRestart[iVar] = Restart_Data[index + iVar]; - nodes->SetSolution(iPoint_Local, SolutionRestart); + + /*--- Load Momentum components (indices 1..nDim) ---*/ + for (auto iDim = 0u; iDim < nDim; iDim++) { + if (fieldIdx_Momentum[iDim] >= 0) { + targetBuffer[iDim + 1] = Restart_Data[baseIndex + fieldIdx_Momentum[iDim]]; + } else { + if (rank == MASTER_NODE && counter == 0) { + cout << "WARNING: " << RestartFieldNames::GetMomentumName(iDim) + << " field not found in restart file, using existing solution value." << endl; + } + } } - /*--- For dynamic meshes, read in and store the - grid coordinates and grid velocities for each node. ---*/ + /*--- Load Energy (index nDim+1) ---*/ + if (fieldIdx_Energy >= 0) { + targetBuffer[nDim + 1] = Restart_Data[baseIndex + fieldIdx_Energy]; + } else { + if (rank == MASTER_NODE && counter == 0) { + cout << "WARNING: " << RestartFieldNames::ENERGY + << " field not found in restart file, using existing solution value." << endl; + } + } - if (dynamic_grid && update_geo) { + if (SolutionRestart != nullptr) { + nodes->SetSolution(iPoint_Local, targetBuffer); + } - /*--- Read in the next 2 or 3 variables which are the grid velocities ---*/ - /*--- If we are restarting the solution from a previously computed static calculation (no grid movement) ---*/ - /*--- the grid velocities are set to 0. This is useful for FSI computations ---*/ + /*--- Handle grid velocities and coordinates ---*/ + LoadRestartGridData(geometry, iPoint_Local, baseIndex, update_geo, static_fsi, steady_restart); - /*--- Rewind the index to retrieve the Coords. ---*/ - index = counter * Restart_Vars[1]; - const auto* Coord = &Restart_Data[index]; + counter++; + } + } + } else { + /*--- Enforce specialization expectation ---*/ + static_assert(R == ENUM_REGIME::INCOMPRESSIBLE || R == ENUM_REGIME::COMPRESSIBLE, + "LoadRestartSolutionFields must be specialized for the flow regime"); + } - su2double GridVel[MAXNDIM] = {0.0}; - if (!steady_restart) { - /*--- Move the index forward to get the grid velocities. ---*/ - index += skipVars + nVar_Restart + config->GetnTurbVar(); - for (auto iDim = 0u; iDim < nDim; iDim++) { GridVel[iDim] = Restart_Data[index+iDim]; } - } + if (counter != nPointDomain) { + SU2_MPI::Error(string("The solution file ") + restart_filename + string(" does not match with the mesh file.\n") + + string("This can be caused by empty lines at the end of the file."), + CURRENT_FUNCTION); + } +} - for (auto iDim = 0u; iDim < nDim; iDim++) { - geometry[MESH_0]->nodes->SetCoord(iPoint_Local, iDim, Coord[iDim]); - geometry[MESH_0]->nodes->SetGridVel(iPoint_Local, iDim, GridVel[iDim]); - } +/*! + * \brief Helper function to load grid data (coordinates and velocities) from restart file. + */ +template +void CFVMFlowSolverBase::LoadRestartGridData(CGeometry** geometry, long iPoint_Local, unsigned long baseIndex, + bool update_geo, bool static_fsi, bool steady_restart) { + auto index = baseIndex; + const auto* Coord = &Restart_Data[index]; + + if (dynamic_grid && update_geo) { + su2double GridVel[MAXNDIM] = {0.0}; + if (!steady_restart) { + for (auto iDim = 0u; iDim < nDim; iDim++) { + int fieldIdx_GridVel = FindFieldIndex(RestartFieldNames::GetGridVelocityName(iDim)); + if (fieldIdx_GridVel >= 0) { + GridVel[iDim] = Restart_Data[baseIndex + fieldIdx_GridVel]; } + } + } - /*--- For static FSI problems, grid_movement is 0 but we need to read in and store the - grid coordinates for each node (but not the grid velocities, as there are none). ---*/ + for (auto iDim = 0u; iDim < nDim; iDim++) { + geometry[MESH_0]->nodes->SetCoord(iPoint_Local, iDim, Coord[iDim]); + geometry[MESH_0]->nodes->SetGridVel(iPoint_Local, iDim, GridVel[iDim]); + } + } - if (static_fsi && update_geo) { - /*--- Rewind the index to retrieve the Coords. ---*/ - index = counter*Restart_Vars[1]; - const auto* Coord = &Restart_Data[index]; + if (static_fsi && update_geo) { + for (auto iDim = 0u; iDim < nDim; iDim++) { + geometry[MESH_0]->nodes->SetCoord(iPoint_Local, iDim, Coord[iDim]); + } + } +} - for (auto iDim = 0u; iDim < nDim; iDim++) { - geometry[MESH_0]->nodes->SetCoord(iPoint_Local, iDim, Coord[iDim]); - } - } +template +void CFVMFlowSolverBase::LoadRestart_impl(CGeometry** geometry, CSolver*** solver, CConfig* config, int iter, + bool update_geo, su2double* SolutionRestart, + unsigned short nVar_Restart) { + /*--- Restart the solution from file information ---*/ - /*--- Increment the overall counter for how many points have been loaded. ---*/ - counter++; - } - } + string restart_filename = config->GetSolution_FileName(); + const bool static_fsi = ((config->GetTime_Marching() == TIME_MARCHING::STEADY) && config->GetFSI_Simulation()); + + /*--- To make this routine safe to call in parallel most of it can only be executed by one thread. ---*/ + BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS { + if (nVar_Restart == 0) nVar_Restart = nVar; - /*--- Detect a wrong solution file ---*/ - if (counter != nPointDomain) { - SU2_MPI::Error(string("The solution file ") + restart_filename + string(" does not match with the mesh file.\n") + - string("This can be caused by empty lines at the end of the file."), CURRENT_FUNCTION); + + /*--- Read the restart data from either an ASCII or binary SU2 file. ---*/ + if (config->GetRead_Binary_Restart()) { + restart_filename = config->GetFilename(restart_filename, ".dat", iter); + Read_SU2_Restart_Binary(geometry[MESH_0], config, restart_filename); + } else { + restart_filename = config->GetFilename(restart_filename, ".csv", iter); + Read_SU2_Restart_ASCII(geometry[MESH_0], config, restart_filename); } + + bool steady_restart = config->GetSteadyRestart(); + if (update_geo && dynamic_grid) { + auto notFound = fields.end(); + string gridVelName = "\"" + string(RestartFieldNames::GRID_VELOCITY_X) + "\""; + if (find(fields.begin(), notFound, gridVelName) == notFound) { + if (rank == MASTER_NODE) + cout << "\nWARNING: The restart file does not contain grid velocities, these will be set to zero.\n" << endl; + steady_restart = true; + } + } + + /*--- Use template specialization to load solution based on flow regime ---*/ + LoadRestartSolutionFields(geometry, solver, config, update_geo, static_fsi, steady_restart, SolutionRestart, + nVar_Restart, restart_filename); } END_SU2_OMP_SAFE_GLOBAL_ACCESS /*--- Update the geometry for flows on deforming meshes. ---*/ if ((dynamic_grid || static_fsi) && update_geo) { - CGeometry::UpdateGeometry(geometry, config); if (dynamic_grid) { for (auto iMesh = 0u; iMesh <= config->GetnMGLevels(); iMesh++) { - /*--- Compute the grid velocities on the coarser levels. ---*/ - if (iMesh) geometry[iMesh]->SetRestricted_GridVelocity(geometry[iMesh - 1]); + if (iMesh) + geometry[iMesh]->SetRestricted_GridVelocity(geometry[iMesh - 1]); else { geometry[MESH_0]->InitiateComms(geometry[MESH_0], config, MPI_QUANTITIES::GRID_VELOCITY); geometry[MESH_0]->CompleteComms(geometry[MESH_0], config, MPI_QUANTITIES::GRID_VELOCITY); @@ -968,22 +1113,22 @@ void CFVMFlowSolverBase::LoadRestart_impl(CGeometry **geometry, CSolver ** /*--- For turbulent/species simulations the flow preprocessing is done by the turbulence/species solver * after it loads its variables (they are needed to compute flow primitives). In case turbulence and species, the * species solver does all the Pre-/Postprocessing. ---*/ - if (config->GetKind_Turb_Model() == TURB_MODEL::NONE && - config->GetKind_Species_Model() == SPECIES_MODEL::NONE) { - solver[MESH_0][FLOW_SOL]->Preprocessing(geometry[MESH_0], solver[MESH_0], config, MESH_0, NO_RK_ITER, RUNTIME_FLOW_SYS, false); + if (config->GetKind_Turb_Model() == TURB_MODEL::NONE && config->GetKind_Species_Model() == SPECIES_MODEL::NONE) { + solver[MESH_0][FLOW_SOL]->Preprocessing(geometry[MESH_0], solver[MESH_0], config, MESH_0, NO_RK_ITER, + RUNTIME_FLOW_SYS, false); } /*--- Interpolate the solution down to the coarse multigrid levels ---*/ for (auto iMesh = 1u; iMesh <= config->GetnMGLevels(); iMesh++) { - MultigridRestriction(*geometry[iMesh - 1], solver[iMesh - 1][FLOW_SOL]->GetNodes()->GetSolution(), - *geometry[iMesh], solver[iMesh][FLOW_SOL]->GetNodes()->GetSolution()); + MultigridRestriction(*geometry[iMesh - 1], solver[iMesh - 1][FLOW_SOL]->GetNodes()->GetSolution(), *geometry[iMesh], + solver[iMesh][FLOW_SOL]->GetNodes()->GetSolution()); solver[iMesh][FLOW_SOL]->InitiateComms(geometry[iMesh], config, MPI_QUANTITIES::SOLUTION); solver[iMesh][FLOW_SOL]->CompleteComms(geometry[iMesh], config, MPI_QUANTITIES::SOLUTION); - if (config->GetKind_Turb_Model() == TURB_MODEL::NONE && - config->GetKind_Species_Model() == SPECIES_MODEL::NONE) { - solver[iMesh][FLOW_SOL]->Preprocessing(geometry[iMesh], solver[iMesh], config, iMesh, NO_RK_ITER, RUNTIME_FLOW_SYS, false); + if (config->GetKind_Turb_Model() == TURB_MODEL::NONE && config->GetKind_Species_Model() == SPECIES_MODEL::NONE) { + solver[iMesh][FLOW_SOL]->Preprocessing(geometry[iMesh], solver[iMesh], config, iMesh, NO_RK_ITER, + RUNTIME_FLOW_SYS, false); } } @@ -996,9 +1141,8 @@ void CFVMFlowSolverBase::LoadRestart_impl(CGeometry **geometry, CSolver ** } /*--- Go back to single threaded execution. ---*/ - BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS - { - /*--- Delete the class memory that is used to load the restart. ---*/ + BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS { + /*--- Delete the class memory that is used to load the restart. ---*/ Restart_Vars = decltype(Restart_Vars){}; Restart_Data = decltype(Restart_Data){}; diff --git a/SU2_CFD/include/solvers/CRestartFieldNames.hpp b/SU2_CFD/include/solvers/CRestartFieldNames.hpp new file mode 100644 index 000000000000..8ad3b7344e45 --- /dev/null +++ b/SU2_CFD/include/solvers/CRestartFieldNames.hpp @@ -0,0 +1,101 @@ +/*! + * \file CRestartFieldNames.hpp + * \brief Constants for restart file field names. + * \note These names must match those used in output classes (CFlowIncOutput, CFlowCompOutput, etc.) + * \version 8.3.0 "Harrier" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2025, SU2 Contributors (cf. AUTHORS.md) + * + * SU2 is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * SU2 is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with SU2. If not, see . + */ + +#pragma once + +#include + +/*! + * \namespace RestartFieldNames + * \brief Namespace containing constants for restart file field names. + * \note Field names in restart files match output field names (without quotes in output, with quotes in restart). + */ +namespace RestartFieldNames { + /*--- Solution field names (must match output class field names) ---*/ + constexpr const char* PRESSURE = "Pressure"; + constexpr const char* TEMPERATURE = "Temperature"; + constexpr const char* DENSITY = "Density"; + constexpr const char* ENERGY = "Energy"; + + /*--- Velocity component field names ---*/ + constexpr const char* VELOCITY_X = "Velocity_x"; + constexpr const char* VELOCITY_Y = "Velocity_y"; + constexpr const char* VELOCITY_Z = "Velocity_z"; + + /*--- Momentum component field names ---*/ + constexpr const char* MOMENTUM_X = "Momentum_x"; + constexpr const char* MOMENTUM_Y = "Momentum_y"; + constexpr const char* MOMENTUM_Z = "Momentum_z"; + + /*--- Grid velocity component field names ---*/ + constexpr const char* GRID_VELOCITY_X = "Grid_Velocity_x"; + constexpr const char* GRID_VELOCITY_Y = "Grid_Velocity_y"; + constexpr const char* GRID_VELOCITY_Z = "Grid_Velocity_z"; + + /*--- Turbulence variable field names ---*/ + constexpr const char* NU_TILDE = "Nu_Tilde"; + constexpr const char* TURB_KIN_ENERGY = "Turb_Kin_Energy"; + constexpr const char* OMEGA = "Omega"; + constexpr const char* LM_GAMMA = "LM_gamma"; + constexpr const char* LM_RET = "LM_Re_t"; + constexpr const char* LM_GAMMA_SEP = "LM_gamma_sep"; + constexpr const char* LM_GAMMA_EFF = "LM_gamma_eff"; + + /*! + * \brief Get velocity field name for a given dimension index. + * \param[in] iDim - Dimension index (0=x, 1=y, 2=z) + * \return Field name string + */ + inline std::string GetVelocityName(unsigned short iDim) { + if (iDim == 0) return VELOCITY_X; + if (iDim == 1) return VELOCITY_Y; + return VELOCITY_Z; + } + + /*! + * \brief Get momentum field name for a given dimension index. + * \param[in] iDim - Dimension index (0=x, 1=y, 2=z) + * \return Field name string + */ + inline std::string GetMomentumName(unsigned short iDim) { + if (iDim == 0) return MOMENTUM_X; + if (iDim == 1) return MOMENTUM_Y; + return MOMENTUM_Z; + } + + /*! + * \brief Get grid velocity field name for a given dimension index. + * \param[in] iDim - Dimension index (0=x, 1=y, 2=z) + * \return Field name string + */ + inline std::string GetGridVelocityName(unsigned short iDim) { + if (iDim == 0) return GRID_VELOCITY_X; + if (iDim == 1) return GRID_VELOCITY_Y; + return GRID_VELOCITY_Z; + } +} + diff --git a/SU2_CFD/include/solvers/CSolver.hpp b/SU2_CFD/include/solvers/CSolver.hpp index 32b7bf43c08e..9f21a1b4d6c0 100644 --- a/SU2_CFD/include/solvers/CSolver.hpp +++ b/SU2_CFD/include/solvers/CSolver.hpp @@ -3420,6 +3420,13 @@ class CSolver { const CConfig *config, string val_filename); + /*! + * \brief Find the index of a field in the restart file fields vector. + * \param[in] fieldName - Name of the field to find (with or without quotes). + * \return Index of the field (0-based, excluding PointID), or -1 if not found. + */ + int FindFieldIndex(const string& fieldName) const; + /*! * \brief Read the metadata from a native SU2 restart file (ASCII or binary). * \param[in] geometry - Geometrical definition of the problem. diff --git a/SU2_CFD/src/solvers/CSolver.cpp b/SU2_CFD/src/solvers/CSolver.cpp index 99be4b04cc40..0b2f61b56118 100644 --- a/SU2_CFD/src/solvers/CSolver.cpp +++ b/SU2_CFD/src/solvers/CSolver.cpp @@ -2913,6 +2913,30 @@ void CSolver::Read_SU2_Restart_ASCII(CGeometry *geometry, const CConfig *config, } +int CSolver::FindFieldIndex(const string& fieldName) const { + /*--- Search for the field name in the fields vector. + Fields are stored with quotes, e.g., "Temperature", so we need to check + both with and without quotes. The first field is "PointID" which we skip. ---*/ + + string fieldNameWithQuotes = "\"" + fieldName + "\""; + string fieldNameNoQuotes = fieldName; + + /*--- Search starting from index 1 (skip PointID) ---*/ + for (size_t i = 1; i < fields.size(); i++) { + string field = fields[i]; + PrintingToolbox::trim(field); + + /*--- Check if field matches (with or without quotes) ---*/ + if (field == fieldNameWithQuotes || field == fieldNameNoQuotes) { + /*--- Return index excluding PointID (0-based) ---*/ + return static_cast(i - 1); + } + } + + /*--- Field not found ---*/ + return -1; +} + void CSolver::Read_SU2_Restart_Binary(CGeometry *geometry, const CConfig *config, string val_filename) { char str_buf[CGNS_STRING_SIZE], fname[100]; diff --git a/SU2_CFD/src/solvers/CTransLMSolver.cpp b/SU2_CFD/src/solvers/CTransLMSolver.cpp index 8e9b94d7f458..0e5ebd2f103d 100644 --- a/SU2_CFD/src/solvers/CTransLMSolver.cpp +++ b/SU2_CFD/src/solvers/CTransLMSolver.cpp @@ -31,6 +31,7 @@ #include "../../include/variables/CTurbSAVariable.hpp" #include "../../../Common/include/parallelization/omp_structure.hpp" #include "../../../Common/include/toolboxes/geometry_toolbox.hpp" +#include "../../include/solvers/CRestartFieldNames.hpp" /*--- This is the implementation of the Langtry-Menter transition model. The main reference for this model is:Langtry, Menter, AIAA J. 47(12) 2009 @@ -515,21 +516,16 @@ void CTransLMSolver::LoadRestart(CGeometry** geometry, CSolver*** solver, CConfi Read_SU2_Restart_ASCII(geometry[MESH_0], config, restart_filename); } - /*--- Skip flow variables ---*/ + /*--- Identify indices for LM transition variables ---*/ + long idx_gamma = solver[MESH_0][TRANS_SOL]->FindFieldIndex(RestartFieldNames::LM_GAMMA); + long idx_ret = solver[MESH_0][TRANS_SOL]->FindFieldIndex(RestartFieldNames::LM_RET); + long idx_sep = solver[MESH_0][TRANS_SOL]->FindFieldIndex(RestartFieldNames::LM_GAMMA_SEP); + long idx_eff = solver[MESH_0][TRANS_SOL]->FindFieldIndex(RestartFieldNames::LM_GAMMA_EFF); - unsigned short skipVars = nDim + solver[MESH_0][FLOW_SOL]->GetnVar() + solver[MESH_0][TURB_SOL] ->GetnVar(); - - /*--- Adjust the number of solution variables in the incompressible - restart. We always carry a space in nVar for the energy equation in the - mean flow solver, but we only write it to the restart if it is active. - Therefore, we must reduce skipVars here if energy is inactive so that - the turbulent variables are read correctly. ---*/ - - const bool incompressible = (config->GetKind_Regime() == ENUM_REGIME::INCOMPRESSIBLE); - const bool energy = config->GetEnergy_Equation(); - const bool weakly_coupled_heat = config->GetWeakly_Coupled_Heat(); - - if (incompressible && ((!energy) && (!weakly_coupled_heat))) skipVars--; + if (rank == MASTER_NODE) { + if (idx_gamma == -1) cout << "WARNING: Variable '" << RestartFieldNames::LM_GAMMA << "' not found. Initializing with default.\n"; + if (idx_ret == -1) cout << "WARNING: Variable '" << RestartFieldNames::LM_RET << "' not found. Initializing with default.\n"; + } /*--- Load data from the restart into correct containers. ---*/ @@ -541,13 +537,21 @@ void CTransLMSolver::LoadRestart(CGeometry** geometry, CSolver*** solver, CConfi const auto iPoint_Local = geometry[MESH_0]->GetGlobal_to_Local_Point(iPoint_Global); if (iPoint_Local > -1) { - /*--- We need to store this point's data, so jump to the correct - offset in the buffer of data from the restart file and load it. ---*/ + // Load Gamma + if (idx_gamma != -1) nodes->SetSolution(iPoint_Local, 0, Restart_Data[counter * Restart_Vars[1] + idx_gamma]); + else nodes->SetSolution(iPoint_Local, 0, Solution_Inf[0]); + + // Load Re_Theta_t + if (idx_ret != -1) nodes->SetSolution(iPoint_Local, 1, Restart_Data[counter * Restart_Vars[1] + idx_ret]); + else nodes->SetSolution(iPoint_Local, 1, Solution_Inf[1]); + + // Load Intermittency Sep + if (idx_sep != -1) nodes->SetIntermittencySep(iPoint_Local, Restart_Data[counter * Restart_Vars[1] + idx_sep]); + else nodes->SetIntermittencySep(iPoint_Local, 0.0); - const auto index = counter * Restart_Vars[1] + skipVars; - for (auto iVar = 0u; iVar < nVar; iVar++) nodes->SetSolution(iPoint_Local, iVar, Restart_Data[index + iVar]); - nodes->SetIntermittencySep(iPoint_Local, Restart_Data[index + 2]); - nodes->SetIntermittencyEff(iPoint_Local, Restart_Data[index + 3]); + // Load Intermittency Eff + if (idx_eff != -1) nodes->SetIntermittencyEff(iPoint_Local, Restart_Data[counter * Restart_Vars[1] + idx_eff]); + else nodes->SetIntermittencyEff(iPoint_Local, 0.0); /*--- Increment the overall counter for how many points have been loaded. ---*/ counter++; diff --git a/SU2_CFD/src/solvers/CTurbSolver.cpp b/SU2_CFD/src/solvers/CTurbSolver.cpp index 05f456883fb2..9f77c8aa48a1 100644 --- a/SU2_CFD/src/solvers/CTurbSolver.cpp +++ b/SU2_CFD/src/solvers/CTurbSolver.cpp @@ -29,6 +29,7 @@ #include "../../../Common/include/parallelization/omp_structure.hpp" #include "../../../Common/include/toolboxes/geometry_toolbox.hpp" #include "../../include/solvers/CScalarSolver.inl" +#include "../../include/solvers/CRestartFieldNames.hpp" /*--- Explicit instantiation of the parent class of CTurbSolver. ---*/ template class CScalarSolver; @@ -118,21 +119,55 @@ void CTurbSolver::LoadRestart(CGeometry** geometry, CSolver*** solver, CConfig* Read_SU2_Restart_ASCII(geometry[MESH_0], config, restart_filename); } - /*--- Skip flow variables ---*/ + /*--- Identify which turbulence variables to look for based on the model. ---*/ + struct TurbFieldInfo { + string name; + unsigned short model_var_index; + }; + vector target_fields; + + TURB_MODEL kind_turb = config->GetKind_Turb_Model(); + if (kind_turb == TURB_MODEL::SA) { + target_fields.push_back({RestartFieldNames::NU_TILDE, 0}); + } else if (kind_turb == TURB_MODEL::SST) { + target_fields.push_back({RestartFieldNames::TURB_KIN_ENERGY, 0}); + target_fields.push_back({RestartFieldNames::OMEGA, 1}); + } + + /*--- Fallback variables for legacy loading ---*/ + unsigned short skipVars = 0; + unsigned short nVar_Restart_Legacy = 0; - unsigned short skipVars = nDim + solver[MESH_0][FLOW_SOL]->GetnVar(); + if (target_fields.empty()) { + /*--- Skip flow variables ---*/ + skipVars = nDim + solver[MESH_0][FLOW_SOL]->GetnVar(); - /*--- Adjust the number of solution variables in the incompressible - restart. We always carry a space in nVar for the energy equation in the - mean flow solver, but we only write it to the restart if it is active. - Therefore, we must reduce skipVars here if energy is inactive so that - the turbulent variables are read correctly. ---*/ + const bool incompressible = (config->GetKind_Regime() == ENUM_REGIME::INCOMPRESSIBLE); + const bool energy = config->GetEnergy_Equation(); + const bool weakly_coupled_heat = config->GetWeakly_Coupled_Heat(); + const bool flamelet = (config->GetKind_FluidModel() == FLUID_FLAMELET); - const bool incompressible = (config->GetKind_Regime() == ENUM_REGIME::INCOMPRESSIBLE); - const bool energy = config->GetEnergy_Equation(); - const bool weakly_coupled_heat = config->GetWeakly_Coupled_Heat(); + if (incompressible && ((!energy) && (!weakly_coupled_heat) && (!flamelet))) { + skipVars--; + } - if (incompressible && ((!energy) && (!weakly_coupled_heat))) skipVars--; + /*--- Determine number of turbulence variables in restart file ---*/ + if (Restart_Vars[1] > skipVars) { + nVar_Restart_Legacy = min(static_cast(Restart_Vars[1] - skipVars), nVar); + } + } + + /*--- Find indices for named fields ---*/ + vector field_indices; + if (!target_fields.empty()) { + for (const auto& field : target_fields) { + long idx = solver[MESH_0][TURB_SOL]->FindFieldIndex(field.name); + field_indices.push_back(idx); + if (idx == -1 && rank == MASTER_NODE) { + cout << "WARNING: Turbulence variable '" << field.name << "' not found in restart file. Initializing with default.\n"; + } + } + } /*--- Load data from the restart into correct containers. ---*/ @@ -146,9 +181,31 @@ void CTurbSolver::LoadRestart(CGeometry** geometry, CSolver*** solver, CConfig* if (iPoint_Local > -1) { /*--- We need to store this point's data, so jump to the correct offset in the buffer of data from the restart file and load it. ---*/ - - const auto index = counter * Restart_Vars[1] + skipVars; - for (auto iVar = 0u; iVar < nVar; iVar++) nodes->SetSolution(iPoint_Local, iVar, Restart_Data[index + iVar]); + const auto base_idx = counter * Restart_Vars[1]; + + if (!target_fields.empty()) { + /*--- Named field loading ---*/ + for (size_t i = 0; i < target_fields.size(); ++i) { + long r_idx = field_indices[i]; + unsigned short v_idx = target_fields[i].model_var_index; + + if (r_idx != -1 && r_idx < Restart_Vars[1]) { + nodes->SetSolution(iPoint_Local, v_idx, Restart_Data[base_idx + r_idx]); + } else { + nodes->SetSolution(iPoint_Local, v_idx, Solution_Inf[v_idx]); + } + } + } else { + /*--- Positional loading (Legacy) ---*/ + const auto index = base_idx + skipVars; + for (auto iVar = 0u; iVar < nVar_Restart_Legacy; iVar++) { + nodes->SetSolution(iPoint_Local, iVar, Restart_Data[index + iVar]); + } + /*--- Initialize missing turbulence variables from defaults (legacy only) ---*/ + for (auto iVar = nVar_Restart_Legacy; iVar < nVar; iVar++) { + nodes->SetSolution(iPoint_Local, iVar, Solution_Inf[iVar]); + } + } /*--- Increment the overall counter for how many points have been loaded. ---*/ counter++; diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index d809b3f5ba56..26ff198c4cbd 100755 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -384,6 +384,24 @@ def main(): turb_flatplate.test_vals = [-4.297192, -6.731227, -0.187632, 0.057700] test_list.append(turb_flatplate) + test_list.append(turb_flatplate) + + # Euler precursor for restart test + euler_NACA0012_precursor = TestCase('euler_NACA0012_precursor') + euler_NACA0012_precursor.cfg_dir = "rans/naca0012" + euler_NACA0012_precursor.cfg_file = "euler_NACA0012_precursor.cfg" + euler_NACA0012_precursor.test_iter = 10 + euler_NACA0012_precursor.test_vals = [0.132688, 0.056260, -1.210312, 1.164664] # CL, CD, rho, rho*v + test_list.append(euler_NACA0012_precursor) + + # SA Restart from Euler (tests missing field initialization) + turb_SA_restart_test = TestCase('turb_SA_restart_test') + turb_SA_restart_test.cfg_dir = "rans/naca0012" + turb_SA_restart_test.cfg_file = "turb_SA_restart_test.cfg" + turb_SA_restart_test.test_iter = 10 + turb_SA_restart_test.test_vals = [0.169344, -0.030184, 0.051485, 1.474687] # CL, CD, rho, nu + test_list.append(turb_SA_restart_test) + # Flat plate (compressible) with species inlet turb_flatplate_species = TestCase('turb_flatplate_species') turb_flatplate_species.cfg_dir = "rans/flatplate" diff --git a/TestCases/rans/naca0012/euler_NACA0012_precursor.cfg b/TestCases/rans/naca0012/euler_NACA0012_precursor.cfg new file mode 100644 index 000000000000..c7dcdd6273fd --- /dev/null +++ b/TestCases/rans/naca0012/euler_NACA0012_precursor.cfg @@ -0,0 +1,54 @@ +SOLVER= EULER +MATH_PROBLEM= DIRECT +RESTART_SOL= NO + +MACH_NUMBER= 0.8 +AOA= 1.25 +FREESTREAM_PRESSURE= 101325.0 +FREESTREAM_TEMPERATURE= 273.15 +GAMMA_VALUE= 1.4 +GAS_CONSTANT= 287.87 + +REF_ORIGIN_MOMENT_X = 0.25 +REF_ORIGIN_MOMENT_Y = 0.00 +REF_ORIGIN_MOMENT_Z = 0.00 +REF_LENGTH= 1.0 +REF_AREA= 1.0 +REF_DIMENSIONALIZATION= DIMENSIONAL + +MARKER_EULER= ( airfoil ) +MARKER_FAR= ( farfield ) +MARKER_PLOTTING = ( airfoil ) +MARKER_MONITORING = ( airfoil ) + +NUM_METHOD_GRAD= WEIGHTED_LEAST_SQUARES +CFL_NUMBER= 10.0 +CFL_ADAPT= NO +ITER= 10 + +LINEAR_SOLVER= FGMRES +LINEAR_SOLVER_PREC= ILU +LINEAR_SOLVER_ERROR= 1E-10 +LINEAR_SOLVER_ITER= 10 + +MGLEVEL= 2 +MGCYCLE= V_CYCLE +MG_PRE_SMOOTH= ( 1, 1, 1, 1 ) +MG_POST_SMOOTH= ( 0, 0, 0, 0 ) +MG_CORRECTION_SMOOTH= ( 0, 0, 0, 0 ) +MG_DAMP_RESTRICTION= 1.0 +MG_DAMP_PROLONGATION= 1.0 + +CONV_NUM_METHOD_FLOW= JST +JST_SENSOR_COEFF= ( 0.5, 0.02 ) +TIME_DISCRE_FLOW= EULER_IMPLICIT + +MESH_FILENAME= ../../../QuickStart/mesh_NACA0012_inv.su2 +MESH_FORMAT= SU2 +MESH_OUT_FILENAME= mesh_out +TABULAR_FORMAT= CSV +CONV_FILENAME= history +RESTART_FILENAME= restart_flow_euler + +OUTPUT_WRT_FREQ= 10 +SCREEN_OUTPUT=(INNER_ITER, WALL_TIME, RMS_RES, LIFT, DRAG) diff --git a/TestCases/rans/naca0012/turb_SA_restart_test.cfg b/TestCases/rans/naca0012/turb_SA_restart_test.cfg new file mode 100644 index 000000000000..bed5ef79eb8d --- /dev/null +++ b/TestCases/rans/naca0012/turb_SA_restart_test.cfg @@ -0,0 +1,64 @@ +SOLVER= RANS +KIND_TURB_MODEL= SA +MATH_PROBLEM= DIRECT +RESTART_SOL= YES + +MACH_NUMBER= 0.8 +AOA= 1.25 +FREESTREAM_PRESSURE= 101325.0 +FREESTREAM_TEMPERATURE= 273.15 +REYNOLDS_NUMBER= 6.5e6 +REYNOLDS_LENGTH= 1.0 +GAMMA_VALUE= 1.4 +GAS_CONSTANT= 287.87 + +REF_ORIGIN_MOMENT_X = 0.25 +REF_ORIGIN_MOMENT_Y = 0.00 +REF_ORIGIN_MOMENT_Z = 0.00 +REF_LENGTH= 1.0 +REF_AREA= 1.0 +REF_DIMENSIONALIZATION= DIMENSIONAL + +MARKER_HEATFLUX= ( airfoil, 0.0 ) +MARKER_FAR= ( farfield ) +MARKER_PLOTTING = ( airfoil ) +MARKER_MONITORING = ( airfoil ) + +NUM_METHOD_GRAD= WEIGHTED_LEAST_SQUARES +CFL_NUMBER= 10.0 +CFL_ADAPT= NO +ITER= 10 + +LINEAR_SOLVER= FGMRES +LINEAR_SOLVER_PREC= ILU +LINEAR_SOLVER_ERROR= 1E-10 +LINEAR_SOLVER_ITER= 10 + +MGLEVEL= 2 +MGCYCLE= V_CYCLE +MG_PRE_SMOOTH= ( 1, 1, 1, 1 ) +MG_POST_SMOOTH= ( 0, 0, 0, 0 ) +MG_CORRECTION_SMOOTH= ( 0, 0, 0, 0 ) +MG_DAMP_RESTRICTION= 1.0 +MG_DAMP_PROLONGATION= 1.0 + +CONV_NUM_METHOD_FLOW= ROE +JST_SENSOR_COEFF= ( 0.5, 0.02 ) +TIME_DISCRE_FLOW= EULER_IMPLICIT + +CONV_NUM_METHOD_TURB= SCALAR_UPWIND +MUSCL_TURB= NO +SLOPE_LIMITER_TURB= VENKATAKRISHNAN +TIME_DISCRE_TURB= EULER_IMPLICIT + +MESH_FILENAME= ../../../QuickStart/mesh_NACA0012_inv.su2 +MESH_FORMAT= SU2 +MESH_OUT_FILENAME= mesh_out +TABULAR_FORMAT= CSV +CONV_FILENAME= history +RESTART_FILENAME= restart_flow_sa +RESTART_ITER= 10 +SOLUTION_FILENAME= restart_flow_euler.dat + +OUTPUT_WRT_FREQ= 10 +SCREEN_OUTPUT=(INNER_ITER, WALL_TIME, RMS_RES, LIFT, DRAG, RMS_NU_TILDE)