Skip to content

Commit 527b64a

Browse files
bmunguiabmunguia
committed
Apply relaxation before limiter
1 parent 176fdc7 commit 527b64a

File tree

4 files changed

+72
-61
lines changed

4 files changed

+72
-61
lines changed

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

Lines changed: 60 additions & 49 deletions
Original file line numberDiff line numberDiff line change
@@ -34,17 +34,17 @@
3434

3535
/*!
3636
* \brief Blended difference for U-MUSCL reconstruction.
37-
* \param[in] grad_proj - Gradient projection at point i: dot(grad_i, vector_ij)
38-
* \param[in] delta - Centered difference: V_j - V_i
39-
* \param[in] kappa - Blending parameter
40-
* \return Blended difference for reconstruction
37+
* \param[in] grad_proj - Gradient projection at point i: dot(grad_i, vector_ij).
38+
* \param[in] delta - Centered difference: V_j - V_i.
39+
* \param[in] kappa - Blending parameter.
40+
* \return Blended difference for reconstruction from point i.
4141
*/
4242
FORCEINLINE Double umusclProjection(Double grad_proj,
4343
Double delta,
4444
Double kappa) {
4545
/*-------------------------------------------------------------------*/
4646
/*--- 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 ---*/
47+
/*--- Scale dV_ij^cent by 0.5 here and return 0.5 * dV_ij^kappa ---*/
4848
/*--- for backward compatibility (i.e. when kappa==0 and this ---*/
4949
/*--- function isn't called, the reconstruction will be ---*/
5050
/*--- V_L = V_i + 0.5 * proj, otherwise V_L = V_i + 0.5 * blend). ---*/
@@ -58,6 +58,36 @@ FORCEINLINE Double umusclProjection(Double grad_proj,
5858
return (1.0 - kappa) * upw + (1.0 + kappa) * cent;
5959
}
6060

61+
/*!
62+
* \brief MUSCL reconstruction of the specified variable.
63+
* \param[in] grad_i - Gradient vector at point i.
64+
* \param[in] delta - Centered difference: V_j - V_i.
65+
* \param[in] vector_ij - Distance vector from i to j.
66+
* \param[in] iVar - Variable index.
67+
* \param[in] umuscl - If true, use U-MUSCL blended difference.
68+
* \param[in] kappa - Blending coefficient.
69+
* \param[in] relax - Newton-Krylov relaxation.
70+
* \return Variable reconstructed from point i.
71+
*/
72+
template<class GradType, size_t nDim>
73+
FORCEINLINE Double musclReconstruction(const GradType& grad,
74+
const Double delta,
75+
const VectorDbl<nDim>& vector_ij,
76+
size_t iVar,
77+
bool umuscl,
78+
Double kappa,
79+
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;
89+
}
90+
6191
/*!
6292
* \brief Unlimited reconstruction.
6393
*/
@@ -72,28 +102,22 @@ FORCEINLINE void musclUnlimited(Int iPoint,
72102
Double relax) {
73103
constexpr auto nVarGrad = nVarGrad_ > 0 ? nVarGrad_ : VarType::nVar;
74104

75-
const auto scale = 0.5 * relax;
76-
77105
auto grad_i = gatherVariables<nVarGrad,nDim>(iPoint, gradient);
78106
auto grad_j = gatherVariables<nVarGrad,nDim>(jPoint, gradient);
79107

80108
for (size_t iVar = 0; iVar < nVarGrad; ++iVar) {
81-
/*--- Gradient projections ---*/
82-
Double proj_i = dot(grad_i[iVar], vector_ij);
83-
Double proj_j = dot(grad_j[iVar], vector_ij);
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);
84113

85-
if (umuscl) {
86-
/*--- Centered difference: dV_ij^cent = V_j - V_i ---*/
87-
const Double cent = V.j.all(iVar) - V.i.all(iVar);
88-
89-
/*--- Blended difference: dV_ij^kappa = (1-kappa)dV_ij^upw + (1+kappa)dV_ij^cent ---*/
90-
proj_i = umusclProjection(proj_i, cent, kappa);
91-
proj_j = umusclProjection(proj_j, cent, kappa);
92-
}
114+
/*--- 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);
93117

94118
/*--- Apply reconstruction: V_L = V_i + 0.25 * dV_ij^kappa ---*/
95-
V.i.all(iVar) += scale * proj_i;
96-
V.j.all(iVar) -= scale * proj_j;
119+
V.i.all(iVar) += 0.5 * proj_i;
120+
V.j.all(iVar) -= 0.5 * proj_j;
97121
}
98122
}
99123

@@ -112,31 +136,25 @@ FORCEINLINE void musclPointLimited(Int iPoint,
112136
Double relax) {
113137
constexpr auto nVarGrad = nVarGrad_ > 0 ? nVarGrad_ : VarType::nVar;
114138

115-
const auto scale = 0.5 * relax;
116-
117139
auto lim_i = gatherVariables<nVarGrad>(iPoint, limiter);
118140
auto lim_j = gatherVariables<nVarGrad>(jPoint, limiter);
119141

120142
auto grad_i = gatherVariables<nVarGrad,nDim>(iPoint, gradient);
121143
auto grad_j = gatherVariables<nVarGrad,nDim>(jPoint, gradient);
122144

123145
for (size_t iVar = 0; iVar < nVarGrad; ++iVar) {
124-
/*--- Gradient projections ---*/
125-
Double proj_i = dot(grad_i[iVar], vector_ij);
126-
Double proj_j = dot(grad_j[iVar], vector_ij);
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);
127150

128-
if (umuscl) {
129-
/*--- Centered difference: dV_ij^cent = Vj - Vi ---*/
130-
const Double cent = V.j.all(iVar) - V.i.all(iVar);
131-
132-
/*--- Blended difference: dV_ij^kappa = (1-kappa)dV_ij^upw + (1+kappa)dV_ij^cent ---*/
133-
proj_i = umusclProjection(proj_i, cent, kappa);
134-
proj_j = umusclProjection(proj_j, cent, kappa);
135-
}
151+
/*--- 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);
136154

137155
/*--- Apply reconstruction: V_L = V_i + 0.25 * lim * dV_ij^kappa ---*/
138-
V.i.all(iVar) += scale * lim_i(iVar) * proj_i;
139-
V.j.all(iVar) -= scale * lim_j(iVar) * proj_j;
156+
V.i.all(iVar) += 0.5 * lim_i(iVar) * proj_i;
157+
V.j.all(iVar) -= 0.5 * lim_j(iVar) * proj_j;
140158
}
141159
}
142160

@@ -154,31 +172,24 @@ FORCEINLINE void musclEdgeLimited(Int iPoint,
154172
Double relax) {
155173
constexpr auto nVarGrad = nVarGrad_ > 0 ? nVarGrad_ : VarType::nVar;
156174

157-
const auto scale = 0.5 * relax;
158-
159175
auto grad_i = gatherVariables<nVarGrad,nDim>(iPoint, gradient);
160176
auto grad_j = gatherVariables<nVarGrad,nDim>(jPoint, gradient);
161177

162178
for (size_t iVar = 0; iVar < nVarGrad; ++iVar) {
163-
Double proj_i = dot(grad_i[iVar], vector_ij);
164-
Double proj_j = dot(grad_j[iVar], vector_ij);
179+
/*--- Centered difference, needed for U-MUSCL projection and limiter ---*/
165180
const Double delta_ij = V.j.all(iVar) - V.i.all(iVar);
166181
const Double delta_ij_2 = pow(delta_ij, 2) + 1e-6;
167182

168-
if (umuscl) {
169-
/*--- Centered difference: dV_ij^cent = V_j - V_i ---*/
170-
const Double cent = V.j.all(iVar) - V.i.all(iVar);
171-
172-
/*--- Blended difference: dV_ij^kappa = (1-kappa)dV_ij^upw + (1+kappa)dV_ij^cent ---*/
173-
proj_i = umusclProjection(proj_i, cent, kappa);
174-
proj_j = umusclProjection(proj_j, cent, kappa);
175-
}
183+
/*--- 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);
176186

177187
/// TODO: Customize the limiter function.
178188
const Double lim_i = (delta_ij_2 + proj_i*delta_ij) / (pow(proj_i,2) + delta_ij_2);
179189
const Double lim_j = (delta_ij_2 + proj_j*delta_ij) / (pow(proj_j,2) + delta_ij_2);
180-
V.i.all(iVar) += scale * lim_i * proj_i;
181-
V.j.all(iVar) -= scale * lim_j * proj_j;
190+
191+
V.i.all(iVar) += 0.5 * lim_i * proj_i;
192+
V.j.all(iVar) -= 0.5 * lim_j * proj_j;
182193
}
183194
}
184195

‎SU2_CFD/src/solvers/CEulerSolver.cpp‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -1885,8 +1885,8 @@ void CEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_contain
18851885
if (umuscl || van_albada)
18861886
V_ij = V_j[iVar] - V_i[iVar];
18871887

1888-
const su2double Project_Grad_i = MUSCL_Reconstruction(Gradient_i[iVar], V_ij, Vector_ij, umuscl, kappa);
1889-
const su2double Project_Grad_j = MUSCL_Reconstruction(Gradient_j[iVar], V_ij, Vector_ij, umuscl, kappa);
1888+
const su2double Project_Grad_i = nkRelax * MUSCL_Reconstruction(Gradient_i[iVar], V_ij, Vector_ij, umuscl, kappa);
1889+
const su2double Project_Grad_j = nkRelax * MUSCL_Reconstruction(Gradient_j[iVar], V_ij, Vector_ij, umuscl, kappa);
18901890

18911891
su2double lim_i = 1.0;
18921892
su2double lim_j = 1.0;
@@ -1899,8 +1899,8 @@ void CEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_contain
18991899
lim_j = nodes->GetLimiter_Primitive(jPoint, iVar);
19001900
}
19011901

1902-
Primitive_i[iVar] = V_i[iVar] + 0.5 * nkRelax * lim_i * Project_Grad_i;
1903-
Primitive_j[iVar] = V_j[iVar] - 0.5 * nkRelax * lim_j * Project_Grad_j;
1902+
Primitive_i[iVar] = V_i[iVar] + 0.5 * lim_i * Project_Grad_i;
1903+
Primitive_j[iVar] = V_j[iVar] - 0.5 * lim_j * Project_Grad_j;
19041904

19051905
}
19061906

‎SU2_CFD/src/solvers/CIncEulerSolver.cpp‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -1283,8 +1283,8 @@ void CIncEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_cont
12831283
if (umuscl || van_albada)
12841284
V_ij = V_j[iVar] - V_i[iVar];
12851285

1286-
const su2double Project_Grad_i = MUSCL_Reconstruction(Gradient_i[iVar], V_ij, Vector_ij, umuscl, kappa);
1287-
const su2double Project_Grad_j = MUSCL_Reconstruction(Gradient_j[iVar], V_ij, Vector_ij, umuscl, kappa);
1286+
const su2double Project_Grad_i = nkRelax * MUSCL_Reconstruction(Gradient_i[iVar], V_ij, Vector_ij, umuscl, kappa);
1287+
const su2double Project_Grad_j = nkRelax * MUSCL_Reconstruction(Gradient_j[iVar], V_ij, Vector_ij, umuscl, kappa);
12881288

12891289
su2double lim_i = 1.0;
12901290
su2double lim_j = 1.0;
@@ -1297,8 +1297,8 @@ void CIncEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_cont
12971297
lim_j = nodes->GetLimiter_Primitive(jPoint, iVar);
12981298
}
12991299

1300-
Primitive_i[iVar] = V_i[iVar] + 0.5 * nkRelax * lim_i * Project_Grad_i;
1301-
Primitive_j[iVar] = V_j[iVar] - 0.5 * nkRelax * lim_j * Project_Grad_j;
1300+
Primitive_i[iVar] = V_i[iVar] + 0.5 * lim_i * Project_Grad_i;
1301+
Primitive_j[iVar] = V_j[iVar] - 0.5 * lim_j * Project_Grad_j;
13021302
}
13031303

13041304
for (auto iVar = nPrimVarGrad; iVar < nPrimVar; iVar++) {

‎SU2_CFD/src/solvers/CNEMOEulerSolver.cpp‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -546,8 +546,8 @@ void CNEMOEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_con
546546
if (umuscl || van_albada)
547547
V_ij = V_j[iVar] - V_i[iVar];
548548

549-
Project_Grad_i[iVar] = MUSCL_Reconstruction(Gradient_i[iVar], V_ij, Vector_ij, umuscl, kappa);
550-
Project_Grad_j[iVar] = MUSCL_Reconstruction(Gradient_j[iVar], V_ij, Vector_ij, umuscl, kappa);
549+
Project_Grad_i[iVar] = nkRelax * MUSCL_Reconstruction(Gradient_i[iVar], V_ij, Vector_ij, umuscl, kappa);
550+
Project_Grad_j[iVar] = nkRelax * MUSCL_Reconstruction(Gradient_j[iVar], V_ij, Vector_ij, umuscl, kappa);
551551

552552
if (limiter) {
553553
if (van_albada) {
@@ -566,8 +566,8 @@ void CNEMOEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_con
566566
su2double lim_ij = min(lim_i, lim_j);
567567

568568
for (auto iVar = 0ul; iVar < nPrimVarGrad; iVar++) {
569-
Primitive_i[iVar] = V_i[iVar] + 0.5 * nkRelax * lim_ij * Project_Grad_i[iVar];
570-
Primitive_j[iVar] = V_j[iVar] - 0.5 * nkRelax * lim_ij * Project_Grad_j[iVar];
569+
Primitive_i[iVar] = V_i[iVar] + 0.5 * lim_ij * Project_Grad_i[iVar];
570+
Primitive_j[iVar] = V_j[iVar] - 0.5 * lim_ij * Project_Grad_j[iVar];
571571
}
572572

573573
/*--- Check for non-physical solutions after reconstruction. If found, use the

0 commit comments

Comments
 (0)