From ffc5bafc7a489f00dbe4be78e02e0eaf874acf1b Mon Sep 17 00:00:00 2001 From: rois1995 Date: Tue, 29 Sep 2026 14:09:19 +0200 Subject: [PATCH 1/4] Clear the whole Jacobian in the non-physical guard of the general-gas Roe scheme Only the diagonal was zeroed (index [iVar][iVar] inside the jVar loop), so the off-diagonal entries kept the values of the previous edge. The guard of CUpwRoeBase_Flow already clears all entries. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01LEL91DW5WPbPwgFtCvHga6 --- SU2_CFD/src/numerics/flow/convection/roe.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/SU2_CFD/src/numerics/flow/convection/roe.cpp b/SU2_CFD/src/numerics/flow/convection/roe.cpp index 73114576d409..371e23059d20 100644 --- a/SU2_CFD/src/numerics/flow/convection/roe.cpp +++ b/SU2_CFD/src/numerics/flow/convection/roe.cpp @@ -811,8 +811,8 @@ CNumerics::ResidualType<> CUpwGeneralRoe_Flow::ComputeResidual(const CConfig* co for (iVar = 0; iVar < nVar; iVar++) { Flux[iVar] = 0.0; for (jVar = 0; jVar < nVar; jVar++) { - Jacobian_i[iVar][iVar] = 0.0; - Jacobian_j[iVar][iVar] = 0.0; + Jacobian_i[iVar][jVar] = 0.0; + Jacobian_j[iVar][jVar] = 0.0; } } AD::SetPreaccOut(Flux, nVar); From a277f59157ab2142ca06a47e5e8a046f270a0da0 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Tue, 29 Sep 2026 14:09:55 +0200 Subject: [PATCH 2/4] Use wave strengths that match the 3D P matrix in L2Roe, LMRoe and the explicit general-gas Roe In 3D the first three columns of the P matrix combine the entropy wave, weighted by n_j, with the shear waves. Their strengths are n_j (Delta rho - Delta p / c^2) + (Delta u x n)_j; the entropy jump was used unweighted in the first one only, so a contact or shear layer on a face not aligned with x got a wrong flux. The standard Roe scheme already projects the jump with P^-1 since #299. Unit test roe_contact_tests.cpp: a contact and a shear layer give the upwind flux in 2D and 3D. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01LEL91DW5WPbPwgFtCvHga6 --- SU2_CFD/src/numerics/flow/convection/roe.cpp | 27 ++-- .../SU2_CFD/numerics/roe_contact_tests.cpp | 125 ++++++++++++++++++ UnitTests/meson.build | 1 + 3 files changed, 144 insertions(+), 9 deletions(-) create mode 100644 UnitTests/SU2_CFD/numerics/roe_contact_tests.cpp diff --git a/SU2_CFD/src/numerics/flow/convection/roe.cpp b/SU2_CFD/src/numerics/flow/convection/roe.cpp index 371e23059d20..37c697515cfb 100644 --- a/SU2_CFD/src/numerics/flow/convection/roe.cpp +++ b/SU2_CFD/src/numerics/flow/convection/roe.cpp @@ -334,9 +334,12 @@ void CUpwL2Roe_Flow::FinalizeResidual(su2double *val_residual, su2double **val_J delta_wave[2] = proj_delta_vel + delta_p/(RoeDensity*RoeSoundSpeed); delta_wave[3] = -proj_delta_vel + delta_p/(RoeDensity*RoeSoundSpeed); } else { - delta_wave[0] = delta_rho - delta_p/RoeSoundSpeed2; - delta_wave[1] = (UnitNormal[0]*delta_vel[2]-UnitNormal[2]*delta_vel[0])*zeta; - delta_wave[2] = (UnitNormal[1]*delta_vel[0]-UnitNormal[0]*delta_vel[1])*zeta; + /*--- The first three columns of the 3D P matrix combine the entropy wave (weighted by the normal) with the + shear waves, so their strengths are n_j (Delta rho - Delta p / c^2) + (Delta u x n)_j. ---*/ + const su2double delta_entropy = delta_rho - delta_p/RoeSoundSpeed2; + delta_wave[0] = UnitNormal[0]*delta_entropy + (delta_vel[1]*UnitNormal[2]-delta_vel[2]*UnitNormal[1])*zeta; + delta_wave[1] = UnitNormal[1]*delta_entropy + (delta_vel[2]*UnitNormal[0]-delta_vel[0]*UnitNormal[2])*zeta; + delta_wave[2] = UnitNormal[2]*delta_entropy + (delta_vel[0]*UnitNormal[1]-delta_vel[1]*UnitNormal[0])*zeta; delta_wave[3] = proj_delta_vel + delta_p/(RoeDensity*RoeSoundSpeed); delta_wave[4] = -proj_delta_vel + delta_p/(RoeDensity*RoeSoundSpeed); } @@ -407,9 +410,12 @@ void CUpwLMRoe_Flow::FinalizeResidual(su2double *val_residual, su2double **val_J delta_wave[2] = proj_delta_vel + delta_p/(RoeDensity*RoeSoundSpeed); delta_wave[3] = -proj_delta_vel + delta_p/(RoeDensity*RoeSoundSpeed); } else { - delta_wave[0] = delta_rho - delta_p/RoeSoundSpeed2; - delta_wave[1] = (UnitNormal[0]*delta_vel[2]-UnitNormal[2]*delta_vel[0]); - delta_wave[2] = (UnitNormal[1]*delta_vel[0]-UnitNormal[0]*delta_vel[1]); + /*--- The first three columns of the 3D P matrix combine the entropy wave (weighted by the normal) with the + shear waves, so their strengths are n_j (Delta rho - Delta p / c^2) + (Delta u x n)_j. ---*/ + const su2double delta_entropy = delta_rho - delta_p/RoeSoundSpeed2; + delta_wave[0] = UnitNormal[0]*delta_entropy + (delta_vel[1]*UnitNormal[2]-delta_vel[2]*UnitNormal[1]); + delta_wave[1] = UnitNormal[1]*delta_entropy + (delta_vel[2]*UnitNormal[0]-delta_vel[0]*UnitNormal[2]); + delta_wave[2] = UnitNormal[2]*delta_entropy + (delta_vel[0]*UnitNormal[1]-delta_vel[1]*UnitNormal[0]); delta_wave[3] = proj_delta_vel + delta_p/(RoeDensity*RoeSoundSpeed); delta_wave[4] = -proj_delta_vel + delta_p/(RoeDensity*RoeSoundSpeed); } @@ -901,9 +907,12 @@ CNumerics::ResidualType<> CUpwGeneralRoe_Flow::ComputeResidual(const CConfig* co delta_wave[2] = proj_delta_vel + delta_p/(RoeDensity*RoeSoundSpeed); delta_wave[3] = -proj_delta_vel + delta_p/(RoeDensity*RoeSoundSpeed); } else { - delta_wave[0] = delta_rho - delta_p/(RoeSoundSpeed*RoeSoundSpeed); - delta_wave[1] = UnitNormal[0]*delta_vel[2]-UnitNormal[2]*delta_vel[0]; - delta_wave[2] = UnitNormal[1]*delta_vel[0]-UnitNormal[0]*delta_vel[1]; + /*--- The first three columns of the 3D P matrix combine the entropy wave (weighted by the normal) with the + shear waves, so their strengths are n_j (Delta rho - Delta p / c^2) + (Delta u x n)_j. ---*/ + const su2double delta_entropy = delta_rho - delta_p/(RoeSoundSpeed*RoeSoundSpeed); + delta_wave[0] = UnitNormal[0]*delta_entropy + delta_vel[1]*UnitNormal[2]-delta_vel[2]*UnitNormal[1]; + delta_wave[1] = UnitNormal[1]*delta_entropy + delta_vel[2]*UnitNormal[0]-delta_vel[0]*UnitNormal[2]; + delta_wave[2] = UnitNormal[2]*delta_entropy + delta_vel[0]*UnitNormal[1]-delta_vel[1]*UnitNormal[0]; delta_wave[3] = proj_delta_vel + delta_p/(RoeDensity*RoeSoundSpeed); delta_wave[4] = -proj_delta_vel + delta_p/(RoeDensity*RoeSoundSpeed); } diff --git a/UnitTests/SU2_CFD/numerics/roe_contact_tests.cpp b/UnitTests/SU2_CFD/numerics/roe_contact_tests.cpp new file mode 100644 index 000000000000..9c99a461f218 --- /dev/null +++ b/UnitTests/SU2_CFD/numerics/roe_contact_tests.cpp @@ -0,0 +1,125 @@ +/*! + * \file roe_contact_tests.cpp + * \brief Roe-type schemes must resolve a contact discontinuity and a shear layer exactly (upwind flux), in 2D + * and 3D, on faces not aligned with the axes. + * \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 +#include +#include +#include "../../../SU2_CFD/include/numerics/flow/convection/roe.hpp" + +namespace { + +constexpr su2double gamma = 1.4; + +std::unique_ptr MakeConfig(bool implicit) { + std::stringstream options; + options << "SOLVER= EULER\nCONV_NUM_METHOD_FLOW= ROE\nENTROPY_FIX_COEFF= 0.0\nROE_KAPPA= 0.5\n"; + options << "TIME_DISCRE_FLOW= " << (implicit ? "EULER_IMPLICIT" : "EULER_EXPLICIT") << "\n"; + return std::make_unique(options, SU2_COMPONENT::SU2_CFD, false); +} + +/*--- Primitive (T, u, v, [w], p, rho, h, c) and secondary (dp/drho_e, dp/de_rho) variables of an ideal gas. ---*/ +void State(unsigned short nDim, su2double rho, const su2double* vel, su2double p, su2double* V, su2double* S) { + su2double q2 = 0.0; + for (unsigned short iDim = 0; iDim < nDim; iDim++) { + V[iDim + 1] = vel[iDim]; + q2 += vel[iDim] * vel[iDim]; + } + V[0] = p / (rho * 287.058); + V[nDim + 1] = p; + V[nDim + 2] = rho; + V[nDim + 3] = gamma / (gamma - 1) * p / rho + 0.5 * q2; + V[nDim + 4] = sqrt(gamma * p / rho); + S[0] = p / rho; + S[1] = (gamma - 1) * rho; +} + +template +std::unique_ptr MakeScheme(unsigned short nDim, const CConfig* config) { + return std::make_unique(nDim, nDim + 2, config); +} +template <> +std::unique_ptr MakeScheme(unsigned short nDim, const CConfig* config) { + return std::make_unique(nDim, nDim + 2, config, false); +} + +/*--- Two states with the same pressure and normal velocity, different density and (with shear) tangential + velocity: a contact discontinuity, whose exact (upwind) flux is the flux of the upstream state. ---*/ +template +void CheckContact(const char* name, unsigned short nDim, const su2double* vel, bool implicit, bool shear) { + auto config = MakeConfig(implicit); + auto numerics = MakeScheme(nDim, config.get()); + + const su2double p = 1.0, rho_i = 1.4, rho_j = 0.7; + const su2double n2[2] = {0.6, 0.8}, n3[3] = {2.0 / 7, 3.0 / 7, 6.0 / 7}; + const su2double* n = (nDim == 2) ? n2 : n3; + const su2double t2[2] = {0.16, -0.12}, t3[3] = {0.15, -0.1, 0.0}; // tangential velocity jumps + const su2double* t = (nDim == 2) ? t2 : t3; + su2double vel_j[3] = {0.0}; + for (unsigned short iDim = 0; iDim < nDim; iDim++) vel_j[iDim] = vel[iDim] + (shear ? t[iDim] : 0.0); + + su2double Vi[10] = {0.0}, Vj[10] = {0.0}, Si[2], Sj[2]; + State(nDim, rho_i, vel, p, Vi, Si); + State(nDim, rho_j, vel_j, p, Vj, Sj); + numerics->SetPrimitive(Vi, Vj); + numerics->SetSecondary(Si, Sj); + numerics->SetNormal(n); + const auto res = numerics->ComputeResidual(config.get()); + + su2double un = 0.0; + for (unsigned short iDim = 0; iDim < nDim; iDim++) un += vel[iDim] * n[iDim]; + const su2double rho_up = (un >= 0.0) ? rho_i : rho_j; + const su2double* V_up = (un >= 0.0) ? Vi : Vj; + + CAPTURE(name, nDim, un, implicit, shear); + CHECK(res.residual[0] == Approx(rho_up * un).margin(1e-12)); + for (unsigned short iDim = 0; iDim < nDim; iDim++) + CHECK(res.residual[iDim + 1] == Approx(rho_up * V_up[iDim + 1] * un + p * n[iDim]).margin(1e-12)); + CHECK(res.residual[nDim + 1] == Approx(rho_up * V_up[nDim + 3] * un).margin(1e-12)); +} + +template +void CheckContacts(const char* name, bool implicit = true, bool shear = true) { + const su2double zero[3] = {0.0, 0.0, 0.0}, moving[3] = {0.3, -0.1, 0.2}, reverse[3] = {-0.3, 0.1, -0.2}; + for (const unsigned short nDim : {2, 3}) { + for (const auto* vel : {zero, moving, reverse}) { + CheckContact(name, nDim, vel, implicit, false); + if (shear) CheckContact(name, nDim, vel, implicit, true); + } + } +} + +} // namespace + +TEST_CASE("Roe schemes resolve a contact and a shear layer exactly", "[Roe]") { + CheckContacts("Roe"); + CheckContacts("general Roe", true); + CheckContacts("general Roe", false); + CheckContacts("LMRoe"); + /*--- L2Roe scales the tangential velocity jump at low Mach number (not upwind by design). ---*/ + CheckContacts("L2Roe", true, false); +} diff --git a/UnitTests/meson.build b/UnitTests/meson.build index c2974c4f11e6..8c2e61614e74 100644 --- a/UnitTests/meson.build +++ b/UnitTests/meson.build @@ -15,6 +15,7 @@ su2_cfd_tests = files(['Common/CConfig_tests.cpp', 'Common/containers/CLookupTable_tests.cpp', 'Common/toolboxes/multilayer_perceptron/CLookUp_ANN_tests.cpp', 'SU2_CFD/numerics/CNumerics_tests.cpp', + 'SU2_CFD/numerics/roe_contact_tests.cpp', 'SU2_CFD/numerics/hllc_jacobian_tests.cpp', 'SU2_CFD/fluid/CFluidModel_tests.cpp', 'SU2_CFD/gradients.cpp', From 07760e64c2fe2407148341615dbfb183551f8419 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Mon, 12 Oct 2026 22:28:45 +0200 Subject: [PATCH 3/4] Use deduced types for Roe entropy wave strengths --- SU2_CFD/src/numerics/flow/convection/roe.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/SU2_CFD/src/numerics/flow/convection/roe.cpp b/SU2_CFD/src/numerics/flow/convection/roe.cpp index 37c697515cfb..3fd23d644584 100644 --- a/SU2_CFD/src/numerics/flow/convection/roe.cpp +++ b/SU2_CFD/src/numerics/flow/convection/roe.cpp @@ -336,7 +336,7 @@ void CUpwL2Roe_Flow::FinalizeResidual(su2double *val_residual, su2double **val_J } else { /*--- The first three columns of the 3D P matrix combine the entropy wave (weighted by the normal) with the shear waves, so their strengths are n_j (Delta rho - Delta p / c^2) + (Delta u x n)_j. ---*/ - const su2double delta_entropy = delta_rho - delta_p/RoeSoundSpeed2; + const auto delta_entropy = delta_rho - delta_p/RoeSoundSpeed2; delta_wave[0] = UnitNormal[0]*delta_entropy + (delta_vel[1]*UnitNormal[2]-delta_vel[2]*UnitNormal[1])*zeta; delta_wave[1] = UnitNormal[1]*delta_entropy + (delta_vel[2]*UnitNormal[0]-delta_vel[0]*UnitNormal[2])*zeta; delta_wave[2] = UnitNormal[2]*delta_entropy + (delta_vel[0]*UnitNormal[1]-delta_vel[1]*UnitNormal[0])*zeta; @@ -412,7 +412,7 @@ void CUpwLMRoe_Flow::FinalizeResidual(su2double *val_residual, su2double **val_J } else { /*--- The first three columns of the 3D P matrix combine the entropy wave (weighted by the normal) with the shear waves, so their strengths are n_j (Delta rho - Delta p / c^2) + (Delta u x n)_j. ---*/ - const su2double delta_entropy = delta_rho - delta_p/RoeSoundSpeed2; + const auto delta_entropy = delta_rho - delta_p/RoeSoundSpeed2; delta_wave[0] = UnitNormal[0]*delta_entropy + (delta_vel[1]*UnitNormal[2]-delta_vel[2]*UnitNormal[1]); delta_wave[1] = UnitNormal[1]*delta_entropy + (delta_vel[2]*UnitNormal[0]-delta_vel[0]*UnitNormal[2]); delta_wave[2] = UnitNormal[2]*delta_entropy + (delta_vel[0]*UnitNormal[1]-delta_vel[1]*UnitNormal[0]); @@ -909,7 +909,7 @@ CNumerics::ResidualType<> CUpwGeneralRoe_Flow::ComputeResidual(const CConfig* co } else { /*--- The first three columns of the 3D P matrix combine the entropy wave (weighted by the normal) with the shear waves, so their strengths are n_j (Delta rho - Delta p / c^2) + (Delta u x n)_j. ---*/ - const su2double delta_entropy = delta_rho - delta_p/(RoeSoundSpeed*RoeSoundSpeed); + const auto delta_entropy = delta_rho - delta_p/(RoeSoundSpeed*RoeSoundSpeed); delta_wave[0] = UnitNormal[0]*delta_entropy + delta_vel[1]*UnitNormal[2]-delta_vel[2]*UnitNormal[1]; delta_wave[1] = UnitNormal[1]*delta_entropy + delta_vel[2]*UnitNormal[0]-delta_vel[0]*UnitNormal[2]; delta_wave[2] = UnitNormal[2]*delta_entropy + delta_vel[0]*UnitNormal[1]-delta_vel[1]*UnitNormal[0]; From 6c571f3d47a221d7270c45e7cb2a4997eb2228fd Mon Sep 17 00:00:00 2001 From: rois1995 Date: Mon, 12 Oct 2026 22:28:44 +0200 Subject: [PATCH 4/4] Separate the HLLC Jacobian option and keep verification access private --- Common/include/CConfig.hpp | 4 +++ Common/src/CConfig.cpp | 4 ++- .../include/numerics/flow/convection/hllc.hpp | 30 +++---------------- SU2_CFD/src/numerics/flow/convection/hllc.cpp | 12 +++----- UnitTests/Common/CConfig_tests.cpp | 16 ++++++++++ .../SU2_CFD/numerics/hllc_jacobian_tests.cpp | 19 ++++++++---- config_template.cfg | 5 +++- 7 files changed, 49 insertions(+), 41 deletions(-) diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp index 9d4920c871b0..183ecde303af 100644 --- a/Common/include/CConfig.hpp +++ b/Common/include/CConfig.hpp @@ -607,6 +607,7 @@ class CConfig { MUSCL_Kappa_Heat, /*!< \brief Blending coefficient for U-MUSCL scheme for the (fvm) heat equation.*/ MUSCL_Kappa_AdjFlow, /*!< \brief Blending coefficient for U-MUSCL scheme for the adj flow equations.*/ MUSCL_Kappa_Species; /*!< \brief Blending coefficient for U-MUSCL scheme for the species equations.*/ + bool Use_Exact_HLLC_Jacobians; /*!< \brief Exact HLLC Jacobians with wave speeds held fixed. */ bool Use_Accurate_Jacobians; /*!< \brief Use numerically computed Jacobians for AUSM+up(2) and SLAU(2). */ bool Use_Accurate_Turb_Jacobians; /*!< \brief Use numerically computed Jacobians for standard SA turbulence model. */ bool EulerPersson; /*!< \brief Boolean to determine whether this is an Euler simulation with Persson shock capturing. */ @@ -4800,6 +4801,9 @@ class CConfig { */ bool GetUse_Accurate_Jacobians(void) const { return Use_Accurate_Jacobians; } + /*!\brief Get whether to use exact HLLC Jacobians with wave speeds held fixed. */ + bool GetUse_Exact_HLLC_Jacobians(void) const { return Use_Exact_HLLC_Jacobians; } + /*! * \brief Get whether to "Use Accurate Jacobians" for Standard SA turbulence model. * \return yes/no. diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index 4cabe516d7c7..550729877fba 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -2176,8 +2176,10 @@ void CConfig::SetConfig_Options() { addDoubleArrayOption("JST_SENSOR_COEFF", 2, false, jst_coeff); /*!\brief LAX_SENSOR_COEFF \n DESCRIPTION: 1st order artificial dissipation coefficients for the Lax-Friedrichs method. \ingroup Config*/ addDoubleOption("LAX_SENSOR_COEFF", Kappa_1st_Flow, 0.15); - /*!\brief USE_ACCURATE_FLUX_JACOBIANS \n DESCRIPTION: Use numerically computed Jacobians for AUSM+up(2) and SLAU(2), exact Jacobians (for fixed wave speeds) for HLLC \ingroup Config*/ + /*!\brief USE_ACCURATE_FLUX_JACOBIANS \n DESCRIPTION: Use numerically computed Jacobians for AUSM+up(2) and SLAU(2) \ingroup Config*/ addBoolOption("USE_ACCURATE_FLUX_JACOBIANS", Use_Accurate_Jacobians, false); + /*!\brief USE_EXACT_HLLC_JACOBIANS \n DESCRIPTION: Exact HLLC Jacobians with wave speeds held fixed \ingroup Config*/ + addBoolOption("USE_EXACT_HLLC_JACOBIANS", Use_Exact_HLLC_Jacobians, false); /*!\brief CENTRAL_JACOBIAN_FIX_FACTOR \n DESCRIPTION: Improve the numerical properties (diagonal dominance) of the global Jacobian matrix, 3 to 4 is "optimum" (central schemes) \ingroup Config*/ addDoubleOption("CENTRAL_JACOBIAN_FIX_FACTOR", Cent_Jac_Fix_Factor, 4.0); /*!\brief CENTRAL_JACOBIAN_FIX_FACTOR \n DESCRIPTION: Control numerical properties of the global Jacobian matrix using a multiplication factor for incompressible central schemes \ingroup Config*/ diff --git a/SU2_CFD/include/numerics/flow/convection/hllc.hpp b/SU2_CFD/include/numerics/flow/convection/hllc.hpp index 29f7ce727a34..8684862fa520 100644 --- a/SU2_CFD/include/numerics/flow/convection/hllc.hpp +++ b/SU2_CFD/include/numerics/flow/convection/hllc.hpp @@ -38,8 +38,9 @@ */ class CUpwHLLC_Flow final : public CNumerics { private: + friend struct HLLCJacobianTestAccess; bool implicit, dynamic_grid; - bool useAccurateJacobian; /*!< \brief Jacobians computed exactly for fixed wave speeds (USE_ACCURATE_FLUX_JACOBIANS). */ + bool useAccurateJacobian; /*!< \brief Jacobians computed exactly for fixed wave speeds (USE_EXACT_HLLC_JACOBIANS). */ unsigned short iDim, jDim, iVar, jVar; su2double *IntermediateState; @@ -81,19 +82,6 @@ class CUpwHLLC_Flow final : public CNumerics { * \return A lightweight const-view (read-only) of the residual/flux and Jacobians. */ ResidualType<> ComputeResidual(const CConfig* config) override; - - /*! - * \brief Fix the left and right wave speeds (relative to the face), for the verification of the Jacobians, which - * treat them as constants. - * \param[in] active - Use the fixed values. - * \param[in] val_sL, val_sR - Left and right wave speeds. - */ - inline void SetFixedWaveSpeeds(bool active, su2double val_sL, su2double val_sR) { - fixedWaveSpeeds = active; - fixed_sL = val_sL; - fixed_sR = val_sR; - } - }; /*! @@ -105,8 +93,9 @@ class CUpwHLLC_Flow final : public CNumerics { */ class CUpwGeneralHLLC_Flow final : public CNumerics { private: + friend struct HLLCJacobianTestAccess; bool implicit, dynamic_grid; - bool useAccurateJacobian; /*!< \brief Jacobians computed exactly for fixed wave speeds (USE_ACCURATE_FLUX_JACOBIANS). */ + bool useAccurateJacobian; /*!< \brief Jacobians computed exactly for fixed wave speeds (USE_EXACT_HLLC_JACOBIANS). */ unsigned short iDim, jDim, iVar, jVar; su2double *IntermediateState; @@ -150,17 +139,6 @@ class CUpwGeneralHLLC_Flow final : public CNumerics { */ ResidualType<> ComputeResidual(const CConfig* config) override; - /*! - * \brief Fix the left and right wave speeds (relative to the face), for the verification of the Jacobians, which - * treat them as constants. - * \param[in] active - Use the fixed values. - * \param[in] val_sL, val_sR - Left and right wave speeds. - */ - inline void SetFixedWaveSpeeds(bool active, su2double val_sL, su2double val_sR) { - fixedWaveSpeeds = active; - fixed_sL = val_sL; - fixed_sR = val_sR; - } /*! * \brief Compute the Average quantities for a general fluid flux between two nodes i and j. diff --git a/SU2_CFD/src/numerics/flow/convection/hllc.cpp b/SU2_CFD/src/numerics/flow/convection/hllc.cpp index 2e315d6a1d91..131ee5b65fa7 100644 --- a/SU2_CFD/src/numerics/flow/convection/hllc.cpp +++ b/SU2_CFD/src/numerics/flow/convection/hllc.cpp @@ -80,7 +80,7 @@ CUpwHLLC_Flow::CUpwHLLC_Flow(unsigned short val_nDim, unsigned short val_nVar, c /* A grid is defined as dynamic if there's rigid grid movement or grid deformation AND the problem is time domain */ dynamic_grid = config->GetDynamic_Grid(); - useAccurateJacobian = config->GetUse_Accurate_Jacobians(); + useAccurateJacobian = config->GetUse_Exact_HLLC_Jacobians(); Gamma = config->GetGamma(); @@ -292,13 +292,11 @@ CNumerics::ResidualType<> CUpwHLLC_Flow::ComputeResidual(const CConfig* config) /*--- Jacobians for the wave speeds held fixed (in the frame of the face). p* = rho_j (u_j - sR)(u_j - sM) + p_j = rho_i (u_i - sL)(u_i - sM) + p_i, so its derivative contains rho_j (sR - u_j) dsM/dU (or rho_i (sL - u_i) dsM/dU). By default the density of the other state is kept, as in the previous implementation (an inexact transcription - of the reference, see #2198). USE_ACCURATE_FLUX_JACOBIANS= YES gives the exact derivative. ---*/ + of the reference, see #2198). USE_EXACT_HLLC_JACOBIANS= YES gives the exact derivative. ---*/ const su2double rhoPStar_R = useAccurateJacobian ? Density_j : Density_i; const su2double rhoPStar_L = useAccurateJacobian ? Density_i : Density_j; - - if (sM > 0.0) { if (sL > 0.0) { @@ -597,7 +595,7 @@ CUpwGeneralHLLC_Flow::CUpwGeneralHLLC_Flow(unsigned short val_nDim, unsigned sho /* A grid is defined as dynamic if there's rigid grid movement or grid deformation AND the problem is time domain */ dynamic_grid = config->GetDynamic_Grid(); - useAccurateJacobian = config->GetUse_Accurate_Jacobians(); + useAccurateJacobian = config->GetUse_Exact_HLLC_Jacobians(); Gamma = config->GetGamma(); @@ -837,13 +835,11 @@ CNumerics::ResidualType<> CUpwGeneralHLLC_Flow::ComputeResidual(const CConfig* c /*--- Jacobians for the wave speeds held fixed (in the frame of the face). p* = rho_j (u_j - sR)(u_j - sM) + p_j = rho_i (u_i - sL)(u_i - sM) + p_i, so its derivative contains rho_j (sR - u_j) dsM/dU (or rho_i (sL - u_i) dsM/dU). By default the density of the other state is kept, as in the previous implementation (an inexact transcription - of the reference, see #2198). USE_ACCURATE_FLUX_JACOBIANS= YES gives the exact derivative. ---*/ + of the reference, see #2198). USE_EXACT_HLLC_JACOBIANS= YES gives the exact derivative. ---*/ const su2double rhoPStar_R = useAccurateJacobian ? Density_j : Density_i; const su2double rhoPStar_L = useAccurateJacobian ? Density_i : Density_j; - - if (sM > 0.0) { if (sL > 0.0) { diff --git a/UnitTests/Common/CConfig_tests.cpp b/UnitTests/Common/CConfig_tests.cpp index 111918785545..a499a075fa1d 100644 --- a/UnitTests/Common/CConfig_tests.cpp +++ b/UnitTests/Common/CConfig_tests.cpp @@ -87,3 +87,19 @@ TEST_CASE("INIT_OPTION_INC defaults", "[Config]") { INIT_OPTION_INC::OPERATING_PRESSURE); CHECK(GetInitOptionInc(ideal_gas_options + "INIT_OPTION_INC= DENSITY_INIT\n") == INIT_OPTION_INC::DENSITY_INIT); } + +TEST_CASE("HLLC and AUSM SLAU Jacobian options are independent", "[Config][HLLCOptions]") { + for (const auto accurate : {false, true}) { + for (const auto exactHLLC : {false, true}) { + std::stringstream options(base_options + "USE_ACCURATE_FLUX_JACOBIANS= " + (accurate ? "YES" : "NO") + + "\nUSE_EXACT_HLLC_JACOBIANS= " + (exactHLLC ? "YES" : "NO") + "\n"); + CConfig config(options, SU2_COMPONENT::SU2_CFD, false); + CHECK(config.GetUse_Accurate_Jacobians() == accurate); + CHECK(config.GetUse_Exact_HLLC_Jacobians() == exactHLLC); + } + } + std::stringstream options(base_options); + CConfig config(options, SU2_COMPONENT::SU2_CFD, false); + CHECK_FALSE(config.GetUse_Accurate_Jacobians()); + CHECK_FALSE(config.GetUse_Exact_HLLC_Jacobians()); +} diff --git a/UnitTests/SU2_CFD/numerics/hllc_jacobian_tests.cpp b/UnitTests/SU2_CFD/numerics/hllc_jacobian_tests.cpp index b03edb3f1414..fc36c6cc9eba 100644 --- a/UnitTests/SU2_CFD/numerics/hllc_jacobian_tests.cpp +++ b/UnitTests/SU2_CFD/numerics/hllc_jacobian_tests.cpp @@ -1,7 +1,7 @@ /*! * \file hllc_jacobian_tests.cpp * \brief Finite-difference verification of the HLLC Jacobians (ideal and general gas, 2D and 3D), on fixed and moving - * faces, exact (USE_ACCURATE_FLUX_JACOBIANS) or default. + * faces, exact (USE_EXACT_HLLC_JACOBIANS) or default. * \version 8.5.0 "Harrier" * * SU2 Project Website: https://su2code.github.io @@ -32,6 +32,15 @@ #include #include "../../../SU2_CFD/include/numerics/flow/convection/hllc.hpp" +struct HLLCJacobianTestAccess { + template + static void SetFixedWaveSpeeds(Numerics& numerics, su2double sL, su2double sR) { + numerics.fixedWaveSpeeds = true; + numerics.fixed_sL = sL; + numerics.fixed_sR = sR; + } +}; + namespace { constexpr su2double gamma = 1.4; @@ -142,12 +151,12 @@ su2double ContactSpeed(unsigned short nDim, const State& Ui, const State& Uj, su std::unique_ptr MakeConfig(bool moving, bool accurate) { std::stringstream options; options << "SOLVER= EULER\nCONV_NUM_METHOD_FLOW= HLLC\nTIME_DISCRE_FLOW= EULER_IMPLICIT\n"; - options << "USE_ACCURATE_FLUX_JACOBIANS= " << (accurate ? "YES" : "NO") << "\n"; + options << "USE_EXACT_HLLC_JACOBIANS= " << (accurate ? "YES" : "NO") << "\n"; if (moving) options << "GRID_MOVEMENT= STEADY_TRANSLATION\nTRANSLATION_RATE= 0.0 0.0 0.0\n"; return std::make_unique(options, SU2_COMPONENT::SU2_CFD, false); } -/*--- With USE_ACCURATE_FLUX_JACOBIANS= YES the Jacobians are exact for fixed wave speeds in all branches. By default +/*--- With USE_EXACT_HLLC_JACOBIANS= YES the Jacobians are exact for fixed wave speeds in all branches. By default the derivative of p* is approximated (star states), the supersonic branches remain exact. ---*/ template void CheckAllBranches(unsigned short nDim, bool moving, su2double k, bool accurate) { @@ -186,7 +195,7 @@ void CheckAllBranches(unsigned short nDim, bool moving, su2double k, bool accura const su2double supersonicSpeeds[2] = {dir > 0 ? 0.02 : -2.0, dir > 0 ? 2.0 : -0.02}; const su2double starSpeeds[2] = {-2.0, 2.0}; for (const auto* s : {supersonicSpeeds, starSpeeds}) { - numerics.SetFixedWaveSpeeds(true, s[0], s[1]); + HLLCJacobianTestAccess::SetFixedWaveSpeeds(numerics, s[0], s[1]); const bool supersonic = (s[0] > 0.0) || (s[1] < 0.0); const su2double sM = ContactSpeed(nDim, Ui, Uj, k, 0.8 * k, n, w, s[0], s[1]); CAPTURE(nDim, moving, k, accurate, dir, s[0], s[1], sM); @@ -235,7 +244,7 @@ TEST_CASE("HLLC supersonic Jacobian uses reconstructed face enthalpy", "[HLLC-fa } numerics.SetNormal(normal); numerics.SetGridVel(gridVelocity, gridVelocity); - numerics.SetFixedWaveSpeeds(true, direction > 0 ? 2.0 : -4.0, direction > 0 ? 4.0 : -2.0); + HLLCJacobianTestAccess::SetFixedWaveSpeeds(numerics, direction > 0 ? 2.0 : -4.0, direction > 0 ? 4.0 : -2.0); const State state = Conservative(nDim, 1.2, velocity, 1.1, 0.0); for (su2double pressureOffset : {0.0, 0.1}) { auto evaluate = [&](const State& Ui, const State& Uj) { diff --git a/config_template.cfg b/config_template.cfg index 0d58bc62727a..8795186ab3d9 100644 --- a/config_template.cfg +++ b/config_template.cfg @@ -1856,9 +1856,12 @@ LOW_MACH_PREC= NO % % Use numerically computed Jacobians for AUSM+up(2) and SLAU(2) % Slower per iteration but potentialy more stable and capable of higher CFL -% For HLLC: exact Jacobians for fixed wave speeds (the default approximation is more robust with inexact linear solves) USE_ACCURATE_FLUX_JACOBIANS= NO % +% Use exact HLLC Jacobians with wave speeds held fixed (NO, YES). +% The default retains the previous approximation; wave-speed derivatives are omitted in both modes. +USE_EXACT_HLLC_JACOBIANS= NO +% % Use the vectorized version of the selected numerical method (available for JST family and Roe). % SU2 should be compiled for an AVX or AVX512 architecture for best performance. % NOTE: Currently vectorization always used for schemes that support it.