Fix solver initialization and multigrid logic bugs

- CTransLMSolver: Add EdgeFluxesDiff init, VelocityMag protection
- CTurbSSTSolver: Add wall_dist and FrictionVelocity protection
- CSpeciesSolver: Fix multigrid loop using wrong mesh level

Signed-off-by: shbhmexe <shubhushukla586@gmail.com>

Author: shbhmexe

SU2_CFD/src/solvers/CSpeciesSolver.cpp

@@ -275,9 +275,9 @@ void CSpeciesSolver::LoadRestart(CGeometry** geometry, CSolver*** solver, CConfi
                                            false);
 
     if (config->GetKind_Turb_Model() != TURB_MODEL::NONE)
-      solver[iMesh][TURB_SOL]->Postprocessing(geometry[MESH_0], solver[MESH_0], config, MESH_0);
+      solver[iMesh][TURB_SOL]->Postprocessing(geometry[iMesh], solver[iMesh], config, iMesh);
 
-    solver[iMesh][SPECIES_SOL]->Preprocessing(geometry[MESH_0], solver[MESH_0], config, MESH_0, NO_RK_ITER,
+    solver[iMesh][SPECIES_SOL]->Preprocessing(geometry[iMesh], solver[iMesh], config, iMesh, NO_RK_ITER,
                                               RUNTIME_SPECIES_SYS, false);
   }
 

SU2_CFD/src/solvers/CTransLMSolver.cpp

@@ -84,8 +84,10 @@ CTransLMSolver::CTransLMSolver(CGeometry *geometry, CConfig *config, unsigned sh
     LinSysRes.Initialize(nPoint, nPointDomain, nVar, 0.0);
     System.SetxIsZero(true);
 
-    if (ReducerStrategy)
+    if (ReducerStrategy) {
       EdgeFluxes.Initialize(geometry->GetnEdge(), geometry->GetnEdge(), nVar, nullptr);
+      EdgeFluxesDiff.Initialize(geometry->GetnEdge(), geometry->GetnEdge(), nVar, nullptr);
+    }
 
     /*--- Initialize the BGS residuals in multizone problems. ---*/
     if (multizone){
@@ -214,7 +216,7 @@ void CTransLMSolver::Postprocessing(CGeometry *geometry, CSolver **solver_contai
     const su2double vel_u = flowNodes->GetVelocity(iPoint, 0);
     const su2double vel_v = flowNodes->GetVelocity(iPoint, 1);
     const su2double vel_w = (nDim ==3) ? flowNodes->GetVelocity(iPoint, 2) : 0.0;
-    const su2double VelocityMag = sqrt(vel_u*vel_u + vel_v*vel_v + vel_w*vel_w);
+    const su2double VelocityMag = max(sqrt(vel_u*vel_u + vel_v*vel_v + vel_w*vel_w), 1e-10);
     su2double omega = 0.0;
     su2double k = 0.0;
     if(TurbFamily == TURB_FAMILY::KW){
@@ -253,8 +255,8 @@ void CTransLMSolver::Postprocessing(CGeometry *geometry, CSolver **solver_contai
     su2double Intermittency_Sep = 2.0*max(0.0, Re_v/(3.235*Corr_Rec)-1.0)*f_reattach;
     Intermittency_Sep = min(Intermittency_Sep,2.0)*f_theta;
     Intermittency_Sep = min(max(0.0, Intermittency_Sep), 2.0);
-    nodes -> SetIntermittencySep(iPoint, Intermittency_Sep);
-    nodes -> SetIntermittencyEff(iPoint, Intermittency_Sep);
+    nodes->SetIntermittencySep(iPoint, Intermittency_Sep);
+    nodes->SetIntermittencyEff(iPoint, Intermittency_Sep);
 
   }
   END_SU2_OMP_FOR
@@ -282,7 +284,6 @@ void CTransLMSolver::Source_Residual(CGeometry *geometry, CSolver **solver_conta
   const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT);
 
   auto* flowNodes = su2staticcast_p<CFlowVariable*>(solver_container[FLOW_SOL]->GetNodes());
-  //auto* turbNodes = su2staticcast_p<CFlowVariable*>(solver_container[TURB_SOL]->GetNodes());
   CVariable* turbNodes = solver_container[TURB_SOL]->GetNodes();
 
   /*--- Pick one numerics object per thread. ---*/
@@ -547,8 +548,8 @@ void CTransLMSolver::LoadRestart(CGeometry** geometry, CSolver*** solver, CConfi
 
         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]);
+        nodes->SetIntermittencySep(iPoint_Local,  Restart_Data[index + 2]);
+        nodes->SetIntermittencyEff(iPoint_Local,  Restart_Data[index + 3]);
 
         /*--- Increment the overall counter for how many points have been loaded. ---*/
         counter++;

SU2_CFD/src/solvers/CTurbSSTSolver.cpp

@@ -272,10 +272,10 @@ void CTurbSSTSolver::Postprocessing(CGeometry *geometry, CSolver **solver_contai
         }
         shearStress = sqrt(shearStress);
 
-        const su2double FrictionVelocity = sqrt(shearStress/flowNodes->GetDensity(iPoint));
+        const su2double FrictionVelocity = max(sqrt(shearStress/flowNodes->GetDensity(iPoint)), 1e-10);
         const su2double wall_dist = geometry->vertex[iMarker][iVertex]->GetNearestNeighborDistance();
 
-        const su2double Derivative = flowNodes->GetLaminarViscosity(jPoint) * pow(nodes->GetSolution(jPoint, 0), 0.673) / wall_dist;
+        const su2double Derivative = flowNodes->GetLaminarViscosity(jPoint) * pow(nodes->GetSolution(jPoint, 0), 0.673) / max(wall_dist, EPS);
         const su2double turbulence_index = 6.1 * Derivative / pow(FrictionVelocity, 2.346);
 
         nodes->SetTurbIndex(iPoint, turbulence_index);