From 595ce6b98d9830830eb8faa96bedb02a790ab718 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Sun, 14 Jun 2026 21:06:08 -0700 Subject: [PATCH 1/2] detect temperature outliers --- Common/include/CConfig.hpp | 6 + Common/src/CConfig.cpp | 7 +- .../numerics_simd/flow/convection/common.hpp | 17 +- .../numerics_simd/flow/convection/upwind.hpp | 17 +- SU2_CFD/include/solvers/CEulerSolver.hpp | 53 +++--- .../include/solvers/CFVMFlowSolverBase.hpp | 44 ++--- .../include/solvers/CFVMFlowSolverBase.inl | 45 ++--- SU2_CFD/include/solvers/CIncEulerSolver.hpp | 7 - SU2_CFD/include/solvers/CNEMOEulerSolver.hpp | 7 - SU2_CFD/include/solvers/CSolver.hpp | 8 + SU2_CFD/include/variables/CEulerVariable.hpp | 10 ++ SU2_CFD/src/iteration/CFluidIteration.cpp | 7 +- SU2_CFD/src/solvers/CEulerSolver.cpp | 156 +++++++++++++++++- SU2_CFD/src/solvers/CIncEulerSolver.cpp | 6 - SU2_CFD/src/solvers/CNEMOEulerSolver.cpp | 6 - SU2_CFD/src/solvers/CSolver.cpp | 7 +- SU2_CFD/src/solvers/CTurbSASolver.cpp | 4 + SU2_CFD/src/variables/CEulerVariable.cpp | 2 + config_template.cfg | 8 + 19 files changed, 289 insertions(+), 128 deletions(-) diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp index e44ad6a10a7..55bb1cfed04 100644 --- a/Common/include/CConfig.hpp +++ b/Common/include/CConfig.hpp @@ -815,6 +815,7 @@ class CConfig { unsigned short ActDisk_Jump; /*!< \brief Format of the output files. */ unsigned long StartWindowIteration; /*!< \brief Starting Iteration for long time Windowing apporach . */ unsigned short nCFL_AdaptParam; /*!< \brief Number of CFL parameters provided in config. */ + unsigned long outlierMitigationParam[4]; /*!< \brief Parameters of outlier mitigation strategy. */ bool CFL_Adapt; /*!< \brief Use adaptive CFL number. */ bool HB_Precondition; /*!< \brief Flag to turn on harmonic balance source term preconditioning */ su2double RefArea, /*!< \brief Reference area for coefficient computation. */ @@ -1712,6 +1713,11 @@ class CConfig { */ bool GetCFL_Adapt(void) const { return CFL_Adapt; } + /*! + * \brief Get the outlier mitigation parameters. + */ + const unsigned long* GetOutlierMitigationParam() const { return outlierMitigationParam; } + /*! * \brief Get the value of the limits for the sections. * \return Value of the limits for the sections. diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index c168ced0830..6bdb8feb99c 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -1819,6 +1819,11 @@ void CConfig::SetConfig_Options() { addDoubleOption("CFL_NUMBER", CFLFineGrid, 1.25); /* DESCRIPTION: Max time step in local time stepping simulations */ addDoubleOption("MAX_DELTA_TIME", Max_DeltaTime, 1000000); + /* !\brief OUTLIER_MITIGATION_PARAM + * DESCRIPTION: Parameters of the outlier mitigation strategy: start iteration, update frequency, print frequency, + * and number of standard deviations (N * sigma) to identify outliers statistically. \ingroup Config*/ + outlierMitigationParam[0] = 999999; outlierMitigationParam[1] = 5; outlierMitigationParam[2] = 2; outlierMitigationParam[3] = 5; + addULongArrayOption("OUTLIER_MITIGATION_PARAM", 4, true, outlierMitigationParam); /* DESCRIPTION: Activate The adaptive CFL number. */ addBoolOption("CFL_ADAPT", CFL_Adapt, false); /* !\brief CFL_ADAPT_PARAM @@ -2029,7 +2034,7 @@ void CConfig::SetConfig_Options() { addDoubleOption("MUSCL_KAPPA_FLOW", MUSCL_Kappa_Flow, 0.0); /*!\brief RAMP_MUSCL \n DESCRIPTION: Enable ramping of the MUSCL scheme from 1st to 2nd order using specified method*/ addBoolOption("RAMP_MUSCL", RampMUSCL, false); - /*! brief RAMP_OUTLET_COEFF \n DESCRIPTION: the 1st coeff is the ramp start iteration, + /*! brief RAMP_MUSCL_COEFF \n DESCRIPTION: the 1st coeff is the ramp start iteration, * the 2nd coeff is the iteration update frequenct, 3rd coeff is the total number of iterations */ RampMUSCLParam.rampMUSCLCoeff[0] = 0.0; RampMUSCLParam.rampMUSCLCoeff[1] = 1.0; RampMUSCLParam.rampMUSCLCoeff[2] = 500.0; addULongArrayOption("RAMP_MUSCL_COEFF", 3, false, RampMUSCLParam.rampMUSCLCoeff); diff --git a/SU2_CFD/include/numerics_simd/flow/convection/common.hpp b/SU2_CFD/include/numerics_simd/flow/convection/common.hpp index 4796f80ce03..b6c4c23c775 100644 --- a/SU2_CFD/include/numerics_simd/flow/convection/common.hpp +++ b/SU2_CFD/include/numerics_simd/flow/convection/common.hpp @@ -188,6 +188,7 @@ FORCEINLINE void musclEdgeLimited(const Int& iPoint, * \param[in] V1st - Pair of compressible flow primitives for nodes i,j. * \param[in] vector_ij - Distance vector from i to j. * \param[in] solution - Entire solution container (a derived CVariable). + * \param[out] nonPhysical - Signals that the edge is treated as non-physical. * \return Pair of primitive variables. */ template @@ -201,7 +202,8 @@ FORCEINLINE CPair reconstructPrimitives(const Int& iEdge, const LIMITER limiterType, const CPair& V1st, const VectorDbl& vector_ij, - const VariableType& solution) { + const VariableType& solution, + Double& nonPhysical) { static_assert(ReconVarType::nVar <= PrimVarType::nVar); const auto& gradients = solution.GetGradient_Reconstruction(); @@ -262,15 +264,20 @@ FORCEINLINE CPair reconstructPrimitives(const Int& iEdge, const Double neg_sound_speed = enthalpy * (R+1) < 0.5 * v_squared; /*--- Revert to first order if the state is non-physical. ---*/ - Double bad_recon = fmax(neg_p_or_rho, neg_sound_speed); + Double nonPhysical = fmax(neg_p_or_rho, neg_sound_speed); /*--- Handle SIMD dimensions 1 by 1. ---*/ for (size_t k = 0; k < Double::Size; ++k) { - bad_recon[k] = solution.UpdateNonPhysicalEdgeCounter(iEdge[k], bad_recon[k]); + nonPhysical[k] = solution.UpdateNonPhysicalEdgeCounter(iEdge[k], nonPhysical[k]); + nonPhysical[k] = fmax(nonPhysical[k], + fmax(solution.OutlierMitigation(iPoint[k]), + solution.OutlierMitigation(jPoint[k])) / VariableType::MAX_OUTLIER_MITIGATION); } for (size_t iVar = 0; iVar < ReconVarType::nVar; ++iVar) { - V.i.all(iVar) = bad_recon * V1st.i.all(iVar) + (1-bad_recon) * V.i.all(iVar); - V.j.all(iVar) = bad_recon * V1st.j.all(iVar) + (1-bad_recon) * V.j.all(iVar); + V.i.all(iVar) = nonPhysical * V1st.i.all(iVar) + (1-nonPhysical) * V.i.all(iVar); + V.j.all(iVar) = nonPhysical * V1st.j.all(iVar) + (1-nonPhysical) * V.j.all(iVar); } + } else { + nonPhysical = 0; } return V; } diff --git a/SU2_CFD/include/numerics_simd/flow/convection/upwind.hpp b/SU2_CFD/include/numerics_simd/flow/convection/upwind.hpp index ad958a7e2f0..090fb111332 100644 --- a/SU2_CFD/include/numerics_simd/flow/convection/upwind.hpp +++ b/SU2_CFD/include/numerics_simd/flow/convection/upwind.hpp @@ -118,8 +118,10 @@ class CUpwindBase : public Base { V1st.j.all = gatherVariables(jPoint, solution.GetPrimitive()); /*--- Recompute density and enthalpy instead of reconstructing. ---*/ + Double nonPhysical; auto V = reconstructPrimitives >( - iEdge, iPoint, jPoint, gamma, gasConst, muscl, umusclKappa, umusclRamp, typeLimiter, V1st, vector_ij, solution); + iEdge, iPoint, jPoint, gamma, gasConst, muscl, umusclKappa, umusclRamp, + typeLimiter, V1st, vector_ij, solution, nonPhysical); /*--- Compute conservative variables. ---*/ @@ -132,8 +134,8 @@ class CUpwindBase : public Base { const auto derived = static_cast(this); VectorDbl flux; MatrixDbl jac_i, jac_j; - derived->finalizeFlux(flux, jac_i, jac_j, implicit, area, unitNormal, - normal, V, U, iPoint, jPoint, solution, geometry); + derived->finalizeFlux(flux, jac_i, jac_j, implicit, area, unitNormal, normal, + V, U, iPoint, jPoint, nonPhysical, solution, geometry); /*--- Add the contributions from the base class (static decorator). ---*/ @@ -199,6 +201,7 @@ class CRoeScheme : public CUpwindBase, Decorator> { const CPair& U, const Int& iPoint, const Int& jPoint, + const Double& nonPhysical, const CEulerVariable& solution, const CGeometry& geometry, Ts&...) const { @@ -227,10 +230,9 @@ class CRoeScheme : public CUpwindBase, Decorator> { /*--- Apply Mavriplis' entropy correction to eigenvalues. ---*/ - Double maxLambda = abs(projVel) + roeAvg.speedSound; - + Double lambdaMin = fmax(entropyFix, nonPhysical) * (abs(projVel) + roeAvg.speedSound); for (size_t iVar = 0; iVar < nVar; ++iVar) { - lambda(iVar) = fmax(abs(lambda(iVar)), entropyFix*maxLambda); + lambda(iVar) = fmax(abs(lambda(iVar)), lambdaMin); } /*--- Inviscid fluxes and Jacobians. ---*/ @@ -348,6 +350,7 @@ class CMSWScheme : public CUpwindBase, Decorator> { const CPair& U, const Int& iPoint, const Int& jPoint, + const Double& nonPhysical, const CEulerVariable& solution, const CGeometry& geometry, Ts&...) const { @@ -358,7 +361,7 @@ class CMSWScheme : public CUpwindBase, Decorator> { const auto sj = gatherVariables(jPoint, solution.GetSensor()); const Double dp = fmax(si, sj) - alpha * 0.06; - const Double w = 0.25 * (1 - sign(dp)) * (1 - exp(-100 * abs(dp))); + const Double w = 0.25 * (1 - sign(dp) * (1 - exp(-100 * abs(dp)))) * (1 - nonPhysical); const Double onemw = 1 - w; CPair> Vweighted; diff --git a/SU2_CFD/include/solvers/CEulerSolver.hpp b/SU2_CFD/include/solvers/CEulerSolver.hpp index 46eed006b04..dc932095bc4 100644 --- a/SU2_CFD/include/solvers/CEulerSolver.hpp +++ b/SU2_CFD/include/solvers/CEulerSolver.hpp @@ -116,6 +116,9 @@ class CEulerSolver : public CFVMFlowSolverBase FluidModel; /*!< \brief fluid model used in the solver. */ + /*!< \brief Variables for outlier detection. */ + su2double MeanTemperature, StdDevTemperature; + /*--- Turbomachinery Solver Variables ---*/ vector AverageFlux; @@ -782,11 +785,17 @@ class CEulerSolver : public CFVMFlowSolverBase - void CompleteImplicitIteration_impl(CGeometry *geometry, CConfig *config) { - - if (compute_ur) ComputeUnderRelaxationFactor(config); - - /*--- Update solution with under-relaxation and communicate it. ---*/ - - if (!config->GetContinuous_Adjoint()) { - SU2_OMP_FOR_STAT(omp_chunk_size) - for (unsigned long iPoint = 0; iPoint < nPointDomain; iPoint++) { - for (unsigned short iVar = 0; iVar < nVar; iVar++) { - nodes->AddSolution(iPoint, iVar, nodes->GetUnderRelaxation(iPoint)*LinSysSol[iPoint*nVar+iVar]); - } - } - END_SU2_OMP_FOR - } - - for (unsigned short iPeriodic = 1; iPeriodic <= config->GetnMarker_Periodic()/2; iPeriodic++) { - InitiatePeriodicComms(geometry, config, iPeriodic, PERIODIC_IMPLICIT); - CompletePeriodicComms(geometry, config, iPeriodic, PERIODIC_IMPLICIT); - } - - InitiateComms(geometry, config, MPI_QUANTITIES::SOLUTION); - CompleteComms(geometry, config, MPI_QUANTITIES::SOLUTION); - - /*--- For verification cases, compute the global error metrics. ---*/ - ComputeVerificationError(geometry, config); - } - /*! * \brief Evaluate the vorticity and strain rate magnitude. */ @@ -1075,6 +1044,13 @@ class CFVMFlowSolverBase : public CSolver { */ void ImplicitEuler_Iteration(CGeometry *geometry, CSolver **solver_container, CConfig *config) final; + /*! + * \brief Complete an implicit iteration. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] config - Definition of the particular problem. + */ + void CompleteImplicitIteration(CGeometry *geometry, CSolver**, CConfig *config) final; + /*! * \brief Set the total residual adding the term that comes from the Dual Time Strategy. * \param[in] geometry - Geometrical definition of the problem. diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl index 05ec25a0e56..f651db091fe 100644 --- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl +++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl @@ -596,42 +596,33 @@ void CFVMFlowSolverBase::ComputeVerificationError(CGeometry* geometry, CCo } template -void CFVMFlowSolverBase::ComputeUnderRelaxationFactor(const CConfig* config) { +void CFVMFlowSolverBase::CompleteImplicitIteration(CGeometry *geometry, CSolver**, CConfig *config) { SU2_ZONE_SCOPED - /* Loop over the solution update given by relaxing the linear - system for this nonlinear iteration. */ + if constexpr (R == ENUM_REGIME::COMPRESSIBLE) ComputeUnderRelaxationFactor(config); - const su2double allowableRatio = config->GetMaxUpdateFractionFlow(); + /*--- Update solution with under-relaxation and communicate it. ---*/ - SU2_OMP_FOR_STAT(omp_chunk_size) - for (unsigned long iPoint = 0; iPoint < nPointDomain; iPoint++) { - su2double localUnderRelaxation = 1.0; - - for (unsigned short iVar = 0; iVar < nVar; iVar++) { - /* We impose a limit on the maximum percentage that the - density and energy can change over a nonlinear iteration. */ - - if ((iVar == 0) || (iVar == nVar - 1)) { - const unsigned long index = iPoint * nVar + iVar; - su2double ratio = fabs(LinSysSol[index]) / (fabs(nodes->GetSolution(iPoint, iVar)) + EPS); - if (ratio > allowableRatio) { - localUnderRelaxation = min(allowableRatio / ratio, localUnderRelaxation); - } + if (!config->GetContinuous_Adjoint()) { + SU2_OMP_FOR_STAT(omp_chunk_size) + for (unsigned long iPoint = 0; iPoint < nPointDomain; iPoint++) { + for (unsigned short iVar = 0; iVar < nVar; iVar++) { + nodes->AddSolution(iPoint, iVar, nodes->GetUnderRelaxation(iPoint) * LinSysSol(iPoint, iVar)); } } + END_SU2_OMP_FOR + } - /* Threshold the relaxation factor in the event that there is - a very small value. This helps avoid catastrophic crashes due - to non-realizable states by canceling the update. */ - - if (localUnderRelaxation < 1e-10) localUnderRelaxation = 0.0; + for (unsigned short iPeriodic = 1; iPeriodic <= config->GetnMarker_Periodic()/2; iPeriodic++) { + InitiatePeriodicComms(geometry, config, iPeriodic, PERIODIC_IMPLICIT); + CompletePeriodicComms(geometry, config, iPeriodic, PERIODIC_IMPLICIT); + } - /* Store the under-relaxation factor for this point. */ + InitiateComms(geometry, config, MPI_QUANTITIES::SOLUTION); + CompleteComms(geometry, config, MPI_QUANTITIES::SOLUTION); - nodes->SetUnderRelaxation(iPoint, localUnderRelaxation); - } - END_SU2_OMP_FOR + /*--- For verification cases, compute the global error metrics. ---*/ + ComputeVerificationError(geometry, config); } template diff --git a/SU2_CFD/include/solvers/CIncEulerSolver.hpp b/SU2_CFD/include/solvers/CIncEulerSolver.hpp index 98c3d21a8f2..9c18691cf56 100644 --- a/SU2_CFD/include/solvers/CIncEulerSolver.hpp +++ b/SU2_CFD/include/solvers/CIncEulerSolver.hpp @@ -365,13 +365,6 @@ class CIncEulerSolver : public CFVMFlowSolverBase::max(); + + /*! + * \brief Marks outliers (0 ok, MAX_OUTLIER_MITIGATION maximum mitigation). + */ + su2vector OutlierMitigation; + }; diff --git a/SU2_CFD/src/iteration/CFluidIteration.cpp b/SU2_CFD/src/iteration/CFluidIteration.cpp index 35e2905a7c5..ca06e928c6c 100644 --- a/SU2_CFD/src/iteration/CFluidIteration.cpp +++ b/SU2_CFD/src/iteration/CFluidIteration.cpp @@ -132,12 +132,13 @@ void CFluidIteration::Iterate(COutput* output, CIntegration**** integration, CGe /*--- Adapt the CFL number using an exponential progression with under-relaxation approach. ---*/ - if ((config[val_iZone]->GetCFL_Adapt() == YES) && (!disc_adj)) { - SU2_OMP_PARALLEL + SU2_OMP_PARALLEL + if (!disc_adj) { solver[val_iZone][val_iInst][MESH_0][FLOW_SOL]->AdaptCFLNumber(geometry[val_iZone][val_iInst], solver[val_iZone][val_iInst], config[val_iZone]); - END_SU2_OMP_PARALLEL + solver[val_iZone][val_iInst][MESH_0][FLOW_SOL]->IdentifySolutionOutliers(config[val_iZone], InnerIter); } + END_SU2_OMP_PARALLEL /*--- Call Dynamic mesh update if AEROELASTIC motion was specified ---*/ diff --git a/SU2_CFD/src/solvers/CEulerSolver.cpp b/SU2_CFD/src/solvers/CEulerSolver.cpp index e5737ceebb1..79aa8fe29ca 100644 --- a/SU2_CFD/src/solvers/CEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CEulerSolver.cpp @@ -2520,10 +2520,162 @@ void CEulerSolver::PrepareImplicitIteration(CGeometry *geometry, CSolver**, CCon PrepareImplicitIteration_impl(precond, geometry, config); } -void CEulerSolver::CompleteImplicitIteration(CGeometry *geometry, CSolver**, CConfig *config) { +void CEulerSolver::ComputeUnderRelaxationFactor(const CConfig* config) { SU2_ZONE_SCOPED - CompleteImplicitIteration_impl(geometry, config); + /*--- Loop over the solution update given by relaxing the linear system for this + * nonlinear iteration and impose a limit on the maximum percentage that the + * density and static energy can change. */ + + const su2double allowableRatio = config->GetMaxUpdateFractionFlow(); + + SU2_OMP_FOR_STAT(omp_chunk_size) + for (unsigned long iPoint = 0; iPoint < nPointDomain; iPoint++) { + su2double ratio = fabs(LinSysSol(iPoint, 0)) / max(nodes->GetSolution(iPoint, 0), EPS); + su2double e_old = nodes->GetSolution(iPoint, nVar - 1); + su2double e_new = e_old + LinSysSol(iPoint, nVar - 1); + for (unsigned short jVar = 1; jVar <= nDim; jVar++) { + e_old -= 0.5 * pow(nodes->GetSolution(iPoint, jVar), 2); + e_new -= 0.5 * pow(nodes->GetSolution(iPoint, jVar) + LinSysSol(iPoint, jVar), 2); + } + ratio = fmax(ratio, fabs(e_new - e_old) / max(e_old, EPS)); + + su2double localUnderRelaxation = fmin(allowableRatio / fmax(ratio, EPS), 1); + + /* Threshold the relaxation factor in the event that there is + a very small value. This helps avoid catastrophic crashes due + to non-realizable states by canceling the update. */ + + if (localUnderRelaxation < 1e-10) localUnderRelaxation = 0.0; + + nodes->SetUnderRelaxation(iPoint, localUnderRelaxation); + } + END_SU2_OMP_FOR +} + +void CEulerSolver::IdentifySolutionOutliers(const CConfig *config, unsigned long iter) { + SU2_ZONE_SCOPED + + const unsigned long startIteration = config->GetOutlierMitigationParam()[0]; + const unsigned long updateFrequency = config->GetOutlierMitigationParam()[1]; + const unsigned long printFrequency = config->GetOutlierMitigationParam()[2]; + const int nSigma = config->GetOutlierMitigationParam()[3]; + + if (iter < startIteration) return; + + auto DetermineBinAndUpdatePoint = [&](const unsigned long iPoint) { + const su2double t = nodes->GetTemperature(iPoint); + const auto i = nSigma + min(max(-nSigma, SU2_TYPE::Int((t - MeanTemperature) / max(StdDevTemperature, EPS))), nSigma); + if (i == 0 || i == 2 * nSigma) { + if (nodes->OutlierMitigation(iPoint) == 0) { + /*--- Start mitigating with maximum strength if the point was just identified as outlier. ---*/ + nodes->OutlierMitigation(iPoint) = CEulerVariable::MAX_OUTLIER_MITIGATION; + } else if (nodes->OutlierMitigation(iPoint) < CEulerVariable::MAX_OUTLIER_MITIGATION) { + /*--- The point became an outlier again after reducing mitigations (below), increase them slowly. ---*/ + ++nodes->OutlierMitigation(iPoint); + } + } else if (nodes->OutlierMitigation(iPoint) > 0) { + /*--- Not an outlier anymore, try to slowly reduce the mitigations. ---*/ + --nodes->OutlierMitigation(iPoint); + } + return i; + }; + + /*--- Recompute mean and std deviation of temperature or use the stored values. ---*/ + if (iter == startIteration || iter % updateFrequency == 0) { + su2double localSum = 0; + static su2double nPointGlobal; + SU2_OMP_MASTER + MeanTemperature = 0; + + SU2_OMP_FOR_STAT(omp_chunk_size) + for (unsigned long iPoint = 0; iPoint < nPointDomain; ++iPoint) { + localSum += nodes->GetTemperature(iPoint); + } + END_SU2_OMP_FOR + + atomicAdd(localSum, MeanTemperature); + + BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS { + su2double tmp[2] = {MeanTemperature, static_cast(nPointDomain)}, global[2]; + SU2_MPI::Allreduce(tmp, global, 2, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm()); + nPointGlobal = global[1]; + MeanTemperature = global[0] / nPointGlobal; + } + END_SU2_OMP_SAFE_GLOBAL_ACCESS + + localSum = 0; + SU2_OMP_MASTER + StdDevTemperature = 0; + + SU2_OMP_FOR_STAT(omp_chunk_size) + for (unsigned long iPoint = 0; iPoint < nPointDomain; iPoint ++) { + localSum += pow(MeanTemperature - nodes->GetTemperature(iPoint), 2); + } + END_SU2_OMP_FOR + + atomicAdd(localSum, StdDevTemperature); + + BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS { + SU2_MPI::Allreduce(&StdDevTemperature, &localSum, 1, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm()); + StdDevTemperature = sqrt(localSum / nPointGlobal); + } + END_SU2_OMP_SAFE_GLOBAL_ACCESS + + const auto nBins = 2 * nSigma + 2; + std::vector bins(nBins, 0); + unsigned long nPointLocal = 0; + SU2_OMP_MASTER + nPointGlobal = 0; + + SU2_OMP_FOR_STAT(omp_chunk_size) + for (unsigned long iPoint = 0; iPoint < nPoint; ++iPoint) { + const auto i = DetermineBinAndUpdatePoint(iPoint); + if (iPoint < nPointDomain) { + nPointLocal += static_cast(nodes->OutlierMitigation(iPoint) > 0); + + SU2_OMP_ATOMIC + ++bins[i]; + } + } + END_SU2_OMP_FOR + + SU2_OMP_ATOMIC + nPointGlobal += nPointLocal; + + BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS + if (iter == 0 || iter % (updateFrequency * printFrequency) == 0) { + bins.back() = nPointGlobal; + std::vector global(nBins); + SU2_MPI::Reduce(bins.data(), global.data(), nBins, MPI_UNSIGNED_LONG, MPI_SUM, MASTER_NODE, SU2_MPI::GetComm()); + if (rank == MASTER_NODE) { + PrintingToolbox::CTablePrinter outlierTable(&std::cout); + outlierTable.AddColumn("Mean T", 12); + outlierTable.AddColumn("StdDev", 12); + outlierTable.AddColumn("< -" + std::to_string(nSigma), 12); + for (int i = -nSigma; i < nSigma; ++i) { + const int jump = (i == -1); + outlierTable.AddColumn(std::to_string(i) + " to " + std::to_string(i + 1 + jump), 12); + i += jump; + } + outlierTable.AddColumn("> " + std::to_string(nSigma), 12); + outlierTable.AddColumn("#Mitig.", 12); + outlierTable.SetAlign(PrintingToolbox::CTablePrinter::RIGHT); + outlierTable.PrintHeader(); + outlierTable << MeanTemperature << StdDevTemperature; + for (const auto n : global) outlierTable << n; + outlierTable.PrintFooter(); + } + } + END_SU2_OMP_SAFE_GLOBAL_ACCESS + + } else { + SU2_OMP_FOR_STAT(omp_chunk_size) + for (unsigned long iPoint = 0; iPoint < nPoint; ++iPoint) { + (void)DetermineBinAndUpdatePoint(iPoint); + } + END_SU2_OMP_FOR + } } void CEulerSolver::SetPreconditioner(const CConfig *config, unsigned long iPoint, diff --git a/SU2_CFD/src/solvers/CIncEulerSolver.cpp b/SU2_CFD/src/solvers/CIncEulerSolver.cpp index f138422f7ea..2c4c9b05e1c 100644 --- a/SU2_CFD/src/solvers/CIncEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CIncEulerSolver.cpp @@ -2014,12 +2014,6 @@ void CIncEulerSolver::PrepareImplicitIteration(CGeometry *geometry, CSolver**, C PrepareImplicitIteration_impl(precond, geometry, config); } -void CIncEulerSolver::CompleteImplicitIteration(CGeometry *geometry, CSolver**, CConfig *config) { - SU2_ZONE_SCOPED - - CompleteImplicitIteration_impl(geometry, config); -} - void CIncEulerSolver::SetBeta_Parameter(CGeometry *geometry, CSolver **solver_container, CConfig *config, unsigned short iMesh) { SU2_ZONE_SCOPED diff --git a/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp b/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp index 1f3090a67cc..37de9a83207 100644 --- a/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp @@ -936,12 +936,6 @@ void CNEMOEulerSolver::PrepareImplicitIteration(CGeometry *geometry, CSolver**, PrepareImplicitIteration_impl(precond, geometry, config); } -void CNEMOEulerSolver::CompleteImplicitIteration(CGeometry *geometry, CSolver**, CConfig *config) { - SU2_ZONE_SCOPED - - CompleteImplicitIteration_impl(geometry, config); -} - void CNEMOEulerSolver::ComputeUnderRelaxationFactor(const CConfig *config) { SU2_ZONE_SCOPED diff --git a/SU2_CFD/src/solvers/CSolver.cpp b/SU2_CFD/src/solvers/CSolver.cpp index 0efafa2b880..db463454cef 100644 --- a/SU2_CFD/src/solvers/CSolver.cpp +++ b/SU2_CFD/src/solvers/CSolver.cpp @@ -1747,6 +1747,8 @@ void CSolver::AdaptCFLNumber(CGeometry **geometry, CConfig *config) { SU2_ZONE_SCOPED + if (config->GetCFL_Adapt() != YES) return; + /* Adapt the CFL number on all multigrid levels using an exponential progression with under-relaxation approach. */ @@ -1811,7 +1813,10 @@ void CSolver::AdaptCFLNumber(CGeometry **geometry, canIncrease = (linRes < linTol) && (iter >= startingIter); - if ((iMesh == MESH_0) && (Res_Count > 0)) { + /* Do not use the residual flip-flop criteria when we are mitigating outliers + * because the former was never very reliable for large cases where monotonic + * residual reduction is impossible to achieve. */ + if (!config->OptionIsSet("OUTLIER_MITIGATION_PARAM") && iMesh == MESH_0 && Res_Count > 0) { Old_Func = New_Func; if (NonLinRes_Series.empty()) NonLinRes_Series.resize(Res_Count,0.0); diff --git a/SU2_CFD/src/solvers/CTurbSASolver.cpp b/SU2_CFD/src/solvers/CTurbSASolver.cpp index 7bc05787be1..7f9e05975df 100644 --- a/SU2_CFD/src/solvers/CTurbSASolver.cpp +++ b/SU2_CFD/src/solvers/CTurbSASolver.cpp @@ -143,6 +143,10 @@ CTurbSASolver::CTurbSASolver(CGeometry *geometry, CConfig *config, unsigned shor fv1 = Ji_3/(Ji_3+cv1_3); muT_Inf = Density_Inf*fv1*nu_tilde_Inf; + if (config->GetSAParsedOptions().version != SA_OPTIONS::NEG) { + lowerlimit[0] = EPS; + } + /*--- Initialize the solution to the far-field state everywhere. ---*/ nodes = new CTurbSAVariable(nu_tilde_Inf, muT_Inf, nPoint, nDim, nVar, config); diff --git a/SU2_CFD/src/variables/CEulerVariable.cpp b/SU2_CFD/src/variables/CEulerVariable.cpp index 78671f780dc..f45c60060b8 100644 --- a/SU2_CFD/src/variables/CEulerVariable.cpp +++ b/SU2_CFD/src/variables/CEulerVariable.cpp @@ -105,6 +105,8 @@ CEulerVariable::CEulerVariable(su2double density, const su2double *velocity, su2 NIterNewtonsolver.resize(nPoint) = 0; FluidEntropy.resize(nPoint) = su2double(0.0); } + + OutlierMitigation.resize(nPoint) = 0; } bool CEulerVariable::SetPrimVar(unsigned long iPoint, CFluidModel *FluidModel) { diff --git a/config_template.cfg b/config_template.cfg index 6f6d5ee9bcc..891cf8bfb3c 100644 --- a/config_template.cfg +++ b/config_template.cfg @@ -1497,6 +1497,14 @@ CFL_ADAPT= NO % It is reset back to min when linear solvers diverge, or if nonlinear residuals increase too much. CFL_ADAPT_PARAM= ( 0.1, 2.0, 10.0, 1e10, 0.001, 0) % +% For the compressible solver only (EULER, NS, RANS). Detect non-physical solutions statistically instead of +% relying just on zero or negative pressure or temperature, treat all points where temperature is further than N +% standard deviations from the mean as non-physical (MUSCL off and potentially higher dissipation). +% Strategy starts at inner iteration I, mean and sigma update every F, and distribution is printed every P updates. +% OUTLIER_MITIGATION_PARAM= (I, F, P, N). With the default start iteration this is essentially off. +% When starting from freestream, use RAMP_MUSCL and set the start iteration to approximately the end of the MUSCL ramp. +OUTLIER_MITIGATION_PARAM= (999999, 5, 2, 5) +% % Maximum Delta Time in local time stepping simulations MAX_DELTA_TIME= 1E6 % From 4449b7e07252aa2901b917e68259da1976d7cdf7 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Sun, 14 Jun 2026 21:10:03 -0700 Subject: [PATCH 2/2] non simd too --- SU2_CFD/src/solvers/CEulerSolver.cpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/SU2_CFD/src/solvers/CEulerSolver.cpp b/SU2_CFD/src/solvers/CEulerSolver.cpp index 79aa8fe29ca..fb332a765d9 100644 --- a/SU2_CFD/src/solvers/CEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CEulerSolver.cpp @@ -1916,11 +1916,12 @@ void CEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_contain if (van_albada) { lim_i = LimiterHelpers<>::vanAlbadaFunction(Project_Grad_i, V_ij, EPS); lim_j = LimiterHelpers<>::vanAlbadaFunction(Project_Grad_j, V_ij, EPS); - } - else if (limiter) { + } else if (limiter) { lim_i = nodes->GetLimiter_Primitive(iPoint, iVar); lim_j = nodes->GetLimiter_Primitive(jPoint, iVar); } + lim_i *= (1 - static_cast(nodes->OutlierMitigation(iPoint)) / CEulerVariable::MAX_OUTLIER_MITIGATION); + lim_j *= (1 - static_cast(nodes->OutlierMitigation(jPoint)) / CEulerVariable::MAX_OUTLIER_MITIGATION); Primitive_i[iVar] = V_i[iVar] + 0.5 * lim_i * Project_Grad_i; Primitive_j[iVar] = V_j[iVar] - 0.5 * lim_j * Project_Grad_j;