Skip to content

Commit 1335c29

Browse files
bmunguiabmunguia
committed
Address numerics comments
- Single code path regardless of kappa - Remove repeated formulation, and clarify scaling
1 parent c9910dd commit 1335c29

File tree

3 files changed

+37
-58
lines changed

3 files changed

+37
-58
lines changed

‎SU2_CFD/include/limiters/CLimiterDetails.hpp‎

Lines changed: 7 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -75,18 +75,16 @@ struct LimiterHelpers
7575
FORCEINLINE static Type umusclProjection(const Type& grad_proj, const Type& delta, const Type& kappa)
7676
{
7777
/*-------------------------------------------------------------------*/
78-
/*--- Reconstruction will be V_L = V_i + 0.25 * dV_ij^kappa. ---*/
79-
/*--- Scale dV_ij^cent by 0.5 here and return 0.5 * dV_ij^kappa ---*/
80-
/*--- for backward compatibility (i.e. when kappa==0 and this ---*/
81-
/*--- function isn't called, the reconstruction will be ---*/
82-
/*--- V_L = V_i + 0.5 * proj, otherwise V_L = V_i + 0.5 * blend). ---*/
78+
/*--- The MUSCL kappa-scheme reconstruction is typically written: ---*/
79+
/*--- V_L = V_i + 0.25 * dV_ij^kap, where ---*/
80+
/*--- dV_ij^kap = (1-kappa) dV_ij^upw + (1+kappa) dV_ij^cen, ---*/
81+
/*--- dV_ij^cen = V_j - V_i, ---*/
82+
/*--- dV_ij^upw = 2 grad(Vi) dot vector_ij - dV_ij^cen. ---*/
83+
/*--- To maintain proper scaling for edge limiters, the result of ---*/
84+
/*--- this function is 0.5 * dV_ij^kap. ---*/
8385
/*-------------------------------------------------------------------*/
8486
const Type cent = 0.5 * delta;
85-
86-
/*--- Upwind difference: dV_ij^upw = 2 grad(Vi) dot vector_ij - dV_ij^cent ---*/
8787
const Type upw = grad_proj - cent;
88-
89-
/*--- Blended difference: dV_ij^kappa = (1-kappa)dV_ij^upw + (1+kappa)dV_ij^cent ---*/
9088
return (1.0 - kappa) * upw + (1.0 + kappa) * cent;
9189
}
9290

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

Lines changed: 29 additions & 46 deletions
Original file line numberDiff line numberDiff line change
@@ -43,49 +43,39 @@ FORCEINLINE Double umusclProjection(Double grad_proj,
4343
Double delta,
4444
Double kappa) {
4545
/*-------------------------------------------------------------------*/
46-
/*--- Reconstruction will be V_L = V_i + 0.25 * dV_ij^kappa. ---*/
47-
/*--- Scale dV_ij^cent by 0.5 here and return 0.5 * dV_ij^kappa ---*/
48-
/*--- for backward compatibility (i.e. when kappa==0 and this ---*/
49-
/*--- function isn't called, the reconstruction will be ---*/
50-
/*--- V_L = V_i + 0.5 * proj, otherwise V_L = V_i + 0.5 * blend). ---*/
46+
/*--- The MUSCL kappa-scheme reconstruction is typically written: ---*/
47+
/*--- V_L = V_i + 0.25 * dV_ij^kap, where ---*/
48+
/*--- dV_ij^kap = (1-kappa) dV_ij^upw + (1+kappa) dV_ij^cen, ---*/
49+
/*--- dV_ij^cen = V_j - V_i, ---*/
50+
/*--- dV_ij^upw = 2 grad(Vi) dot vector_ij - dV_ij^cen. ---*/
51+
/*--- To maintain proper scaling for edge limiters, the result of ---*/
52+
/*--- this function is 0.5 * dV_ij^kap. ---*/
5153
/*-------------------------------------------------------------------*/
5254
const Double cent = 0.5 * delta;
53-
54-
/*--- Upwind difference: dV_ij^upw = 2 grad(Vi) dot vector_ij - dV_ij^cent ---*/
5555
const Double upw = grad_proj - cent;
56-
57-
/*--- Blended difference: dV_ij^kappa = (1-kappa)dV_ij^upw + (1+kappa)dV_ij^cent ---*/
5856
return (1.0 - kappa) * upw + (1.0 + kappa) * cent;
5957
}
6058

6159
/*!
6260
* \brief MUSCL reconstruction of the specified variable.
61+
* \note The result should be halved when added to i (or subtracted from j).
6362
* \param[in] grad_i - Gradient vector at point i.
64-
* \param[in] delta - Centered difference: V_j - V_i.
6563
* \param[in] vector_ij - Distance vector from i to j.
64+
* \param[in] delta - Centered difference: V_j - V_i.
6665
* \param[in] iVar - Variable index.
67-
* \param[in] umuscl - If true, use U-MUSCL blended difference.
6866
* \param[in] kappa - Blending coefficient.
6967
* \param[in] relax - Newton-Krylov relaxation.
7068
* \return Variable reconstructed from point i.
7169
*/
7270
template<class GradType, size_t nDim>
7371
FORCEINLINE Double musclReconstruction(const GradType& grad,
74-
const Double delta,
7572
const VectorDbl<nDim>& vector_ij,
73+
const Double delta,
7674
size_t iVar,
77-
bool umuscl,
7875
Double kappa,
7976
Double relax) {
80-
/*--- Gradient projection ---*/
81-
Double proj = dot(grad[iVar], vector_ij);
82-
83-
if (umuscl) {
84-
/*--- Blended difference: dV_ij^kappa = (1-kappa)dV_ij^upw + (1+kappa)dV_ij^cent ---*/
85-
proj = umusclProjection(proj, delta, kappa);
86-
}
87-
88-
return relax * proj;
77+
const Double proj = dot(grad[iVar], vector_ij);
78+
return relax * umusclProjection(proj, delta, kappa);
8979
}
9080

9181
/*!
@@ -97,7 +87,6 @@ FORCEINLINE void musclUnlimited(Int iPoint,
9787
const VectorDbl<nDim>& vector_ij,
9888
const Gradient_t& gradient,
9989
CPair<VarType>& V,
100-
bool umuscl,
10190
Double kappa,
10291
Double relax) {
10392
constexpr auto nVarGrad = nVarGrad_ > 0 ? nVarGrad_ : VarType::nVar;
@@ -106,16 +95,14 @@ FORCEINLINE void musclUnlimited(Int iPoint,
10695
auto grad_j = gatherVariables<nVarGrad,nDim>(jPoint, gradient);
10796

10897
for (size_t iVar = 0; iVar < nVarGrad; ++iVar) {
109-
/*--- Centered difference, only needed for U-MUSCL projection ---*/
110-
Double delta_ij = 0.0;
111-
if (umuscl)
112-
delta_ij = V.j.all(iVar) - V.i.all(iVar);
98+
/*--- Centered difference, needed for U-MUSCL projection ---*/
99+
const Double delta_ij = V.j.all(iVar) - V.i.all(iVar);
113100

114101
/*--- U-MUSCL reconstructed variables ---*/
115-
const Double proj_i = musclReconstruction(grad_i, delta_ij, vector_ij, iVar, umuscl, kappa, relax);
116-
const Double proj_j = musclReconstruction(grad_j, delta_ij, vector_ij, iVar, umuscl, kappa, relax);
102+
const Double proj_i = musclReconstruction(grad_i, vector_ij, delta_ij, iVar, kappa, relax);
103+
const Double proj_j = musclReconstruction(grad_j, vector_ij, delta_ij, iVar, kappa, relax);
117104

118-
/*--- Apply reconstruction: V_L = V_i + 0.25 * dV_ij^kappa ---*/
105+
/*--- Apply reconstruction: V_L = V_i + 0.5 * dV_ij^kap ---*/
119106
V.i.all(iVar) += 0.5 * proj_i;
120107
V.j.all(iVar) -= 0.5 * proj_j;
121108
}
@@ -131,7 +118,6 @@ FORCEINLINE void musclPointLimited(Int iPoint,
131118
const Limiter_t& limiter,
132119
const Gradient_t& gradient,
133120
CPair<VarType>& V,
134-
bool umuscl,
135121
Double kappa,
136122
Double relax) {
137123
constexpr auto nVarGrad = nVarGrad_ > 0 ? nVarGrad_ : VarType::nVar;
@@ -143,16 +129,14 @@ FORCEINLINE void musclPointLimited(Int iPoint,
143129
auto grad_j = gatherVariables<nVarGrad,nDim>(jPoint, gradient);
144130

145131
for (size_t iVar = 0; iVar < nVarGrad; ++iVar) {
146-
/*--- Centered difference, only needed for U-MUSCL projection ---*/
147-
Double delta_ij = 0.0;
148-
if (umuscl)
149-
delta_ij = V.j.all(iVar) - V.i.all(iVar);
132+
/*--- Centered difference, needed for U-MUSCL projection ---*/
133+
const Double delta_ij = V.j.all(iVar) - V.i.all(iVar);
150134

151135
/*--- U-MUSCL reconstructed variables ---*/
152-
const Double proj_i = musclReconstruction(grad_i, delta_ij, vector_ij, iVar, umuscl, kappa, relax);
153-
const Double proj_j = musclReconstruction(grad_j, delta_ij, vector_ij, iVar, umuscl, kappa, relax);
136+
const Double proj_i = musclReconstruction(grad_i, vector_ij, delta_ij, iVar, kappa, relax);
137+
const Double proj_j = musclReconstruction(grad_j, vector_ij, delta_ij, iVar, kappa, relax);
154138

155-
/*--- Apply reconstruction: V_L = V_i + 0.25 * lim * dV_ij^kappa ---*/
139+
/*--- Apply reconstruction: V_L = V_i + 0.5 * lim * dV_ij^kap ---*/
156140
V.i.all(iVar) += 0.5 * lim_i(iVar) * proj_i;
157141
V.j.all(iVar) -= 0.5 * lim_j(iVar) * proj_j;
158142
}
@@ -167,7 +151,6 @@ FORCEINLINE void musclEdgeLimited(Int iPoint,
167151
const VectorDbl<nDim>& vector_ij,
168152
const Gradient_t& gradient,
169153
CPair<VarType>& V,
170-
bool umuscl,
171154
Double kappa,
172155
Double relax) {
173156
constexpr auto nVarGrad = nVarGrad_ > 0 ? nVarGrad_ : VarType::nVar;
@@ -181,13 +164,14 @@ FORCEINLINE void musclEdgeLimited(Int iPoint,
181164
const Double delta_ij_2 = pow(delta_ij, 2) + 1e-6;
182165

183166
/*--- U-MUSCL reconstructed variables ---*/
184-
const Double proj_i = musclReconstruction(grad_i, delta_ij, vector_ij, iVar, umuscl, kappa, relax);
185-
const Double proj_j = musclReconstruction(grad_j, delta_ij, vector_ij, iVar, umuscl, kappa, relax);
167+
const Double proj_i = musclReconstruction(grad_i, vector_ij, delta_ij, iVar, kappa, relax);
168+
const Double proj_j = musclReconstruction(grad_j, vector_ij, delta_ij, iVar, kappa, relax);
186169

187170
/// TODO: Customize the limiter function.
188171
const Double lim_i = (delta_ij_2 + proj_i*delta_ij) / (pow(proj_i,2) + delta_ij_2);
189172
const Double lim_j = (delta_ij_2 + proj_j*delta_ij) / (pow(proj_j,2) + delta_ij_2);
190173

174+
/*--- Apply reconstruction: V_L = V_i + 0.5 * lim * dV_ij^kap ---*/
191175
V.i.all(iVar) += 0.5 * lim_i * proj_i;
192176
V.j.all(iVar) -= 0.5 * lim_j * proj_j;
193177
}
@@ -200,7 +184,6 @@ FORCEINLINE void musclEdgeLimited(Int iPoint,
200184
* \param[in] gamma - Heat capacity ratio.
201185
* \param[in] gasConst - Specific gas constant.
202186
* \param[in] muscl - If true, reconstruct, else simply fetch.
203-
* \param[in] umuscl - If true, reconstruct with a blended difference.
204187
* \param[in] kappa - Blending coefficient for MUSCL reconstruction.
205188
* \param[in] relax - Newton-Krylov relaxation.
206189
* \param[in] limiterType - Type of flux limiter.
@@ -213,7 +196,7 @@ template<class ReconVarType, class PrimVarType, size_t nDim, class VariableType>
213196
FORCEINLINE CPair<ReconVarType> reconstructPrimitives(Int iEdge, Int iPoint, Int jPoint,
214197
const su2double& gamma,
215198
const su2double& gasConst,
216-
bool muscl, bool umuscl,
199+
bool muscl,
217200
const su2double& kappa,
218201
const su2double& relax,
219202
LIMITER limiterType,
@@ -237,13 +220,13 @@ FORCEINLINE CPair<ReconVarType> reconstructPrimitives(Int iEdge, Int iPoint, Int
237220
constexpr auto nVarGrad = ReconVarType::nVar - 2;
238221
switch (limiterType) {
239222
case LIMITER::NONE:
240-
musclUnlimited<nVarGrad>(iPoint, jPoint, vector_ij, gradients, V, umuscl, kappa, relax);
223+
musclUnlimited<nVarGrad>(iPoint, jPoint, vector_ij, gradients, V, kappa, relax);
241224
break;
242225
case LIMITER::VAN_ALBADA_EDGE:
243-
musclEdgeLimited<nVarGrad>(iPoint, jPoint, vector_ij, gradients, V, umuscl, kappa, relax);
226+
musclEdgeLimited<nVarGrad>(iPoint, jPoint, vector_ij, gradients, V, kappa, relax);
244227
break;
245228
default:
246-
musclPointLimited<nVarGrad>(iPoint, jPoint, vector_ij, limiters, gradients, V, umuscl, kappa, relax);
229+
musclPointLimited<nVarGrad>(iPoint, jPoint, vector_ij, limiters, gradients, V, kappa, relax);
247230
break;
248231
}
249232
V.i.density() = V.i.pressure() / (gasConst * V.i.temperature());

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

Lines changed: 1 addition & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -61,7 +61,6 @@ class CRoeBase : public Base {
6161
const bool finestGrid;
6262
const bool dynamicGrid;
6363
const bool muscl;
64-
const bool umuscl;
6564
const su2double umusclKappa;
6665
const LIMITER typeLimiter;
6766

@@ -77,7 +76,6 @@ class CRoeBase : public Base {
7776
finestGrid(iMesh == MESH_0),
7877
dynamicGrid(config.GetDynamic_Grid()),
7978
muscl(finestGrid && config.GetMUSCL_Flow()),
80-
umuscl(muscl && config.GetMUSCL_Kappa_Flow() != 0.0),
8179
umusclKappa(config.GetMUSCL_Kappa_Flow()),
8280
typeLimiter(config.GetKind_SlopeLimit_Flow()) {
8381
}
@@ -125,7 +123,7 @@ class CRoeBase : public Base {
125123

126124
/*--- Recompute density and enthalpy instead of reconstructing. ---*/
127125
auto V = reconstructPrimitives<CCompressiblePrimitives<nDim,nPrimVarGrad> >(
128-
iEdge, iPoint, jPoint, gamma, gasConst, muscl, umuscl, umusclKappa, nkRelax, typeLimiter, V1st, vector_ij, solution);
126+
iEdge, iPoint, jPoint, gamma, gasConst, muscl, umusclKappa, nkRelax, typeLimiter, V1st, vector_ij, solution);
129127

130128
/*--- Compute conservative variables. ---*/
131129

0 commit comments

Comments
 (0)