diff --git a/SU2_CFD/include/limiters/computeLimiters.hpp b/SU2_CFD/include/limiters/computeLimiters.hpp index b7fe304a4f5..cd2a9ad5ea9 100644 --- a/SU2_CFD/include/limiters/computeLimiters.hpp +++ b/SU2_CFD/include/limiters/computeLimiters.hpp @@ -50,7 +50,8 @@ 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); @@ -58,10 +59,10 @@ void computeLimiters(LIMITER LimiterKind, #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: diff --git a/SU2_CFD/include/limiters/computeLimiters_impl.hpp b/SU2_CFD/include/limiters/computeLimiters_impl.hpp index b53b49018e0..3eb3b94d443 100644 --- a/SU2_CFD/include/limiters/computeLimiters_impl.hpp +++ b/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); diff --git a/SU2_CFD/include/solvers/CSolver.hpp b/SU2_CFD/include/solvers/CSolver.hpp index e5cb9881077..3c0bdc68975 100644 --- a/SU2_CFD/include/solvers/CSolver.hpp +++ b/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. diff --git a/SU2_CFD/include/solvers/CTurbSolver.hpp b/SU2_CFD/include/solvers/CTurbSolver.hpp index 2b978f08b85..600c3fd988b 100644 --- a/SU2_CFD/include/solvers/CTurbSolver.hpp +++ b/SU2_CFD/include/solvers/CTurbSolver.hpp @@ -56,6 +56,12 @@ class CTurbSolver : public CScalarSolver { */ 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. diff --git a/SU2_CFD/src/solvers/CSolver.cpp b/SU2_CFD/src/solvers/CSolver.cpp index 407e1c9bc7a..cd1f0f41ca3 100644 --- a/SU2_CFD/src/solvers/CSolver.cpp +++ b/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) { diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index 97eaace03df..f3846751aee 100755 --- a/TestCases/parallel_regression.py +++ b/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" diff --git a/TestCases/rans/naca0012/turb_NACA0012_sa_limiter_eps.cfg b/TestCases/rans/naca0012/turb_NACA0012_sa_limiter_eps.cfg new file mode 100644 index 00000000000..0fd5e24972e --- /dev/null +++ b/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) diff --git a/TestCases/rans/naca0012/turb_NACA0012_sst_limiter_eps.cfg b/TestCases/rans/naca0012/turb_NACA0012_sst_limiter_eps.cfg new file mode 100644 index 00000000000..df1bc6f7ae8 --- /dev/null +++ b/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) diff --git a/TestCases/rans/scalar_limiter_eps_convergence/README.md b/TestCases/rans/scalar_limiter_eps_convergence/README.md new file mode 100644 index 00000000000..61d82840ba3 --- /dev/null +++ b/TestCases/rans/scalar_limiter_eps_convergence/README.md @@ -0,0 +1,110 @@ +# Turbulence limiter epsilon: convergence comparisons + +These configurations reproduce the full convergence comparisons for PR #2978. +They complement the short `turb_naca0012_{sa,sst}_limiter_eps` MPI regressions +in `parallel_regression.py`, which retain `LIMITER_ITER=3` to exercise frozen +point limiting within ten iterations. These longer cases are manual benchmarks; +they do not add 10,000-iteration jobs to the regression suite. + +The [recorded configurations, compressed full histories and convergence plot](https://github.com/rois1995/SU2/tree/pr-images/scalarLimiterEps) +include the three compared heads and the SST NACA0012 controls. + +## Inputs and reproduction + +Use the existing inputs from the [SU2 TestCases repository](https://github.com/su2code/TestCases): + +| Configuration | Input directory | Required files | Initial condition | +|---|---|---|---| +| `sa_naca.cfg` | `rans/naca0012` | `n0012_225-65.su2`, `solution_flow_sa.dat` | Existing SA restart | +| `sst_flatplate.cfg` | `rans/flatplate` | `mesh_flatplate_turb_137x97.su2` | Free stream | +| `sst_naca.cfg` | `rans/naca0012` | `n0012_225-65.su2` | Free stream; limiter updated throughout | +| `sst_naca_frozen500.cfg` | `rans/naca0012` | `n0012_225-65.su2` | Free stream; limiter frozen at iteration 500 | + +SHA-256 checksums of the inputs used: + +```text +20dfc96875b1d1a8377c2e64248b6ee42a1ecf525c50d822fb5f539d06d68d15 n0012_225-65.su2 +3e665e2d239f71c80c8625235be7e03a8fd5bd9440ef2d930e997291ed412fbc solution_flow_sa.dat +8c720c21c0fc5a84f33565df9c0b0b361e4d3ac67d97fa96ab928783a52c6707 mesh_flatplate_turb_137x97.su2 +``` + +Run each configuration in a separate directory, with its input files linked or +copied there. For example: + +```bash +mkdir -p /tmp/su2-sst-frozen500 +cd /tmp/su2-sst-frozen500 +cp /path/to/SU2/TestCases/rans/scalar_limiter_eps_convergence/sst_naca_frozen500.cfg run.cfg +ln -s /path/to/TestCases/rans/naca0012/n0012_225-65.su2 . +OMP_NUM_THREADS=1 OPENBLAS_NUM_THREADS=1 mpirun -n 2 /path/to/SU2_CFD run.cfg +``` + +For the baseline comparisons, run the same configuration and inputs with each +head's executable. Recorded builds use GCC 9.4.0, release/skylake, MPI and OpenMP +enabled, two MPI ranks, and one OpenMP thread per rank. CGNS and TecIO are +disabled. Compare iteration counts; concurrent runs make their wall times +unsuitable for speedup claims. + +## Results + +Compared numerical source revisions: + +- `develop`: `3aa1a95558f0af7ff43ccc3ec2fe403aff655811` +- PR #2956: `a3968c94950c18314a61e76703022883071cf655` +- PR #2978: `f1b5967a57cba9d1aa1eea9d33d33ca07531d6f6` + +The PR #2956 executable was built at `eef97f4ae5`; its `Common`/`SU2_CFD` +sources are identical to `a3968c9495`. Subsequent changes updated test references +and unit-test attribution. The configurations added here preserve the numerical +settings of the recorded runs. + +Counts include iteration zero. A stop below the **density residual** target is +not a claim that every flow and turbulence residual has converged. + +| Configuration | Density target (log10 RMS) | develop | PR #2956 | PR #2978 | +|---|---:|---:|---:|---:| +| SA NACA0012 | -12 | 272 | 272 | 268 | +| SST flat plate | -12 | Not reached in 10,000 | Not reached in 10,000 | 758 | +| SST NACA0012, limiter updated throughout | -11 | Not reached in 10,000 | 467 | Not reached in 10,000 | +| SST NACA0012, limiter frozen at 500 | -11 | Not run | Not run | 620 | + +The improvement is case dependent: the flat plate reaches its target, SA changes +little, and continuously updating the SST NACA0012 limiter stalls near a density +residual of -10 on the epsilon branch. + +## SST NACA0012 and the current workaround + +Both SST NACA0012 configurations have `MUSCL_TURB=YES`, +`SLOPE_LIMITER_TURB=VENKATAKRISHNAN`, `VENKAT_LIMITER_COEFF=0.03` and fixed +`CFL_NUMBER=1000`. Their only numerical difference is `LIMITER_ITER`: +99999 versus 500. + +Diagnostics point to nonlinear cycling of the turbulence reconstruction in the +wake. In two consecutive stalled steps, the largest change of the k point +limiter is about 0.103 at x=2.284, y=0.063, beside the largest density residuals. +Nine reconstructed wake faces cross the lower k bound. On edge 7122--7123, the +reconstructed k changes from 2.57e-8 to -8.67e-9; the existing positivity guard +then switches that face to its first-order nodal value (about 4.50e-7). + +This switching is an observed contributor to a discontinuous reconstruction, +not a proven complete explanation. The implicit convective Jacobian omits the +gradient/limiter dependence of the reconstructed states, which is a plausible +contributor to the cycle. Lowering CFL to 100 did not reach the target in 10,000 +iterations. A more accurate linear solve and isolated continuous face-clipping +experiments also failed to remove the stalled cycle within their tested budgets; +no general production correction has been validated. + +For this case, the current workaround is to **freeze the stored point limiter**: + +```text +MUSCL_TURB= YES +SLOPE_LIMITER_TURB= VENKATAKRISHNAN +LIMITER_ITER= 500 +``` + +The stored limiter continues to be applied after freezing, and turbulence MUSCL +remains enabled. This run reaches the density target in 620 iterations +(log10 RMS density = -11.0072, CL = 1.070484, CD = 0.0157461). The freeze +iteration is a demonstrated setting for this configuration, rather than a +validated universal choice. `VAN_ALBADA_EDGE` is evaluated at edges and is not +frozen by `LIMITER_ITER`. diff --git a/TestCases/rans/scalar_limiter_eps_convergence/sa_naca.cfg b/TestCases/rans/scalar_limiter_eps_convergence/sa_naca.cfg new file mode 100644 index 00000000000..1223df35a7e --- /dev/null +++ b/TestCases/rans/scalar_limiter_eps_convergence/sa_naca.cfg @@ -0,0 +1,111 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% % +% SU2 configuration file % +% Case description: 2D NACA 0012 Airfoil Validation Case (compressible) % +% 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 ) + +% ----------------------- SLOPE LIMITER DEFINITION ----------------------------% +VENKAT_LIMITER_COEFF= 0.03 + +% ------------------------ 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 +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, LINSOL_ITER, LINSOL_RESIDUAL, LINSOL_ITER_TURB, LINSOL_RESIDUAL_TURB, TOTAL_HEATFLUX) +ITER= 3000 +OUTPUT_FILES= (RESTART, PARAVIEW) +VOLUME_OUTPUT= (COORDINATES, SOLUTION, PRIMITIVE, LIMITER) +SLOPE_LIMITER_FLOW= NONE +HISTORY_OUTPUT= (ITER, RMS_RES, MAX_RES, AERO_COEFF, LINSOL) +SLOPE_LIMITER_TURB= VENKATAKRISHNAN +LIMITER_ITER= 99999 diff --git a/TestCases/rans/scalar_limiter_eps_convergence/sst_flatplate.cfg b/TestCases/rans/scalar_limiter_eps_convergence/sst_flatplate.cfg new file mode 100644 index 00000000000..f05b98591a4 --- /dev/null +++ b/TestCases/rans/scalar_limiter_eps_convergence/sst_flatplate.cfg @@ -0,0 +1,115 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% % +% SU2 configuration file % +% Case description: Turbulent flow over flat plate with zero pressure gradient % +% Author: Thomas D. Economon % +% Co-author: rois1995 +% Institution: Stanford University % +% Date: 2011.11.10 % +% File Version 8.5.0 "Harrier" % +% % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% ------------- DIRECT, ADJOINT, AND LINEARIZED PROBLEM DEFINITION ------------% +% +SOLVER= RANS +KIND_TURB_MODEL= SST +MATH_PROBLEM= DIRECT +RESTART_SOL= NO + +% ----------- COMPRESSIBLE AND INCOMPRESSIBLE FREE-STREAM DEFINITION ----------% +% +MACH_NUMBER= 0.2 +AOA= 0.0 +SIDESLIP_ANGLE= 0.0 +FREESTREAM_TEMPERATURE= 300.0 +REYNOLDS_NUMBER= 5000000.0 +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= 2.0 + +% -------------------- BOUNDARY CONDITION DEFINITION --------------------------% +% +MARKER_HEATFLUX= ( wall, 0.0 ) +MARKER_INLET= ( inlet, 302.4, 118309.784, 1.0, 0.0, 0.0 ) +MARKER_OUTLET= ( outlet, 115056.0, farfield, 115056.0 ) +MARKER_SYM= ( symmetry ) +MARKER_PLOTTING= ( wall ) +MARKER_MONITORING= ( wall ) + +% ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------% +% +NUM_METHOD_GRAD= WEIGHTED_LEAST_SQUARES +CFL_NUMBER= 1000.0 +CFL_ADAPT= NO +CFL_ADAPT_PARAM= ( 1.5, 0.5, 1.0, 100.0 ) +RK_ALPHA_COEFF= ( 0.66667, 0.66667, 1.000000 ) + +% ----------------------- SLOPE LIMITER DEFINITION ----------------------------% +% +VENKAT_LIMITER_COEFF= 0.1 +ADJ_SHARP_LIMITER_COEFF= 3.0 +REF_SHARP_EDGES= 3.0 +SENS_REMOVE_SHARP= NO + +% -------------------------- MULTIGRID PARAMETERS -----------------------------% +% +MGLEVEL= 3 +MGCYCLE= V_CYCLE +MG_PRE_SMOOTH= ( 2, 2, 2, 2 ) +MG_POST_SMOOTH= ( 2, 2, 2, 2) +MG_CORRECTION_SMOOTH= ( 1, 1, 1, 1 ) +MG_DAMP_RESTRICTION= 0.5 +MG_DAMP_PROLONGATION= 0.5 +MG_IMPLICIT_LINES= YES + +% -------------------- FLOW NUMERICAL METHOD DEFINITION -----------------------% +% +CONV_NUM_METHOD_FLOW= ROE +MUSCL_FLOW= YES +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 +TIME_DISCRE_TURB= EULER_IMPLICIT + +% --------------------------- CONVERGENCE PARAMETERS --------------------------% +CONV_RESIDUAL_MINVAL= -12 +CONV_STARTITER= 10 +CONV_CAUCHY_ELEMS= 100 +CONV_CAUCHY_EPS= 1E-6 + +% ------------------------- INPUT/OUTPUT INFORMATION --------------------------% +% +MESH_FILENAME= mesh_flatplate_turb_137x97.su2 +MESH_FORMAT= SU2 +MESH_OUT_FILENAME= mesh_out +SOLUTION_FILENAME= solution_flow +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= 1000 +ITER= 10000 +OUTPUT_FILES= (RESTART, PARAVIEW, SURFACE_PARAVIEW) +VOLUME_OUTPUT= (COORDINATES, SOLUTION, PRIMITIVE, LIMITER) +HISTORY_OUTPUT= (ITER, RMS_RES, MAX_RES, AERO_COEFF, LINSOL) +FREESTREAM_TURBULENCEINTENSITY= 0.0003873 +FREESTREAM_TURB2LAMVISCRATIO= 0.009 +MUSCL_TURB= YES +SLOPE_LIMITER_TURB= VENKATAKRISHNAN diff --git a/TestCases/rans/scalar_limiter_eps_convergence/sst_naca.cfg b/TestCases/rans/scalar_limiter_eps_convergence/sst_naca.cfg new file mode 100644 index 00000000000..540d41eef52 --- /dev/null +++ b/TestCases/rans/scalar_limiter_eps_convergence/sst_naca.cfg @@ -0,0 +1,113 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% % +% SU2 configuration file % +% Case description: 2D NACA 0012 Airfoil Validation Case (compressible) % +% 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 + +% -------------------- 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 ) + +% ----------------------- SLOPE LIMITER DEFINITION ----------------------------% +VENKAT_LIMITER_COEFF= 0.03 +LIMITER_ITER= 99999 + +% ------------------------ 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_KAPPA_TURB= 0.0 +TIME_DISCRE_TURB= EULER_IMPLICIT +CFL_REDUCTION_TURB= 1.0 + +% --------------------------- CONVERGENCE PARAMETERS --------------------------% +% +% Convergence field +CONV_FIELD= RMS_DENSITY +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, LINSOL_RESIDUAL, TOTAL_HEATFLUX) +OUTPUT_FILES= (RESTART, PARAVIEW, SURFACE_PARAVIEW) +ITER= 10000 +HISTORY_OUTPUT= (ITER, RMS_RES, MAX_RES, AERO_COEFF, LINSOL) +RESTART_SOL= NO +CONV_RESIDUAL_MINVAL= -11 +MUSCL_TURB= YES +SLOPE_LIMITER_TURB= VENKATAKRISHNAN diff --git a/TestCases/rans/scalar_limiter_eps_convergence/sst_naca_frozen500.cfg b/TestCases/rans/scalar_limiter_eps_convergence/sst_naca_frozen500.cfg new file mode 100644 index 00000000000..24f70e6238c --- /dev/null +++ b/TestCases/rans/scalar_limiter_eps_convergence/sst_naca_frozen500.cfg @@ -0,0 +1,114 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% % +% SU2 configuration file % +% Case description: 2D NACA 0012 Airfoil Validation Case (compressible) % +% 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 + +% -------------------- 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 ) + +% ----------------------- SLOPE LIMITER DEFINITION ----------------------------% +VENKAT_LIMITER_COEFF= 0.03 +% Freeze the stored point limiter to suppress the observed wake cycling. +LIMITER_ITER= 500 + +% ------------------------ 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_KAPPA_TURB= 0.0 +TIME_DISCRE_TURB= EULER_IMPLICIT +CFL_REDUCTION_TURB= 1.0 + +% --------------------------- CONVERGENCE PARAMETERS --------------------------% +% +% Convergence field +CONV_FIELD= RMS_DENSITY +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, LINSOL_RESIDUAL, TOTAL_HEATFLUX) +OUTPUT_FILES= (RESTART, PARAVIEW, SURFACE_PARAVIEW) +ITER= 10000 +HISTORY_OUTPUT= (ITER, RMS_RES, MAX_RES, AERO_COEFF, LINSOL) +RESTART_SOL= NO +CONV_RESIDUAL_MINVAL= -11 +MUSCL_TURB= YES +SLOPE_LIMITER_TURB= VENKATAKRISHNAN diff --git a/UnitTests/SU2_CFD/scalar_limiter_eps.cpp b/UnitTests/SU2_CFD/scalar_limiter_eps.cpp new file mode 100644 index 00000000000..402e55a83fb --- /dev/null +++ b/UnitTests/SU2_CFD/scalar_limiter_eps.cpp @@ -0,0 +1,171 @@ +/*! + * \file scalar_limiter_eps.cpp + * \brief Tests for turbulence limiter units, zero references and derivatives. + * \author rois1995 + * \version 8.5.0 "Harrier" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2026, SU2 Contributors (cf. AUTHORS.md) + * + * SU2 is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * SU2 is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with SU2. If not, see . + */ + +#include "catch.hpp" +#include "../UnitQuadTestCase.hpp" +#include "../../SU2_CFD/include/solvers/CTurbSolver.hpp" +#include + +namespace { +struct ScalarLimiterTestCase : UnitQuadTestCase { + ScalarLimiterTestCase(const std::string& model, const std::string& limiter, unsigned short dimension) { + config_options = + "SOLVER= RANS\nKIND_TURB_MODEL= " + model + "\nSLOPE_LIMITER_TURB= " + limiter + + "\nMUSCL_TURB= YES\nVENKAT_LIMITER_COEFF= 0.2\nMACH_NUMBER= 0.5\nREYNOLDS_NUMBER= 1e6\nMESH_FORMAT= " + + (dimension == 2 ? "RECTANGLE" : "BOX") + + "\nMARKER_HEATFLUX= (y_minus,0,y_plus,0)\nMESH_BOX_SIZE= 5,5,5\nMARKER_FAR= (x_minus,x_plus" + + (dimension == 2 ? ")\n" : ",z_minus,z_plus)\n"); + InitConfig(); + InitGeometry(); + InitSolver(); + config->SetGlobalParam(config->GetKind_Solver(), RUNTIME_TURB_SYS); + } + ~ScalarLimiterTestCase() { + delete solver[TURB_SOL]; + solver[TURB_SOL] = nullptr; + } + CTurbSolver& Turb() { return *static_cast(solver[TURB_SOL]); } + su2double Evaluate(const su2double& amplitude, const su2double& reference) { + auto* nodes = Turb().GetNodes(); + auto& gradient = nodes->GetGradient_Reconstruction(); + auto* ref = const_cast(Turb().GetLimiterReference()); + for (auto var = 0u; var < Turb().GetnVar(); ++var) ref[var] = reference * (var + 1.0); + for (auto point = 0ul; point < geometry->GetnPoint(); ++point) { + const auto x = geometry->nodes->GetCoord(point, 0); + for (auto var = 0u; var < Turb().GetnVar(); ++var) { + nodes->GetSolution(point)[var] = amplitude * (1.0 + var + x * x); + for (auto dim = 0u; dim < geometry->GetnDim(); ++dim) { + su2double slope = 0.0; + if (dim == 0) slope = 4.0 + 2.0 * x; + gradient(point, var, dim) = amplitude * slope; + } + } + } + SU2_OMP_PARALLEL { Turb().SetSolution_Limiter(geometry.get(), config.get()); } + su2double result = 0.0; + for (auto point = 0ul; point < geometry->GetnPointDomain(); ++point) + for (auto var = 0u; var < Turb().GetnVar(); ++var) result += nodes->GetLimiter()(point, var); + return result; + } + std::vector Snapshot() { + std::vector values; + for (auto point = 0ul; point < geometry->GetnPointDomain(); ++point) + for (auto var = 0u; var < Turb().GetnVar(); ++var) + values.push_back(SU2_TYPE::GetValue(Turb().GetNodes()->GetLimiter()(point, var))); + return values; + } +}; +const std::vector scaledLimiters = {"VENKATAKRISHNAN", "NISHIKAWA_R3", "NISHIKAWA_R4", "NISHIKAWA_R5"}; +} // namespace + +TEST_CASE("Turbulence point limiters preserve their values across field units", "[ScalarLimiterEps]") { + for (const std::string model : {"SA", "SST"}) { + for (const unsigned short dim : {2, 3}) { + for (const auto& kind : scaledLimiters) { + ScalarLimiterTestCase test(model, kind, dim); + CHECK(test.solver[FLOW_SOL]->GetLimiterReference() == nullptr); + test.Evaluate(1.0, 10.0); + const auto reference = test.Snapshot(); + REQUIRE_FALSE(reference.empty()); + for (const double units : {1e-12, 1e-6, 1.0, 1e6, 1e12}) { + for (const double sign : {-1.0, 1.0}) { + CAPTURE(model, dim, kind, units, sign); + test.Evaluate(units, sign * 10.0 * units); + const auto values = test.Snapshot(); + REQUIRE(values.size() == reference.size()); + for (size_t i = 0; i < values.size(); ++i) CHECK(values[i] == Approx(reference[i]).margin(1e-12)); + } + } + test.Evaluate(1.0, 0.0); + const auto zeroReference = test.Snapshot(); + test.Evaluate(1.0, 1e-300); + const auto tinyReference = test.Snapshot(); + for (size_t i = 0; i < zeroReference.size(); ++i) { + CHECK(std::isfinite(zeroReference[i])); + CHECK(tinyReference[i] == Approx(zeroReference[i]).margin(1e-12)); + } + test.Evaluate(0.0, 0.0); + for (const auto value : test.Snapshot()) CHECK(value == Approx(1.0)); + } + } + } +} + +TEST_CASE("Wang's range-relative turbulence limiter ignores the reference magnitude", "[ScalarLimiterEps]") { + for (const unsigned short dim : {2, 3}) { + ScalarLimiterTestCase test("SST", "VENKATAKRISHNAN_WANG", dim); + test.Evaluate(1.0, 10.0); + const auto reference = test.Snapshot(); + for (const double value : {0.0, -10.0, 1e-300, 1e12}) { + test.Evaluate(1.0, value); + const auto values = test.Snapshot(); + for (size_t i = 0; i < values.size(); ++i) CHECK(values[i] == Approx(reference[i]).margin(1e-12)); + } + } +} + +#if defined(CODI_REVERSE_TYPE) || defined(CODI_FORWARD_TYPE) +TEST_CASE("Scalar limiter reference derivatives match finite differences", "[ScalarLimiterEps][AD]") { + for (const auto& kind : scaledLimiters) { + AD::Reset(); + { + ScalarLimiterTestCase test("SST", kind, 2); + su2double amplitude = 1.0, reference = 10.0; + double dAmplitude, dReference; +#ifdef CODI_REVERSE_TYPE + AD::StartRecording(); + AD::RegisterInput(amplitude); + AD::RegisterInput(reference); + su2double result = test.Evaluate(amplitude, reference); + AD::RegisterOutput(result); + AD::StopRecording(); + SU2_TYPE::SetDerivative(result, 1.0); + AD::ComputeAdjoint(); + dAmplitude = SU2_TYPE::GetDerivative(amplitude); + dReference = SU2_TYPE::GetDerivative(reference); +#else + SU2_TYPE::SetDerivative(amplitude, 1.0); + SU2_TYPE::SetDerivative(reference, 0.0); + dAmplitude = SU2_TYPE::GetDerivative(test.Evaluate(amplitude, reference)); + SU2_TYPE::SetDerivative(amplitude, 0.0); + SU2_TYPE::SetDerivative(reference, 1.0); + dReference = SU2_TYPE::GetDerivative(test.Evaluate(amplitude, reference)); +#endif + const double h = 1e-4; + const double plusAmplitude = SU2_TYPE::GetValue(test.Evaluate(su2double(1.0 + h), su2double(10.0))); + const double minusAmplitude = SU2_TYPE::GetValue(test.Evaluate(su2double(1.0 - h), su2double(10.0))); + const double plusReference = SU2_TYPE::GetValue(test.Evaluate(su2double(1.0), su2double(10.0 + h))); + const double minusReference = SU2_TYPE::GetValue(test.Evaluate(su2double(1.0), su2double(10.0 - h))); + CHECK(dAmplitude == Approx((plusAmplitude - minusAmplitude) / (2 * h)).epsilon(1e-5).margin(1e-7)); + CHECK(dReference == Approx((plusReference - minusReference) / (2 * h)).epsilon(1e-5).margin(1e-7)); + CHECK(std::abs(dReference) > 1e-8); + CHECK(dAmplitude + 10.0 * dReference == Approx(0.0).margin(1e-8)); + } + AD::Reset(); + } +} +#endif diff --git a/UnitTests/meson.build b/UnitTests/meson.build index 1883237365e..baf6dfa7ecd 100644 --- a/UnitTests/meson.build +++ b/UnitTests/meson.build @@ -22,6 +22,7 @@ su2_cfd_tests = files(['Common/CConfig_tests.cpp', 'SU2_CFD/gradients.cpp', 'SU2_CFD/periodic_limiters.cpp', 'SU2_CFD/scalar_muscl.cpp', + 'SU2_CFD/scalar_limiter_eps.cpp', 'SU2_CFD/solvers/CNEMOEulerSolver_tests.cpp', 'SU2_CFD/nemo_viscous_assembly.cpp', 'SU2_CFD/windowing.cpp', @@ -31,12 +32,14 @@ su2_cfd_tests = files(['Common/CConfig_tests.cpp', # Reverse-mode (algorithmic differentiation) tests: su2_cfd_tests_ad = files(['Common/simple_ad_test.cpp', + 'SU2_CFD/scalar_limiter_eps.cpp', 'SU2_CFD/solvers/under_relaxation_ad_tests.cpp', 'SU2_CFD/fluid/CFluidModel_tests_AD.cpp', 'Common/toolboxes/multilayer_perceptron/MLP_Jacobian_tests.cpp']) # Forward-mode (direct differentiation) tests: su2_cfd_tests_dd = files(['Common/simple_directdiff_test.cpp', + 'SU2_CFD/scalar_limiter_eps.cpp', 'SU2_CFD/solvers/under_relaxation_ad_tests.cpp']) # -------------------------------------------------------------------------