From 52be6cade66f4ebca2bcef4067c368b07b447103 Mon Sep 17 00:00:00 2001 From: Josh Kelly Date: Fri, 17 Apr 2026 15:04:13 +0200 Subject: [PATCH 1/4] Added blended back-flow for thermodynamically consistent flamelets --- Common/include/CConfig.hpp | 7 + Common/src/CConfig.cpp | 2 + .../solvers/CSpeciesFlameletSolver.hpp | 16 ++ SU2_CFD/include/solvers/CSpeciesSolver.hpp | 2 +- SU2_CFD/src/solvers/CIncEulerSolver.cpp | 46 ++++++ .../src/solvers/CSpeciesFlameletSolver.cpp | 137 ++++++++++++++++++ 6 files changed, 209 insertions(+), 1 deletion(-) diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp index bbb2cf12bf8d..f454199e246b 100644 --- a/Common/include/CConfig.hpp +++ b/Common/include/CConfig.hpp @@ -635,6 +635,7 @@ class CConfig { unsigned short nInc_Outlet; /*!< \brief Number of inlet boundary treatment types listed. */ su2double Inc_Inlet_Damping; /*!< \brief Damping factor applied to the iterative updates to the velocity at a pressure inlet in incompressible flow. */ su2double Inc_Outlet_Damping; /*!< \brief Damping factor applied to the iterative updates to the pressure at a mass flow outlet in incompressible flow. */ + bool Inc_Outlet_BackflowPrevention; /*!< \brief Enable removal of the reversed normal velocity component at outlet faces where backflow is detected. */ bool InletUseNormal; /*!< \brief Flag for whether to use the local normal as the flow direction for a pressure inlet. */ su2double Linear_Solver_Error; /*!< \brief Min error of the linear solver for the implicit formulation. */ su2double Deform_Linear_Solver_Error; /*!< \brief Min error of the linear solver for the implicit formulation. */ @@ -5191,6 +5192,12 @@ class CConfig { */ su2double GetInc_Outlet_Damping(void) const { return Inc_Outlet_Damping; } + /*! + * \brief Check whether backflow prevention is enabled at incompressible outlets. + * \return True if the reversed normal velocity component should be zeroed on backflow faces. + */ + bool GetInc_Outlet_BackflowPrevention(void) const { return Inc_Outlet_BackflowPrevention; } + /*! * \brief Get the kind of mixing process for averaging quantities at the boundaries. * \return Kind of mixing process. diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index 2ba453b197b0..8b0886fceb6d 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -1396,6 +1396,8 @@ void CConfig::SetConfig_Options() { addDoubleOption("INC_INLET_DAMPING", Inc_Inlet_Damping, 0.1); /*!\brief INC_OUTLET_DAMPING \n DESCRIPTION: Damping factor applied to the iterative updates to the pressure at a mass flow outlet in incompressible flow (0.1 by default). \ingroup Config*/ addDoubleOption("INC_OUTLET_DAMPING", Inc_Outlet_Damping, 0.1); + /*!\brief INC_OUTLET_BACKFLOW_PREVENTION \n DESCRIPTION: Remove the reversed normal velocity component from the outlet ghost state when backflow is detected, preventing spurious inflow through the outlet. \ingroup Config*/ + addBoolOption("INC_OUTLET_BACKFLOW_PREVENTION", Inc_Outlet_BackflowPrevention, false); /*--- Options related to the species solver. ---*/ diff --git a/SU2_CFD/include/solvers/CSpeciesFlameletSolver.hpp b/SU2_CFD/include/solvers/CSpeciesFlameletSolver.hpp index df61813a3e27..9b26db0e5c29 100644 --- a/SU2_CFD/include/solvers/CSpeciesFlameletSolver.hpp +++ b/SU2_CFD/include/solvers/CSpeciesFlameletSolver.hpp @@ -173,6 +173,22 @@ class CSpeciesFlameletSolver final : public CSpeciesSolver { void BC_Inlet(CGeometry* geometry, CSolver** solver_container, CNumerics* conv_numerics, CNumerics* visc_numerics, CConfig* config, unsigned short val_marker) override; + /*! + * \brief Impose the outlet boundary condition. + * When backflow prevention is enabled and reversed flow is detected at a face, + * the scalar ghost state is set to the inlet prescribed values (fresh mixture) + * instead of the Neumann interior values, preventing burned-gas scalars from + * re-entering the domain and pushing FGM look-ups outside the manifold. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] solver_container - Container vector with all the solutions. + * \param[in] conv_numerics - Description of the numerical method. + * \param[in] visc_numerics - Description of the numerical method. + * \param[in] config - Definition of the particular problem. + * \param[in] val_marker - Surface marker where the boundary condition is applied. + */ + void BC_Outlet(CGeometry* geometry, CSolver** solver_container, CNumerics* conv_numerics, CNumerics* visc_numerics, + CConfig* config, unsigned short val_marker) override; + /*! * \brief Impose the (received) conjugate heat variables. * \param[in] geometry - Geometrical definition of the problem. diff --git a/SU2_CFD/include/solvers/CSpeciesSolver.hpp b/SU2_CFD/include/solvers/CSpeciesSolver.hpp index 82898f83b83f..b9c09f06ad01 100644 --- a/SU2_CFD/include/solvers/CSpeciesSolver.hpp +++ b/SU2_CFD/include/solvers/CSpeciesSolver.hpp @@ -146,7 +146,7 @@ class CSpeciesSolver : public CScalarSolver { * \param[in] val_marker - Surface marker where the boundary condition is applied. */ void BC_Outlet(CGeometry* geometry, CSolver** solver_container, CNumerics* conv_numerics, CNumerics* visc_numerics, - CConfig* config, unsigned short val_marker) final; + CConfig* config, unsigned short val_marker) override; /*! * \brief Impose the isothermal wall Dirichlet boundary condition (value). diff --git a/SU2_CFD/src/solvers/CIncEulerSolver.cpp b/SU2_CFD/src/solvers/CIncEulerSolver.cpp index b7896ae4380d..8e68de647e4d 100644 --- a/SU2_CFD/src/solvers/CIncEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CIncEulerSolver.cpp @@ -2620,7 +2620,9 @@ void CIncEulerSolver::BC_Outlet(CGeometry *geometry, CSolver **solver_container, string Marker_Tag = config->GetMarker_All_TagBound(val_marker); su2double Normal[MAXNDIM] = {0.0}; + passivedouble nBackflow_loc = 0.0; + const bool backflow_prevention = config->GetInc_Outlet_BackflowPrevention(); INC_OUTLET_TYPE Kind_Outlet = config->GetKind_Inc_Outlet(Marker_Tag); /*--- Loop over all the vertices on this boundary marker ---*/ @@ -2655,6 +2657,16 @@ void CIncEulerSolver::BC_Outlet(CGeometry *geometry, CSolver **solver_container, P_domain = nodes->GetPressure(iPoint); + /*--- Compute the face area and normal velocity (positive = outflow, negative = backflow). ---*/ + + const su2double Area = GeometryToolbox::Norm(nDim, Normal); + const su2double Vn = GeometryToolbox::DotProduct(nDim, &V_domain[prim_idx.Velocity()], Normal) / Area; + + if (Vn < 0.0) { + SU2_OMP_ATOMIC + nBackflow_loc += 1.0; + } + /*--- Compute a boundary value for the pressure depending on whether we are prescribing a back pressure or a mass flow target. ---*/ @@ -2725,6 +2737,28 @@ void CIncEulerSolver::BC_Outlet(CGeometry *geometry, CSolver **solver_container, } + /*--- Backflow prevention via velocity reflection. + * + * When Vn < 0 (reversed flow entering through the outlet), the ghost velocity + * is reflected about the face: Vn_ghost = -Vn > 0. In vector form: + * + * V_ghost = V_domain - 2 * Vn * n_hat + * + * The factor of 2 (vs. the factor-1 zero-flux projection) gives the upwind + * Riemann solver a genuine outflow state on the ghost side, creating an active + * restoring force proportional to the backflow magnitude. + * + * A pressure penalty (0.5*rho*Vn^2) is intentionally NOT applied here. For + * variable-density / FGM cases the penalty is O(0.01-0.1 Pa) — too small to + * overcome expansion-driven backflow — yet it creates a thermodynamically + * inconsistent ghost state (modified pressure with Neumann enthalpy) that + * pushes FGM controlling variables outside the manifold. ---*/ + + if (Vn < 0.0 && backflow_prevention) { + for (iDim = 0; iDim < nDim; iDim++) + V_outlet[iDim + prim_idx.Velocity()] -= 2.0 * Vn * Normal[iDim] / Area; + } + /*--- Neumann condition for the temperature. ---*/ V_outlet[prim_idx.Temperature()] = nodes->GetTemperature(iPoint); @@ -2811,6 +2845,18 @@ void CIncEulerSolver::BC_Outlet(CGeometry *geometry, CSolver **solver_container, } END_SU2_OMP_FOR + /*--- Print a warning if backflow was detected on this outlet marker. + * Note: no MPI_Allreduce here — BC_Outlet is not called symmetrically + * across all ranks, so collective calls would corrupt the message queue. ---*/ + BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS { + if (nBackflow_loc > 0.5) { + cout << "WARNING [Rank " << rank << "]: Backflow detected at outlet marker \"" + << Marker_Tag << "\": " << static_cast(nBackflow_loc) + << " face(s) have reversed normal velocity." << endl; + } + } + END_SU2_OMP_SAFE_GLOBAL_ACCESS + } void CIncEulerSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver_container, CConfig *config, diff --git a/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp b/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp index de8d5e618531..22babc89342f 100644 --- a/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp +++ b/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp @@ -424,6 +424,143 @@ void CSpeciesFlameletSolver::BC_Inlet(CGeometry* geometry, CSolver** solver_cont CSpeciesSolver::BC_Inlet(geometry, solver_container, conv_numerics, visc_numerics, config, val_marker); } +void CSpeciesFlameletSolver::BC_Outlet(CGeometry* geometry, CSolver** solver_container, CNumerics* conv_numerics, + CNumerics* visc_numerics, CConfig* config, unsigned short val_marker) { + SU2_ZONE_SCOPED + + /*--- When backflow prevention is disabled, fall back to the base Neumann outlet. ---*/ + if (!config->GetInc_Outlet_BackflowPrevention()) { + CSpeciesSolver::BC_Outlet(geometry, solver_container, conv_numerics, visc_numerics, config, val_marker); + return; + } + + const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); + const string Marker_Tag = config->GetMarker_All_TagBound(val_marker); + + passivedouble nBackflow_loc = 0.0; + passivedouble Vn_min_loc = 0.0; /*--- Most negative Vn seen this call. ---*/ + + // Compute a manifold-consistent enthalpy state for backflow faces. + // V_ref is the inlet bulk velocity; used to normalise the blending weight so that + // alpha -> 0 for marginal backflow and alpha -> 1 for strong backflow. + string Inlet_Tag; + bool found_inlet = false; + for (auto iMarker = 0u; iMarker < config->GetnMarker_CfgFile(); iMarker++) { + const string tag = config->GetMarker_CfgFile_TagBound(iMarker); + if (config->GetMarker_CfgFile_KindBC(tag) == INLET_FLOW) { + Inlet_Tag = tag; + found_inlet = true; + break; + } + } + if (!found_inlet) { + /*--- No inlet boundary defined — nothing to blend against; fall back to Neumann. ---*/ + CSpeciesSolver::BC_Outlet(geometry, solver_container, conv_numerics, visc_numerics, config, val_marker); + return; + } + + const su2double T_inlet = config->GetInletTtotal(Inlet_Tag); + const su2double V_ref = config->GetInletPtotal(Inlet_Tag) / config->GetVelocity_Ref(); + + /*--- Determine thermodynamically consistent inlet enthalpy via Newton method. ---*/ + su2double H_inlet; + GetEnthFromTemp(solver_container[FLOW_SOL]->GetFluidModel(), + T_inlet, config->GetInlet_SpeciesVal(Inlet_Tag), &H_inlet); + + SU2_OMP_FOR_STAT(OMP_MIN_SIZE) + for (auto iVertex = 0u; iVertex < geometry->nVertex[val_marker]; iVertex++) { + const auto iPoint = geometry->vertex[val_marker][iVertex]->GetNode(); + + if (!geometry->nodes->GetDomain(iPoint)) continue; + + /*--- Outward face normal and normal velocity — needed for both strong and weak paths. ---*/ + su2double Normal[MAXNDIM] = {0.0}; + for (auto iDim = 0u; iDim < nDim; iDim++) Normal[iDim] = -geometry->vertex[val_marker][iVertex]->GetNormal(iDim); + + auto V_domain = solver_container[FLOW_SOL]->GetNodes()->GetPrimitive(iPoint); + const su2double Area = GeometryToolbox::Norm(nDim, Normal); + const su2double Vn = GeometryToolbox::DotProduct(nDim, &V_domain[prim_idx.Velocity()], Normal) / Area; + + if (Vn < 0.0) { + SU2_OMP_ATOMIC + nBackflow_loc += 1.0; + const passivedouble Vn_passive = SU2_TYPE::GetValue(Vn); + /*--- Bug fix: SU2_OMP_ATOMIC (#pragma omp atomic) only supports simple op= expressions, + * not min(). Use atomicMin which falls back to a critical section when needed. ---*/ + atomicMin(Vn_passive, Vn_min_loc); + } + + /*--- Strong BC path: prescribe inlet scalars on backflow faces; otherwise copy + from the interior neighbour as usual. ---*/ + if (config->GetMarker_StrongBC(Marker_Tag)) { + if (Vn < 0.0) { + nodes->SetSolution_Old(iPoint, Inlet_SpeciesVars[val_marker][iVertex]); + } else { + auto Point_Normal = geometry->vertex[val_marker][iVertex]->GetNormal_Neighbor(); + nodes->SetSolution_Old(iPoint, nodes->GetSolution(Point_Normal)); + } + LinSysRes.SetBlock_Zero(iPoint); + for (auto iVar = 0u; iVar < nVar; iVar++) Jacobian.DeleteValsRowi(iPoint, iVar); + continue; + } + + /*--- Weak BC: flow ghost state set by CIncEulerSolver::BC_Outlet. ---*/ + auto V_outlet = solver_container[FLOW_SOL]->GetCharacPrimVar(val_marker, iVertex); + + /*--- Choose scalar ghost state: + * Outflow (Vn >= 0): Neumann — interior values propagate outward (standard). + * Backflow (Vn < 0): Dirichlet — Use blended approach to create thermodynamically consistent state. ---*/ + su2double backflow_scalars[MAXNVAR]; + const su2double* scalar_ghost; + + if (Vn < 0.0) { + const su2double alpha = min(-Vn / V_ref, 1.0); + const su2double H_interior = nodes->GetSolution(iPoint, I_ENTH); + + /*--- Start from the interior solution, then blend the enthalpy toward the + * thermodynamically consistent inlet value and clamp PV to zero. ---*/ + for (auto iVar = 0u; iVar < nVar; iVar++) + backflow_scalars[iVar] = nodes->GetSolution(iPoint, iVar); + + backflow_scalars[I_ENTH] = (1.0 - alpha) * H_interior + alpha * H_inlet; + backflow_scalars[I_PROGVAR] = config->GetSpecies_Clipping_Min(I_PROGVAR); + + scalar_ghost = backflow_scalars; + } else { + scalar_ghost = nodes->GetSolution(iPoint); + } + + conv_numerics->SetPrimitive(V_domain, V_outlet); + conv_numerics->SetScalarVar(nodes->GetSolution(iPoint), scalar_ghost); + conv_numerics->SetNormal(Normal); + + if (dynamic_grid) + conv_numerics->SetGridVel(geometry->nodes->GetGridVel(iPoint), geometry->nodes->GetGridVel(iPoint)); + + if (conv_numerics->GetBoundedScalar()) { + const su2double* velocity = &V_outlet[prim_idx.Velocity()]; + const su2double density = solver_container[FLOW_SOL]->GetNodes()->GetDensity(iPoint); + conv_numerics->SetMassFlux(BoundedScalarBCFlux(iPoint, implicit, density, velocity, Normal)); + } + + auto residual = conv_numerics->ComputeResidual(config); + LinSysRes.AddBlock(iPoint, residual); + if (implicit) Jacobian.AddBlock2Diag(iPoint, residual.jacobian_i); + } + END_SU2_OMP_FOR + + /*--- Print backflow diagnostics once per BC call (no MPI reduction — safe). ---*/ + BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS { + if (nBackflow_loc > 0.5 && rank == MASTER_NODE) { + cout << " Flamelet BC_Outlet [" << Marker_Tag << "]: " + << static_cast(nBackflow_loc) + << " backflow face(s), Vn_min = " << Vn_min_loc << " m/s" + << " --> applying inlet scalars (PV=0, H=H_inlet) on backflow faces." << endl; + } + } + END_SU2_OMP_SAFE_GLOBAL_ACCESS +} + void CSpeciesFlameletSolver::BC_Isothermal_Wall_Generic(CGeometry* geometry, CSolver** solver_container, CNumerics* conv_numerics, CNumerics* visc_numerics, CConfig* config, unsigned short val_marker, bool cht_mode) { From 4e55810b9aa9eb6170ca98b25b7511eaf5089518 Mon Sep 17 00:00:00 2001 From: Josh Kelly Date: Wed, 29 Apr 2026 10:52:56 +0200 Subject: [PATCH 2/4] Add iteration switch for backflow treatment --- Common/include/CConfig.hpp | 10 +++++++- .../include/parallelization/mpi_structure.cpp | 24 +++++++++---------- Common/src/CConfig.cpp | 2 ++ SU2_CFD/src/solvers/CIncEulerSolver.cpp | 3 ++- .../src/solvers/CSpeciesFlameletSolver.cpp | 7 ++++-- 5 files changed, 30 insertions(+), 16 deletions(-) diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp index f454199e246b..272b9c6ddd64 100644 --- a/Common/include/CConfig.hpp +++ b/Common/include/CConfig.hpp @@ -635,7 +635,8 @@ class CConfig { unsigned short nInc_Outlet; /*!< \brief Number of inlet boundary treatment types listed. */ su2double Inc_Inlet_Damping; /*!< \brief Damping factor applied to the iterative updates to the velocity at a pressure inlet in incompressible flow. */ su2double Inc_Outlet_Damping; /*!< \brief Damping factor applied to the iterative updates to the pressure at a mass flow outlet in incompressible flow. */ - bool Inc_Outlet_BackflowPrevention; /*!< \brief Enable removal of the reversed normal velocity component at outlet faces where backflow is detected. */ + bool Inc_Outlet_BackflowPrevention; /*!< \brief Enable removal of the reversed normal velocity component at outlet faces where backflow is detected. */ + unsigned long Inc_Outlet_BackflowPrevention_Iter; /*!< \brief Number of outer iterations for which backflow prevention is active; after this the BC reverts to pure Neumann. Default ULONG_MAX (always active). */ bool InletUseNormal; /*!< \brief Flag for whether to use the local normal as the flow direction for a pressure inlet. */ su2double Linear_Solver_Error; /*!< \brief Min error of the linear solver for the implicit formulation. */ su2double Deform_Linear_Solver_Error; /*!< \brief Min error of the linear solver for the implicit formulation. */ @@ -5198,6 +5199,13 @@ class CConfig { */ bool GetInc_Outlet_BackflowPrevention(void) const { return Inc_Outlet_BackflowPrevention; } + /*! + * \brief Get the outer-iteration limit for backflow prevention at incompressible outlets. + * Backflow prevention is active while GetOuterIter() < this value. + * Returns ULONG_MAX (always active) when INC_OUTLET_BACKFLOW_PREVENTION_ITER is not set. + */ + unsigned long GetInc_Outlet_BackflowPrevention_Iter(void) const { return Inc_Outlet_BackflowPrevention_Iter; } + /*! * \brief Get the kind of mixing process for averaging quantities at the boundaries. * \return Kind of mixing process. diff --git a/Common/include/parallelization/mpi_structure.cpp b/Common/include/parallelization/mpi_structure.cpp index ae9a0ee25aef..fe2895c33d8c 100644 --- a/Common/include/parallelization/mpi_structure.cpp +++ b/Common/include/parallelization/mpi_structure.cpp @@ -101,12 +101,12 @@ void CBaseMPIWrapper::Error(const std::string& ErrorMsg, const std::string& Func /* Check if this rank must write the error message and do so. */ if (Rank == MinRankError) { - std::cout << std::endl << std::endl; - std::cout << "Error in \"" << FunctionName << "\": " << std::endl; - std::cout << "-------------------------------------------------------------------------" << std::endl; - std::cout << ErrorMsg << std::endl; - std::cout << "------------------------------ Error Exit -------------------------------" << std::endl; - std::cout << std::endl << std::endl; + std::cerr << std::endl << std::endl; + std::cerr << "Error in \"" << FunctionName << "\": " << std::endl; + std::cerr << "-------------------------------------------------------------------------" << std::endl; + std::cerr << ErrorMsg << std::endl; + std::cerr << "------------------------------ Error Exit -------------------------------" << std::endl; + std::cerr << std::endl << std::endl; } Abort(currentComm, EXIT_FAILURE); } @@ -130,12 +130,12 @@ void CBaseMPIWrapper::CopyData(const void* sendbuf, void* recvbuf, int size, Dat template void CBaseMPIWrapper::Error(const std::string& ErrorMsg, const std::string& FunctionName) { if (Rank == 0) { - std::cout << std::endl << std::endl; - std::cout << "Error in \"" << FunctionName << "\": " << std::endl; - std::cout << "-------------------------------------------------------------------------" << std::endl; - std::cout << ErrorMsg << std::endl; - std::cout << "------------------------------ Error Exit -------------------------------" << std::endl; - std::cout << std::endl << std::endl; + std::cerr << std::endl << std::endl; + std::cerr << "Error in \"" << FunctionName << "\": " << std::endl; + std::cerr << "-------------------------------------------------------------------------" << std::endl; + std::cerr << ErrorMsg << std::endl; + std::cerr << "------------------------------ Error Exit -------------------------------" << std::endl; + std::cerr << std::endl << std::endl; } Abort(currentComm, 0); } diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index 8b0886fceb6d..c681cac80ada 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -1398,6 +1398,8 @@ void CConfig::SetConfig_Options() { addDoubleOption("INC_OUTLET_DAMPING", Inc_Outlet_Damping, 0.1); /*!\brief INC_OUTLET_BACKFLOW_PREVENTION \n DESCRIPTION: Remove the reversed normal velocity component from the outlet ghost state when backflow is detected, preventing spurious inflow through the outlet. \ingroup Config*/ addBoolOption("INC_OUTLET_BACKFLOW_PREVENTION", Inc_Outlet_BackflowPrevention, false); + /*!\brief INC_OUTLET_BACKFLOW_PREVENTION_ITER \n DESCRIPTION: Number of outer iterations for which backflow prevention is active. After this iteration count the outlet BC reverts to pure Neumann, allowing physically real backflow to develop. Defaults to ULONG_MAX (always active). \ingroup Config*/ + addUnsignedLongOption("INC_OUTLET_BACKFLOW_PREVENTION_ITER", Inc_Outlet_BackflowPrevention_Iter, numeric_limits::max()); /*--- Options related to the species solver. ---*/ diff --git a/SU2_CFD/src/solvers/CIncEulerSolver.cpp b/SU2_CFD/src/solvers/CIncEulerSolver.cpp index 8e68de647e4d..32717bdd8ab8 100644 --- a/SU2_CFD/src/solvers/CIncEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CIncEulerSolver.cpp @@ -2622,7 +2622,8 @@ void CIncEulerSolver::BC_Outlet(CGeometry *geometry, CSolver **solver_container, su2double Normal[MAXNDIM] = {0.0}; passivedouble nBackflow_loc = 0.0; - const bool backflow_prevention = config->GetInc_Outlet_BackflowPrevention(); + const bool backflow_prevention = config->GetInc_Outlet_BackflowPrevention() && + (config->GetOuterIter() < config->GetInc_Outlet_BackflowPrevention_Iter()); INC_OUTLET_TYPE Kind_Outlet = config->GetKind_Inc_Outlet(Marker_Tag); /*--- Loop over all the vertices on this boundary marker ---*/ diff --git a/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp b/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp index 22babc89342f..fefce339ba0d 100644 --- a/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp +++ b/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp @@ -428,8 +428,11 @@ void CSpeciesFlameletSolver::BC_Outlet(CGeometry* geometry, CSolver** solver_con CNumerics* visc_numerics, CConfig* config, unsigned short val_marker) { SU2_ZONE_SCOPED - /*--- When backflow prevention is disabled, fall back to the base Neumann outlet. ---*/ - if (!config->GetInc_Outlet_BackflowPrevention()) { + /*--- When backflow prevention is disabled, or the outer-iteration window has elapsed, + * fall back to the base Neumann outlet so that physically real backflow (e.g. + * recirculation zones at a converged solution) is not artificially suppressed. ---*/ + if (!config->GetInc_Outlet_BackflowPrevention() || + config->GetOuterIter() >= config->GetInc_Outlet_BackflowPrevention_Iter()) { CSpeciesSolver::BC_Outlet(geometry, solver_container, conv_numerics, visc_numerics, config, val_marker); return; } From a30859fa9621db9243b8838264ebf3f709e2d7f0 Mon Sep 17 00:00:00 2001 From: Josh Kelly Date: Thu, 21 May 2026 14:27:57 +0200 Subject: [PATCH 3/4] Added documentaiton, updated comments --- SU2_CFD/src/solvers/CIncEulerSolver.cpp | 34 ++++++++----------- .../src/solvers/CSpeciesFlameletSolver.cpp | 14 +++----- config_template.cfg | 10 ++++++ 3 files changed, 30 insertions(+), 28 deletions(-) diff --git a/SU2_CFD/src/solvers/CIncEulerSolver.cpp b/SU2_CFD/src/solvers/CIncEulerSolver.cpp index 32717bdd8ab8..7b1fdf9a6b76 100644 --- a/SU2_CFD/src/solvers/CIncEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CIncEulerSolver.cpp @@ -2620,7 +2620,7 @@ void CIncEulerSolver::BC_Outlet(CGeometry *geometry, CSolver **solver_container, string Marker_Tag = config->GetMarker_All_TagBound(val_marker); su2double Normal[MAXNDIM] = {0.0}; - passivedouble nBackflow_loc = 0.0; + int nBackflow_loc = 0; const bool backflow_prevention = config->GetInc_Outlet_BackflowPrevention() && (config->GetOuterIter() < config->GetInc_Outlet_BackflowPrevention_Iter()); @@ -2665,7 +2665,7 @@ void CIncEulerSolver::BC_Outlet(CGeometry *geometry, CSolver **solver_container, if (Vn < 0.0) { SU2_OMP_ATOMIC - nBackflow_loc += 1.0; + nBackflow_loc += 1; } /*--- Compute a boundary value for the pressure depending on whether @@ -2740,20 +2740,18 @@ void CIncEulerSolver::BC_Outlet(CGeometry *geometry, CSolver **solver_container, /*--- Backflow prevention via velocity reflection. * - * When Vn < 0 (reversed flow entering through the outlet), the ghost velocity - * is reflected about the face: Vn_ghost = -Vn > 0. In vector form: + * When reverse flow is detected at the outlet (Vn < 0), + * we reflect the normal velocity component to prevent + * backflow and stabilize the solution. + * + * Factor of 2 is used to ensure that the reflected velocity + * is sufficient to counteract the backflow, by applying a + * restorative force that is proportional to the detected backflow velocity. * - * V_ghost = V_domain - 2 * Vn * n_hat - * - * The factor of 2 (vs. the factor-1 zero-flux projection) gives the upwind - * Riemann solver a genuine outflow state on the ghost side, creating an active - * restoring force proportional to the backflow magnitude. - * - * A pressure penalty (0.5*rho*Vn^2) is intentionally NOT applied here. For - * variable-density / FGM cases the penalty is O(0.01-0.1 Pa) — too small to - * overcome expansion-driven backflow — yet it creates a thermodynamically - * inconsistent ghost state (modified pressure with Neumann enthalpy) that - * pushes FGM controlling variables outside the manifold. ---*/ + * Dynamic pressure is intentionally NOT applied here. For + * variable-density / FGM cases the penalty is too small to + * overcome expansion-driven backflow and creates a thermodynamically + * inconsistent ghost state that pushes FGM controlling variables outside the manifold. ---*/ if (Vn < 0.0 && backflow_prevention) { for (iDim = 0; iDim < nDim; iDim++) @@ -2846,11 +2844,9 @@ void CIncEulerSolver::BC_Outlet(CGeometry *geometry, CSolver **solver_container, } END_SU2_OMP_FOR - /*--- Print a warning if backflow was detected on this outlet marker. - * Note: no MPI_Allreduce here — BC_Outlet is not called symmetrically - * across all ranks, so collective calls would corrupt the message queue. ---*/ + /*--- Print a warning if backflow was detected on this outlet marker. ---*/ BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS { - if (nBackflow_loc > 0.5) { + if (nBackflow_loc > 0) { cout << "WARNING [Rank " << rank << "]: Backflow detected at outlet marker \"" << Marker_Tag << "\": " << static_cast(nBackflow_loc) << " face(s) have reversed normal velocity." << endl; diff --git a/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp b/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp index fefce339ba0d..a6df817eea68 100644 --- a/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp +++ b/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp @@ -440,8 +440,8 @@ void CSpeciesFlameletSolver::BC_Outlet(CGeometry* geometry, CSolver** solver_con const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); const string Marker_Tag = config->GetMarker_All_TagBound(val_marker); - passivedouble nBackflow_loc = 0.0; - passivedouble Vn_min_loc = 0.0; /*--- Most negative Vn seen this call. ---*/ + int nBackflow_loc = 0; /*--- Number of faces with backflow seen this call. ---*/ + su2double Vn_min_loc = 0.0; /*--- Most negative Vn seen this call. ---*/ // Compute a manifold-consistent enthalpy state for backflow faces. // V_ref is the inlet bulk velocity; used to normalise the blending weight so that @@ -486,11 +486,7 @@ void CSpeciesFlameletSolver::BC_Outlet(CGeometry* geometry, CSolver** solver_con if (Vn < 0.0) { SU2_OMP_ATOMIC - nBackflow_loc += 1.0; - const passivedouble Vn_passive = SU2_TYPE::GetValue(Vn); - /*--- Bug fix: SU2_OMP_ATOMIC (#pragma omp atomic) only supports simple op= expressions, - * not min(). Use atomicMin which falls back to a critical section when needed. ---*/ - atomicMin(Vn_passive, Vn_min_loc); + nBackflow_loc += 1; } /*--- Strong BC path: prescribe inlet scalars on backflow faces; otherwise copy @@ -554,9 +550,9 @@ void CSpeciesFlameletSolver::BC_Outlet(CGeometry* geometry, CSolver** solver_con /*--- Print backflow diagnostics once per BC call (no MPI reduction — safe). ---*/ BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS { - if (nBackflow_loc > 0.5 && rank == MASTER_NODE) { + if (nBackflow_loc > 0 && rank == MASTER_NODE) { cout << " Flamelet BC_Outlet [" << Marker_Tag << "]: " - << static_cast(nBackflow_loc) + << nBackflow_loc << " backflow face(s), Vn_min = " << Vn_min_loc << " m/s" << " --> applying inlet scalars (PV=0, H=H_inlet) on backflow faces." << endl; } diff --git a/config_template.cfg b/config_template.cfg index ab8175bb6802..453fa1391803 100644 --- a/config_template.cfg +++ b/config_template.cfg @@ -304,6 +304,16 @@ INC_OUTLET_TYPE= PRESSURE_OUTLET % Damping coefficient for iterative updates at mass flow outlets. (0.1 by default) INC_OUTLET_DAMPING= 0.1 % +% Activate backflow prevention at outlet markers for incompressible flows (NO, YES). +% If YES, the normal velocity at outlet faces will be set to zero if backflow is detected. +% This is only active for the first INC_OUTLET_BACKFLOW_PREVENTION_ITER iterations, +% after which the BC will switch to the specified type in INC_OUTLET_TYPE. +INC_OUTLET_BACKFLOW_PREVENTION= NO +% +% Number of iterations to apply backflow prevention at outlet markers for incompressible flows +% If no value is specified, backflow prevention will be applied for the entire simulation. +INC_OUTLET_BACKFLOW_PREVENTION_ITER= 1000 +% % Bulk Modulus for computing the Mach number BULK_MODULUS= 1.42E5 % Epsilon^2 multipier in Beta calculation for incompressible preconditioner. From 90735cc96648b66b4a7ad922b83897787d442e9c Mon Sep 17 00:00:00 2001 From: Josh Kelly Date: Thu, 21 May 2026 14:39:53 +0200 Subject: [PATCH 4/4] Trailing whitespace --- config_template.cfg | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/config_template.cfg b/config_template.cfg index 453fa1391803..711df7cdb808 100644 --- a/config_template.cfg +++ b/config_template.cfg @@ -305,7 +305,7 @@ INC_OUTLET_TYPE= PRESSURE_OUTLET INC_OUTLET_DAMPING= 0.1 % % Activate backflow prevention at outlet markers for incompressible flows (NO, YES). -% If YES, the normal velocity at outlet faces will be set to zero if backflow is detected. +% If YES, the normal velocity at outlet faces will be set to zero if backflow is detected. % This is only active for the first INC_OUTLET_BACKFLOW_PREVENTION_ITER iterations, % after which the BC will switch to the specified type in INC_OUTLET_TYPE. INC_OUTLET_BACKFLOW_PREVENTION= NO