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 1ab7ced2d354..550729877fba 100644
--- a/Common/src/CConfig.cpp
+++ b/Common/src/CConfig.cpp
@@ -2178,6 +2178,8 @@ void CConfig::SetConfig_Options() {
   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) \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 b525df4bf9cf..8684862fa520 100644
--- a/SU2_CFD/include/numerics/flow/convection/hllc.hpp
+++ b/SU2_CFD/include/numerics/flow/convection/hllc.hpp
@@ -38,7 +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_EXACT_HLLC_JACOBIANS). */
   unsigned short iDim, jDim, iVar, jVar;
 
   su2double *IntermediateState;
@@ -55,6 +57,8 @@ class CUpwHLLC_Flow final : public CNumerics {
   su2double *dSm_dU, *dPI_dU, *drhoStar_dU, *dpStar_dU, *dEStar_dU;
 
   su2double* Flux;        /*!< \brief The flux accross the face. */
+  bool fixedWaveSpeeds = false;     /*!< \brief Use fixed_sL, fixed_sR as wave speeds (verification of the Jacobians). */
+  su2double fixed_sL = 0.0, fixed_sR = 0.0;
   su2double** Jacobian_i; /*!< \brief The Jacobian w.r.t. point i after computation. */
   su2double** Jacobian_j; /*!< \brief The Jacobian w.r.t. point j after computation. */
 
@@ -78,7 +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;
-
 };
 
 /*!
@@ -90,7 +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_EXACT_HLLC_JACOBIANS). */
   unsigned short iDim, jDim, iVar, jVar;
 
   su2double *IntermediateState;
@@ -108,6 +113,8 @@ class CUpwGeneralHLLC_Flow final : public CNumerics {
   su2double *dSm_dU, *dPI_dU, *drhoStar_dU, *dpStar_dU, *dEStar_dU;
 
   su2double* Flux;        /*!< \brief The flux accross the face. */
+  bool fixedWaveSpeeds = false;     /*!< \brief Use fixed_sL, fixed_sR as wave speeds (verification of the Jacobians). */
+  su2double fixed_sL = 0.0, fixed_sR = 0.0;
   su2double** Jacobian_i; /*!< \brief The Jacobian w.r.t. point i after computation. */
   su2double** Jacobian_j; /*!< \brief The Jacobian w.r.t. point j after computation. */
 
@@ -132,6 +139,7 @@ class CUpwGeneralHLLC_Flow final : public CNumerics {
    */
   ResidualType<> ComputeResidual(const CConfig* config) override;
 
+
   /*!
    * \brief Compute the Average quantities for a general fluid flux between two nodes i and j.
    * Using the approach of Vinokur and Montagne'
diff --git a/SU2_CFD/src/numerics/flow/convection/hllc.cpp b/SU2_CFD/src/numerics/flow/convection/hllc.cpp
index ef6a2a5a41e9..131ee5b65fa7 100644
--- a/SU2_CFD/src/numerics/flow/convection/hllc.cpp
+++ b/SU2_CFD/src/numerics/flow/convection/hllc.cpp
@@ -28,12 +28,59 @@
 #include "../../../../include/numerics/flow/convection/hllc.hpp"
 #include "../../../../../Common/include/toolboxes/geometry_toolbox.hpp"
 
+namespace {
+/*--- Moving grids. The flux across a face that moves with normal velocity w is the flux of the states seen from the
+ face (velocity u' = u - c, c = w n, same thermodynamic state) transformed back: with U' = T U,
+ F_ALE(U) = B F'(U'), where F' is the flux for a fixed face. Mass: unchanged; momentum: + c F'_0;
+ energy: + c.F'_m + |c|^2/2 F'_0. The Jacobians are B J' T. ---*/
+
+void ShiftToFaceFrame(unsigned short nDim, su2double w, const su2double* unitNormal, su2double* velocity,
+                      su2double& enthalpy) {
+  const su2double projVel = GeometryToolbox::DotProduct(nDim, velocity, unitNormal);
+  for (unsigned short iDim = 0; iDim < nDim; iDim++) velocity[iDim] -= w * unitNormal[iDim];
+  enthalpy += w * (0.5 * w - projVel);  // H' = H - u.c + |c|^2/2
+}
+
+void FaceFrameToMovingFlux(unsigned short nDim, su2double w, const su2double* unitNormal, su2double* flux) {
+  const unsigned short nVar = nDim + 2;
+  su2double cDotMom = 0.0;
+  for (unsigned short iDim = 0; iDim < nDim; iDim++) cDotMom += w * unitNormal[iDim] * flux[iDim+1];
+  flux[nVar-1] += cDotMom + 0.5 * w * w * flux[0];
+  for (unsigned short iDim = 0; iDim < nDim; iDim++) flux[iDim+1] += w * unitNormal[iDim] * flux[0];
+}
+
+void FaceFrameToMovingJacobian(unsigned short nDim, su2double w, const su2double* unitNormal, su2double** jac) {
+  const unsigned short nVar = nDim + 2;
+  constexpr unsigned short maxDim = 3, maxVar = maxDim + 2;
+  su2double c[maxDim] = {0.0}, T[maxVar][maxVar] = {{0.0}}, B[maxVar][maxVar] = {{0.0}}, tmp[maxVar][maxVar];
+  su2double c2 = 0.0;
+  for (unsigned short iDim = 0; iDim < nDim; iDim++) { c[iDim] = w * unitNormal[iDim]; c2 += c[iDim] * c[iDim]; }
+  for (unsigned short iVar = 0; iVar < nVar; iVar++) { T[iVar][iVar] = 1.0; B[iVar][iVar] = 1.0; }
+  for (unsigned short iDim = 0; iDim < nDim; iDim++) {
+    T[iDim+1][0] = -c[iDim];  T[nVar-1][iDim+1] = -c[iDim];
+    B[iDim+1][0] = c[iDim];   B[nVar-1][iDim+1] = c[iDim];
+  }
+  T[nVar-1][0] = 0.5 * c2;  B[nVar-1][0] = 0.5 * c2;
+  for (unsigned short iVar = 0; iVar < nVar; iVar++)
+    for (unsigned short jVar = 0; jVar < nVar; jVar++) {
+      tmp[iVar][jVar] = 0.0;
+      for (unsigned short kVar = 0; kVar < nVar; kVar++) tmp[iVar][jVar] += jac[iVar][kVar] * T[kVar][jVar];
+    }
+  for (unsigned short iVar = 0; iVar < nVar; iVar++)
+    for (unsigned short jVar = 0; jVar < nVar; jVar++) {
+      jac[iVar][jVar] = 0.0;
+      for (unsigned short kVar = 0; kVar < nVar; kVar++) jac[iVar][jVar] += B[iVar][kVar] * tmp[kVar][jVar];
+    }
+}
+}  // namespace
+
 CUpwHLLC_Flow::CUpwHLLC_Flow(unsigned short val_nDim, unsigned short val_nVar, const CConfig* config) : CNumerics(val_nDim, val_nVar, config) {
 
   implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT);
 
   /* 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_Exact_HLLC_Jacobians();
 
   Gamma = config->GetGamma();
 
@@ -114,6 +161,16 @@ CNumerics::ResidualType<> CUpwHLLC_Flow::ComputeResidual(const CConfig* config)
   Density_j  = V_j[nDim+2];
   Enthalpy_j = V_j[nDim+3];
 
+  /*--- Moving grid: states seen from the face (see ShiftToFaceFrame). ---*/
+
+  ProjInterfaceVel = 0.0;
+  if (dynamic_grid) {
+    for (iDim = 0; iDim < nDim; iDim++)
+      ProjInterfaceVel += 0.5 * (GridVel_i[iDim] + GridVel_j[iDim]) * UnitNormal[iDim];
+    ShiftToFaceFrame(nDim, ProjInterfaceVel, UnitNormal, Velocity_i, Enthalpy_i);
+    ShiftToFaceFrame(nDim, ProjInterfaceVel, UnitNormal, Velocity_j, Enthalpy_j);
+  }
+
   sq_vel_i = GeometryToolbox::SquaredNorm(nDim, Velocity_i);
   sq_vel_j = GeometryToolbox::SquaredNorm(nDim, Velocity_j);
 
@@ -128,22 +185,6 @@ CNumerics::ResidualType<> CUpwHLLC_Flow::ComputeResidual(const CConfig* config)
   ProjVelocity_i = GeometryToolbox::DotProduct(nDim, Velocity_i, UnitNormal);
   ProjVelocity_j = GeometryToolbox::DotProduct(nDim, Velocity_j, UnitNormal);
 
-  /*--- Projected Grid Velocity ---*/
-
-  ProjInterfaceVel = 0;
-
-  if (dynamic_grid) {
-
-    for (iDim = 0; iDim < nDim; iDim++)
-      ProjInterfaceVel += 0.5 * (GridVel_i[iDim] + GridVel_j[iDim]) * UnitNormal[iDim];
-
-    SoundSpeed_i -= ProjInterfaceVel;
-    SoundSpeed_j += ProjInterfaceVel;
-
-    ProjVelocity_i -= ProjInterfaceVel;
-    ProjVelocity_j -= ProjInterfaceVel;
-  }
-
   /*--- Roe's averaging ---*/
 
   Rrho = ( sqrt(Density_i) + sqrt(Density_j) );
@@ -152,7 +193,7 @@ CNumerics::ResidualType<> CUpwHLLC_Flow::ComputeResidual(const CConfig* config)
     RoeVelocity[iDim] = ( Velocity_i[iDim] * sqrt(Density_i) + Velocity_j[iDim] * sqrt(Density_j) ) / Rrho;
   }
   sq_velRoe = GeometryToolbox::SquaredNorm(nDim, RoeVelocity);
-  RoeProjVelocity = GeometryToolbox::DotProduct(nDim, RoeVelocity, UnitNormal) - ProjInterfaceVel;
+  RoeProjVelocity = GeometryToolbox::DotProduct(nDim, RoeVelocity, UnitNormal);
 
   /*--- Mean Roe variables iPoint and jPoint ---*/
 
@@ -161,13 +202,15 @@ CNumerics::ResidualType<> CUpwHLLC_Flow::ComputeResidual(const CConfig* config)
 
   /*--- Roe-averaged speed of sound ---*/
 
-  RoeSoundSpeed  = sqrt( Gamma_Minus_One * ( RoeEnthalpy - 0.5 * sq_velRoe  ) ) - ProjInterfaceVel;
+  RoeSoundSpeed  = sqrt( Gamma_Minus_One * ( RoeEnthalpy - 0.5 * sq_velRoe ) );
 
   /*--- Speed of sound at L and R ---*/
 
   sL = min( RoeProjVelocity - RoeSoundSpeed, ProjVelocity_i - SoundSpeed_i);
   sR = max( RoeProjVelocity + RoeSoundSpeed, ProjVelocity_j + SoundSpeed_j);
 
+  if (fixedWaveSpeeds) { sL = fixed_sL; sR = fixed_sR; }
+
   /*--- speed of contact surface ---*/
 
   RHO = Density_j * (sR - ProjVelocity_j) - Density_i * (sL - ProjVelocity_i);
@@ -205,7 +248,7 @@ CNumerics::ResidualType<> CUpwHLLC_Flow::ComputeResidual(const CConfig* config)
       Flux[0] = sM * IntermediateState[0];
       for (iDim = 0; iDim < nDim; iDim++)
         Flux[iDim+1] = sM * IntermediateState[iDim+1] + pStar * UnitNormal[iDim];
-      Flux[nVar-1] = sM * ( IntermediateState[nVar-1] + pStar ) + pStar * ProjInterfaceVel;
+      Flux[nVar-1] = sM * ( IntermediateState[nVar-1] + pStar );
     }
   }
   else {
@@ -234,16 +277,26 @@ CNumerics::ResidualType<> CUpwHLLC_Flow::ComputeResidual(const CConfig* config)
       Flux[0] = sM * IntermediateState[0];
       for (iDim = 0; iDim < nDim; iDim++)
         Flux[iDim+1] = sM * IntermediateState[iDim+1] + pStar * UnitNormal[iDim];
-      Flux[nVar-1] = sM * (IntermediateState[nVar-1] + pStar ) + pStar * ProjInterfaceVel;
+      Flux[nVar-1] = sM * (IntermediateState[nVar-1] + pStar );
     }
   }
 
+  if (dynamic_grid) FaceFrameToMovingFlux(nDim, ProjInterfaceVel, UnitNormal, Flux);
+
   for (iVar = 0; iVar < nVar; iVar++) Flux[iVar] *= Area;
 
   /*--- Return early if the Jacobians do not need to be computed. ---*/
 
   if (!implicit) return ResidualType<>(Flux, Jacobian_i, Jacobian_j);
 
+  /*--- 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_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) {
@@ -254,7 +307,8 @@ CNumerics::ResidualType<> CUpwHLLC_Flow::ComputeResidual(const CConfig* config)
         for (jVar = 0; jVar < nVar; jVar++)
           Jacobian_j[iVar][jVar] = 0;
 
-      GetInviscidProjJac(Velocity_i, &Energy_i, UnitNormal, 1.0, Jacobian_i);
+      const su2double chi = 0.0;
+      GetInviscidProjJac(Velocity_i, &Enthalpy_i, &chi, &Gamma_Minus_One, UnitNormal, 1.0, Jacobian_i);
 
     }
     else {
@@ -295,7 +349,7 @@ CNumerics::ResidualType<> CUpwHLLC_Flow::ComputeResidual(const CConfig* config)
       /*--- Computing d/dU_L (pStar) ---*/
 
       for (iVar = 0; iVar < nVar; iVar++)
-        dpStar_dU[iVar] = Density_i * (sR - ProjVelocity_j) * dSm_dU[iVar];
+        dpStar_dU[iVar] = rhoPStar_R * (sR - ProjVelocity_j) * dSm_dU[iVar];
 
 
       /*--- Computing d/dU_L (EStar) ---*/
@@ -352,7 +406,7 @@ CNumerics::ResidualType<> CUpwHLLC_Flow::ComputeResidual(const CConfig* config)
       /*--- Computing d/dU_R (pStar) ---*/
 
       for (iVar = 0; iVar < nVar; iVar++)
-        dpStar_dU[iVar] = Density_j * (sL - ProjVelocity_i) * dSm_dU[iVar];
+        dpStar_dU[iVar] = rhoPStar_L * (sL - ProjVelocity_i) * dSm_dU[iVar];
 
 
       /*--- Computing d/dU_R (EStar) ---*/
@@ -389,7 +443,8 @@ CNumerics::ResidualType<> CUpwHLLC_Flow::ComputeResidual(const CConfig* config)
         for (jVar = 0; jVar < nVar; jVar++)
           Jacobian_i[iVar][jVar] = 0;
 
-      GetInviscidProjJac(Velocity_j, &Energy_j, UnitNormal, 1.0, Jacobian_j);
+      const su2double chi = 0.0;
+      GetInviscidProjJac(Velocity_j, &Enthalpy_j, &chi, &Gamma_Minus_One, UnitNormal, 1.0, Jacobian_j);
 
     }
     else {
@@ -414,7 +469,7 @@ CNumerics::ResidualType<> CUpwHLLC_Flow::ComputeResidual(const CConfig* config)
       /*--- Computing d/dU_L (pStar) ---*/
 
       for (iVar = 0; iVar < nVar; iVar++)
-        dpStar_dU[iVar] = Density_i * (sR - ProjVelocity_j) * dSm_dU[iVar];
+        dpStar_dU[iVar] = rhoPStar_R * (sR - ProjVelocity_j) * dSm_dU[iVar];
 
 
       /*--- Computing d/dU_L (EStar) ---*/
@@ -466,7 +521,7 @@ CNumerics::ResidualType<> CUpwHLLC_Flow::ComputeResidual(const CConfig* config)
       /*--- Computing d/dU_R (pStar) ---*/
 
       for (iVar = 0; iVar < nVar; iVar++)
-        dpStar_dU[iVar] = Density_j * (sL - ProjVelocity_i) * dSm_dU[iVar];
+        dpStar_dU[iVar] = rhoPStar_L * (sL - ProjVelocity_i) * dSm_dU[iVar];
 
 
       /*--- Computing d/dU_R (rhoStar) ---*/
@@ -519,6 +574,11 @@ CNumerics::ResidualType<> CUpwHLLC_Flow::ComputeResidual(const CConfig* config)
     }
   }
 
+  if (dynamic_grid) {
+    FaceFrameToMovingJacobian(nDim, ProjInterfaceVel, UnitNormal, Jacobian_i);
+    FaceFrameToMovingJacobian(nDim, ProjInterfaceVel, UnitNormal, Jacobian_j);
+  }
+
   /*--- Scale Jacobians by area (from Flux *= Area). ---*/
   for (iVar = 0; iVar < nVar; iVar++) {
     for (jVar = 0; jVar < nVar; jVar++) {
@@ -535,6 +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_Exact_HLLC_Jacobians();
 
   Gamma = config->GetGamma();
 
@@ -613,6 +674,16 @@ CNumerics::ResidualType<> CUpwGeneralHLLC_Flow::ComputeResidual(const CConfig* c
   Density_j  = V_j[nDim+2];
   Enthalpy_j = V_j[nDim+3];
 
+  /*--- Moving grid: states seen from the face (see ShiftToFaceFrame). ---*/
+
+  ProjInterfaceVel = 0.0;
+  if (dynamic_grid) {
+    for (iDim = 0; iDim < nDim; iDim++)
+      ProjInterfaceVel += 0.5 * (GridVel_i[iDim] + GridVel_j[iDim]) * UnitNormal[iDim];
+    ShiftToFaceFrame(nDim, ProjInterfaceVel, UnitNormal, Velocity_i, Enthalpy_i);
+    ShiftToFaceFrame(nDim, ProjInterfaceVel, UnitNormal, Velocity_j, Enthalpy_j);
+  }
+
 
   sq_vel_i = 0.0;
   sq_vel_j = 0.0;
@@ -650,28 +721,12 @@ CNumerics::ResidualType<> CUpwGeneralHLLC_Flow::ComputeResidual(const CConfig* c
   }
 
 
-  /*--- Projected Grid Velocity ---*/
-
-  ProjInterfaceVel = 0;
-
-  if (dynamic_grid) {
-
-    for (iDim = 0; iDim < nDim; iDim++)
-      ProjInterfaceVel += 0.5 * ( GridVel_i[iDim] + GridVel_j[iDim] )*UnitNormal[iDim];
-
-    SoundSpeed_i -= ProjInterfaceVel;
-    SoundSpeed_j += ProjInterfaceVel;
-
-    ProjVelocity_i -= ProjInterfaceVel;
-    ProjVelocity_j -= ProjInterfaceVel;
-  }
-
   /*--- Roe's averaging ---*/
 
   Rrho = ( sqrt(Density_i) + sqrt(Density_j) );
 
   sq_velRoe        = 0.0;
-  RoeProjVelocity  = - ProjInterfaceVel;
+  RoeProjVelocity  = 0.0;
 
   for (iDim = 0; iDim < nDim; iDim++) {
     RoeVelocity[iDim] = ( Velocity_i[iDim] * sqrt(Density_i) + Velocity_j[iDim] * sqrt(Density_j) ) / Rrho;
@@ -691,13 +746,15 @@ CNumerics::ResidualType<> CUpwGeneralHLLC_Flow::ComputeResidual(const CConfig* c
   /*--- Roe-averaged speed of sound ---*/
 
   //RoeSoundSpeed2 = RoeChi + RoeKappa * ( RoeEnthalpy - 0.5 * sq_velRoe );
-  RoeSoundSpeed  = sqrt( RoeChi + RoeKappa * ( RoeEnthalpy - 0.5 * sq_velRoe ) ) - ProjInterfaceVel;
+  RoeSoundSpeed  = sqrt( RoeChi + RoeKappa * ( RoeEnthalpy - 0.5 * sq_velRoe ) );
 
   /*--- Speed of sound at L and R ---*/
 
   sL = min( RoeProjVelocity - RoeSoundSpeed, ProjVelocity_i - SoundSpeed_i );
   sR = max( RoeProjVelocity + RoeSoundSpeed, ProjVelocity_j + SoundSpeed_j );
 
+  if (fixedWaveSpeeds) { sL = fixed_sL; sR = fixed_sR; }
+
   /*--- speed of contact surface ---*/
 
   RHO = Density_j * (sR - ProjVelocity_j) - Density_i * (sL - ProjVelocity_i);
@@ -734,7 +791,7 @@ CNumerics::ResidualType<> CUpwGeneralHLLC_Flow::ComputeResidual(const CConfig* c
       Flux[0] = sM * IntermediateState[0];
       for (iDim = 0; iDim < nDim; iDim++)
         Flux[iDim+1] = sM * IntermediateState[iDim+1] + pStar * UnitNormal[iDim];
-      Flux[nVar-1] = sM * ( IntermediateState[nVar-1] + pStar )  + pStar * ProjInterfaceVel;
+      Flux[nVar-1] = sM * ( IntermediateState[nVar-1] + pStar );
     }
   }
   else {
@@ -763,16 +820,26 @@ CNumerics::ResidualType<> CUpwGeneralHLLC_Flow::ComputeResidual(const CConfig* c
       Flux[0] = sM * IntermediateState[0];
       for (iDim = 0; iDim < nDim; iDim++)
         Flux[iDim+1] = sM * IntermediateState[iDim+1] + pStar * UnitNormal[iDim];
-      Flux[nVar-1] = sM * (IntermediateState[nVar-1] + pStar )  + pStar * ProjInterfaceVel;
+      Flux[nVar-1] = sM * (IntermediateState[nVar-1] + pStar );
     }
   }
 
+  if (dynamic_grid) FaceFrameToMovingFlux(nDim, ProjInterfaceVel, UnitNormal, Flux);
+
   for (iVar = 0; iVar < nVar; iVar++) Flux[iVar] *= Area;
 
   /*--- Return early if the Jacobians do not need to be computed. ---*/
 
   if (!implicit) return ResidualType<>(Flux, Jacobian_i, Jacobian_j);
 
+  /*--- 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_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) {
@@ -800,7 +867,7 @@ CNumerics::ResidualType<> CUpwGeneralHLLC_Flow::ComputeResidual(const CConfig* c
 
       /*--- Computing pressure derivatives d/dU_L (PI) ---*/
 
-      dPI_dU[0] = Chi_i - 0.5 * Kappa_i * sq_vel_i;
+      dPI_dU[0] = Chi_i + 0.5 * Kappa_i * sq_vel_i;
       for (iDim = 0; iDim < nDim; iDim++)
         dPI_dU[iDim+1] = - Kappa_i * Velocity_i[iDim];
       dPI_dU[nVar-1] = Kappa_i;
@@ -825,7 +892,7 @@ CNumerics::ResidualType<> CUpwGeneralHLLC_Flow::ComputeResidual(const CConfig* c
       /*--- Computing d/dU_L (pStar) ---*/
 
       for (iVar = 0; iVar < nVar; iVar++)
-        dpStar_dU[iVar] = Density_i * (sR - ProjVelocity_j) * dSm_dU[iVar];
+        dpStar_dU[iVar] = rhoPStar_R * (sR - ProjVelocity_j) * dSm_dU[iVar];
 
 
       /*--- Computing d/dU_L (EStar) ---*/
@@ -874,7 +941,7 @@ CNumerics::ResidualType<> CUpwGeneralHLLC_Flow::ComputeResidual(const CConfig* c
 
       /*--- Computing pressure derivatives d/dU_R (PI) ---*/
 
-      dPI_dU[0] = Chi_j - 0.5 * Kappa_j * sq_vel_j;
+      dPI_dU[0] = Chi_j + 0.5 * Kappa_j * sq_vel_j;
       for (iDim = 0; iDim < nDim; iDim++)
         dPI_dU[iDim+1] = - Kappa_j * Velocity_j[iDim];
       dPI_dU[nVar-1] = Kappa_j;
@@ -891,7 +958,7 @@ CNumerics::ResidualType<> CUpwGeneralHLLC_Flow::ComputeResidual(const CConfig* c
       /*--- Computing d/dU_R (pStar) ---*/
 
       for (iVar = 0; iVar < nVar; iVar++)
-        dpStar_dU[iVar] = Density_j * (sL - ProjVelocity_i) * dSm_dU[iVar];
+        dpStar_dU[iVar] = rhoPStar_L * (sL - ProjVelocity_i) * dSm_dU[iVar];
 
 
       /*--- Computing d/dU_R (EStar) ---*/
@@ -944,7 +1011,7 @@ CNumerics::ResidualType<> CUpwGeneralHLLC_Flow::ComputeResidual(const CConfig* c
 
       /*--- Computing pressure derivatives d/dU_L (PI) ---*/
 
-      dPI_dU[0] = Chi_i - 0.5 * Kappa_i * sq_vel_i;
+      dPI_dU[0] = Chi_i + 0.5 * Kappa_i * sq_vel_i;
       for (iDim = 0; iDim < nDim; iDim++)
         dPI_dU[iDim+1] = - Kappa_i * Velocity_i[iDim];
       dPI_dU[nVar-1] = Kappa_i;
@@ -961,7 +1028,7 @@ CNumerics::ResidualType<> CUpwGeneralHLLC_Flow::ComputeResidual(const CConfig* c
       /*--- Computing d/dU_L (pStar) ---*/
 
       for (iVar = 0; iVar < nVar; iVar++)
-        dpStar_dU[iVar] = Density_i * (sR - ProjVelocity_j) * dSm_dU[iVar];
+        dpStar_dU[iVar] = rhoPStar_R * (sR - ProjVelocity_j) * dSm_dU[iVar];
 
 
       /*--- Computing d/dU_L (EStar) ---*/
@@ -995,7 +1062,7 @@ CNumerics::ResidualType<> CUpwGeneralHLLC_Flow::ComputeResidual(const CConfig* c
 
       /*--- Computing pressure derivatives d/dU_R (PI) ---*/
 
-      dPI_dU[0] = Chi_j - 0.5 * Kappa_j * sq_vel_j;
+      dPI_dU[0] = Chi_j + 0.5 * Kappa_j * sq_vel_j;
       for (iDim = 0; iDim < nDim; iDim++)
         dPI_dU[iDim+1] = - Kappa_j * Velocity_j[iDim];
       dPI_dU[nVar-1] = Kappa_j;
@@ -1012,7 +1079,7 @@ CNumerics::ResidualType<> CUpwGeneralHLLC_Flow::ComputeResidual(const CConfig* c
       /*--- Computing d/dU_R (pStar) ---*/
 
       for (iVar = 0; iVar < nVar; iVar++)
-        dpStar_dU[iVar] = Density_j * (sL - ProjVelocity_i) * dSm_dU[iVar];
+        dpStar_dU[iVar] = rhoPStar_L * (sL - ProjVelocity_i) * dSm_dU[iVar];
 
 
       /*--- Computing d/dU_R (rhoStar) ---*/
@@ -1064,6 +1131,11 @@ CNumerics::ResidualType<> CUpwGeneralHLLC_Flow::ComputeResidual(const CConfig* c
     }
   }
 
+  if (dynamic_grid) {
+    FaceFrameToMovingJacobian(nDim, ProjInterfaceVel, UnitNormal, Jacobian_i);
+    FaceFrameToMovingJacobian(nDim, ProjInterfaceVel, UnitNormal, Jacobian_j);
+  }
+
   /*--- Scale Jacobians by area (from Flux *= Area). ---*/
   for (iVar = 0; iVar < nVar; iVar++) {
     for (jVar = 0; jVar < nVar; jVar++) {
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
new file mode 100644
index 000000000000..fc36c6cc9eba
--- /dev/null
+++ b/UnitTests/SU2_CFD/numerics/hllc_jacobian_tests.cpp
@@ -0,0 +1,292 @@
+/*!
+ * \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_EXACT_HLLC_JACOBIANS) or default.
+ * \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 <http://www.gnu.org/licenses/>.
+ */
+
+#include "catch.hpp"
+#include <array>
+#include <cmath>
+#include <memory>
+#include <sstream>
+#include "../../../SU2_CFD/include/numerics/flow/convection/hllc.hpp"
+
+struct HLLCJacobianTestAccess {
+  template <class Numerics>
+  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;
+constexpr unsigned short maxVar = 5;
+
+using State = std::array<su2double, maxVar>;  // conservative variables (rho, rho u, rho v, [rho w], rho E)
+
+/*--- Primitive (T, u, ..., p, rho, h, c) and secondary (dp/drho_e, dp/de_rho) variables of an ideal gas whose total
+ energy contains the turbulent kinetic energy k. ---*/
+void Primitives(unsigned short nDim, const State& U, su2double k, su2double* V, su2double* S) {
+  const su2double rho = U[0];
+  su2double q2 = 0.0;
+  for (unsigned short iDim = 0; iDim < nDim; iDim++) {
+    V[iDim + 1] = U[iDim + 1] / rho;
+    q2 += V[iDim + 1] * V[iDim + 1];
+  }
+  const su2double e = U[nDim + 1] / rho - 0.5 * q2 - k;
+  const su2double p = (gamma - 1) * rho * e;
+  V[0] = p / (rho * 287.058);
+  V[nDim + 1] = p;
+  V[nDim + 2] = rho;
+  V[nDim + 3] = (U[nDim + 1] + p) / rho;
+  V[nDim + 4] = sqrt(gamma * p / rho);
+  S[0] = (gamma - 1) * e;
+  S[1] = (gamma - 1) * rho;
+}
+
+/*--- Conservative state from density, velocity, pressure and k. ---*/
+State Conservative(unsigned short nDim, su2double rho, const su2double* vel, su2double p, su2double k) {
+  State U{};
+  su2double q2 = 0.0;
+  U[0] = rho;
+  for (unsigned short iDim = 0; iDim < nDim; iDim++) {
+    U[iDim + 1] = rho * vel[iDim];
+    q2 += vel[iDim] * vel[iDim];
+  }
+  U[nDim + 1] = p / (gamma - 1) + rho * (0.5 * q2 + k);
+  return U;
+}
+
+struct Result {
+  State flux;
+  su2double jac_i[maxVar][maxVar], jac_j[maxVar][maxVar];
+};
+
+Result Evaluate(unsigned short nDim, CNumerics& numerics, const CConfig* config, const State& Ui, const State& Uj,
+                su2double k_i, su2double k_j) {
+  su2double Vi[20] = {0.0}, Vj[20] = {0.0}, Si[2], Sj[2];
+  Primitives(nDim, Ui, k_i, Vi, Si);
+  Primitives(nDim, Uj, k_j, Vj, Sj);
+  numerics.SetPrimitive(Vi, Vj);
+  numerics.SetSecondary(Si, Sj);
+  numerics.SetTurbKineticEnergy(k_i, k_j);
+  const auto res = numerics.ComputeResidual(config);
+  Result r;
+  for (unsigned short iVar = 0; iVar < nDim + 2; iVar++) {
+    r.flux[iVar] = res.residual[iVar];
+    for (unsigned short jVar = 0; jVar < nDim + 2; jVar++) {
+      r.jac_i[iVar][jVar] = res.jacobian_i[iVar][jVar];
+      r.jac_j[iVar][jVar] = res.jacobian_j[iVar][jVar];
+    }
+  }
+  return r;
+}
+
+/*--- Largest difference between the Jacobians and central differences of the flux, relative to the largest entry. ---*/
+su2double JacobianError(unsigned short nDim, CNumerics& numerics, const CConfig* config, const State& Ui,
+                        const State& Uj, su2double k_i, su2double k_j) {
+  const unsigned short nVar = nDim + 2;
+  const auto ref = Evaluate(nDim, numerics, config, Ui, Uj, k_i, k_j);
+  su2double maxErr = 0.0, maxJac = 0.0;
+  for (int side = 0; side < 2; side++) {
+    for (unsigned short jVar = 0; jVar < nVar; jVar++) {
+      State Up = side ? Uj : Ui, Um = Up;
+      const su2double h = 1e-6 * std::max(std::abs(Up[jVar]), su2double(1.0));
+      Up[jVar] += h;
+      Um[jVar] -= h;
+      const auto fp = side ? Evaluate(nDim, numerics, config, Ui, Up, k_i, k_j)
+                           : Evaluate(nDim, numerics, config, Up, Uj, k_i, k_j);
+      const auto fm = side ? Evaluate(nDim, numerics, config, Ui, Um, k_i, k_j)
+                           : Evaluate(nDim, numerics, config, Um, Uj, k_i, k_j);
+      for (unsigned short iVar = 0; iVar < nVar; iVar++) {
+        const su2double fd = (fp.flux[iVar] - fm.flux[iVar]) / (2 * h);
+        const su2double jac = side ? ref.jac_j[iVar][jVar] : ref.jac_i[iVar][jVar];
+        maxErr = std::max(maxErr, std::abs(jac - fd));
+        maxJac = std::max(maxJac, std::abs(jac));
+      }
+    }
+  }
+  return maxErr / maxJac;
+}
+
+/*--- Contact speed sM of the scheme, in the frame of the face (w = grid velocity along the unit normal n). ---*/
+su2double ContactSpeed(unsigned short nDim, const State& Ui, const State& Uj, su2double k_i, su2double k_j,
+                       const su2double* n, su2double w, su2double sL, su2double sR) {
+  su2double Vi[20] = {0.0}, Vj[20] = {0.0}, S[2];
+  Primitives(nDim, Ui, k_i, Vi, S);
+  Primitives(nDim, Uj, k_j, Vj, S);
+  su2double u_i = -w, u_j = -w;
+  for (unsigned short iDim = 0; iDim < nDim; iDim++) {
+    u_i += Vi[iDim + 1] * n[iDim];
+    u_j += Vj[iDim + 1] * n[iDim];
+  }
+  const su2double p_i = Vi[nDim + 1], p_j = Vj[nDim + 1], rho_i = Vi[nDim + 2], rho_j = Vj[nDim + 2];
+  return (p_i - p_j - rho_i * u_i * (sL - u_i) + rho_j * u_j * (sR - u_j)) / (rho_j * (sR - u_j) - rho_i * (sL - u_i));
+}
+
+std::unique_ptr<CConfig> MakeConfig(bool moving, bool accurate) {
+  std::stringstream options;
+  options << "SOLVER= EULER\nCONV_NUM_METHOD_FLOW= HLLC\nTIME_DISCRE_FLOW= EULER_IMPLICIT\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<CConfig>(options, SU2_COMPONENT::SU2_CFD, false);
+}
+
+/*--- 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 <class HLLC>
+void CheckAllBranches(unsigned short nDim, bool moving, su2double k, bool accurate) {
+  auto config = MakeConfig(moving, accurate);
+  HLLC numerics(nDim, nDim + 2, config.get());
+
+  /*--- Oblique unit normal (face area 1.3) and, on moving faces, a grid velocity with a normal component. ---*/
+  const su2double n2[3] = {0.8, 0.6, 0.0}, n3[3] = {2.0 / 7, 3.0 / 7, 6.0 / 7};
+  const su2double* n = (nDim == 2) ? n2 : n3;
+  su2double normal[3] = {0.0}, gridVel[3] = {0.0};
+  for (unsigned short iDim = 0; iDim < nDim; iDim++) normal[iDim] = 1.3 * n[iDim];
+  if (moving) {
+    gridVel[0] = 0.3;
+    gridVel[1] = -0.1;
+    gridVel[2] = 0.2;
+  }
+  numerics.SetNormal(normal);
+  numerics.SetGridVel(gridVel, gridVel);
+  su2double w = 0.0;
+  for (unsigned short iDim = 0; iDim < nDim; iDim++) w += gridVel[iDim] * n[iDim];
+
+  /*--- Two state pairs with different density, pressure and velocity (also tangential): flow along +n (sM > 0,
+   left branches) and along -n (sM < 0, right branches). ---*/
+  const su2double tangent[3] = {0.1, -0.2, 0.15};
+  for (const su2double dir : {1.0, -1.0}) {
+    su2double vel_i[3], vel_j[3];
+    for (unsigned short iDim = 0; iDim < nDim; iDim++) {
+      vel_i[iDim] = dir * 0.5 * n[iDim] + tangent[iDim];
+      vel_j[iDim] = dir * 0.25 * n[iDim] - tangent[iDim];
+    }
+    const State Ui = Conservative(nDim, 1.2, vel_i, 1.1, k);
+    const State Uj = Conservative(nDim, 0.9, vel_j, 0.8, 0.8 * k);
+
+    /*--- Wave speeds relative to the face, held fixed (the Jacobians treat them as constants): the supersonic and
+     the star branch on the side of the flow. ---*/
+    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}) {
+      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);
+
+      /*--- Make sure the intended branch is the one selected. ---*/
+      if (supersonic)
+        REQUIRE(((dir > 0) ? s[0] : -s[1]) > 0.0);
+      else
+        REQUIRE(sM * dir > 0.0);
+
+      const su2double error = JacobianError(nDim, numerics, config.get(), Ui, Uj, k, 0.8 * k);
+      if (accurate || supersonic)
+        CHECK(error < 1e-6);
+      else
+        CHECK(error > 1e-3);  // the approximation is active
+    }
+  }
+}
+
+}  // namespace
+
+TEST_CASE("HLLC Jacobians match finite differences", "[HLLC]") {
+  for (const unsigned short nDim : {2, 3}) {
+    for (const bool moving : {false, true}) {
+      for (const bool accurate : {false, true}) {
+        CheckAllBranches<CUpwHLLC_Flow>(nDim, moving, 0.0, accurate);
+        CheckAllBranches<CUpwGeneralHLLC_Flow>(nDim, moving, 0.0, accurate);
+      }
+    }
+  }
+}
+
+TEST_CASE("HLLC supersonic Jacobian uses reconstructed face enthalpy", "[HLLC-face]") {
+  for (unsigned short nDim : {2, 3}) {
+    for (bool moving : {false, true}) {
+      for (const bool accurate : {false, true}) {
+        for (su2double direction : {-1.0, 1.0}) {
+          auto config = MakeConfig(moving, accurate);
+          CUpwHLLC_Flow numerics(nDim, nDim + 2, config.get());
+          const su2double n2[3] = {0.8, 0.6, 0.0}, n3[3] = {2.0 / 7, 3.0 / 7, 6.0 / 7};
+          const auto* normal = nDim == 2 ? n2 : n3;
+          su2double velocity[3] = {}, gridVelocity[3] = {};
+          for (auto dim = 0u; dim < nDim; ++dim) {
+            velocity[dim] = direction * 2.5 * normal[dim];
+            gridVelocity[dim] = moving ? 0.3 * normal[dim] : 0.0;
+          }
+          numerics.SetNormal(normal);
+          numerics.SetGridVel(gridVelocity, gridVelocity);
+          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) {
+              su2double Vi[20] = {}, Vj[20] = {}, Si[2], Sj[2];
+              Primitives(nDim, Ui, 0.0, Vi, Si);
+              Primitives(nDim, Uj, 0.0, Vj, Sj);
+              // A fixed pressure offset models independently reconstructed p and h.
+              Vi[nDim + 1] += pressureOffset;
+              Vj[nDim + 1] += pressureOffset;
+              Vi[nDim + 3] += pressureOffset / Ui[0];
+              Vj[nDim + 3] += pressureOffset / Uj[0];
+              numerics.SetPrimitive(Vi, Vj);
+              numerics.SetTurbKineticEnergy(0.0, 0.0);
+              const auto residual = numerics.ComputeResidual(config.get());
+              Result result{};
+              for (auto i = 0u; i < nDim + 2; ++i) {
+                result.flux[i] = residual.residual[i];
+                for (auto j = 0u; j < nDim + 2; ++j) {
+                  result.jac_i[i][j] = residual.jacobian_i[i][j];
+                  result.jac_j[i][j] = residual.jacobian_j[i][j];
+                }
+              }
+              return result;
+            };
+            CAPTURE(nDim, moving, accurate, direction, pressureOffset);
+            const auto reference = evaluate(state, state);
+            for (auto column = 0u; column < nDim + 2; ++column) {
+              State plus = state, minus = state;
+              const su2double step = 1e-6;
+              plus[column] += step;
+              minus[column] -= step;
+              const auto fp = direction > 0 ? evaluate(plus, state) : evaluate(state, plus);
+              const auto fm = direction > 0 ? evaluate(minus, state) : evaluate(state, minus);
+              for (auto row = 0u; row < nDim + 2; ++row) {
+                const su2double fd = (fp.flux[row] - fm.flux[row]) / (2 * step);
+                const su2double jac = direction > 0 ? reference.jac_i[row][column] : reference.jac_j[row][column];
+                CHECK(jac == Approx(fd).margin(1e-6));
+              }
+            }
+          }
+        }
+      }
+    }
+  }
+}
diff --git a/UnitTests/meson.build b/UnitTests/meson.build
index 76d85e251574..176abed8dba6 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/hllc_jacobian_tests.cpp',
                        'SU2_CFD/fluid/CFluidModel_tests.cpp',
                        'SU2_CFD/output/CNEMOCompOutput_tests.cpp',
                        'SU2_CFD/output/COutput_convergence_tests.cpp',
diff --git a/config_template.cfg b/config_template.cfg
index 5e611f22ed48..8795186ab3d9 100644
--- a/config_template.cfg
+++ b/config_template.cfg
@@ -1858,6 +1858,10 @@ LOW_MACH_PREC= NO
 % Slower per iteration but potentialy more stable and capable of higher CFL
 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.