Detect non-physical solutions in a statistic way instead of just P or T <= 0

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.

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);

SU2_CFD/include/numerics_simd/flow/convection/common.hpp

@@ -188,6 +188,7 @@
  * \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<class ReconVarType, class PrimVarType, size_t nDim, class VariableType>
@@ -201,7 +202,8 @@
                                                       const LIMITER limiterType,
                                                       const CPair<PrimVarType>& V1st,
                                                       const VectorDbl<nDim>& 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 @@
     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;
 }

SU2_CFD/include/numerics_simd/flow/convection/upwind.hpp

@@ -118,8 +118,10 @@ class CUpwindBase : public Base {
     V1st.j.all = gatherVariables<nPrimVar>(jPoint, solution.GetPrimitive());
 
     /*--- Recompute density and enthalpy instead of reconstructing. ---*/
+    Double nonPhysical;
     auto V = reconstructPrimitives<CCompressiblePrimitives<nDim,nPrimVarGrad> >(
-        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<const Derived*>(this);
     VectorDbl<nVar> flux;
     MatrixDbl<nVar> 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<CRoeScheme<Decorator>, Decorator> {
                                 const CPair<ConsVarType>& 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<CRoeScheme<Decorator>, 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<CMSWScheme<Decorator>, Decorator> {
                                 const CPair<ConsVarType>& 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<CMSWScheme<Decorator>, 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<CCompressiblePrimitives<nDim, nPrimVarGrad>> Vweighted;

SU2_CFD/include/solvers/CEulerSolver.hpp

@@ -116,6 +116,9 @@ class CEulerSolver : public CFVMFlowSolverBase<CEulerVariable, ENUM_REGIME::COMP
 
   vector<CFluidModel*> FluidModel;   /*!< \brief fluid model used in the solver. */
 
+  /*!< \brief Variables for outlier detection. */
+  su2double MeanTemperature, StdDevTemperature;
+
   /*--- Turbomachinery Solver Variables ---*/
 
   vector<su2activematrix> AverageFlux;
@@ -782,11 +785,17 @@ class CEulerSolver : public CFVMFlowSolverBase<CEulerVariable, ENUM_REGIME::COMP
   void PrepareImplicitIteration(CGeometry *geometry, CSolver**, CConfig *config) final;
 
   /*!
-   * \brief Complete an implicit iteration.
-   * \param[in] geometry - Geometrical definition of the problem.
+   * \brief Compute a suitable under-relaxation parameter to limit the change in the solution variables over
+   * a nonlinear iteration for stability.
    * \param[in] config - Definition of the particular problem.
    */
-  void CompleteImplicitIteration(CGeometry *geometry, CSolver**, CConfig *config) final;
+  void ComputeUnderRelaxationFactor(const CConfig *config) final;
+
+  /*!
+   * \brief Identify points where the static temperature is an outlier to treat them
+   *        as non-physical, e.g. force the reconstruction to be 1st order.
+   */
+  void IdentifySolutionOutliers(const CConfig *config, unsigned long iter) final;
 
   /*!
    * \brief Provide the mass flow rate.
@@ -1127,20 +1136,21 @@ class CEulerSolver : public CFVMFlowSolverBase<CEulerVariable, ENUM_REGIME::COMP
                                    unsigned short marker_flag,
                                    TURBOMACHINERY_TYPE kind_turb){
 
-    if ((kind_turb == TURBOMACHINERY_TYPE::AXIAL && nDim == 3) || (kind_turb == TURBOMACHINERY_TYPE::CENTRIPETAL_AXIAL && marker_flag == OUTFLOW) || (kind_turb == TURBOMACHINERY_TYPE::AXIAL_CENTRIFUGAL && marker_flag == INFLOW) ){
-      turboVelocity[2] =  turboNormal[0]*cartesianVelocity[0] + cartesianVelocity[1]*turboNormal[1];
-      turboVelocity[1] =  turboNormal[0]*cartesianVelocity[1] - turboNormal[1]*cartesianVelocity[0];
+    if ((kind_turb == TURBOMACHINERY_TYPE::AXIAL && nDim == 3) ||
+        (kind_turb == TURBOMACHINERY_TYPE::CENTRIPETAL_AXIAL && marker_flag == OUTFLOW) ||
+        (kind_turb == TURBOMACHINERY_TYPE::AXIAL_CENTRIFUGAL && marker_flag == INFLOW)) {
+      turboVelocity[2] = turboNormal[0]*cartesianVelocity[0] + cartesianVelocity[1]*turboNormal[1];
+      turboVelocity[1] = turboNormal[0]*cartesianVelocity[1] - turboNormal[1]*cartesianVelocity[0];
       turboVelocity[0] = cartesianVelocity[2];
-    }
-    else{
-      turboVelocity[0] =  turboNormal[0]*cartesianVelocity[0] + cartesianVelocity[1]*turboNormal[1];
-      turboVelocity[1] =  turboNormal[0]*cartesianVelocity[1] - turboNormal[1]*cartesianVelocity[0];
-      if (marker_flag == INFLOW){
+    } else {
+      turboVelocity[0] = turboNormal[0]*cartesianVelocity[0] + cartesianVelocity[1]*turboNormal[1];
+      turboVelocity[1] = turboNormal[0]*cartesianVelocity[1] - turboNormal[1]*cartesianVelocity[0];
+
+      if (marker_flag == INFLOW) {
         turboVelocity[0] *= -1.0;
         turboVelocity[1] *= -1.0;
       }
-      if(nDim == 3)
-        turboVelocity[2] = cartesianVelocity[2];
+      if (nDim == 3) turboVelocity[2] = cartesianVelocity[2];
     }
   }
 
@@ -1156,22 +1166,21 @@ class CEulerSolver : public CFVMFlowSolverBase<CEulerVariable, ENUM_REGIME::COMP
                                   unsigned short marker_flag,
                                   TURBOMACHINERY_TYPE kind_turb){
 
-    if ((kind_turb == TURBOMACHINERY_TYPE::AXIAL && nDim == 3) || (kind_turb == TURBOMACHINERY_TYPE::CENTRIPETAL_AXIAL && marker_flag == OUTFLOW) || (kind_turb == TURBOMACHINERY_TYPE::AXIAL_CENTRIFUGAL && marker_flag == INFLOW)){
+    if ((kind_turb == TURBOMACHINERY_TYPE::AXIAL && nDim == 3) ||
+        (kind_turb == TURBOMACHINERY_TYPE::CENTRIPETAL_AXIAL && marker_flag == OUTFLOW) ||
+        (kind_turb == TURBOMACHINERY_TYPE::AXIAL_CENTRIFUGAL && marker_flag == INFLOW)) {
       cartesianVelocity[0] = turboVelocity[2]*turboNormal[0] - turboVelocity[1]*turboNormal[1];
       cartesianVelocity[1] = turboVelocity[2]*turboNormal[1] + turboVelocity[1]*turboNormal[0];
       cartesianVelocity[2] = turboVelocity[0];
-    }
-    else{
-      cartesianVelocity[0] =  turboVelocity[0]*turboNormal[0] - turboVelocity[1]*turboNormal[1];
-      cartesianVelocity[1] =  turboVelocity[0]*turboNormal[1] + turboVelocity[1]*turboNormal[0];
+    } else {
+      cartesianVelocity[0] = turboVelocity[0]*turboNormal[0] - turboVelocity[1]*turboNormal[1];
+      cartesianVelocity[1] = turboVelocity[0]*turboNormal[1] + turboVelocity[1]*turboNormal[0];
 
-      if (marker_flag == INFLOW){
+      if (marker_flag == INFLOW) {
         cartesianVelocity[0] *= -1.0;
         cartesianVelocity[1] *= -1.0;
       }
-
-      if(nDim == 3)
-        cartesianVelocity[2] = turboVelocity[2];
+      if (nDim == 3) cartesianVelocity[2] = turboVelocity[2];
     }
   }
 

SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp

@@ -345,7 +345,9 @@ class CFVMFlowSolverBase : public CSolver {
    * \brief Compute a suitable under-relaxation parameter to limit the change in the solution variables over a nonlinear
    * iteration for stability.
    */
-  virtual void ComputeUnderRelaxationFactor(const CConfig* config);
+  virtual void ComputeUnderRelaxationFactor(const CConfig* config) {
+    SU2_MPI::Error("Not implemented for this solver.", CURRENT_FUNCTION);
+  }
 
   /*!
    * \brief General implementation to load a flow solution from a restart file.
@@ -976,39 +978,6 @@ class CFVMFlowSolverBase : public CSolver {
 
   }
 
-  /*!
-   * \brief Generic implementation to complete an implicit iteration, i.e. update the solution.
-   * \tparam compute_ur - Whether to use automatic under-relaxation for the update.
-   */
-  template<bool compute_ur>
-  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.

SU2_CFD/include/solvers/CFVMFlowSolverBase.inl

@@ -596,42 +596,33 @@ void CFVMFlowSolverBase<V, R>::ComputeVerificationError(CGeometry* geometry, CCo
 }
 
 template <class V, ENUM_REGIME R>
-void CFVMFlowSolverBase<V, R>::ComputeUnderRelaxationFactor(const CConfig* config) {
+void CFVMFlowSolverBase<V, R>::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 <class V, ENUM_REGIME R>

SU2_CFD/include/solvers/CIncEulerSolver.hpp

@@ -365,13 +365,6 @@ class CIncEulerSolver : public CFVMFlowSolverBase<CIncEulerVariable, ENUM_REGIME
    */
   void PrepareImplicitIteration(CGeometry *geometry, CSolver**, 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.

SU2_CFD/include/solvers/CNEMOEulerSolver.hpp

@@ -329,13 +329,6 @@ class CNEMOEulerSolver : public CFVMFlowSolverBase<CNEMOEulerVariable, ENUM_REGI
    */
   void PrepareImplicitIteration(CGeometry *geometry, CSolver**, 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 Print verification error to screen.
    * \param[in] config - Definition of the particular problem.