Skip to content

Make reconstructed rhoE more stable #2608

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
pcarruscag merged 17 commits into from
Nov 30, 2025
Merged

Make reconstructed rhoE more stable #2608

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
17 commits
Select commit Hold shift + click to select a range
a2125ea
reconstruct but still without gradient
pcarruscag Nov 1, 2025
68106c2
Merge remote-tracking branch 'upstream/develop' into pedro/fix_muscl_…
pcarruscag Nov 15, 2025
c18d8f6
some tweaks and docs
pcarruscag Nov 15, 2025
dc731ca
make some cases run better
pcarruscag Nov 16, 2025
83db958
re-test
pcarruscag Nov 24, 2025
f26cf7e
Merge remote-tracking branch 'origin/develop' into pedro/fix_muscl_en…
pcarruscag Nov 29, 2025
7aff386
update vandv
pcarruscag Nov 29, 2025
1656474
hybrid
pcarruscag Nov 29, 2025
d9c1bfb
hybrid AD
pcarruscag Nov 29, 2025
a3c116b
parallel
pcarruscag Nov 29, 2025
45d940f
parallel AD
pcarruscag Nov 29, 2025
d4708b8
serial
pcarruscag Nov 29, 2025
7302c37
serial AD
pcarruscag Nov 29, 2025
e6d49d5
[skip ci]
pcarruscag Nov 29, 2025
2553761
update diff
pcarruscag Nov 29, 2025
4508d88
revert opdi [skip ci]
pcarruscag Nov 30, 2025
610c66c
push
pcarruscag Nov 30, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
21 changes: 15 additions & 6 deletions SU2_CFD/include/numerics_simd/flow/convection/common.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -51,9 +51,7 @@ FORCEINLINE Double umusclProjection(Double grad_proj,
/*--- To maintain proper scaling for edge limiters, the result of ---*/
/*--- this function is 0.5 * dV_ij^kap. ---*/
/*-------------------------------------------------------------------*/
const Double cent = 0.5 * delta;
const Double upw = grad_proj - cent;
return (1.0 - kappa) * upw + (1.0 + kappa) * cent;
return (1.0 - kappa) * grad_proj + kappa * delta;
}

/*!
Expand Down Expand Up @@ -216,7 +214,7 @@ FORCEINLINE CPair<ReconVarType> reconstructPrimitives(Int iEdge, Int iPoint, Int
}

if (muscl) {
/*--- Recompute density and enthalpy instead of reconstructing. ---*/
/*--- Reconstruct density and enthalpy without using their gradients. ---*/
constexpr auto nVarGrad = ReconVarType::nVar - 2;
switch (limiterType) {
case LIMITER::NONE:
Expand All @@ -229,12 +227,23 @@ FORCEINLINE CPair<ReconVarType> reconstructPrimitives(Int iEdge, Int iPoint, Int
musclPointLimited<nVarGrad>(iPoint, jPoint, vector_ij, limiters, gradients, V, kappa, relax);
break;
}
/*--- Recompute density using the reconstructed pressure and temperature. ---*/
V.i.density() = V.i.pressure() / (gasConst * V.i.temperature());
V.j.density() = V.j.pressure() / (gasConst * V.j.temperature());

/*--- Reconstruct enthalpy using dH/dT = Cp and dH/dv = v. Recomputing enthalpy would cause
* stability issues because we use rho E = rho H - P, which loses its relation to temperature
* if both rho and H are recomputed. This only seems to be an issue for wall-function meshes.
* NOTE: This "one-sided" reconstruction does not lose much of the U-MUSCL benefit, because
* the static enthalpy is linear, and for the KE term we do the equivalent of using the
* average dH/dv, which is exact for quadratic functions. ---*/
const su2double cp = gasConst * gamma / (gamma - 1);
V.i.enthalpy() = cp * V.i.temperature() + 0.5 * squaredNorm<nDim>(V.i.velocity());
V.j.enthalpy() = cp * V.j.temperature() + 0.5 * squaredNorm<nDim>(V.j.velocity());
V.i.enthalpy() += cp * (V.i.temperature() - V1st.i.temperature());
V.j.enthalpy() += cp * (V.j.temperature() - V1st.j.temperature());
for (size_t iDim = 0; iDim < nDim; ++iDim) {
V.i.enthalpy() += 0.5 * (pow(V.i.velocity(iDim), 2) - pow(V1st.i.velocity(iDim), 2));
V.j.enthalpy() += 0.5 * (pow(V.j.velocity(iDim), 2) - pow(V1st.j.velocity(iDim), 2));
}

/*--- Detect a non-physical reconstruction based on negative pressure or density. ---*/
const Double neg_p_or_rho = fmax(fmin(V.i.pressure(), V.j.pressure()) < 0.0,
Expand Down
16 changes: 8 additions & 8 deletions TestCases/disc_adj_fsi/dyn_fsi/grad_dv.opt.ref
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,9 +1,9 @@
INDEX GRAD
0 -4.570869180445024e-04
1 -2.401466726069483e-04
2 -9.134698219534478e-05
3 -1.628086953114332e-05
4 -1.741052432426715e-05
5 -9.462489852612668e-05
6 -2.452466275723706e-04
7 -4.635483655821636e-04
0 -4.570869126273463e-04
1 -2.401466681907557e-04
2 -9.134697929286435e-05
3 -1.628086861828660e-05
4 -1.741052581631559e-05
5 -9.462490278520256e-05
6 -2.452466349703724e-04
7 -4.635483765277457e-04
58 changes: 29 additions & 29 deletions TestCases/hybrid_regression.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -103,7 +103,7 @@ def main():
flatplate.cfg_dir = "navierstokes/flatplate"
flatplate.cfg_file = "lam_flatplate.cfg"
flatplate.test_iter = 100
flatplate.test_vals = [-7.680034, -2.207912, 0.001084, 0.036233, 2.361500, -2.325300, 0.000000, 0.000000]
flatplate.test_vals = [-7.680046, -2.207922, 0.001084, 0.036233, 2.361500, -2.325300, 0.000000, 0.000000]
test_list.append(flatplate)

# Laminar cylinder (steady)
Expand All @@ -119,7 +119,7 @@ def main():
cylinder_lowmach.cfg_dir = "navierstokes/cylinder"
cylinder_lowmach.cfg_file = "cylinder_lowmach.cfg"
cylinder_lowmach.test_iter = 25
cylinder_lowmach.test_vals = [-6.830996, -1.368850, -0.143956, 73.963354, 0]
cylinder_lowmach.test_vals = [-6.830997, -1.368850, -0.143986, 73.963289, 0.000000]
cylinder_lowmach.test_vals_aarch64 = [-6.830996, -1.368850, -0.143956, 73.963354, 0]
test_list.append(cylinder_lowmach)

Expand All @@ -128,15 +128,15 @@ def main():
poiseuille.cfg_dir = "navierstokes/poiseuille"
poiseuille.cfg_file = "lam_poiseuille.cfg"
poiseuille.test_iter = 10
poiseuille.test_vals = [-5.046131, 0.652984, 0.008355, 13.735818, 0]
poiseuille.test_vals = [-5.046131, 0.652984, 0.008355, 13.735816, 0.000000]
test_list.append(poiseuille)

# 2D Poiseuille flow (inlet profile file)
poiseuille_profile = TestCase('poiseuille_profile')
poiseuille_profile.cfg_dir = "navierstokes/poiseuille"
poiseuille_profile.cfg_file = "profile_poiseuille.cfg"
poiseuille_profile.test_iter = 10
poiseuille_profile.test_vals = [-12.009008, -7.262217, -0.000000, 2.089953]
poiseuille_profile.test_vals = [-12.008997, -7.262292, -0.000000, 2.089953]
poiseuille_profile.test_vals_aarch64 = [-12.009012, -7.262530, -0.000000, 2.089953]
test_list.append(poiseuille_profile)

Expand All @@ -145,7 +145,7 @@ def main():
periodic2d.cfg_dir = "navierstokes/periodic2D"
periodic2d.cfg_file = "config.cfg"
periodic2d.test_iter = 1400
periodic2d.test_vals = [-10.817611, -8.363544, -8.287460, -5.334104, -1.088411, -2945.2]
periodic2d.test_vals = [-10.817611, -8.363544, -8.287461, -5.334104, -1.088411, -2945.200000]
test_list.append(periodic2d)

##########################
Expand All @@ -157,7 +157,7 @@ def main():
rae2822_sa.cfg_dir = "rans/rae2822"
rae2822_sa.cfg_file = "turb_SA_RAE2822.cfg"
rae2822_sa.test_iter = 20
rae2822_sa.test_vals = [-2.020123, -5.269339, 0.807147, 0.060499, 0]
rae2822_sa.test_vals = [-2.020123, -5.269264, 0.807147, 0.060494, 0.000000]
test_list.append(rae2822_sa)

# RAE2822 SST
Expand Down Expand Up @@ -197,7 +197,7 @@ def main():
turb_naca0012_sa.cfg_dir = "rans/naca0012"
turb_naca0012_sa.cfg_file = "turb_NACA0012_sa.cfg"
turb_naca0012_sa.test_iter = 5
turb_naca0012_sa.test_vals = [-12.038087, -16.332090, 1.080346, 0.018385, 20.000000, -2.873540, 0.000000, -14.250271, 0.000000]
turb_naca0012_sa.test_vals = [-12.038066, -16.332090, 1.080346, 0.018385, 20.000000, -2.873343, 0.000000, -14.250271, 0.000000]
turb_naca0012_sa.test_vals_aarch64 = [-12.038091, -16.332090, 1.080346, 0.018385, 20.000000, -2.873236, 0.000000, -14.250271, 0.000000]
test_list.append(turb_naca0012_sa)

Expand All @@ -206,7 +206,7 @@ def main():
turb_naca0012_sst.cfg_dir = "rans/naca0012"
turb_naca0012_sst.cfg_file = "turb_NACA0012_sst.cfg"
turb_naca0012_sst.test_iter = 10
turb_naca0012_sst.test_vals = [-12.075727, -15.246732, -5.861249, 1.070036, 0.015841, -2.835749, 0.000000]
turb_naca0012_sst.test_vals = [-12.093953, -15.251079, -5.906327, 1.070413, 0.015775, -2.855686, 0.000000]
turb_naca0012_sst.test_vals_aarch64 = [-12.075928, -15.246732, -5.861249, 1.070036, 0.015841, -2.835263, 0.000000]
test_list.append(turb_naca0012_sst)

Expand All @@ -215,7 +215,7 @@ def main():
turb_naca0012_sst_sust.cfg_dir = "rans/naca0012"
turb_naca0012_sst_sust.cfg_file = "turb_NACA0012_sst_sust.cfg"
turb_naca0012_sst_sust.test_iter = 10
turb_naca0012_sst_sust.test_vals = [-12.073281, -14.836724, -5.732627, 1.000050, 0.019144, -2.629409]
turb_naca0012_sst_sust.test_vals = [-12.080740, -14.837176, -5.732917, 1.000893, 0.019109, -2.120257]
turb_naca0012_sst_sust.test_vals_aarch64 = [-12.073210, -14.836724, -5.732627, 1.000050, 0.019144, -2.629689]
test_list.append(turb_naca0012_sst_sust)

Expand All @@ -224,15 +224,15 @@ def main():
turb_naca0012_sst_fixedvalues.cfg_dir = "rans/naca0012"
turb_naca0012_sst_fixedvalues.cfg_file = "turb_NACA0012_sst_fixedvalues.cfg"
turb_naca0012_sst_fixedvalues.test_iter = 10
turb_naca0012_sst_fixedvalues.test_vals = [-5.192397, -10.445236, 0.774100, 1.022534, 0.040529, -2.383406]
turb_naca0012_sst_fixedvalues.test_vals = [-5.192389, -10.445226, 0.774100, 1.022534, 0.040529, -2.383436]
test_list.append(turb_naca0012_sst_fixedvalues)

# NACA0012 (SST, explicit Euler for flow and turbulence equations)
turb_naca0012_sst_expliciteuler = TestCase('turb_naca0012_sst_expliciteuler')
turb_naca0012_sst_expliciteuler.cfg_dir = "rans/naca0012"
turb_naca0012_sst_expliciteuler.cfg_file = "turb_NACA0012_sst_expliciteuler.cfg"
turb_naca0012_sst_expliciteuler.test_iter = 10
turb_naca0012_sst_expliciteuler.test_vals = [-3.533766, -3.157224, 3.743381, 1.124757, 0.501700, -float("inf")]
turb_naca0012_sst_expliciteuler.test_vals = [-3.532365, -3.157224, 3.743381, 1.124798, 0.501715, -float("inf")]
test_list.append(turb_naca0012_sst_expliciteuler)

# PROPELLER
Expand All @@ -252,7 +252,7 @@ def main():
axi_rans_air_nozzle_restart.cfg_dir = "axisymmetric_rans/air_nozzle"
axi_rans_air_nozzle_restart.cfg_file = "air_nozzle_restart.cfg"
axi_rans_air_nozzle_restart.test_iter = 10
axi_rans_air_nozzle_restart.test_vals = [-14.137306, -9.104064, -10.882700, -5.808219, 0.000000]
axi_rans_air_nozzle_restart.test_vals = [-12.066672, -7.446472, -8.813335, -3.730660, 0.000000]
axi_rans_air_nozzle_restart.test_vals_aarch64 = [-14.140441, -9.154674, -10.886121, -5.806594, 0.000000]
test_list.append(axi_rans_air_nozzle_restart)

Expand All @@ -278,7 +278,7 @@ def main():
turb_naca0012_1c.cfg_dir = "rans_uq/naca0012"
turb_naca0012_1c.cfg_file = "turb_NACA0012_uq_1c.cfg"
turb_naca0012_1c.test_iter = 10
turb_naca0012_1c.test_vals = [-4.976620, 1.345983, 0.433171, -0.033685]
turb_naca0012_1c.test_vals = [-4.979757, 1.345555, 0.450656, -0.030217]
turb_naca0012_1c.test_vals_aarch64 = [-4.976620, 1.345983, 0.433171, -0.033685]
test_list.append(turb_naca0012_1c)

Expand All @@ -287,7 +287,7 @@ def main():
turb_naca0012_2c.cfg_dir = "rans_uq/naca0012"
turb_naca0012_2c.cfg_file = "turb_NACA0012_uq_2c.cfg"
turb_naca0012_2c.test_iter = 10
turb_naca0012_2c.test_vals = [-5.485484, 1.263406, 0.411442, -0.040859]
turb_naca0012_2c.test_vals = [-5.482860, 1.263778, 0.403145, -0.043182]
turb_naca0012_2c.test_vals_aarch64 = [-5.485484, 1.263406, 0.411442, -0.040859]
test_list.append(turb_naca0012_2c)

Expand All @@ -296,7 +296,7 @@ def main():
turb_naca0012_3c.cfg_dir = "rans_uq/naca0012"
turb_naca0012_3c.cfg_file = "turb_NACA0012_uq_3c.cfg"
turb_naca0012_3c.test_iter = 10
turb_naca0012_3c.test_vals = [-5.583737, 1.232005, 0.390258, -0.046305]
turb_naca0012_3c.test_vals = [-5.583729, 1.232137, 0.384897, -0.047679]
turb_naca0012_3c.test_vals_aarch64 = [-5.583737, 1.232005, 0.390258, -0.046305]
test_list.append(turb_naca0012_3c)

Expand All @@ -305,7 +305,7 @@ def main():
turb_naca0012_p1c1.cfg_dir = "rans_uq/naca0012"
turb_naca0012_p1c1.cfg_file = "turb_NACA0012_uq_p1c1.cfg"
turb_naca0012_p1c1.test_iter = 10
turb_naca0012_p1c1.test_vals = [-5.114189, 1.285037, 0.406851, -0.043003]
turb_naca0012_p1c1.test_vals = [-5.133790, 1.285469, 0.551790, 0.009703]
turb_naca0012_p1c1.test_vals_aarch64 = [-5.114189, 1.285037, 0.406851, -0.043003]
test_list.append(turb_naca0012_p1c1)

Expand All @@ -314,7 +314,7 @@ def main():
turb_naca0012_p1c2.cfg_dir = "rans_uq/naca0012"
turb_naca0012_p1c2.cfg_file = "turb_NACA0012_uq_p1c2.cfg"
turb_naca0012_p1c2.test_iter = 10
turb_naca0012_p1c2.test_vals = [-5.548245, 1.236384, 0.381823, -0.050336]
turb_naca0012_p1c2.test_vals = [-5.553940, 1.237339, 0.427689, -0.034727]
turb_naca0012_p1c2.test_vals_aarch64 = [-5.548245, 1.236384, 0.381821, -0.050337]
test_list.append(turb_naca0012_p1c2)

Expand Down Expand Up @@ -432,15 +432,15 @@ def main():
cavity.cfg_dir = "moving_wall/cavity"
cavity.cfg_file = "lam_cavity.cfg"
cavity.test_iter = 25
cavity.test_vals = [-5.627869, -0.164403, 0.054734, 2.545856]
cavity.test_vals = [-5.627869, -0.164403, 0.054734, 2.545857]
test_list.append(cavity)

# Spinning cylinder
spinning_cylinder = TestCase('spinning_cylinder')
spinning_cylinder.cfg_dir = "moving_wall/spinning_cylinder"
spinning_cylinder.cfg_file = "spinning_cylinder.cfg"
spinning_cylinder.test_iter = 25
spinning_cylinder.test_vals = [-8.008048, -2.611074, 1.497289, 1.487468]
spinning_cylinder.test_vals = [-8.008023, -2.611064, 1.497308, 1.487483]
spinning_cylinder.test_vals_aarch64 = [-8.008023, -2.611064, 1.497308, 1.487483]
test_list.append(spinning_cylinder)

Expand All @@ -453,7 +453,7 @@ def main():
square_cylinder.cfg_dir = "unsteady/square_cylinder"
square_cylinder.cfg_file = "turb_square.cfg"
square_cylinder.test_iter = 3
square_cylinder.test_vals = [-2.560665, -1.175915, 0.062109, 1.399401, 2.220361, 1.399349, 2.218600, 0.000000]
square_cylinder.test_vals = [-2.560678, -1.175981, 0.062203, 1.399351, 2.219223, 1.399297, 2.217470, 0.000000]
square_cylinder.unsteady = True
test_list.append(square_cylinder)

Expand Down Expand Up @@ -512,7 +512,7 @@ def main():
unst_inc_turb_naca0015_sa.cfg_dir = "unsteady/pitching_naca0015_rans_inc"
unst_inc_turb_naca0015_sa.cfg_file = "config_incomp_turb_sa.cfg"
unst_inc_turb_naca0015_sa.test_iter = 1
unst_inc_turb_naca0015_sa.test_vals = [-3.008629, -6.888996, 1.435193, 0.433537]
unst_inc_turb_naca0015_sa.test_vals = [-3.008629, -6.889003, 1.435193, 0.433537]
unst_inc_turb_naca0015_sa.unsteady = True
test_list.append(unst_inc_turb_naca0015_sa)

Expand Down Expand Up @@ -564,15 +564,15 @@ def main():
axial_stage2D.cfg_dir = "turbomachinery/axial_stage_2D"
axial_stage2D.cfg_file = "Axial_stage2D.cfg"
axial_stage2D.test_iter = 20
axial_stage2D.test_vals = [1.084446, 1.526931, -2.895084, 2.607567, -2.479664, 3.063778, 106380.000000, 106380.000000, 5.733600, 64.749000]
axial_stage2D.test_vals = [1.084448, 1.526930, -2.895083, 2.607569, -2.479664, 3.063779, 106380.000000, 106380.000000, 5.733600, 64.737000]
test_list.append(axial_stage2D)

# 2D transonic stator restart
transonic_stator_restart = TestCase('transonic_stator_restart')
transonic_stator_restart.cfg_dir = "turbomachinery/transonic_stator_2D"
transonic_stator_restart.cfg_file = "transonic_stator_restart.cfg"
transonic_stator_restart.test_iter = 20
transonic_stator_restart.test_vals = [-4.442508, -2.561367, -2.165776, 1.652753, -1.355494, 3.172717, -471620.000000, 94.843000, -0.043825]
transonic_stator_restart.test_vals = [-4.442510, -2.561369, -2.165778, 1.652750, -1.355494, 3.172711, -471620.000000, 94.843000, -0.043825]
transonic_stator_restart.test_vals_aarch64 = [-4.442510, -2.561369, -2.165778, 1.652750, -1.355494, 3.172712, -471620.000000, 94.843000, -0.043825]
test_list.append(transonic_stator_restart)

Expand All @@ -594,7 +594,7 @@ def main():
uniform_flow.cfg_dir = "sliding_interface/uniform_flow"
uniform_flow.cfg_file = "uniform_NN.cfg"
uniform_flow.test_iter = 5
uniform_flow.test_vals = [5.000000, 0.000000, -0.195002, -10.624444]
uniform_flow.test_vals = [5.000000, 0.000000, -0.195002, -10.624451]
uniform_flow.unsteady = True
uniform_flow.multizone = True
test_list.append(uniform_flow)
Expand All @@ -604,7 +604,7 @@ def main():
channel_2D.cfg_dir = "sliding_interface/channel_2D"
channel_2D.cfg_file = "channel_2D_WA.cfg"
channel_2D.test_iter = 2
channel_2D.test_vals = [2.000000, 0.000000, 0.464907, 0.348052, 0.397452]
channel_2D.test_vals = [2.000000, 0.000000, 0.464906, 0.348048, 0.397471]
channel_2D.unsteady = True
channel_2D.multizone = True
test_list.append(channel_2D)
Expand All @@ -614,7 +614,7 @@ def main():
channel_3D.cfg_dir = "sliding_interface/channel_3D"
channel_3D.cfg_file = "channel_3D_WA.cfg"
channel_3D.test_iter = 2
channel_3D.test_vals = [2.000000, 0.000000, 0.629106, 0.524901, 0.422380]
channel_3D.test_vals = [2.000000, 0.000000, 0.629106, 0.524897, 0.422366]
channel_3D.test_vals_aarch64 = [2.000000, 0.000000, 0.629112, 0.524948, 0.422396]
channel_3D.unsteady = True
channel_3D.multizone = True
Expand Down Expand Up @@ -697,7 +697,7 @@ def main():
fsi2d.cfg_dir = "fea_fsi/WallChannel_2d"
fsi2d.cfg_file = "configFSI.cfg"
fsi2d.test_iter = 4
fsi2d.test_vals = [4.000000, 0.000000, -3.726029, -4.277769]
fsi2d.test_vals = [4.000000, 0.000000, -3.726028, -4.277769]
fsi2d.multizone= True
fsi2d.unsteady = True
fsi2d.enabled_with_tsan = False
Expand All @@ -708,7 +708,7 @@ def main():
dyn_fsi.cfg_dir = "fea_fsi/dyn_fsi"
dyn_fsi.cfg_file = "config.cfg"
dyn_fsi.test_iter = 4
dyn_fsi.test_vals = [-4.330725, -4.152983, 0.000000, 103.000000]
dyn_fsi.test_vals = [-4.330725, -4.152983, 0.000000, 102.000000]
dyn_fsi.multizone = True
dyn_fsi.unsteady = True
test_list.append(dyn_fsi)
Expand All @@ -732,7 +732,7 @@ def main():
mms_fvm_ns.cfg_dir = "mms/fvm_navierstokes"
mms_fvm_ns.cfg_file = "lam_mms_roe.cfg"
mms_fvm_ns.test_iter = 20
mms_fvm_ns.test_vals = [-2.808514, 2.152655, 0.000000, 0.000000]
mms_fvm_ns.test_vals = [-2.808514, 2.152654, 0.000000, 0.000000]
test_list.append(mms_fvm_ns)

# FVM, incompressible, euler
Expand Down
Loading
Loading