Skip to content

Commit 3cdc4f7

Browse files
bmunguiabmunguia
committed
Extend U-MUSCL extrapolation to non-SIMD upwind residuals
1 parent ea18471 commit 3cdc4f7

File tree

9 files changed

+104
-69
lines changed

9 files changed

+104
-69
lines changed

‎Common/include/CConfig.hpp‎

Lines changed: 14 additions & 20 deletions
Original file line numberDiff line numberDiff line change
@@ -593,12 +593,12 @@ class CConfig {
593593
MUSCL_AdjFlow, /*!< \brief MUSCL scheme for the adj flow equations.*/
594594
MUSCL_AdjTurb, /*!< \brief MUSCL scheme for the adj turbulence equations.*/
595595
MUSCL_Species; /*!< \brief MUSCL scheme for the species equations.*/
596-
su2double MUSCL_Kappa, /*!< \brief Blending coefficient for MUSCL scheme. */
597-
MUSCL_Kappa_Flow, /*!< \brief Blending coefficient for MUSCL scheme for the flow equations.*/
598-
MUSCL_Kappa_Turb, /*!< \brief Blending coefficient for MUSCL scheme for the turbulence equations.*/
599-
MUSCL_Kappa_Heat, /*!< \brief Blending coefficient for MUSCL scheme for the (fvm) heat equation.*/
600-
MUSCL_Kappa_AdjFlow, /*!< \brief Blending coefficient for MUSCL scheme for the adj flow equations.*/
601-
MUSCL_Kappa_Species; /*!< \brief Blending coefficient for MUSCL scheme for the species equations.*/
596+
su2double MUSCL_Kappa, /*!< \brief Blending coefficient for U-MUSCL scheme. */
597+
MUSCL_Kappa_Flow, /*!< \brief Blending coefficient for U-MUSCL scheme for the flow equations.*/
598+
MUSCL_Kappa_Turb, /*!< \brief Blending coefficient for U-MUSCL scheme for the turbulence equations.*/
599+
MUSCL_Kappa_Heat, /*!< \brief Blending coefficient for U-MUSCL scheme for the (fvm) heat equation.*/
600+
MUSCL_Kappa_AdjFlow, /*!< \brief Blending coefficient for U-MUSCL scheme for the adj flow equations.*/
601+
MUSCL_Kappa_Species; /*!< \brief Blending coefficient for U-MUSCL scheme for the species equations.*/
602602
bool Use_Accurate_Jacobians; /*!< \brief Use numerically computed Jacobians for AUSM+up(2) and SLAU(2). */
603603
bool Use_Accurate_Turb_Jacobians; /*!< \brief Use numerically computed Jacobians for standard SA turbulence model. */
604604
bool EulerPersson; /*!< \brief Boolean to determine whether this is an Euler simulation with Persson shock capturing. */
@@ -4624,27 +4624,18 @@ class CConfig {
46244624

46254625
/*!
46264626
* \brief Get if the upwind scheme used MUSCL or not.
4627-
* \note This is the information that the code will use, the method will
4628-
* change in runtime depending of the specific equation (direct, adjoint,
4629-
* linearized) that is being solved.
46304627
* \return MUSCL scheme.
46314628
*/
46324629
bool GetMUSCL_Flow(void) const { return MUSCL_Flow; }
46334630

46344631
/*!
46354632
* \brief Get if the upwind scheme used MUSCL or not.
4636-
* \note This is the information that the code will use, the method will
4637-
* change in runtime depending of the specific equation (direct, adjoint,
4638-
* linearized) that is being solved.
46394633
* \return MUSCL scheme.
46404634
*/
46414635
bool GetMUSCL_Heat(void) const { return MUSCL_Heat; }
46424636

46434637
/*!
46444638
* \brief Get if the upwind scheme used MUSCL or not.
4645-
* \note This is the information that the code will use, the method will
4646-
* change in runtime depending of the specific equation (direct, adjoint,
4647-
* linearized) that is being solved.
46484639
* \return MUSCL scheme.
46494640
*/
46504641
bool GetMUSCL_Turb(void) const { return MUSCL_Turb; }
@@ -4657,21 +4648,24 @@ class CConfig {
46574648

46584649
/*!
46594650
* \brief Get if the upwind scheme used MUSCL or not.
4660-
* \note This is the information that the code will use, the method will
4661-
* change in runtime depending of the specific equation (direct, adjoint,
4662-
* linearized) that is being solved.
46634651
* \return MUSCL scheme.
46644652
*/
46654653
bool GetMUSCL_AdjFlow(void) const { return MUSCL_AdjFlow; }
46664654

46674655
/*!
46684656
* \brief Get if the upwind scheme used MUSCL or not.
4657+
* \return MUSCL scheme.
4658+
*/
4659+
bool GetMUSCL_AdjTurb(void) const { return MUSCL_AdjTurb; }
4660+
4661+
/*!
4662+
* \brief Get the blending coefficient for the U-MUSCL scheme.
46694663
* \note This is the information that the code will use, the method will
46704664
* change in runtime depending of the specific equation (direct, adjoint,
46714665
* linearized) that is being solved.
4672-
* \return MUSCL scheme.
4666+
* \return Blending coefficient for the U-MUSCL scheme.
46734667
*/
4674-
bool GetMUSCL_AdjTurb(void) const { return MUSCL_AdjTurb; }
4668+
su2double GetMUSCL_Kappa(void) const { return MUSCL_Kappa; }
46754669

46764670
/*!
46774671
* \brief Get the blending coefficient for the MUSCL scheme.

‎Common/src/CConfig.cpp‎

Lines changed: 5 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -1984,7 +1984,7 @@ void CConfig::SetConfig_Options() {
19841984

19851985
/*!\brief MUSCL_FLOW \n DESCRIPTION: Check if the MUSCL scheme should be used \ingroup Config*/
19861986
addBoolOption("MUSCL_FLOW", MUSCL_Flow, true);
1987-
/*!\brief MUSCL_KAPPA_FLOW \n DESCRIPTION: Blending coefficient for the MUSCL scheme \ingroup Config*/
1987+
/*!\brief MUSCL_KAPPA_FLOW \n DESCRIPTION: Blending coefficient for the U-MUSCL scheme \ingroup Config*/
19881988
addDoubleOption("MUSCL_KAPPA_FLOW", MUSCL_Kappa_Flow, 0.0);
19891989
/*!\brief SLOPE_LIMITER_FLOW
19901990
* DESCRIPTION: Slope limiter for the direct solution. \n OPTIONS: See \link Limiter_Map \endlink \n DEFAULT VENKATAKRISHNAN \ingroup Config*/
@@ -2007,7 +2007,7 @@ void CConfig::SetConfig_Options() {
20072007
addConvectOption("CONV_NUM_METHOD_ADJFLOW", Kind_ConvNumScheme_AdjFlow, Kind_Centered_AdjFlow, Kind_Upwind_AdjFlow);
20082008
/*!\brief MUSCL_ADJFLOW \n DESCRIPTION: Check if the MUSCL scheme should be used \ingroup Config*/
20092009
addBoolOption("MUSCL_ADJFLOW", MUSCL_AdjFlow, true);
2010-
/*!\brief MUSCL_KAPPA_ADJFLOW \n DESCRIPTION: Blending coefficient for the MUSCL scheme \ingroup Config*/
2010+
/*!\brief MUSCL_KAPPA_ADJFLOW \n DESCRIPTION: Blending coefficient for the U-MUSCL scheme \ingroup Config*/
20112011
addDoubleOption("MUSCL_KAPPA_ADJFLOW", MUSCL_Kappa_AdjFlow, 0.0);
20122012
/*!\brief SLOPE_LIMITER_ADJFLOW
20132013
* DESCRIPTION: Slope limiter for the adjoint solution. \n OPTIONS: See \link Limiter_Map \endlink \n DEFAULT VENKATAKRISHNAN \ingroup Config*/
@@ -2020,7 +2020,7 @@ void CConfig::SetConfig_Options() {
20202020

20212021
/*!\brief MUSCL_TURB \n DESCRIPTION: Check if the MUSCL scheme should be used \ingroup Config*/
20222022
addBoolOption("MUSCL_TURB", MUSCL_Turb, false);
2023-
/*!\brief MUSCL_KAPPA_TURB \n DESCRIPTION: Blending coefficient for the MUSCL scheme \ingroup Config*/
2023+
/*!\brief MUSCL_KAPPA_TURB \n DESCRIPTION: Blending coefficient for the U-MUSCL scheme \ingroup Config*/
20242024
addDoubleOption("MUSCL_KAPPA_TURB", MUSCL_Kappa_Turb, 0.0);
20252025
/*!\brief SLOPE_LIMITER_TURB
20262026
* \n DESCRIPTION: Slope limiter \n OPTIONS: See \link Limiter_Map \endlink \n DEFAULT VENKATAKRISHNAN \ingroup Config*/
@@ -2041,7 +2041,7 @@ void CConfig::SetConfig_Options() {
20412041

20422042
/*!\brief MUSCL_SPECIES \n DESCRIPTION: Check if the MUSCL scheme should be used \n DEFAULT false \ingroup Config*/
20432043
addBoolOption("MUSCL_SPECIES", MUSCL_Species, false);
2044-
/*!\brief MUSCL_KAPPA_SPECIES \n DESCRIPTION: Blending coefficient for the MUSCL scheme \ingroup Config*/
2044+
/*!\brief MUSCL_KAPPA_SPECIES \n DESCRIPTION: Blending coefficient for the U-MUSCL scheme \ingroup Config*/
20452045
addDoubleOption("MUSCL_KAPPA_SPECIES", MUSCL_Kappa_Species, 0.0);
20462046
/*!\brief SLOPE_LIMITER_SPECIES \n DESCRIPTION: Slope limiter \n OPTIONS: See \link Limiter_Map \endlink \n DEFAULT NONE \ingroup Config*/
20472047
addEnumOption("SLOPE_LIMITER_SPECIES", Kind_SlopeLimit_Species, Limiter_Map, LIMITER::NONE);
@@ -2050,7 +2050,7 @@ void CConfig::SetConfig_Options() {
20502050

20512051
/*!\brief MUSCL_HEAT \n DESCRIPTION: Check if the MUSCL scheme should be used \ingroup Config*/
20522052
addBoolOption("MUSCL_HEAT", MUSCL_Heat, false);
2053-
/*!\brief MUSCL_KAPPA_HEAT \n DESCRIPTION: Blending coefficient for the MUSCL scheme \ingroup Config*/
2053+
/*!\brief MUSCL_KAPPA_HEAT \n DESCRIPTION: Blending coefficient for the U-MUSCL scheme \ingroup Config*/
20542054
addDoubleOption("MUSCL_KAPPA_HEAT", MUSCL_Kappa_Heat, 0.0);
20552055
/*!\brief SLOPE_LIMITER_HEAT \n DESCRIPTION: Slope limiter \n OPTIONS: See \link Limiter_Map \endlink \n DEFAULT NONE \ingroup Config*/
20562056
addEnumOption("SLOPE_LIMITER_HEAT", Kind_SlopeLimit_Heat, Limiter_Map, LIMITER::NONE);

‎SU2_CFD/include/limiters/CLimiterDetails.hpp‎

Lines changed: 9 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -71,6 +71,15 @@ struct LimiterHelpers
7171
{
7272
FORCEINLINE static Type epsilon() {return std::numeric_limits<passivedouble>::epsilon();}
7373

74+
FORCEINLINE static Type umusclProjection(const Type& grad_proj, const Type& cent, const Type& kappa)
75+
{
76+
/*--- Upwind difference: dV_ij^upw = 2 grad(Vi) dot vector_ij - dV_ij^cent ---*/
77+
const Type upw = grad_proj - 0.5 * cent;
78+
79+
/*--- Blended difference: dV_ij^kappa = (1-kappa)dV_ij^upw + (1+kappa)dV_ij^cent ---*/
80+
return (1.0 - kappa) * upw + (1.0 + kappa) * cent;
81+
}
82+
7483
FORCEINLINE static Type venkatFunction(const Type& proj, const Type& delta, const Type& eps2)
7584
{
7685
Type y = delta*(delta+proj) + eps2;

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

Lines changed: 12 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -33,20 +33,20 @@
3333
#include "../../../variables/CNSVariable.hpp"
3434

3535
/*!
36-
* \brief Compute blended difference for MUSCL reconstruction.
36+
* \brief Blended difference for U-MUSCL reconstruction.
3737
* \param[in] grad_proj - Gradient projection at point i: dot(grad_i, vector_ij)
3838
* \param[in] cent - Centered difference: V_j - V_i
3939
* \param[in] kappa - Blending parameter
4040
* \return Blended difference for reconstruction
4141
*/
42-
FORCEINLINE Double blendedDifference(Double grad_proj,
43-
Double cent,
44-
Double kappa) {
42+
FORCEINLINE Double umusclProjection(Double grad_proj,
43+
Double cent,
44+
Double kappa) {
4545
/*--- Upwind difference: dV_ij^upw = 2 grad(Vi) dot vector_ij - dV_ij^cent ---*/
46-
const Double delta_upwind = grad_proj - 0.5 * cent;
46+
const Double upw = grad_proj - 0.5 * cent;
4747

4848
/*--- Blended difference: dV_ij^kappa = (1-kappa)dV_ij^upw + (1+kappa)dV_ij^cent ---*/
49-
return (1.0 - kappa) * delta_upwind + (1.0 + kappa) * cent;
49+
return (1.0 - kappa) * upw + (1.0 + kappa) * cent;
5050
}
5151

5252
/*!
@@ -75,8 +75,8 @@ FORCEINLINE void musclUnlimited(Int iPoint,
7575
const Double cent = V.j.all(iVar) - V.i.all(iVar);
7676

7777
/*--- Blended difference: dV_ij^kappa = (1-kappa)dV_ij^upw + (1+kappa)dV_ij^cent ---*/
78-
proj_i = blendedDifference(proj_i, cent, kappa);
79-
proj_j = blendedDifference(proj_j, cent, kappa);
78+
proj_i = umusclProjection(proj_i, cent, kappa);
79+
proj_j = umusclProjection(proj_j, cent, kappa);
8080
}
8181

8282
/*--- Apply reconstruction: V_ij = Vi + 0.25 * dV_ij^kappa ---*/
@@ -115,8 +115,8 @@ FORCEINLINE void musclPointLimited(Int iPoint,
115115
const Double cent = V.j.all(iVar) - V.i.all(iVar);
116116

117117
/*--- Blended difference: dV_ij^kappa = (1-kappa)dV_ij^upw + (1+kappa)dV_ij^cent ---*/
118-
proj_i = blendedDifference(proj_i, cent, kappa);
119-
proj_j = blendedDifference(proj_j, cent, kappa);
118+
proj_i = umusclProjection(proj_i, cent, kappa);
119+
proj_j = umusclProjection(proj_j, cent, kappa);
120120
}
121121

122122
/*--- Apply reconstruction: V_ij = Vi + 0.25 * lim * dV_ij^kappa ---*/
@@ -152,8 +152,8 @@ FORCEINLINE void musclEdgeLimited(Int iPoint,
152152
const Double cent = V.j.all(iVar) - V.i.all(iVar);
153153

154154
/*--- Blended difference: dV_ij^kappa = (1-kappa)dV_ij^upw + (1+kappa)dV_ij^cent ---*/
155-
proj_i = blendedDifference(proj_i, cent, kappa);
156-
proj_j = blendedDifference(proj_j, cent, kappa);
155+
proj_i = umusclProjection(proj_i, cent, kappa);
156+
proj_j = umusclProjection(proj_j, cent, kappa);
157157
}
158158

159159
/// TODO: Customize the limiter function.

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

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -77,7 +77,7 @@ class CRoeBase : public Base {
7777
finestGrid(iMesh == MESH_0),
7878
dynamicGrid(config.GetDynamic_Grid()),
7979
muscl(finestGrid && config.GetMUSCL_Flow()),
80-
umuscl(config.GetMUSCL_Kappa_Flow() != 0.0),
80+
umuscl(muscl && config.GetMUSCL_Kappa_Flow() != 0.0),
8181
umusclKappa(config.GetMUSCL_Kappa_Flow()),
8282
typeLimiter(config.GetKind_SlopeLimit_Flow()) {
8383
}

‎SU2_CFD/include/solvers/CScalarSolver.inl‎

Lines changed: 29 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -28,6 +28,7 @@
2828
#include "../../../Common/include/toolboxes/geometry_toolbox.hpp"
2929
#include "../../include/solvers/CScalarSolver.hpp"
3030
#include "../../include/variables/CFlowVariable.hpp"
31+
#include "../../include/limiters/CLimiterDetails.hpp"
3132

3233
template <class VariableType>
3334
CScalarSolver<VariableType>::CScalarSolver(CGeometry* geometry, CConfig* config, bool conservative)
@@ -142,6 +143,12 @@ void CScalarSolver<VariableType>::Upwind_Residual(CGeometry* geometry, CSolver**
142143
const bool limiterFlow =
143144
(config->GetKind_SlopeLimit_Flow() != LIMITER::NONE) && (config->GetKind_SlopeLimit_Flow() != LIMITER::VAN_ALBADA_EDGE);
144145

146+
/*--- U-MUSCL reconstruction ---*/
147+
const su2double kappa = config->GetMUSCL_Kappa();
148+
const su2double kappaFlow = config->GetMUSCL_Kappa_Flow();
149+
const bool umuscl = muscl && (kappa != 0.0);
150+
const bool umusclFlow = musclFlow && (kappaFlow != 0.0);
151+
145152
auto* flowNodes = su2staticcast_p<CFlowVariable*>(solver_container[FLOW_SOL]->GetNodes());
146153
const auto& edgeMassFluxes = *(solver_container[FLOW_SOL]->GetEdgeMassFluxes());
147154

@@ -220,18 +227,23 @@ void CScalarSolver<VariableType>::Upwind_Residual(CGeometry* geometry, CSolver**
220227
}
221228

222229
for (iVar = 0; iVar < solver_container[FLOW_SOL]->GetnPrimVarGrad(); iVar++) {
223-
su2double Project_Grad_i = 0.0;
224-
su2double Project_Grad_j = 0.0;
225-
for (iDim = 0; iDim < nDim; iDim++) {
226-
Project_Grad_i += Vector_ij[iDim] * Gradient_i[iVar][iDim];
227-
Project_Grad_j -= Vector_ij[iDim] * Gradient_j[iVar][iDim];
230+
su2double Project_Grad_i = GeometryToolbox::DotProduct(nDim, Gradient_i[iVar], Vector_ij);
231+
su2double Project_Grad_j = GeometryToolbox::DotProduct(nDim, Gradient_j[iVar], Vector_ij);
232+
233+
const su2double V_ij = (umusclFlow)? V_j[iVar] - V_i[iVar] : 0.0;
234+
235+
if (umusclFlow) {
236+
Project_Grad_i = LimiterHelpers<>::umusclProjection(Project_Grad_i, V_ij, kappaFlow);
237+
Project_Grad_j = LimiterHelpers<>::umusclProjection(Project_Grad_j, V_ij, kappaFlow);
228238
}
239+
229240
if (limiterFlow) {
230241
Project_Grad_i *= Limiter_i[iVar];
231242
Project_Grad_j *= Limiter_j[iVar];
232243
}
244+
233245
flowPrimVar_i[iVar] = V_i[iVar] + Project_Grad_i;
234-
flowPrimVar_j[iVar] = V_j[iVar] + Project_Grad_j;
246+
flowPrimVar_j[iVar] = V_j[iVar] - Project_Grad_j;
235247
}
236248

237249
numerics->SetPrimitive(flowPrimVar_i, flowPrimVar_j);
@@ -249,17 +261,23 @@ void CScalarSolver<VariableType>::Upwind_Residual(CGeometry* geometry, CSolver**
249261
}
250262

251263
for (iVar = 0; iVar < nVar; iVar++) {
252-
su2double Project_Grad_i = 0.0, Project_Grad_j = 0.0;
253-
for (iDim = 0; iDim < nDim; iDim++) {
254-
Project_Grad_i += Vector_ij[iDim] * Gradient_i[iVar][iDim];
255-
Project_Grad_j -= Vector_ij[iDim] * Gradient_j[iVar][iDim];
264+
su2double Project_Grad_i = GeometryToolbox::DotProduct(nDim, Gradient_i[iVar], Vector_ij);
265+
su2double Project_Grad_j = GeometryToolbox::DotProduct(nDim, Gradient_j[iVar], Vector_ij);
266+
267+
const su2double U_ij = (umuscl)? Scalar_j[iVar] - Scalar_i[iVar] : 0.0;
268+
269+
if (umuscl) {
270+
Project_Grad_i = LimiterHelpers<>::umusclProjection(Project_Grad_i, U_ij, kappa);
271+
Project_Grad_j = LimiterHelpers<>::umusclProjection(Project_Grad_j, U_ij, kappa);
256272
}
273+
257274
if (limiter) {
258275
Project_Grad_i *= Limiter_i[iVar];
259276
Project_Grad_j *= Limiter_j[iVar];
260277
}
278+
261279
solution_i[iVar] = Scalar_i[iVar] + Project_Grad_i;
262-
solution_j[iVar] = Scalar_j[iVar] + Project_Grad_j;
280+
solution_j[iVar] = Scalar_j[iVar] - Project_Grad_j;
263281
}
264282

265283
numerics->SetScalarVar(solution_i, solution_j);

‎SU2_CFD/src/solvers/CEulerSolver.cpp‎

Lines changed: 11 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -1804,6 +1804,9 @@ void CEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_contain
18041804
const bool limiter = (config->GetKind_SlopeLimit_Flow() != LIMITER::NONE);
18051805
const bool van_albada = (config->GetKind_SlopeLimit_Flow() == LIMITER::VAN_ALBADA_EDGE);
18061806

1807+
const su2double kappa = config->GetMUSCL_Kappa_Flow();
1808+
const bool umuscl = muscl && (kappa != 0.0);
1809+
18071810
/*--- Non-physical counter. ---*/
18081811
unsigned long counter_local = 0;
18091812
SU2_OMP_MASTER
@@ -1882,13 +1885,14 @@ void CEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_contain
18821885
auto Gradient_j = nodes->GetGradient_Reconstruction(jPoint);
18831886

18841887
for (iVar = 0; iVar < nPrimVarGrad; iVar++) {
1888+
su2double Project_Grad_i = GeometryToolbox::DotProduct(nDim, Gradient_i[iVar], Vector_ij);
1889+
su2double Project_Grad_j = GeometryToolbox::DotProduct(nDim, Gradient_j[iVar], Vector_ij);
18851890

1886-
su2double Project_Grad_i = 0.0;
1887-
su2double Project_Grad_j = 0.0;
1891+
const su2double V_ij = (umuscl || van_albada)? V_j[iVar] - V_i[iVar] : 0.0;
18881892

1889-
for (iDim = 0; iDim < nDim; iDim++) {
1890-
Project_Grad_i += Vector_ij[iDim]*Gradient_i[iVar][iDim];
1891-
Project_Grad_j -= Vector_ij[iDim]*Gradient_j[iVar][iDim];
1893+
if (umuscl) {
1894+
Project_Grad_i = LimiterHelpers<>::umusclProjection(Project_Grad_i, V_ij, kappa);
1895+
Project_Grad_j = LimiterHelpers<>::umusclProjection(Project_Grad_j, V_ij, kappa);
18921896
}
18931897

18941898
su2double lim_i = 1.0;
@@ -1897,15 +1901,15 @@ void CEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_contain
18971901
if (van_albada) {
18981902
su2double V_ij = V_j[iVar] - V_i[iVar];
18991903
lim_i = LimiterHelpers<>::vanAlbadaFunction(Project_Grad_i, V_ij, EPS);
1900-
lim_j = LimiterHelpers<>::vanAlbadaFunction(-Project_Grad_j, V_ij, EPS);
1904+
lim_j = LimiterHelpers<>::vanAlbadaFunction(Project_Grad_j, V_ij, EPS);
19011905
}
19021906
else if (limiter) {
19031907
lim_i = nodes->GetLimiter_Primitive(iPoint, iVar);
19041908
lim_j = nodes->GetLimiter_Primitive(jPoint, iVar);
19051909
}
19061910

19071911
Primitive_i[iVar] = V_i[iVar] + lim_i * Project_Grad_i;
1908-
Primitive_j[iVar] = V_j[iVar] + lim_j * Project_Grad_j;
1912+
Primitive_j[iVar] = V_j[iVar] - lim_j * Project_Grad_j;
19091913

19101914
}
19111915

‎SU2_CFD/src/solvers/CIncEulerSolver.cpp‎

Lines changed: 11 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -1233,6 +1233,9 @@ void CIncEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_cont
12331233
const bool bounded_scalar = config->GetBounded_Scalar();
12341234
const su2double nkRelax = config->GetNewtonKrylovRelaxation();
12351235

1236+
const su2double kappa = config->GetMUSCL_Kappa_Flow();
1237+
const bool umuscl = muscl && (kappa != 0.0);
1238+
12361239
/*--- For hybrid parallel AD, pause preaccumulation if there is shared reading of
12371240
* variables, otherwise switch to the faster adjoint evaluation mode. ---*/
12381241
bool pausePreacc = false;
@@ -1279,13 +1282,14 @@ void CIncEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_cont
12791282
auto Gradient_j = nodes->GetGradient_Reconstruction(jPoint);
12801283

12811284
for (iVar = 0; iVar < nPrimVarGrad; iVar++) {
1285+
su2double Project_Grad_i = GeometryToolbox::DotProduct(nDim, Gradient_i[iVar], Vector_ij);
1286+
su2double Project_Grad_j = GeometryToolbox::DotProduct(nDim, Gradient_j[iVar], Vector_ij);
12821287

1283-
su2double Project_Grad_i = 0.0;
1284-
su2double Project_Grad_j = 0.0;
1288+
const su2double V_ij = (umuscl || van_albada)? V_j[iVar] - V_i[iVar] : 0.0;
12851289

1286-
for (iDim = 0; iDim < nDim; iDim++) {
1287-
Project_Grad_i += Vector_ij[iDim]*Gradient_i[iVar][iDim];
1288-
Project_Grad_j -= Vector_ij[iDim]*Gradient_j[iVar][iDim];
1290+
if (umuscl) {
1291+
Project_Grad_i = LimiterHelpers<>::umusclProjection(Project_Grad_i, V_ij, kappa);
1292+
Project_Grad_j = LimiterHelpers<>::umusclProjection(Project_Grad_j, V_ij, kappa);
12891293
}
12901294

12911295
su2double lim_i = 1.0;
@@ -1294,15 +1298,15 @@ void CIncEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_cont
12941298
if (van_albada) {
12951299
su2double V_ij = V_j[iVar] - V_i[iVar];
12961300
lim_i = LimiterHelpers<>::vanAlbadaFunction(Project_Grad_i, V_ij, EPS);
1297-
lim_j = LimiterHelpers<>::vanAlbadaFunction(-Project_Grad_j, V_ij, EPS);
1301+
lim_j = LimiterHelpers<>::vanAlbadaFunction(Project_Grad_j, V_ij, EPS);
12981302
}
12991303
else if (limiter) {
13001304
lim_i = nodes->GetLimiter_Primitive(iPoint, iVar);
13011305
lim_j = nodes->GetLimiter_Primitive(jPoint, iVar);
13021306
}
13031307

13041308
Primitive_i[iVar] = V_i[iVar] + lim_i * Project_Grad_i;
1305-
Primitive_j[iVar] = V_j[iVar] + lim_j * Project_Grad_j;
1309+
Primitive_j[iVar] = V_j[iVar] - lim_j * Project_Grad_j;
13061310
}
13071311

13081312
for (iVar = nPrimVarGrad; iVar < nPrimVar; iVar++) {

0 commit comments

Comments
 (0)