Skip to content

Commit 74e30cf

Browse files
pcarruscagpcarruscag
authored
Merge pull request #2608 from su2code/pedro/fix_muscl_enthalpy
Make reconstructed rhoE more stable
2 parents 948a5ea + 610c66c commit 74e30cf

File tree

11 files changed

+128
-119
lines changed

11 files changed

+128
-119
lines changed

‎SU2_CFD/include/numerics_simd/flow/convection/common.hpp‎

Lines changed: 15 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -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
/*!
@@ -216,7 +214,7 @@ FORCEINLINE CPair<ReconVarType> reconstructPrimitives(Int iEdge, Int iPoint, Int
216214
}
217215

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

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

239248
/*--- Detect a non-physical reconstruction based on negative pressure or density. ---*/
240249
const Double neg_p_or_rho = fmax(fmin(V.i.pressure(), V.j.pressure()) < 0.0,

‎TestCases/disc_adj_fsi/dyn_fsi/grad_dv.opt.ref‎

Lines changed: 8 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -1,9 +1,9 @@
11
INDEX GRAD
2-
0 -4.570869180445024e-04
3-
1 -2.401466726069483e-04
4-
2 -9.134698219534478e-05
5-
3 -1.628086953114332e-05
6-
4 -1.741052432426715e-05
7-
5 -9.462489852612668e-05
8-
6 -2.452466275723706e-04
9-
7 -4.635483655821636e-04
2+
0 -4.570869126273463e-04
3+
1 -2.401466681907557e-04
4+
2 -9.134697929286435e-05
5+
3 -1.628086861828660e-05
6+
4 -1.741052581631559e-05
7+
5 -9.462490278520256e-05
8+
6 -2.452466349703724e-04
9+
7 -4.635483765277457e-04

‎TestCases/hybrid_regression.py‎

Lines changed: 29 additions & 29 deletions
Original file line numberDiff line numberDiff line change
@@ -103,7 +103,7 @@ def main():
103103
flatplate.cfg_dir = "navierstokes/flatplate"
104104
flatplate.cfg_file = "lam_flatplate.cfg"
105105
flatplate.test_iter = 100
106-
flatplate.test_vals = [-7.680034, -2.207912, 0.001084, 0.036233, 2.361500, -2.325300, 0.000000, 0.000000]
106+
flatplate.test_vals = [-7.680046, -2.207922, 0.001084, 0.036233, 2.361500, -2.325300, 0.000000, 0.000000]
107107
test_list.append(flatplate)
108108

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

@@ -432,15 +432,15 @@ def main():
432432
cavity.cfg_dir = "moving_wall/cavity"
433433
cavity.cfg_file = "lam_cavity.cfg"
434434
cavity.test_iter = 25
435-
cavity.test_vals = [-5.627869, -0.164403, 0.054734, 2.545856]
435+
cavity.test_vals = [-5.627869, -0.164403, 0.054734, 2.545857]
436436
test_list.append(cavity)
437437

438438
# Spinning cylinder
439439
spinning_cylinder = TestCase('spinning_cylinder')
440440
spinning_cylinder.cfg_dir = "moving_wall/spinning_cylinder"
441441
spinning_cylinder.cfg_file = "spinning_cylinder.cfg"
442442
spinning_cylinder.test_iter = 25
443-
spinning_cylinder.test_vals = [-8.008048, -2.611074, 1.497289, 1.487468]
443+
spinning_cylinder.test_vals = [-8.008023, -2.611064, 1.497308, 1.487483]
444444
spinning_cylinder.test_vals_aarch64 = [-8.008023, -2.611064, 1.497308, 1.487483]
445445
test_list.append(spinning_cylinder)
446446

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

@@ -512,7 +512,7 @@ def main():
512512
unst_inc_turb_naca0015_sa.cfg_dir = "unsteady/pitching_naca0015_rans_inc"
513513
unst_inc_turb_naca0015_sa.cfg_file = "config_incomp_turb_sa.cfg"
514514
unst_inc_turb_naca0015_sa.test_iter = 1
515-
unst_inc_turb_naca0015_sa.test_vals = [-3.008629, -6.888996, 1.435193, 0.433537]
515+
unst_inc_turb_naca0015_sa.test_vals = [-3.008629, -6.889003, 1.435193, 0.433537]
516516
unst_inc_turb_naca0015_sa.unsteady = True
517517
test_list.append(unst_inc_turb_naca0015_sa)
518518

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

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

@@ -594,7 +594,7 @@ def main():
594594
uniform_flow.cfg_dir = "sliding_interface/uniform_flow"
595595
uniform_flow.cfg_file = "uniform_NN.cfg"
596596
uniform_flow.test_iter = 5
597-
uniform_flow.test_vals = [5.000000, 0.000000, -0.195002, -10.624444]
597+
uniform_flow.test_vals = [5.000000, 0.000000, -0.195002, -10.624451]
598598
uniform_flow.unsteady = True
599599
uniform_flow.multizone = True
600600
test_list.append(uniform_flow)
@@ -604,7 +604,7 @@ def main():
604604
channel_2D.cfg_dir = "sliding_interface/channel_2D"
605605
channel_2D.cfg_file = "channel_2D_WA.cfg"
606606
channel_2D.test_iter = 2
607-
channel_2D.test_vals = [2.000000, 0.000000, 0.464907, 0.348052, 0.397452]
607+
channel_2D.test_vals = [2.000000, 0.000000, 0.464906, 0.348048, 0.397471]
608608
channel_2D.unsteady = True
609609
channel_2D.multizone = True
610610
test_list.append(channel_2D)
@@ -614,7 +614,7 @@ def main():
614614
channel_3D.cfg_dir = "sliding_interface/channel_3D"
615615
channel_3D.cfg_file = "channel_3D_WA.cfg"
616616
channel_3D.test_iter = 2
617-
channel_3D.test_vals = [2.000000, 0.000000, 0.629106, 0.524901, 0.422380]
617+
channel_3D.test_vals = [2.000000, 0.000000, 0.629106, 0.524897, 0.422366]
618618
channel_3D.test_vals_aarch64 = [2.000000, 0.000000, 0.629112, 0.524948, 0.422396]
619619
channel_3D.unsteady = True
620620
channel_3D.multizone = True
@@ -697,7 +697,7 @@ def main():
697697
fsi2d.cfg_dir = "fea_fsi/WallChannel_2d"
698698
fsi2d.cfg_file = "configFSI.cfg"
699699
fsi2d.test_iter = 4
700-
fsi2d.test_vals = [4.000000, 0.000000, -3.726029, -4.277769]
700+
fsi2d.test_vals = [4.000000, 0.000000, -3.726028, -4.277769]
701701
fsi2d.multizone= True
702702
fsi2d.unsteady = True
703703
fsi2d.enabled_with_tsan = False
@@ -708,7 +708,7 @@ def main():
708708
dyn_fsi.cfg_dir = "fea_fsi/dyn_fsi"
709709
dyn_fsi.cfg_file = "config.cfg"
710710
dyn_fsi.test_iter = 4
711-
dyn_fsi.test_vals = [-4.330725, -4.152983, 0.000000, 103.000000]
711+
dyn_fsi.test_vals = [-4.330725, -4.152983, 0.000000, 102.000000]
712712
dyn_fsi.multizone = True
713713
dyn_fsi.unsteady = True
714714
test_list.append(dyn_fsi)
@@ -732,7 +732,7 @@ def main():
732732
mms_fvm_ns.cfg_dir = "mms/fvm_navierstokes"
733733
mms_fvm_ns.cfg_file = "lam_mms_roe.cfg"
734734
mms_fvm_ns.test_iter = 20
735-
mms_fvm_ns.test_vals = [-2.808514, 2.152655, 0.000000, 0.000000]
735+
mms_fvm_ns.test_vals = [-2.808514, 2.152654, 0.000000, 0.000000]
736736
test_list.append(mms_fvm_ns)
737737

738738
# FVM, incompressible, euler

0 commit comments

Comments
 (0)