@@ -51,9 +51,7 @@ FORCEINLINE Double umusclProjection(Double grad_proj,
5151 /* --- To maintain proper scaling for edge limiters, the result of ---*/
5252 /* --- this function is 0.5 * dV_ij^kap. ---*/
5353 /* -------------------------------------------------------------------*/
54- const Double cent = 0.5 * delta;
55- const Double upw = grad_proj - cent;
56- return (1.0 - kappa) * upw + (1.0 + kappa) * cent;
54+ return (1.0 - kappa) * grad_proj + kappa * delta;
5755}
5856
5957/* !
@@ -221,7 +219,7 @@ FORCEINLINE CPair<ReconVarType> reconstructPrimitives(Int iEdge, Int iPoint, Int
221219 }
222220
223221 if (muscl) {
224- /* --- Recompute density and enthalpy instead of reconstructing . ---*/
222+ /* --- Reconstruct density and enthalpy without using their gradients . ---*/
225223 constexpr auto nVarGrad = ReconVarType::nVar - 2 ;
226224 switch (limiterType) {
227225 case LIMITER::NONE:
@@ -234,12 +232,23 @@ FORCEINLINE CPair<ReconVarType> reconstructPrimitives(Int iEdge, Int iPoint, Int
234232 musclPointLimited<nVarGrad>(iPoint, jPoint, vector_ij, limiters, gradients, V, kappa, relax, ramp_val);
235233 break ;
236234 }
235+ /* --- Recompute density using the reconstructed pressure and temperature. ---*/
237236 V.i .density () = V.i .pressure () / (gasConst * V.i .temperature ());
238237 V.j .density () = V.j .pressure () / (gasConst * V.j .temperature ());
239238
239+ /* --- Reconstruct enthalpy using dH/dT = Cp and dH/dv = v. Recomputing enthalpy would cause
240+ * stability issues because we use rho E = rho H - P, which loses its relation to temperature
241+ * if both rho and H are recomputed. This only seems to be an issue for wall-function meshes.
242+ * NOTE: This "one-sided" reconstruction does not lose much of the U-MUSCL benefit, because
243+ * the static enthalpy is linear, and for the KE term we do the equivalent of using the
244+ * average dH/dv, which is exact for quadratic functions. ---*/
240245 const su2double cp = gasConst * gamma / (gamma - 1 );
241- V.i .enthalpy () = cp * V.i .temperature () + 0.5 * squaredNorm<nDim>(V.i .velocity ());
242- V.j .enthalpy () = cp * V.j .temperature () + 0.5 * squaredNorm<nDim>(V.j .velocity ());
246+ V.i .enthalpy () += cp * (V.i .temperature () - V1st.i .temperature ());
247+ V.j .enthalpy () += cp * (V.j .temperature () - V1st.j .temperature ());
248+ for (size_t iDim = 0 ; iDim < nDim; ++iDim) {
249+ V.i .enthalpy () += 0.5 * (pow (V.i .velocity (iDim), 2 ) - pow (V1st.i .velocity (iDim), 2 ));
250+ V.j .enthalpy () += 0.5 * (pow (V.j .velocity (iDim), 2 ) - pow (V1st.j .velocity (iDim), 2 ));
251+ }
243252
244253 /* --- Detect a non-physical reconstruction based on negative pressure or density. ---*/
245254 const Double neg_p_or_rho = fmax (fmin (V.i .pressure (), V.j .pressure ()) < 0.0 ,
0 commit comments