Merge pull request #1450 from su2code/clean_loadrestart

Clean LoadRestart for Flow/Turb/Heat

Author: TobiKattmann

Common/include/CConfig.hpp

@@ -554,6 +554,7 @@ class CConfig {
   STRUCT_DEFORMATION Kind_Struct_Solver;   /*!< \brief Determines the geometric condition (small or large deformations) for structural analysis. */
   unsigned short Kind_DV_FEA;              /*!< \brief Kind of Design Variable for FEA problems.*/
 
+  unsigned short nTurbVar;          /*!< \brief Number of Turbulence variables, i.e. 1 for SA-types, 2 for SST. */
   TURB_MODEL Kind_Turb_Model;       /*!< \brief Turbulent model definition. */
   unsigned short Kind_SGS_Model;    /*!< \brief LES SGS model definition. */
   unsigned short Kind_Trans_Model,  /*!< \brief Transition model definition. */
@@ -4179,6 +4180,12 @@ class CConfig {
    */
   void SetKind_SU2(SU2_COMPONENT val_kind_su2) { Kind_SU2 = val_kind_su2 ; }
 
+  /*!
+   * \brief Get the number of Turbulence Variables.
+   * \return Number of Turbulence Variables.
+   */
+  unsigned short GetnTurbVar(void) const { return nTurbVar; }
+
   /*!
    * \brief Get the kind of the turbulence model.
    * \return Kind of the turbulence model.

Common/src/CConfig.cpp

@@ -5075,6 +5075,17 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i
   if (GetGasModel() == "ARGON") {monoatomic = true;}
   else {monoatomic = false;}
 
+  /*--- Set number of Turbulence Variables. ---*/
+  switch(Kind_Turb_Model) {
+    case TURB_MODEL::NONE:
+      nTurbVar = 0; break;
+    case TURB_MODEL::SA: case TURB_MODEL::SA_COMP: case TURB_MODEL::SA_E_COMP: case TURB_MODEL::SA_E:
+    case TURB_MODEL::SA_NEG:
+      nTurbVar = 1; break;
+    case TURB_MODEL::SST: case TURB_MODEL::SST_SUST:
+      nTurbVar = 2; break;
+  }
+
   // This option is deprecated. After a grace period until 7.2.0 the usage warning should become an error.
   if(OptionIsSet("CONV_CRITERIA") && rank == MASTER_NODE) {
     cout << "\n\nWARNING: CONV_CRITERIA is deprecated. SU2 will choose the criteria automatically based on the CONV_FIELD.\n"

SU2_CFD/include/solvers/CFVMFlowSolverBase.inl

@@ -718,139 +718,119 @@ void CFVMFlowSolverBase<V, R>::SetUniformInlet(const CConfig* config, unsigned s
 }
 
 template <class V, ENUM_REGIME R>
-void CFVMFlowSolverBase<V, R>::LoadRestart_impl(CGeometry **geometry, CSolver ***solver, CConfig *config,
-                                                int iter, bool update_geo, su2double* SolutionRestart,
+void CFVMFlowSolverBase<V, R>::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 ---*/
 
-  unsigned short iDim, iVar, iMesh;
-  unsigned long iPoint, index, iChildren, Point_Fine;
-  TURB_MODEL turb_model = config->GetKind_Turb_Model();
-  su2double Area_Children, Area_Parent;
-  const su2double* Solution_Fine = nullptr;
-  const passivedouble* Coord = nullptr;
-  bool dual_time = ((config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_1ST) ||
-                    (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_2ND));
-  bool static_fsi = ((config->GetTime_Marching() == TIME_MARCHING::STEADY) && config->GetFSI_Simulation());
-  bool steady_restart = config->GetSteadyRestart();
-  bool turbulent = (config->GetKind_Turb_Model() != TURB_MODEL::NONE);
-
-  string restart_filename = config->GetFilename(config->GetSolution_FileName(), "", iter);
+  const string restart_filename = config->GetFilename(config->GetSolution_FileName(), "", iter);
+  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. ---*/
   SU2_OMP_MASTER {
 
-  if (nVar_Restart == 0) nVar_Restart = nVar;
+    if (nVar_Restart == 0) nVar_Restart = nVar;
 
-  /*--- Skip coordinates ---*/
+    /*--- Skip coordinates ---*/
 
-  unsigned short skipVars = geometry[MESH_0]->GetnDim();
+    unsigned short skipVars = nDim;
 
-  /*--- Store the number of variables for the turbulence model
-   (that could appear in the restart file before the grid velocities). ---*/
-  unsigned short turbVars = 0;
-  if (turbulent){
-    if ((turb_model == TURB_MODEL::SST) || (turb_model == TURB_MODEL::SST_SUST)) turbVars = 2;
-    else turbVars = 1;
-  }
+    /*--- Read the restart data from either an ASCII or binary SU2 file. ---*/
 
-  /*--- Read the restart data from either an ASCII or binary SU2 file. ---*/
-
-  if (config->GetRead_Binary_Restart()) {
-    Read_SU2_Restart_Binary(geometry[MESH_0], config, restart_filename);
-  } else {
-    Read_SU2_Restart_ASCII(geometry[MESH_0], config, restart_filename);
-  }
+    if (config->GetRead_Binary_Restart()) {
+      Read_SU2_Restart_Binary(geometry[MESH_0], config, restart_filename);
+    } else {
+      Read_SU2_Restart_ASCII(geometry[MESH_0], config, restart_filename);
+    }
 
-  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;
+    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 data from the restart into correct containers. ---*/
+    /*--- Load data from the restart into correct containers. ---*/
 
-  unsigned long counter = 0, iPoint_Global = 0;
-  for (; iPoint_Global < geometry[MESH_0]->GetGlobal_nPointDomain(); iPoint_Global++) {
+    unsigned long counter = 0;
+    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. ---*/
+      /*--- Retrieve local index. If this node from the restart file lives
+      on the current processor, we will load and instantiate the vars. ---*/
 
-    auto iPoint_Local = geometry[MESH_0]->GetGlobal_to_Local_Point(iPoint_Global);
+      const auto iPoint_Local = geometry[MESH_0]->GetGlobal_to_Local_Point(iPoint_Global);
 
-    if (iPoint_Local > -1) {
+      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. ---*/
+        /*--- 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. ---*/
 
-      index = counter*Restart_Vars[1] + skipVars;
+        auto index = counter * Restart_Vars[1] + skipVars;
 
-      if (SolutionRestart == nullptr) {
-        for (iVar = 0; iVar < nVar_Restart; iVar++)
-          nodes->SetSolution(iPoint_Local, iVar, Restart_Data[index+iVar]);
-      }
-      else {
-        /*--- Used as buffer, allows defaults for nVar > nVar_Restart. ---*/
-        for (iVar = 0; iVar < nVar_Restart; iVar++)
-          SolutionRestart[iVar] = Restart_Data[index+iVar];
-        nodes->SetSolution(iPoint_Local, SolutionRestart);
-      }
+        if (SolutionRestart == nullptr) {
+          for (auto iVar = 0u; iVar < nVar_Restart; iVar++)
+            nodes->SetSolution(iPoint_Local, iVar, Restart_Data[index+iVar]);
+        }
+        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);
+        }
 
-      /*--- For dynamic meshes, read in and store the
-       grid coordinates and grid velocities for each node. ---*/
+        /*--- For dynamic meshes, read in and store the
+        grid coordinates and grid velocities for each node. ---*/
 
-      if (dynamic_grid && update_geo) {
+        if (dynamic_grid && update_geo) {
 
-        /*--- 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 ---*/
+          /*--- 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 ---*/
 
-        /*--- Rewind the index to retrieve the Coords. ---*/
-        index = counter*Restart_Vars[1];
-        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 + turbVars;
-          for (iDim = 0; iDim < nDim; iDim++) { GridVel[iDim] = Restart_Data[index+iDim]; }
-        }
+          /*--- Rewind the index to retrieve the Coords. ---*/
+          index = counter * Restart_Vars[1];
+          const auto* Coord = &Restart_Data[index];
 
-        for (iDim = 0; iDim < nDim; iDim++) {
-          geometry[MESH_0]->nodes->SetCoord(iPoint_Local, iDim, Coord[iDim]);
-          geometry[MESH_0]->nodes->SetGridVel(iPoint_Local, iDim, GridVel[iDim]);
+          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]; }
+          }
+
+          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]);
+          }
         }
-      }
 
-      /*--- 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 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). ---*/
 
-      if (static_fsi && update_geo) {
-       /*--- Rewind the index to retrieve the Coords. ---*/
-        index = counter*Restart_Vars[1];
-        Coord = &Restart_Data[index];
+        if (static_fsi && update_geo) {
+        /*--- Rewind the index to retrieve the Coords. ---*/
+          index = counter*Restart_Vars[1];
+          const auto* Coord = &Restart_Data[index];
 
-        for (iDim = 0; 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]);
+          }
         }
-      }
 
-      /*--- Increment the overall counter for how many points have been loaded. ---*/
-      counter++;
+        /*--- Increment the overall counter for how many points have been loaded. ---*/
+        counter++;
+      }
     }
 
-  }
-
-  /*--- Detect a wrong solution file ---*/
+    /*--- Detect a wrong solution file ---*/
 
-  if (counter != nPointDomain) {
-    SU2_MPI::Error(string("The solution file ") + restart_filename + string(" doesn't match with the mesh file!\n") +
-                   string("It could be empty lines at the end of the file."), CURRENT_FUNCTION);
-  }
+    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);
+    }
   }
   END_SU2_OMP_MASTER
   SU2_OMP_BARRIER
@@ -862,10 +842,10 @@ void CFVMFlowSolverBase<V, R>::LoadRestart_impl(CGeometry **geometry, CSolver **
     CGeometry::UpdateGeometry(geometry, config);
 
     if (dynamic_grid) {
-      for (iMesh = 0; iMesh <= config->GetnMGLevels(); iMesh++) {
+      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, GRID_VELOCITY);
           geometry[MESH_0]->CompleteComms(geometry[MESH_0], config, GRID_VELOCITY);
@@ -884,23 +864,26 @@ void CFVMFlowSolverBase<V, R>::LoadRestart_impl(CGeometry **geometry, CSolver **
 
   /*--- For turbulent simulations the flow preprocessing is done by the turbulence solver
    *    after it loads its variables (they are needed to compute flow primitives). ---*/
-  if (!turbulent) {
+  if (config->GetKind_Turb_Model() == TURB_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 (iMesh = 1; iMesh <= config->GetnMGLevels(); iMesh++) {
+  for (auto iMesh = 1u; iMesh <= config->GetnMGLevels(); iMesh++) {
+
     SU2_OMP_FOR_STAT(omp_chunk_size)
-    for (iPoint = 0; iPoint < geometry[iMesh]->GetnPoint(); iPoint++) {
-      Area_Parent = geometry[iMesh]->nodes->GetVolume(iPoint);
+    for (auto iPoint = 0ul; iPoint < geometry[iMesh]->GetnPoint(); iPoint++) {
+      const su2double Area_Parent = geometry[iMesh]->nodes->GetVolume(iPoint);
       su2double Solution_Coarse[MAXNVAR] = {0.0};
-      for (iChildren = 0; iChildren < geometry[iMesh]->nodes->GetnChildren_CV(iPoint); iChildren++) {
-        Point_Fine = geometry[iMesh]->nodes->GetChildren_CV(iPoint, iChildren);
-        Area_Children = geometry[iMesh-1]->nodes->GetVolume(Point_Fine);
-        Solution_Fine = solver[iMesh-1][FLOW_SOL]->GetNodes()->GetSolution(Point_Fine);
-        for (iVar = 0; iVar < nVar; iVar++) {
-          Solution_Coarse[iVar] += Solution_Fine[iVar]*Area_Children/Area_Parent;
+
+      for (auto iChildren = 0ul; iChildren < geometry[iMesh]->nodes->GetnChildren_CV(iPoint); iChildren++) {
+        const auto Point_Fine = geometry[iMesh]->nodes->GetChildren_CV(iPoint, iChildren);
+        const su2double Area_Children = geometry[iMesh - 1]->nodes->GetVolume(Point_Fine);
+        const su2double* Solution_Fine = solver[iMesh - 1][FLOW_SOL]->GetNodes()->GetSolution(Point_Fine);
+
+        for (auto iVar = 0u; iVar < nVar; iVar++) {
+          Solution_Coarse[iVar] += Solution_Fine[iVar] * Area_Children / Area_Parent;
         }
       }
       solver[iMesh][FLOW_SOL]->GetNodes()->SetSolution(iPoint,Solution_Coarse);
@@ -910,12 +893,14 @@ void CFVMFlowSolverBase<V, R>::LoadRestart_impl(CGeometry **geometry, CSolver **
     solver[iMesh][FLOW_SOL]->InitiateComms(geometry[iMesh], config, SOLUTION);
     solver[iMesh][FLOW_SOL]->CompleteComms(geometry[iMesh], config, SOLUTION);
 
-    if (!turbulent) {
+    if (config->GetKind_Turb_Model() == TURB_MODEL::NONE) {
       solver[iMesh][FLOW_SOL]->Preprocessing(geometry[iMesh], solver[iMesh], config, iMesh, NO_RK_ITER, RUNTIME_FLOW_SYS, false);
     }
   }
 
   /*--- Update the old geometry (coordinates n and n-1) in dual time-stepping strategy. ---*/
+  const bool dual_time = ((config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_1ST) ||
+                          (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_2ND));
   if (dual_time && config->GetGrid_Movement() && !config->GetDeform_Mesh() &&
       (config->GetKind_GridMovement() != RIGID_MOTION)) {
     Restart_OldGeometry(geometry[MESH_0], config);
@@ -926,13 +911,13 @@ void CFVMFlowSolverBase<V, R>::LoadRestart_impl(CGeometry **geometry, CSolver **
   {
   /*--- Delete the class memory that is used to load the restart. ---*/
 
-  delete [] Restart_Vars; Restart_Vars = nullptr;
-  delete [] Restart_Data; Restart_Data = nullptr;
-
+    delete [] Restart_Vars;
+    Restart_Vars = nullptr;
+    delete [] Restart_Data;
+    Restart_Data = nullptr;
   }
   END_SU2_OMP_MASTER
   SU2_OMP_BARRIER
-
 }
 
 template <class V, ENUM_REGIME R>