Skip to content

Commit 00d7a2f

Browse files
pcarruscagpcarruscag
committed
use CFlowVariable in scalar and turbulence solvers
1 parent 2736dc1 commit 00d7a2f

10 files changed

Lines changed: 155 additions & 141 deletions

File tree

‎Common/include/code_config.hpp‎

Lines changed: 11 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -67,6 +67,17 @@ template<class T, class F> struct su2conditional<false, T, F> { using type = F;
6767
template<bool B, class T, class F>
6868
using su2conditional_t = typename su2conditional<B,T,F>::type;
6969

70+
/*! \brief Static cast "In" to "Out", in debug builds a dynamic cast is used. */
71+
template<class Out, class In>
72+
FORCEINLINE Out su2staticcast_p(In ptr) {
73+
static_assert(std::is_pointer<In>::value, "This expects a pointer");
74+
#ifndef NDEBUG
75+
return static_cast<Out>(ptr);
76+
#else
77+
return dynamic_cast<Out>(ptr);
78+
#endif
79+
}
80+
7081
/*--- Detect compilation with OpenMP. ---*/
7182
#if defined(_OPENMP)
7283
#define HAVE_OMP

‎SU2_CFD/include/limiters/CLimiterDetails.hpp‎

Lines changed: 25 additions & 19 deletions
Original file line numberDiff line numberDiff line change
@@ -64,22 +64,28 @@ struct CLimiterDetails
6464
/*!
6565
* \brief Common small functions used by limiters.
6666
*/
67-
namespace LimiterHelpers
67+
template<class Type = su2double>
68+
struct LimiterHelpers
6869
{
69-
inline passivedouble epsilon() {return std::numeric_limits<passivedouble>::epsilon();}
70+
FORCEINLINE static Type epsilon() {return std::numeric_limits<passivedouble>::epsilon();}
7071

71-
inline su2double venkatFunction(su2double proj, su2double delta, su2double eps2)
72+
FORCEINLINE static Type venkatFunction(const Type& proj, const Type& delta, const Type& eps2)
7273
{
73-
su2double y = delta*(delta+proj) + eps2;
74+
Type y = delta*(delta+proj) + eps2;
7475
return (y + delta*proj) / (y + 2*proj*proj);
7576
}
7677

77-
inline su2double raisedSine(su2double dist)
78+
FORCEINLINE static Type vanAlbadaFunction(const Type& proj, const Type& delta, const Type& eps)
7879
{
79-
su2double factor = 0.5*(1.0+dist+sin(PI_NUMBER*dist)/PI_NUMBER);
80+
return delta * (2*proj + delta) / (4*pow(proj, 2) + pow(delta, 2) + eps);
81+
}
82+
83+
FORCEINLINE static Type raisedSine(const Type& dist)
84+
{
85+
Type factor = 0.5*(1.0+dist+sin(PI_NUMBER*dist)/PI_NUMBER);
8086
return max(0.0, min(factor, 1.0));
8187
}
82-
}
88+
};
8389

8490

8591
/*!
@@ -94,7 +100,7 @@ struct CLimiterDetails<BARTH_JESPERSEN>
94100
* \brief Set a small epsilon to avoid divisions by 0.
95101
*/
96102
template<class... Ts>
97-
inline void preprocess(Ts&...) {eps2 = LimiterHelpers::epsilon();}
103+
inline void preprocess(Ts&...) {eps2 = LimiterHelpers<>::epsilon();}
98104

99105
/*!
100106
* \brief No geometric modification for this kind of limiter.
@@ -107,7 +113,7 @@ struct CLimiterDetails<BARTH_JESPERSEN>
107113
*/
108114
inline su2double limiterFunction(size_t, su2double proj, su2double delta) const
109115
{
110-
return LimiterHelpers::venkatFunction(proj, delta, eps2);
116+
return LimiterHelpers<>::venkatFunction(proj, delta, eps2);
111117
}
112118
};
113119

@@ -130,7 +136,7 @@ struct CLimiterDetails<VENKATAKRISHNAN>
130136
su2double L = config.GetRefElemLength();
131137
su2double K = config.GetVenkat_LimiterCoeff();
132138
su2double eps1 = fabs(L*K);
133-
eps2 = max(eps1*eps1*eps1, LimiterHelpers::epsilon());
139+
eps2 = max(eps1*eps1*eps1, LimiterHelpers<>::epsilon());
134140
}
135141

136142
/*!
@@ -144,7 +150,7 @@ struct CLimiterDetails<VENKATAKRISHNAN>
144150
*/
145151
inline su2double limiterFunction(size_t, su2double proj, su2double delta) const
146152
{
147-
return LimiterHelpers::venkatFunction(proj, delta, eps2);
153+
return LimiterHelpers<>::venkatFunction(proj, delta, eps2);
148154
}
149155
};
150156

@@ -229,7 +235,7 @@ struct CLimiterDetails<VENKATAKRISHNAN_WANG>
229235
for(size_t iVar = varBegin; iVar < varEnd; ++iVar)
230236
{
231237
su2double range = sharedMax(iVar) - sharedMin(iVar);
232-
eps2(iVar) = max(pow(K*range, 2), LimiterHelpers::epsilon());
238+
eps2(iVar) = max(pow(K*range, 2), LimiterHelpers<>::epsilon());
233239
}
234240
}
235241

@@ -245,7 +251,7 @@ struct CLimiterDetails<VENKATAKRISHNAN_WANG>
245251
inline su2double limiterFunction(size_t iVar, su2double proj, su2double delta) const
246252
{
247253
AD::SetPreaccIn(eps2(iVar));
248-
return LimiterHelpers::venkatFunction(proj, delta, eps2(iVar));
254+
return LimiterHelpers<>::venkatFunction(proj, delta, eps2(iVar));
249255
}
250256
};
251257

@@ -268,7 +274,7 @@ struct CLimiterDetails<SHARP_EDGES>
268274
su2double L = config.GetRefElemLength();
269275
su2double K = config.GetVenkat_LimiterCoeff();
270276
eps1 = fabs(L*K);
271-
eps2 = max(eps1*eps1*eps1, LimiterHelpers::epsilon());
277+
eps2 = max(eps1*eps1*eps1, LimiterHelpers<>::epsilon());
272278
}
273279

274280
/*!
@@ -278,15 +284,15 @@ struct CLimiterDetails<SHARP_EDGES>
278284
{
279285
AD::SetPreaccIn(geometry.nodes->GetSharpEdge_Distance(iPoint));
280286
su2double dist = geometry.nodes->GetSharpEdge_Distance(iPoint)/(sharpCoeff*eps1)-1.0;
281-
return LimiterHelpers::raisedSine(dist);
287+
return LimiterHelpers<>::raisedSine(dist);
282288
}
283289

284290
/*!
285291
* \brief Smooth function that disables limiting in smooth regions.
286292
*/
287293
inline su2double limiterFunction(size_t, su2double proj, su2double delta) const
288294
{
289-
return LimiterHelpers::venkatFunction(proj, delta, eps2);
295+
return LimiterHelpers<>::venkatFunction(proj, delta, eps2);
290296
}
291297
};
292298

@@ -309,7 +315,7 @@ struct CLimiterDetails<WALL_DISTANCE>
309315
su2double L = config.GetRefElemLength();
310316
su2double K = config.GetVenkat_LimiterCoeff();
311317
eps1 = fabs(L*K);
312-
eps2 = max(eps1*eps1*eps1, LimiterHelpers::epsilon());
318+
eps2 = max(eps1*eps1*eps1, LimiterHelpers<>::epsilon());
313319
}
314320

315321
/*!
@@ -319,14 +325,14 @@ struct CLimiterDetails<WALL_DISTANCE>
319325
{
320326
AD::SetPreaccIn(geometry.nodes->GetWall_Distance(iPoint));
321327
su2double dist = geometry.nodes->GetWall_Distance(iPoint)/(sharpCoeff*eps1)-1.0;
322-
return LimiterHelpers::raisedSine(dist);
328+
return LimiterHelpers<>::raisedSine(dist);
323329
}
324330

325331
/*!
326332
* \brief Smooth function that disables limiting in smooth regions.
327333
*/
328334
inline su2double limiterFunction(size_t, su2double proj, su2double delta) const
329335
{
330-
return LimiterHelpers::venkatFunction(proj, delta, eps2);
336+
return LimiterHelpers<>::venkatFunction(proj, delta, eps2);
331337
}
332338
};

‎SU2_CFD/include/solvers/CScalarSolver.inl‎

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -27,6 +27,7 @@
2727
#include "../../../Common/include/parallelization/omp_structure.hpp"
2828
#include "../../../Common/include/toolboxes/geometry_toolbox.hpp"
2929
#include "../../include/solvers/CScalarSolver.hpp"
30+
#include "../../include/variables/CFlowVariable.hpp"
3031

3132
template <class VariableType>
3233
CScalarSolver<VariableType>::CScalarSolver(bool conservative) : CSolver(), Conservative(conservative) {}
@@ -92,10 +93,10 @@ void CScalarSolver<VariableType>::Upwind_Residual(CGeometry* geometry, CSolver**
9293
const bool limiterFlow =
9394
(config->GetKind_SlopeLimit_Flow() != NO_LIMITER) && (config->GetKind_SlopeLimit_Flow() != VAN_ALBADA_EDGE);
9495

95-
CVariable* flowNodes = solver_container[FLOW_SOL]->GetNodes();
96+
auto* flowNodes = su2staticcast_p<CFlowVariable*>(solver_container[FLOW_SOL]->GetNodes());
9697

9798
/*--- Pick one numerics object per thread. ---*/
98-
CNumerics* numerics = numerics_container[CONV_TERM + omp_get_thread_num() * MAX_TERMS];
99+
auto* numerics = numerics_container[CONV_TERM + omp_get_thread_num() * MAX_TERMS];
99100

100101
/*--- Static arrays of MUSCL-reconstructed flow primitives and turbulence variables (thread safety). ---*/
101102
su2double solution_i[MAXNVAR] = {0.0}, flowPrimVar_i[MAXNVARFLOW] = {0.0};

‎SU2_CFD/include/variables/CScalarVariable.hpp‎

Lines changed: 1 addition & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -50,12 +50,7 @@ class CScalarVariable : public CVariable {
5050
* \param[in] nvar - Number of variables of the problem.
5151
* \param[in] config - Definition of the particular problem.
5252
*/
53-
CScalarVariable(unsigned long npoint, unsigned long ndim, unsigned long nvar, CConfig* config);
54-
55-
/*!
56-
* \brief Destructor of the class.
57-
*/
58-
~CScalarVariable() override = default;
53+
CScalarVariable(unsigned long npoint, unsigned long ndim, unsigned long nvar, const CConfig* config);
5954

6055
/*!
6156
* \brief Get the array of the reconstruction variables gradient at a node.

‎SU2_CFD/src/solvers/CEulerSolver.cpp‎

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -33,6 +33,7 @@
3333
#include "../../include/fluid/CVanDerWaalsGas.hpp"
3434
#include "../../include/fluid/CPengRobinson.hpp"
3535
#include "../../include/numerics_simd/CNumericsSIMD.hpp"
36+
#include "../../include/limiters/CLimiterDetails.hpp"
3637

3738

3839
CEulerSolver::CEulerSolver(CGeometry *geometry, CConfig *config,
@@ -1782,8 +1783,8 @@ void CEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_contain
17821783

17831784
if (van_albada) {
17841785
su2double V_ij = V_j[iVar] - V_i[iVar];
1785-
lim_i = V_ij*( 2.0*Project_Grad_i + V_ij) / (4*pow(Project_Grad_i, 2) + pow(V_ij, 2) + EPS);
1786-
lim_j = V_ij*(-2.0*Project_Grad_j + V_ij) / (4*pow(Project_Grad_j, 2) + pow(V_ij, 2) + EPS);
1786+
lim_i = LimiterHelpers<>::vanAlbadaFunction(Project_Grad_i, V_ij, EPS);
1787+
lim_j = LimiterHelpers<>::vanAlbadaFunction(-Project_Grad_j, V_ij, EPS);
17871788
}
17881789
else if (limiter) {
17891790
lim_i = nodes->GetLimiter_Primitive(iPoint, iVar);

‎SU2_CFD/src/solvers/CIncEulerSolver.cpp‎

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -31,6 +31,7 @@
3131
#include "../../include/fluid/CIncIdealGas.hpp"
3232
#include "../../include/fluid/CIncIdealGasPolynomial.hpp"
3333
#include "../../include/variables/CIncNSVariable.hpp"
34+
#include "../../include/limiters/CLimiterDetails.hpp"
3435
#include "../../../Common/include/toolboxes/geometry_toolbox.hpp"
3536

3637

@@ -1184,8 +1185,8 @@ void CIncEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_cont
11841185

11851186
if (van_albada) {
11861187
su2double V_ij = V_j[iVar] - V_i[iVar];
1187-
lim_i = V_ij*( 2.0*Project_Grad_i + V_ij) / (4*pow(Project_Grad_i, 2) + pow(V_ij, 2) + EPS);
1188-
lim_j = V_ij*(-2.0*Project_Grad_j + V_ij) / (4*pow(Project_Grad_j, 2) + pow(V_ij, 2) + EPS);
1188+
lim_i = LimiterHelpers<>::vanAlbadaFunction(Project_Grad_i, V_ij, EPS);
1189+
lim_j = LimiterHelpers<>::vanAlbadaFunction(-Project_Grad_j, V_ij, EPS);
11891190
}
11901191
else if (limiter) {
11911192
lim_i = nodes->GetLimiter_Primitive(iPoint, iVar);

‎SU2_CFD/src/solvers/CNEMOEulerSolver.cpp‎

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -31,6 +31,7 @@
3131
#include "../../../Common/include/toolboxes/printing_toolbox.hpp"
3232
#include "../../include/fluid/CMutationTCLib.hpp"
3333
#include "../../include/fluid/CSU2TCLib.hpp"
34+
#include "../../include/limiters/CLimiterDetails.hpp"
3435

3536
CNEMOEulerSolver::CNEMOEulerSolver(CGeometry *geometry, CConfig *config,
3637
unsigned short iMesh, const bool navier_stokes) :
@@ -562,8 +563,8 @@ void CNEMOEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_con
562563

563564
if (van_albada) {
564565
su2double V_ij = V_j[iVar] - V_i[iVar];
565-
lim_i = V_ij*( 2.0*Project_Grad_i + V_ij) / (4*pow(Project_Grad_i, 2) + pow(V_ij, 2) + EPS);
566-
lim_j = V_ij*(-2.0*Project_Grad_j + V_ij) / (4*pow(Project_Grad_j, 2) + pow(V_ij, 2) + EPS);
566+
lim_i = LimiterHelpers<>::vanAlbadaFunction(Project_Grad_i, V_ij, EPS);
567+
lim_j = LimiterHelpers<>::vanAlbadaFunction(-Project_Grad_j, V_ij, EPS);
567568
}
568569
else if (limiter) {
569570
lim_i = nodes->GetLimiter_Primitive(iPoint, iVar);

0 commit comments

Comments
 (0)