Skip to content

Commit 01103cb

Browse files
WallyMaierWallyMaier
authored
Merge branch 'develop' into feature_NEMO_implicit_viscous
2 parents b70962c + 641d254 commit 01103cb

37 files changed

Lines changed: 818 additions & 800 deletions

‎SU2_CFD/include/drivers/CDriver.hpp‎

Lines changed: 7 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -213,6 +213,13 @@ class CDriver {
213213
*/
214214
void Numerics_Preprocessing(CConfig *config, CGeometry **geometry, CSolver ***solver, CNumerics ****&numerics) const;
215215

216+
/*!
217+
* \brief Helper to instantiate turbulence numerics specialized for different flow solvers.
218+
*/
219+
template <class FlowIndices>
220+
void InstantiateTurbulentNumerics(unsigned short nVar_Turb, int offset, const CConfig *config,
221+
const CSolver* turb_solver, CNumerics ****&numerics) const;
222+
216223
/*!
217224
* \brief Definition and allocation of all solver classes.
218225
* \param[in] numerics_container - Description of the numerical method (the way in which the equations are solved).

‎SU2_CFD/include/numerics/scalar/scalar_convection.hpp‎

Lines changed: 65 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -44,18 +44,24 @@
4444
* \ingroup ConvDiscr
4545
* \author C. Pederson, A. Bueno., and A. Campos.
4646
*/
47+
template <class FlowIndices>
4748
class CUpwScalar : public CNumerics {
4849
protected:
49-
su2double a0 = 0.0; /*!< \brief The maximum of the face-normal velocity and 0 */
50-
su2double a1 = 0.0; /*!< \brief The minimum of the face-normal velocity and 0 */
51-
su2double* Flux = nullptr; /*!< \brief Final result, diffusive flux/residual. */
52-
su2double** Jacobian_i = nullptr; /*!< \brief Flux Jacobian w.r.t. node i. */
53-
su2double** Jacobian_j = nullptr; /*!< \brief Flux Jacobian w.r.t. node j. */
50+
enum : unsigned short {MAXNVAR = 8};
51+
52+
const FlowIndices idx; /*!< \brief Object to manage the access to the flow primitives. */
53+
su2double a0 = 0.0; /*!< \brief The maximum of the face-normal velocity and 0. */
54+
su2double a1 = 0.0; /*!< \brief The minimum of the face-normal velocity and 0. */
55+
su2double Flux[MAXNVAR]; /*!< \brief Final result, diffusive flux/residual. */
56+
su2double* Jacobian_i[MAXNVAR]; /*!< \brief Flux Jacobian w.r.t. node i. */
57+
su2double* Jacobian_j[MAXNVAR]; /*!< \brief Flux Jacobian w.r.t. node j. */
58+
su2double JacobianBuffer[2*MAXNVAR*MAXNVAR]; /*!< \brief Static storage for the two Jacobians. */
5459

5560
const bool implicit = false, incompressible = false, dynamic_grid = false;
5661

5762
/*!
58-
* \brief A pure virtual function; Adds any extra variables to AD
63+
* \brief A pure virtual function. Derived classes must use it to register the additional
64+
* variables they use as preaccumulation inputs, e.g. the density for SST.
5965
*/
6066
virtual void ExtraADPreaccIn() = 0;
6167

@@ -69,21 +75,65 @@ class CUpwScalar : public CNumerics {
6975
public:
7076
/*!
7177
* \brief Constructor of the class.
72-
* \param[in] val_nDim - Number of dimensions of the problem.
73-
* \param[in] val_nVar - Number of variables of the problem.
78+
* \param[in] ndim - Number of dimensions of the problem.
79+
* \param[in] nvar - Number of variables of the problem.
7480
* \param[in] config - Definition of the particular problem.
7581
*/
76-
CUpwScalar(unsigned short val_nDim, unsigned short val_nVar, const CConfig* config);
77-
78-
/*!
79-
* \brief Destructor of the class.
80-
*/
81-
~CUpwScalar(void) override;
82+
CUpwScalar(unsigned short ndim, unsigned short nvar, const CConfig* config)
83+
: CNumerics(ndim, nvar, config),
84+
idx(ndim, config->GetnSpecies()),
85+
implicit(config->GetKind_TimeIntScheme_Turb() == EULER_IMPLICIT),
86+
incompressible(config->GetKind_Regime() == ENUM_REGIME::INCOMPRESSIBLE),
87+
dynamic_grid(config->GetDynamic_Grid()) {
88+
if (nVar > MAXNVAR) {
89+
SU2_MPI::Error("Static arrays are too small.", CURRENT_FUNCTION);
90+
}
91+
for (unsigned short iVar = 0; iVar < nVar; iVar++) {
92+
Jacobian_i[iVar] = &JacobianBuffer[iVar * nVar];
93+
Jacobian_j[iVar] = &JacobianBuffer[iVar * nVar + MAXNVAR * MAXNVAR];
94+
}
95+
}
8296

8397
/*!
8498
* \brief Compute the scalar upwind flux between two nodes i and j.
8599
* \param[in] config - Definition of the particular problem.
86100
* \return A lightweight const-view (read-only) of the residual/flux and Jacobians.
87101
*/
88-
ResidualType<> ComputeResidual(const CConfig* config) override;
102+
CNumerics::ResidualType<> ComputeResidual(const CConfig* config) final {
103+
AD::StartPreacc();
104+
AD::SetPreaccIn(Normal, nDim);
105+
AD::SetPreaccIn(ScalarVar_i, nVar);
106+
AD::SetPreaccIn(ScalarVar_j, nVar);
107+
if (dynamic_grid) {
108+
AD::SetPreaccIn(GridVel_i, nDim);
109+
AD::SetPreaccIn(GridVel_j, nDim);
110+
}
111+
AD::SetPreaccIn(&V_i[idx.Velocity()], nDim);
112+
AD::SetPreaccIn(&V_j[idx.Velocity()], nDim);
113+
114+
ExtraADPreaccIn();
115+
116+
su2double q_ij = 0.0;
117+
if (dynamic_grid) {
118+
for (unsigned short iDim = 0; iDim < nDim; iDim++) {
119+
su2double Velocity_i = V_i[iDim + idx.Velocity()] - GridVel_i[iDim];
120+
su2double Velocity_j = V_j[iDim + idx.Velocity()] - GridVel_j[iDim];
121+
q_ij += 0.5 * (Velocity_i + Velocity_j) * Normal[iDim];
122+
}
123+
} else {
124+
for (unsigned short iDim = 0; iDim < nDim; iDim++) {
125+
q_ij += 0.5 * (V_i[iDim + idx.Velocity()] + V_j[iDim + idx.Velocity()]) * Normal[iDim];
126+
}
127+
}
128+
129+
a0 = 0.5 * (q_ij + fabs(q_ij));
130+
a1 = 0.5 * (q_ij - fabs(q_ij));
131+
132+
FinishResidualCalc(config);
133+
134+
AD::SetPreaccOut(Flux, nVar);
135+
AD::EndPreacc();
136+
137+
return ResidualType<>(Flux, Jacobian_i, Jacobian_j);
138+
}
89139
};

‎SU2_CFD/include/numerics/scalar/scalar_diffusion.hpp‎

Lines changed: 59 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -44,14 +44,18 @@
4444
* \ingroup ViscDiscr
4545
* \author C. Pederson, A. Bueno, and F. Palacios
4646
*/
47+
template <class FlowIndices>
4748
class CAvgGrad_Scalar : public CNumerics {
4849
protected:
49-
su2double* Proj_Mean_GradScalarVar_Normal = nullptr; /*!< \brief Mean_gradScalarVar DOT normal. */
50-
su2double* Proj_Mean_GradScalarVar = nullptr; /*!< \brief Mean_gradScalarVar DOT normal, corrected if required. */
51-
su2double proj_vector_ij = 0.0; /*!< \brief (Edge_Vector DOT normal)/|Edge_Vector|^2 */
52-
su2double* Flux = nullptr; /*!< \brief Final result, diffusive flux/residual. */
53-
su2double** Jacobian_i = nullptr; /*!< \brief Flux Jacobian w.r.t. node i. */
54-
su2double** Jacobian_j = nullptr; /*!< \brief Flux Jacobian w.r.t. node j. */
50+
enum : unsigned short {MAXNVAR = 8};
51+
52+
const FlowIndices idx; /*!< \brief Object to manage the access to the flow primitives. */
53+
su2double Proj_Mean_GradScalarVar[MAXNVAR]; /*!< \brief Mean_gradScalarVar DOT normal, corrected if required. */
54+
su2double proj_vector_ij = 0.0; /*!< \brief (Edge_Vector DOT normal)/|Edge_Vector|^2 */
55+
su2double Flux[MAXNVAR]; /*!< \brief Final result, diffusive flux/residual. */
56+
su2double* Jacobian_i[MAXNVAR]; /*!< \brief Flux Jacobian w.r.t. node i. */
57+
su2double* Jacobian_j[MAXNVAR]; /*!< \brief Flux Jacobian w.r.t. node j. */
58+
su2double JacobianBuffer[2*MAXNVAR*MAXNVAR];/*!< \brief Static storage for the two Jacobians. */
5559

5660
const bool correct_gradient = false, implicit = false, incompressible = false;
5761

@@ -75,17 +79,59 @@ class CAvgGrad_Scalar : public CNumerics {
7579
* \param[in] correct_gradient - Whether to correct gradient for skewness.
7680
* \param[in] config - Definition of the particular problem.
7781
*/
78-
CAvgGrad_Scalar(unsigned short val_nDim, unsigned short val_nVar, bool correct_gradient, const CConfig* config);
79-
80-
/*!
81-
* \brief Destructor of the class.
82-
*/
83-
~CAvgGrad_Scalar(void) override;
82+
CAvgGrad_Scalar(unsigned short val_nDim, unsigned short val_nVar, bool correct_grad,
83+
const CConfig* config)
84+
: CNumerics(val_nDim, val_nVar, config),
85+
idx(val_nDim, config->GetnSpecies()),
86+
correct_gradient(correct_grad),
87+
implicit(config->GetKind_TimeIntScheme_Turb() == EULER_IMPLICIT),
88+
incompressible(config->GetKind_Regime() == ENUM_REGIME::INCOMPRESSIBLE) {
89+
if (nVar > MAXNVAR) {
90+
SU2_MPI::Error("Static arrays are too small.", CURRENT_FUNCTION);
91+
}
92+
for (unsigned short iVar = 0; iVar < nVar; iVar++) {
93+
Jacobian_i[iVar] = &JacobianBuffer[iVar * nVar];
94+
Jacobian_j[iVar] = &JacobianBuffer[iVar * nVar + MAXNVAR * MAXNVAR];
95+
}
96+
}
8497

8598
/*!
8699
* \brief Compute the viscous residual using an average of gradients without correction.
87100
* \param[in] config - Definition of the particular problem.
88101
* \return A lightweight const-view (read-only) of the residual/flux and Jacobians.
89102
*/
90-
ResidualType<> ComputeResidual(const CConfig* config) override;
103+
ResidualType<> ComputeResidual(const CConfig* config) final {
104+
AD::StartPreacc();
105+
AD::SetPreaccIn(Coord_i, nDim);
106+
AD::SetPreaccIn(Coord_j, nDim);
107+
AD::SetPreaccIn(Normal, nDim);
108+
AD::SetPreaccIn(ScalarVar_Grad_i, nVar, nDim);
109+
AD::SetPreaccIn(ScalarVar_Grad_j, nVar, nDim);
110+
if (correct_gradient) {
111+
AD::SetPreaccIn(ScalarVar_i, nVar);
112+
AD::SetPreaccIn(ScalarVar_j, nVar);
113+
}
114+
AD::SetPreaccIn(V_i[idx.Density()], V_i[idx.LaminarViscosity()], V_i[idx.EddyViscosity()]);
115+
AD::SetPreaccIn(V_j[idx.Density()], V_j[idx.LaminarViscosity()], V_j[idx.EddyViscosity()]);
116+
117+
Density_i = V_i[idx.Density()];
118+
Density_j = V_j[idx.Density()];
119+
Laminar_Viscosity_i = V_i[idx.LaminarViscosity()];
120+
Laminar_Viscosity_j = V_j[idx.LaminarViscosity()];
121+
Eddy_Viscosity_i = V_i[idx.EddyViscosity()];
122+
Eddy_Viscosity_j = V_j[idx.EddyViscosity()];
123+
124+
ExtraADPreaccIn();
125+
126+
su2double ProjGradScalarVarNoCorr[MAXNVAR];
127+
proj_vector_ij = ComputeProjectedGradient(nDim, nVar, Normal, Coord_i, Coord_j, ScalarVar_Grad_i, ScalarVar_Grad_j,
128+
correct_gradient, ScalarVar_i, ScalarVar_j, ProjGradScalarVarNoCorr,
129+
Proj_Mean_GradScalarVar);
130+
FinishResidualCalc(config);
131+
132+
AD::SetPreaccOut(Flux, nVar);
133+
AD::EndPreacc();
134+
135+
return ResidualType<>(Flux, Jacobian_i, Jacobian_j);
136+
}
91137
};

‎SU2_CFD/include/numerics/turbulent/turb_convection.hpp‎

Lines changed: 58 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -36,18 +36,36 @@
3636
* \ingroup ConvDiscr
3737
* \author A. Bueno.
3838
*/
39-
class CUpwSca_TurbSA final : public CUpwScalar {
39+
template <class FlowIndices>
40+
class CUpwSca_TurbSA final : public CUpwScalar<FlowIndices> {
4041
private:
42+
using Base = CUpwScalar<FlowIndices>;
43+
using Base::a0;
44+
using Base::a1;
45+
using Base::Flux;
46+
using Base::Jacobian_i;
47+
using Base::Jacobian_j;
48+
using Base::ScalarVar_i;
49+
using Base::ScalarVar_j;
50+
using Base::implicit;
51+
4152
/*!
42-
* \brief Adds any extra variables to AD
53+
* \brief Adds any extra variables to AD.
4354
*/
44-
void ExtraADPreaccIn() override;
55+
void ExtraADPreaccIn() override {}
4556

4657
/*!
4758
* \brief SA specific steps in the ComputeResidual method
4859
* \param[in] config - Definition of the particular problem.
4960
*/
50-
void FinishResidualCalc(const CConfig* config) override;
61+
void FinishResidualCalc(const CConfig* config) override {
62+
Flux[0] = a0*ScalarVar_i[0] + a1*ScalarVar_j[0];
63+
64+
if (implicit) {
65+
Jacobian_i[0][0] = a0;
66+
Jacobian_j[0][0] = a1;
67+
}
68+
}
5169

5270
public:
5371
/*!
@@ -56,8 +74,8 @@ class CUpwSca_TurbSA final : public CUpwScalar {
5674
* \param[in] val_nVar - Number of variables of the problem.
5775
* \param[in] config - Definition of the particular problem.
5876
*/
59-
CUpwSca_TurbSA(unsigned short val_nDim, unsigned short val_nVar, const CConfig* config);
60-
77+
CUpwSca_TurbSA(unsigned short val_nDim, unsigned short val_nVar, const CConfig* config)
78+
: CUpwScalar<FlowIndices>(val_nDim, val_nVar, config) {}
6179
};
6280

6381
/*!
@@ -66,18 +84,47 @@ class CUpwSca_TurbSA final : public CUpwScalar {
6684
* \ingroup ConvDiscr
6785
* \author A. Campos.
6886
*/
69-
class CUpwSca_TurbSST final : public CUpwScalar {
87+
template <class FlowIndices>
88+
class CUpwSca_TurbSST final : public CUpwScalar<FlowIndices> {
7089
private:
90+
using Base = CUpwScalar<FlowIndices>;
91+
using Base::nDim;
92+
using Base::V_i;
93+
using Base::V_j;
94+
using Base::a0;
95+
using Base::a1;
96+
using Base::Flux;
97+
using Base::Jacobian_i;
98+
using Base::Jacobian_j;
99+
using Base::ScalarVar_i;
100+
using Base::ScalarVar_j;
101+
using Base::implicit;
102+
using Base::idx;
103+
71104
/*!
72105
* \brief Adds any extra variables to AD
73106
*/
74-
void ExtraADPreaccIn() override;
107+
void ExtraADPreaccIn() override {
108+
AD::SetPreaccIn(V_i[idx.Density()]);
109+
AD::SetPreaccIn(V_j[idx.Density()]);
110+
}
75111

76112
/*!
77113
* \brief SST specific steps in the ComputeResidual method
78114
* \param[in] config - Definition of the particular problem.
79115
*/
80-
void FinishResidualCalc(const CConfig* config) override;
116+
void FinishResidualCalc(const CConfig* config) override {
117+
Flux[0] = a0*V_i[idx.Density()]*ScalarVar_i[0] + a1*V_j[idx.Density()]*ScalarVar_j[0];
118+
Flux[1] = a0*V_i[idx.Density()]*ScalarVar_i[1] + a1*V_j[idx.Density()]*ScalarVar_j[1];
119+
120+
if (implicit) {
121+
Jacobian_i[0][0] = a0; Jacobian_i[0][1] = 0.0;
122+
Jacobian_i[1][0] = 0.0; Jacobian_i[1][1] = a0;
123+
124+
Jacobian_j[0][0] = a1; Jacobian_j[0][1] = 0.0;
125+
Jacobian_j[1][0] = 0.0; Jacobian_j[1][1] = a1;
126+
}
127+
}
81128

82129
public:
83130
/*!
@@ -86,6 +133,6 @@ class CUpwSca_TurbSST final : public CUpwScalar {
86133
* \param[in] val_nVar - Number of variables of the problem.
87134
* \param[in] config - Definition of the particular problem.
88135
*/
89-
CUpwSca_TurbSST(unsigned short val_nDim, unsigned short val_nVar, const CConfig* config);
90-
136+
CUpwSca_TurbSST(unsigned short val_nDim, unsigned short val_nVar, const CConfig* config)
137+
: CUpwScalar<FlowIndices>(val_nDim, val_nVar, config) {}
91138
};

0 commit comments

Comments
 (0)