From f517cb722f9654524b734dedb457282be848574e Mon Sep 17 00:00:00 2001 From: shbhmexe Date: Sat, 27 Dec 2025 13:48:41 +0530 Subject: [PATCH 1/6] 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 --- SU2_CFD/src/solvers/CSpeciesSolver.cpp | 4 ++-- SU2_CFD/src/solvers/CTransLMSolver.cpp | 15 ++++++++------- SU2_CFD/src/solvers/CTurbSSTSolver.cpp | 4 ++-- 3 files changed, 12 insertions(+), 11 deletions(-) diff --git a/SU2_CFD/src/solvers/CSpeciesSolver.cpp b/SU2_CFD/src/solvers/CSpeciesSolver.cpp index 3f25af7d1183..d8def3dcca34 100644 --- a/SU2_CFD/src/solvers/CSpeciesSolver.cpp +++ b/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); } diff --git a/SU2_CFD/src/solvers/CTransLMSolver.cpp b/SU2_CFD/src/solvers/CTransLMSolver.cpp index e8a6dcf389d0..6abe3b0d208b 100644 --- a/SU2_CFD/src/solvers/CTransLMSolver.cpp +++ b/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(solver_container[FLOW_SOL]->GetNodes()); - //auto* turbNodes = su2staticcast_p(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++; diff --git a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp index b36f5c2b8765..9f12ef76651e 100644 --- a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp +++ b/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); From 3e8292c2de06dbd59cc1cf307afc3d2fc43ad495 Mon Sep 17 00:00:00 2001 From: Shubham shukla Date: Sat, 27 Dec 2025 23:54:07 +0530 Subject: [PATCH 2/6] Update SU2_CFD/src/solvers/CTurbSSTSolver.cpp Co-authored-by: Pedro Gomes <38071223+pcarruscag@users.noreply.github.com> --- SU2_CFD/src/solvers/CTurbSSTSolver.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp index 9f12ef76651e..0dcf4c3cd2d5 100644 --- a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp +++ b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp @@ -272,7 +272,7 @@ void CTurbSSTSolver::Postprocessing(CGeometry *geometry, CSolver **solver_contai } shearStress = sqrt(shearStress); - const su2double FrictionVelocity = max(sqrt(shearStress/flowNodes->GetDensity(iPoint)), 1e-10); + const su2double FrictionVelocity = max(sqrt(shearStress/flowNodes->GetDensity(iPoint)), EPS); const su2double wall_dist = geometry->vertex[iMarker][iVertex]->GetNearestNeighborDistance(); const su2double Derivative = flowNodes->GetLaminarViscosity(jPoint) * pow(nodes->GetSolution(jPoint, 0), 0.673) / max(wall_dist, EPS); From aef9af985662e2d776689884658f8de22360dc80 Mon Sep 17 00:00:00 2001 From: Shubham shukla Date: Sat, 27 Dec 2025 23:54:15 +0530 Subject: [PATCH 3/6] Update SU2_CFD/src/solvers/CTurbSSTSolver.cpp Co-authored-by: Pedro Gomes <38071223+pcarruscag@users.noreply.github.com> --- SU2_CFD/src/solvers/CTurbSSTSolver.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp index 0dcf4c3cd2d5..e5ea6a6634e6 100644 --- a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp +++ b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp @@ -275,7 +275,7 @@ void CTurbSSTSolver::Postprocessing(CGeometry *geometry, CSolver **solver_contai const su2double FrictionVelocity = max(sqrt(shearStress/flowNodes->GetDensity(iPoint)), EPS); const su2double wall_dist = geometry->vertex[iMarker][iVertex]->GetNearestNeighborDistance(); - const su2double Derivative = flowNodes->GetLaminarViscosity(jPoint) * pow(nodes->GetSolution(jPoint, 0), 0.673) / max(wall_dist, EPS); + const su2double Derivative = flowNodes->GetLaminarViscosity(jPoint) * pow(nodes->GetSolution(jPoint, 0), 0.673) / wall_dist; const su2double turbulence_index = 6.1 * Derivative / pow(FrictionVelocity, 2.346); nodes->SetTurbIndex(iPoint, turbulence_index); From ca9a6bc00bf34f4c1bab388c1c6dc12e46da0e73 Mon Sep 17 00:00:00 2001 From: Shubham shukla Date: Sat, 27 Dec 2025 23:56:14 +0530 Subject: [PATCH 4/6] Update SU2_CFD/src/solvers/CTransLMSolver.cpp Co-authored-by: Pedro Gomes <38071223+pcarruscag@users.noreply.github.com> --- SU2_CFD/src/solvers/CTransLMSolver.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SU2_CFD/src/solvers/CTransLMSolver.cpp b/SU2_CFD/src/solvers/CTransLMSolver.cpp index 6abe3b0d208b..0710ec86bfc8 100644 --- a/SU2_CFD/src/solvers/CTransLMSolver.cpp +++ b/SU2_CFD/src/solvers/CTransLMSolver.cpp @@ -216,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 = max(sqrt(vel_u*vel_u + vel_v*vel_v + vel_w*vel_w), 1e-10); + const su2double VelocityMag = max(sqrt(pow(vel_u, 2) + pow(vel_v, 2) + pow(vel_w, 2)), EPS); su2double omega = 0.0; su2double k = 0.0; if(TurbFamily == TURB_FAMILY::KW){ From 55b7bab584f2474760c4ad8ef74cbfae1f197930 Mon Sep 17 00:00:00 2001 From: shbhmexe Date: Sun, 28 Dec 2025 00:11:05 +0530 Subject: [PATCH 5/6] Update CTransLMSolver.cpp --- SU2_CFD/src/solvers/CTransLMSolver.cpp | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/SU2_CFD/src/solvers/CTransLMSolver.cpp b/SU2_CFD/src/solvers/CTransLMSolver.cpp index 0710ec86bfc8..8e9b94d7f458 100644 --- a/SU2_CFD/src/solvers/CTransLMSolver.cpp +++ b/SU2_CFD/src/solvers/CTransLMSolver.cpp @@ -84,10 +84,8 @@ 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){ From 8143f4f6221d67cb11e63a604ffd57ca91b5d130 Mon Sep 17 00:00:00 2001 From: shbhmexe Date: Sun, 28 Dec 2025 11:50:46 +0530 Subject: [PATCH 6/6] Update CSpeciesSolver.cpp --- SU2_CFD/src/solvers/CSpeciesSolver.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/SU2_CFD/src/solvers/CSpeciesSolver.cpp b/SU2_CFD/src/solvers/CSpeciesSolver.cpp index d8def3dcca34..3f25af7d1183 100644 --- a/SU2_CFD/src/solvers/CSpeciesSolver.cpp +++ b/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[iMesh], solver[iMesh], config, iMesh); + solver[iMesh][TURB_SOL]->Postprocessing(geometry[MESH_0], solver[MESH_0], config, MESH_0); - solver[iMesh][SPECIES_SOL]->Preprocessing(geometry[iMesh], solver[iMesh], config, iMesh, NO_RK_ITER, + solver[iMesh][SPECIES_SOL]->Preprocessing(geometry[MESH_0], solver[MESH_0], config, MESH_0, NO_RK_ITER, RUNTIME_SPECIES_SYS, false); }