Fix #2483: Add field-based restart loading with missing field initialization

guptapratykshh · guptapratykshh · commit 3bc9ca21582f · 2026-01-03T18:59:56.000+05:30
diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl
@@ -856,6 +856,47 @@ void CFVMFlowSolverBase<V, R>::LoadRestart_impl(CGeometry **geometry, CSolver **
       }
     }
 
+    /*--- Check for missing fields and prepare field-based lookup ---*/
+    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);
+    
+    /*--- Find field indices for solution variables ---*/
+    int fieldIdx_Pressure = -1, fieldIdx_Temperature = -1;
+    int fieldIdx_Velocity[MAXNDIM] = {-1, -1, -1};
+    int fieldIdx_Density = -1, fieldIdx_Energy = -1;
+    int fieldIdx_Momentum[MAXNDIM] = {-1, -1, -1};
+    
+    if (incompressible) {
+      fieldIdx_Pressure = FindFieldIndex("Pressure");
+      for (auto iDim = 0u; iDim < nDim; iDim++) {
+        string velName = "Velocity_" + string(iDim == 0 ? "x" : (iDim == 1 ? "y" : "z"));
+        fieldIdx_Velocity[iDim] = FindFieldIndex(velName);
+      }
+      if (energy || weakly_coupled_heat || flamelet) {
+        fieldIdx_Temperature = FindFieldIndex("Temperature");
+      }
+    } else {
+      /*--- Compressible flow ---*/
+      fieldIdx_Density = FindFieldIndex("Density");
+      for (auto iDim = 0u; iDim < nDim; iDim++) {
+        string momName = "Momentum_" + string(iDim == 0 ? "x" : (iDim == 1 ? "y" : "z"));
+        fieldIdx_Momentum[iDim] = FindFieldIndex(momName);
+      }
+      fieldIdx_Energy = FindFieldIndex("Energy");
+    }
+    
+    /*--- Warn about missing fields and prepare defaults ---*/
+    su2double default_Temperature = 0.0;
+    if (incompressible && (energy || weakly_coupled_heat || flamelet) && fieldIdx_Temperature < 0) {
+      if (rank == MASTER_NODE) {
+        cout << "\nWARNING: The restart file does not contain Temperature field." << endl;
+        cout << "Temperature will be initialized from config (INC_TEMPERATURE_INIT)." << endl;
+      }
+      default_Temperature = config->GetInc_Temperature_Init() / config->GetInc_Temperature_Ref();
+    }
+
     /*--- Load data from the restart into correct containers. ---*/
 
     unsigned long counter = 0;
@@ -868,19 +909,103 @@ void CFVMFlowSolverBase<V, R>::LoadRestart_impl(CGeometry **geometry, CSolver **
 
       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. ---*/
-
-        auto index = counter * Restart_Vars[1] + skipVars;
-
-        if (SolutionRestart == nullptr) {
-          for (auto iVar = 0u; iVar < nVar_Restart; iVar++)
-            nodes->SetSolution(iPoint_Local, iVar, Restart_Data[index+iVar]);
+        /*--- Base index for this point in the restart data ---*/
+        auto baseIndex = counter * Restart_Vars[1];
+        
+        /*--- Initialize solution buffer with defaults ---*/
+        su2double solutionBuffer[MAXNVAR] = {0.0};
+        
+        /*--- Load solution variables using field-based lookup ---*/
+        if (incompressible) {
+          /*--- Incompressible flow: Pressure, Velocity_x, Velocity_y, [Velocity_z], [Temperature] ---*/
+          
+          /*--- Load Pressure (index 0) ---*/
+          if (fieldIdx_Pressure >= 0) {
+            solutionBuffer[0] = Restart_Data[baseIndex + fieldIdx_Pressure];
+          } else {
+            /*--- Use default from config ---*/
+            solutionBuffer[0] = config->GetPressure_FreeStreamND();
+            if (rank == MASTER_NODE && counter == 0) {
+              cout << "WARNING: Pressure field not found in restart file, using config default." << endl;
+            }
+          }
+          
+          /*--- Load Velocity components (indices 1..nDim) ---*/
+          const su2double* Velocity_Inf = config->GetVelocity_FreeStreamND();
+          for (auto iDim = 0u; iDim < nDim; iDim++) {
+            if (fieldIdx_Velocity[iDim] >= 0) {
+              solutionBuffer[iDim + 1] = Restart_Data[baseIndex + fieldIdx_Velocity[iDim]];
+            } else {
+              /*--- Use default from config ---*/
+              solutionBuffer[iDim + 1] = Velocity_Inf[iDim];
+              if (rank == MASTER_NODE && counter == 0) {
+                cout << "WARNING: Velocity_" << (iDim == 0 ? "x" : (iDim == 1 ? "y" : "z")) 
+                     << " field not found in restart file, using config default." << endl;
+              }
+            }
+          }
+          
+          /*--- Load Temperature (index nDim+1) if energy equation is enabled ---*/
+          if (energy || weakly_coupled_heat || flamelet) {
+            if (fieldIdx_Temperature >= 0) {
+              solutionBuffer[nDim + 1] = Restart_Data[baseIndex + fieldIdx_Temperature];
+            } else {
+              /*--- Use default from config ---*/
+              solutionBuffer[nDim + 1] = default_Temperature;
+            }
+          }
+        } else {
+          /*--- Compressible flow: Density, Momentum_x, Momentum_y, [Momentum_z], Energy ---*/
+          
+          /*--- Load Density (index 0) ---*/
+          if (fieldIdx_Density >= 0) {
+            solutionBuffer[0] = Restart_Data[baseIndex + fieldIdx_Density];
+          } else {
+            solutionBuffer[0] = config->GetDensity_FreeStreamND();
+            if (rank == MASTER_NODE && counter == 0) {
+              cout << "WARNING: Density field not found in restart file, using config default." << endl;
+            }
+          }
+          
+          /*--- Load Momentum components (indices 1..nDim) ---*/
+          const su2double* Velocity_Inf_Comp = config->GetVelocity_FreeStreamND();
+          for (auto iDim = 0u; iDim < nDim; iDim++) {
+            if (fieldIdx_Momentum[iDim] >= 0) {
+              solutionBuffer[iDim + 1] = Restart_Data[baseIndex + fieldIdx_Momentum[iDim]];
+            } else {
+              /*--- Compute from density and velocity defaults ---*/
+              solutionBuffer[iDim + 1] = config->GetDensity_FreeStreamND() * Velocity_Inf_Comp[iDim];
+              if (rank == MASTER_NODE && counter == 0) {
+                cout << "WARNING: Momentum_" << (iDim == 0 ? "x" : (iDim == 1 ? "y" : "z")) 
+                     << " field not found in restart file, using config default." << endl;
+              }
+            }
+          }
+          
+          /*--- Load Energy (index nDim+1) ---*/
+          if (fieldIdx_Energy >= 0) {
+            solutionBuffer[nDim + 1] = Restart_Data[baseIndex + fieldIdx_Energy];
+          } else {
+            /*--- Compute from config defaults ---*/
+            su2double vel2 = 0.0;
+            for (auto iDim = 0u; iDim < nDim; iDim++) {
+              su2double vel = Velocity_Inf_Comp[iDim];
+              vel2 += vel * vel;
+            }
+            solutionBuffer[nDim + 1] = config->GetEnergy_FreeStreamND();
+            if (rank == MASTER_NODE && counter == 0) {
+              cout << "WARNING: Energy field not found in restart file, using config default." << 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];
+        
+        /*--- Set solution from buffer ---*/
+        if (SolutionRestart == nullptr) {
+          for (auto iVar = 0u; iVar < nVar; iVar++)
+            nodes->SetSolution(iPoint_Local, iVar, solutionBuffer[iVar]);
+        } else {
+          for (auto iVar = 0u; iVar < nVar; iVar++)
+            SolutionRestart[iVar] = solutionBuffer[iVar];
           nodes->SetSolution(iPoint_Local, SolutionRestart);
         }
 
@@ -894,14 +1019,22 @@ void CFVMFlowSolverBase<V, R>::LoadRestart_impl(CGeometry **geometry, CSolver **
           /*--- the grid velocities are set to 0. This is useful for FSI computations ---*/
 
           /*--- Rewind the index to retrieve the Coords. ---*/
-          index = counter * Restart_Vars[1];
+          auto index = baseIndex;
           const auto* Coord = &Restart_Data[index];
 
           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]; }
+            /*--- Load grid velocities using field-based lookup ---*/
+            int fieldIdx_GridVel[MAXNDIM] = {-1, -1, -1};
+            for (auto iDim = 0u; iDim < nDim; iDim++) {
+              string gridVelName = "Grid_Velocity_" + string(iDim == 0 ? "x" : (iDim == 1 ? "y" : "z"));
+              fieldIdx_GridVel[iDim] = FindFieldIndex(gridVelName);
+              if (fieldIdx_GridVel[iDim] >= 0) {
+                GridVel[iDim] = Restart_Data[baseIndex + fieldIdx_GridVel[iDim]];
+              } else {
+                GridVel[iDim] = 0.0;  /*--- Default to zero if not found ---*/
+              }
+            }
           }
 
           for (auto iDim = 0u; iDim < nDim; iDim++) {
@@ -915,7 +1048,7 @@ void CFVMFlowSolverBase<V, R>::LoadRestart_impl(CGeometry **geometry, CSolver **
 
         if (static_fsi && update_geo) {
         /*--- Rewind the index to retrieve the Coords. ---*/
-          index = counter*Restart_Vars[1];
+          auto index = baseIndex;
           const auto* Coord = &Restart_Data[index];
 
           for (auto iDim = 0u; iDim < nDim; iDim++) {
diff --git a/SU2_CFD/include/solvers/CSolver.hpp b/SU2_CFD/include/solvers/CSolver.hpp
@@ -3392,6 +3392,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
@@ -2907,6 +2907,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<int>(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/CTurbSolver.cpp b/SU2_CFD/src/solvers/CTurbSolver.cpp
@@ -135,6 +135,26 @@ void CTurbSolver::LoadRestart(CGeometry** geometry, CSolver*** solver, CConfig*
 
     if (incompressible && ((!energy) && (!weakly_coupled_heat) && (!flamelet))) skipVars--;
 
+    /*--- Determine number of turbulence variables in restart file ---*/
+    /*--- Total vars in restart = coordinates + flow vars + turbulence vars + grid velocities ---*/
+    unsigned short nVar_Restart = 0;
+    if (Restart_Vars[1] > skipVars) {
+      /*--- Estimate turbulence vars: total - skipVars, but cap at nVar ---*/
+      nVar_Restart = min(static_cast<unsigned short>(Restart_Vars[1] - skipVars), nVar);
+    }
+
+    /*--- Check if turbulence fields are present in restart file ---*/
+    /*--- Note: This is a simplified check. Full field-based lookup would require
+          knowing the specific turbulence model variable names (Nu_Tilde, Turb_Kin_Energy, etc.) ---*/
+    bool turb_fields_missing = false;
+    if (nVar_Restart < nVar) {
+      turb_fields_missing = true;
+      if (rank == MASTER_NODE) {
+        cout << "\nWARNING: The restart file appears to be missing turbulence variables." << endl;
+        cout << "Turbulence variables will be initialized from config defaults." << endl;
+      }
+    }
+
     /*--- Load data from the restart into correct containers. ---*/
 
     unsigned long counter = 0;
@@ -149,7 +169,18 @@ void CTurbSolver::LoadRestart(CGeometry** geometry, CSolver*** solver, CConfig*
          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]);
+        
+        /*--- Load available turbulence variables ---*/
+        for (auto iVar = 0u; iVar < nVar_Restart; iVar++) {
+          nodes->SetSolution(iPoint_Local, iVar, Restart_Data[index + iVar]);
+        }
+        
+        /*--- Initialize missing turbulence variables from defaults ---*/
+        if (turb_fields_missing) {
+          for (auto iVar = nVar_Restart; iVar < nVar; iVar++) {
+            nodes->SetSolution(iPoint_Local, iVar, Solution_Inf[iVar]);
+          }
+        }
 
         /*--- Increment the overall counter for how many points have been loaded. ---*/
         counter++;
