Use variable-relative epsilon for turbulence point limiters

SU2_CFD/include/limiters/computeLimiters.hpp

@@ -50,18 +50,19 @@ void computeLimiters(LIMITER LimiterKind,
                      const GradientType& gradient,
                      FieldType& fieldMin,
                      FieldType& fieldMax,
-                     FieldType& limiter)
+                     FieldType& limiter,
+                     const su2double* refValue = nullptr)
 {
   if (geometry.GetnDim() != 2 && geometry.GetnDim() != 3)
     SU2_MPI::Error("Too many dimensions to compute limiters.", CURRENT_FUNCTION);
 
 #define INSTANTIATE(KIND)\
 if (geometry.GetnDim() == 2) {\
   computeLimiters_impl<2,KIND>(solver, kindMpiComm, kindPeriodicComm1, kindPeriodicComm2, geometry,\
-                               config, varBegin, varEnd, umusclKappa, field, gradient, fieldMin, fieldMax, limiter);\
+                               config, varBegin, varEnd, umusclKappa, field, gradient, fieldMin, fieldMax, limiter, refValue);\
 } else {\
   computeLimiters_impl<3,KIND>(solver, kindMpiComm, kindPeriodicComm1, kindPeriodicComm2, geometry,\
-                               config, varBegin, varEnd, umusclKappa, field, gradient, fieldMin, fieldMax, limiter);\
+                               config, varBegin, varEnd, umusclKappa, field, gradient, fieldMin, fieldMax, limiter, refValue);\
 }
   switch (LimiterKind) {
     case LIMITER::NONE:

SU2_CFD/include/limiters/computeLimiters_impl.hpp

@@ -55,6 +55,10 @@
  * \param[out] fieldMin - Minimum field values over direct neighbors of each point.
  * \param[out] fieldMax - As above but maximum values.
  * \param[out] limiter - Reconstruction limiter for the field.
+ * \param[in] refValue - Optional reference value of each variable. When given, the limiters whose epsilon is an
+ *            absolute number (e.g. Venkatakrishnan) see the differences divided by max(|field|, |refValue|) of the
+ *            point, so that they act in the same way whatever the scale of the variable (dimensionally consistent
+ *            epsilon, see Nishikawa, White and O'Connell, AIAA 2023, following Luke, AIAA 2007-3956).
  *
  * Template parameters:
  * \param nDim - Number of dimensions.
@@ -76,7 +80,8 @@ void computeLimiters_impl(CSolver* solver,
                           const GradientType& gradient,
                           FieldType& fieldMin,
                           FieldType& fieldMax,
-                          FieldType& limiter)
+                          FieldType& limiter,
+                          const su2double* refValue = nullptr)
 {
   constexpr size_t MAXNVAR = 32;
   static_assert(MAXNVAR >= MAX_TRANSPORTED_SPECIES, "MAXNVAR (size of the limiter work arrays) must be at least MAX_TRANSPORTED_SPECIES.");
@@ -234,16 +239,26 @@ void computeLimiters_impl(CSolver* solver,
 
     su2double geoFactor = limiterDetails.geometricFactor(iPoint, geometry);
 
+    /*--- Wang's epsilon is already relative to the range of the field. ---*/
+    constexpr bool absoluteEpsilon = LimiterKind != LIMITER::VENKATAKRISHNAN_WANG;
+
     /*--- Final limiter computation for each variable, get the min limiter
      *    out of the positive/negative projections and deltas. ---*/
 
     for(size_t iVar = varBegin; iVar < varEnd; ++iVar)
     {
-      su2double limMax = limiterDetails.limiterFunction(iVar, projMax[iVar],
-                         fieldMax(iPoint,iVar) - field(iPoint,iVar));
+      su2double scale = 1.0;
+      if (absoluteEpsilon && refValue != nullptr) {
+        AD::SetPreaccIn(refValue[iVar]);
+        scale = max(fabs(field(iPoint,iVar)), fabs(refValue[iVar]));
+        if (scale <= 0.0) scale = 1.0;
+      }
+
+      su2double limMax = limiterDetails.limiterFunction(iVar, projMax[iVar] / scale,
+                         (fieldMax(iPoint,iVar) - field(iPoint,iVar)) / scale);
 
-      su2double limMin = limiterDetails.limiterFunction(iVar, projMin[iVar],
-                         fieldMin(iPoint,iVar) - field(iPoint,iVar));
+      su2double limMin = limiterDetails.limiterFunction(iVar, projMin[iVar] / scale,
+                         (fieldMin(iPoint,iVar) - field(iPoint,iVar)) / scale);
 
       limiter(iPoint,iVar) = geoFactor * min(limMax, limMin);
 

SU2_CFD/include/solvers/CSolver.hpp

@@ -600,6 +600,12 @@ class CSolver {
    */
   void SetSolution_Limiter(CGeometry *geometry, const CConfig *config);
 
+  /*!
+   * \brief Reference value of each variable for the limiters of SetSolution_Limiter, which then use an epsilon
+   *        relative to max(|solution|, |reference|) (see computeLimiters_impl). Null: absolute epsilon.
+   */
+  inline virtual const su2double* GetLimiterReference() const { return nullptr; }
+
   /*!
    * \brief A virtual member.
    * \param[in] geometry - Geometrical definition of the problem.

SU2_CFD/include/solvers/CTurbSolver.hpp

@@ -56,6 +56,12 @@ class CTurbSolver : public CScalarSolver<CTurbVariable> {
    */
   CTurbSolver(CGeometry* geometry, CConfig *config, const CSolver* flow_solver, bool conservative);
 
+  /*!
+   * \brief The turbulence (and transition) variables span orders of magnitude and have units, so their limiters
+   *        use an epsilon relative to the free-stream values.
+   */
+  inline const su2double* GetLimiterReference() const override { return Solution_Inf; }
+
 
   /*!
    * \brief Impose via the residual the Euler wall boundary condition.

SU2_CFD/src/solvers/CSolver.cpp

@@ -2495,7 +2495,8 @@ void CSolver::SetSolution_Limiter(CGeometry *geometry, const CConfig *config) {
   auto& limiter = base_nodes->GetLimiter();
 
   computeLimiters(kindLimiter, this, MPI_QUANTITIES::SOLUTION_LIMITER, PERIODIC_LIM_SOL_1, PERIODIC_LIM_SOL_2,
-                  *geometry, *config, 0, nVar, umusclKappa, solution, gradient, solMin, solMax, limiter);
+                  *geometry, *config, 0, nVar, umusclKappa, solution, gradient, solMin, solMax, limiter,
+                  GetLimiterReference());
 }
 
 void CSolver::GaussElimination(su2double** A, su2double* rhs, unsigned short nVar) {

TestCases/parallel_regression.py

@@ -619,6 +619,24 @@ def main():
     turb_naca0012_sst_muscl_edge.timeout = 3200
     test_list.append(turb_naca0012_sst_muscl_edge)
 
+    # SA point limiting with epsilon relative to the turbulence variable.
+    turb_naca0012_sa_limiter_eps = TestCase('turb_naca0012_sa_limiter_eps')
+    turb_naca0012_sa_limiter_eps.cfg_dir = "rans/naca0012"
+    turb_naca0012_sa_limiter_eps.cfg_file = "turb_NACA0012_sa_limiter_eps.cfg"
+    turb_naca0012_sa_limiter_eps.test_iter = 10
+    turb_naca0012_sa_limiter_eps.test_vals = [-7.852233, -10.813265, 1.080059, 0.017464]
+    turb_naca0012_sa_limiter_eps.timeout = 3200
+    test_list.append(turb_naca0012_sa_limiter_eps)
+
+    # SST point limiting with epsilon relative to the turbulence variable.
+    turb_naca0012_sst_limiter_eps = TestCase('turb_naca0012_sst_limiter_eps')
+    turb_naca0012_sst_limiter_eps.cfg_dir = "rans/naca0012"
+    turb_naca0012_sst_limiter_eps.cfg_file = "turb_NACA0012_sst_limiter_eps.cfg"
+    turb_naca0012_sst_limiter_eps.test_iter = 10
+    turb_naca0012_sst_limiter_eps.test_vals = [-8.162919, -9.358528, -3.590806, 1.071595, 0.015827]
+    turb_naca0012_sst_limiter_eps.timeout = 3200
+    test_list.append(turb_naca0012_sst_limiter_eps)
+
     # Keep applying the stored turbulence limiter after LIMITER_ITER.
     turb_naca0012_sst_muscl_frozen = TestCase('turb_naca0012_sst_muscl_frozen')
     turb_naca0012_sst_muscl_frozen.cfg_dir = "rans/naca0012"

TestCases/rans/naca0012/turb_NACA0012_sa_limiter_eps.cfg

@@ -0,0 +1,107 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%                                                                              %
+% SU2 configuration file                                                       %
+% Case description: NACA0012 SA with variable-relative point limiting
+%                   http://turbmodels.larc.nasa.gov/naca0012_val_sa.html       %
+% Author: Francisco Palacios                                                   %
+% Co-author: rois1995
+% Institution: Stanford University                                             %
+% Date: Feb 18th, 2013                                                         %
+% File Version 8.5.0 "Harrier"                                                 %
+%                                                                              %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+% ------------- DIRECT, ADJOINT, AND LINEARIZED PROBLEM DEFINITION ------------%
+%
+SOLVER= RANS
+KIND_TURB_MODEL= SA
+SA_OPTIONS= NEGATIVE, EXPERIMENTAL
+MATH_PROBLEM= DIRECT
+RESTART_SOL= YES
+
+% -------------------- COMPRESSIBLE FREE-STREAM DEFINITION --------------------%
+%
+MACH_NUMBER= 0.15
+AOA= 10.0
+FREESTREAM_TEMPERATURE= 300.0
+REYNOLDS_NUMBER= 6.0E6
+REYNOLDS_LENGTH= 1.0
+
+% ---------------------- REFERENCE VALUE DEFINITION ---------------------------%
+%
+REF_ORIGIN_MOMENT_X = 0.25
+REF_ORIGIN_MOMENT_Y = 0.00
+REF_ORIGIN_MOMENT_Z = 0.00
+REF_LENGTH= 1.0
+REF_AREA= 1.0
+REF_DIMENSIONALIZATION= FREESTREAM_PRESS_EQ_ONE
+
+% -------------------- BOUNDARY CONDITION DEFINITION --------------------------%
+%
+MARKER_HEATFLUX= ( airfoil, 0.0 )
+MARKER_FAR= ( farfield )
+MARKER_PLOTTING= ( airfoil )
+MARKER_MONITORING= ( airfoil )
+
+% ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------%
+%
+NUM_METHOD_GRAD= WEIGHTED_LEAST_SQUARES
+NUM_METHOD_GRAD_RECON= LEAST_SQUARES
+CFL_NUMBER= 1000.0
+MAX_DELTA_TIME= 1E10
+CFL_ADAPT= NO
+CFL_ADAPT_PARAM= ( 1.5, 0.5, 1.0, 100.0 )
+ITER= 99999
+
+% ----------------------- SLOPE LIMITER DEFINITION ----------------------------%
+VENKAT_LIMITER_COEFF= 0.05
+LIMITER_ITER= 3
+
+% ------------------------ LINEAR SOLVER DEFINITION ---------------------------%
+%
+LINEAR_SOLVER= FGMRES
+LINEAR_SOLVER_PREC= ILU
+LINEAR_SOLVER_ERROR= 1E-10
+LINEAR_SOLVER_ITER= 20
+
+% -------------------- FLOW NUMERICAL METHOD DEFINITION -----------------------%
+%
+CONV_NUM_METHOD_FLOW= ROE
+MUSCL_FLOW= YES
+MUSCL_KAPPA_FLOW= 0.5
+JST_SENSOR_COEFF= ( 0.5, 0.02 )
+TIME_DISCRE_FLOW= EULER_IMPLICIT
+
+% -------------------- TURBULENT NUMERICAL METHOD DEFINITION ------------------%
+%
+CONV_NUM_METHOD_TURB= SCALAR_UPWIND
+MUSCL_TURB= YES
+MUSCL_KAPPA_TURB= 0.0
+SLOPE_LIMITER_TURB= VENKATAKRISHNAN
+TIME_DISCRE_TURB= EULER_IMPLICIT
+CFL_REDUCTION_TURB= 1.0
+
+% --------------------------- CONVERGENCE PARAMETERS --------------------------%
+CONV_RESIDUAL_MINVAL= -12
+CONV_STARTITER= 10
+CONV_CAUCHY_ELEMS= 100
+CONV_CAUCHY_EPS= 1E-6
+
+% ------------------------- INPUT/OUTPUT INFORMATION --------------------------%
+%
+MESH_FILENAME= n0012_225-65.su2
+MESH_FORMAT= SU2
+MESH_OUT_FILENAME= mesh_out
+SOLUTION_FILENAME= solution_flow_sa
+SOLUTION_ADJ_FILENAME= solution_adj
+TABULAR_FORMAT= CSV
+CONV_FILENAME= history
+RESTART_FILENAME= restart_flow
+RESTART_ADJ_FILENAME= restart_adj
+VOLUME_FILENAME= flow
+VOLUME_ADJ_FILENAME= adjoint
+GRAD_OBJFUNC_FILENAME= of_grad
+SURFACE_FILENAME= surface_flow
+SURFACE_ADJ_FILENAME= surface_adjoint
+OUTPUT_WRT_FREQ= 10000
+SCREEN_OUTPUT= (INNER_ITER, RMS_DENSITY, RMS_NU_TILDE, LIFT, DRAG)

TestCases/rans/naca0012/turb_NACA0012_sst_limiter_eps.cfg

@@ -0,0 +1,113 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%                                                                              %
+% SU2 configuration file                                                       %
+% Case description: NACA0012 SST with variable-relative point limiting
+%                   Relative turbulence limiter epsilon                       %
+%                   http://turbmodels.larc.nasa.gov/naca0012_val_sa.html       %
+% Author: Thomas D. Economon                                                   %
+% Co-author: rois1995                                                          %
+% Institution: Stanford University                                             %
+% Date: Feb 18th, 2013                                                         %
+% File Version 8.5.0 "Harrier"                                                 %
+%                                                                              %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+% ------------- DIRECT, ADJOINT, AND LINEARIZED PROBLEM DEFINITION ------------%
+%
+SOLVER= RANS
+KIND_TURB_MODEL= SST
+SST_OPTIONS= V1994m
+MATH_PROBLEM= DIRECT
+RESTART_SOL= YES
+
+% -------------------- COMPRESSIBLE FREE-STREAM DEFINITION --------------------%
+%
+MACH_NUMBER= 0.15
+AOA= 10.0
+FREESTREAM_TEMPERATURE= 300.0
+REYNOLDS_NUMBER= 6.0E6
+REYNOLDS_LENGTH= 1.0
+
+% ---------------------- REFERENCE VALUE DEFINITION ---------------------------%
+%
+REF_ORIGIN_MOMENT_X = 0.25
+REF_ORIGIN_MOMENT_Y = 0.00
+REF_ORIGIN_MOMENT_Z = 0.00
+REF_LENGTH= 1.0
+REF_AREA= 1.0
+REF_DIMENSIONALIZATION= FREESTREAM_PRESS_EQ_ONE
+
+% -------------------- BOUNDARY CONDITION DEFINITION --------------------------%
+%
+MARKER_HEATFLUX= ( airfoil, 0.0 )
+MARKER_FAR= ( farfield )
+MARKER_PLOTTING= ( airfoil )
+MARKER_MONITORING= ( airfoil )
+
+% ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------%
+%
+NUM_METHOD_GRAD= WEIGHTED_LEAST_SQUARES
+NUM_METHOD_GRAD_RECON= LEAST_SQUARES
+CFL_NUMBER= 1000.0
+MAX_DELTA_TIME= 1E10
+CFL_ADAPT= NO
+CFL_ADAPT_PARAM= ( 1.5, 0.5, 1.0, 100.0 )
+ITER= 99999
+
+% ----------------------- SLOPE LIMITER DEFINITION ----------------------------%
+VENKAT_LIMITER_COEFF= 0.05
+LIMITER_ITER= 3
+
+% ------------------------ LINEAR SOLVER DEFINITION ---------------------------%
+%
+LINEAR_SOLVER= FGMRES
+LINEAR_SOLVER_PREC= ILU
+LINEAR_SOLVER_ERROR= 1E-10
+LINEAR_SOLVER_ITER= 20
+
+% -------------------- FLOW NUMERICAL METHOD DEFINITION -----------------------%
+%
+CONV_NUM_METHOD_FLOW= ROE
+MUSCL_FLOW= YES
+MUSCL_KAPPA_FLOW= 0.5
+SLOPE_LIMITER_FLOW= NONE
+JST_SENSOR_COEFF= ( 0.5, 0.02 )
+TIME_DISCRE_FLOW= EULER_IMPLICIT
+
+% -------------------- TURBULENT NUMERICAL METHOD DEFINITION ------------------%
+%
+CONV_NUM_METHOD_TURB= SCALAR_UPWIND
+MUSCL_TURB= YES
+MUSCL_KAPPA_TURB= 0.0
+SLOPE_LIMITER_TURB= VENKATAKRISHNAN
+TIME_DISCRE_TURB= EULER_IMPLICIT
+CFL_REDUCTION_TURB= 1.0
+
+% --------------------------- CONVERGENCE PARAMETERS --------------------------%
+%
+% Convergence field
+CONV_FIELD= RMS_DENSITY
+CONV_RESIDUAL_MINVAL= -12
+CONV_STARTITER= 10
+CONV_CAUCHY_ELEMS= 100
+CONV_CAUCHY_EPS= 1E-6
+
+% ------------------------- INPUT/OUTPUT INFORMATION --------------------------%
+%
+MESH_FILENAME= n0012_225-65.su2
+MESH_FORMAT= SU2
+MESH_OUT_FILENAME= mesh_out
+SOLUTION_FILENAME= solution_flow_sst
+SOLUTION_ADJ_FILENAME= solution_adj
+TABULAR_FORMAT= CSV
+CONV_FILENAME= history
+RESTART_FILENAME= restart_flow
+RESTART_ADJ_FILENAME= restart_adj
+VOLUME_FILENAME= flow
+VOLUME_ADJ_FILENAME= adjoint
+GRAD_OBJFUNC_FILENAME= of_grad
+SURFACE_FILENAME= surface_flow
+SURFACE_ADJ_FILENAME= surface_adjoint
+OUTPUT_WRT_FREQ= 10000
+SCREEN_OUTPUT= (INNER_ITER, RMS_DENSITY, RMS_TKE, RMS_DISSIPATION, LIFT, DRAG)
+OUTPUT_FILES= (RESTART, PARAVIEW, SURFACE_PARAVIEW)