From 8e561ee771c97c9a6601f52e4ae23c5248aff4b7 Mon Sep 17 00:00:00 2001
From: Cristopher Morales
 <98025159+Cristopher-Morales@users.noreply.github.com>
Date: Mon, 16 Sep 2024 22:33:30 +0200
Subject: [PATCH] Extending CFL adapt to include species transport (#2298)

* extending CFL adapt to include species transport

* add underrelaxation factor for the last species that is not computed during simulations

* removing underrelaxation factor for species transport

* adding CFLSpeciesReduction factor
---
 SU2_CFD/src/solvers/CSolver.cpp | 16 ++++++++++++++--
 1 file changed, 14 insertions(+), 2 deletions(-)

diff --git a/SU2_CFD/src/solvers/CSolver.cpp b/SU2_CFD/src/solvers/CSolver.cpp
index bfbd4edad58..f8327c868bd 100644
--- a/SU2_CFD/src/solvers/CSolver.cpp
+++ b/SU2_CFD/src/solvers/CSolver.cpp
@@ -1730,6 +1730,7 @@ void CSolver::AdaptCFLNumber(CGeometry **geometry,
 
     CSolver *solverFlow = solver_container[iMesh][FLOW_SOL];
     CSolver *solverTurb = solver_container[iMesh][TURB_SOL];
+    CSolver *solverSpecies = solver_container[iMesh][SPECIES_SOL];
 
     /* Compute the reduction factor for CFLs on the coarse levels. */
 
@@ -1744,10 +1745,12 @@ void CSolver::AdaptCFLNumber(CGeometry **geometry,
      solver residual within the specified number of linear iterations. */
 
     su2double linResTurb = 0.0;
+    su2double linResSpecies = 0.0;
     if ((iMesh == MESH_0) && solverTurb) linResTurb = solverTurb->GetResLinSolver();
+    if ((iMesh == MESH_0) && solverSpecies) linResSpecies = solverSpecies->GetResLinSolver();
 
-    /* Max linear residual between flow and turbulence. */
-    const su2double linRes = max(solverFlow->GetResLinSolver(), linResTurb);
+    /* Max linear residual between flow and turbulence/species transport. */
+    const su2double linRes = max(solverFlow->GetResLinSolver(), max(linResTurb, linResSpecies));
 
     /* Tolerance limited to an acceptable value. */
     const su2double linTol = max(acceptableLinTol, config->GetLinear_Solver_Error());
@@ -1779,6 +1782,11 @@ void CSolver::AdaptCFLNumber(CGeometry **geometry,
           New_Func += log10(solverTurb->GetRes_RMS(iVar));
         }
       }
+      if ((iMesh == MESH_0) && solverSpecies) {
+        for (unsigned short iVar = 0; iVar < solverSpecies->GetnVar(); iVar++) {
+          New_Func += log10(solverSpecies->GetRes_RMS(iVar));
+        }
+      }
 
       /* Compute the difference in the nonlinear residuals between the
        current and previous iterations, taking care with very low initial
@@ -1822,6 +1830,7 @@ void CSolver::AdaptCFLNumber(CGeometry **geometry,
 
     su2double myCFLMin = 1e30, myCFLMax = 0.0, myCFLSum = 0.0;
     const su2double CFLTurbReduction = config->GetCFLRedCoeff_Turb();
+    const su2double CFLSpeciesReduction = config->GetCFLRedCoeff_Species();
 
     SU2_OMP_MASTER
     if ((iMesh == MESH_0) && fullComms) {
@@ -1888,6 +1897,9 @@ void CSolver::AdaptCFLNumber(CGeometry **geometry,
       if ((iMesh == MESH_0) && solverTurb) {
         solverTurb->GetNodes()->SetLocalCFL(iPoint, CFL * CFLTurbReduction);
       }
+      if ((iMesh == MESH_0) && solverSpecies) {
+        solverSpecies->GetNodes()->SetLocalCFL(iPoint, CFL * CFLSpeciesReduction);
+      }
 
       /* Store min and max CFL for reporting on the fine grid. */