Skip to content

Commit

Permalink
Merge pull request #1392 from su2code/feature_NEMO_turbulent
Browse files Browse the repository at this point in the history 
Fix bug in SA-neg diffusion term  (and generalize indices of flow variables for use by scalar solvers and numerics)
  • Loading branch information
@pcarruscag
pcarruscag authored Nov 8, 2021
2 parents 147dfac + 57a736a commit 641d254
Show file tree
Hide file tree
Showing 37 changed files with 818 additions and 800 deletions.
7 changes: 7 additions & 0 deletions SU2_CFD/include/drivers/CDriver.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -213,6 +213,13 @@ class CDriver {
*/
void Numerics_Preprocessing(CConfig *config, CGeometry **geometry, CSolver ***solver, CNumerics ****&numerics) const;

/*!
* \brief Helper to instantiate turbulence numerics specialized for different flow solvers.
*/
template <class FlowIndices>
void InstantiateTurbulentNumerics(unsigned short nVar_Turb, int offset, const CConfig *config,
const CSolver* turb_solver, CNumerics ****&numerics) const;

/*!
* \brief Definition and allocation of all solver classes.
* \param[in] numerics_container - Description of the numerical method (the way in which the equations are solved).
Expand Down
80 changes: 65 additions & 15 deletions SU2_CFD/include/numerics/scalar/scalar_convection.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -44,18 +44,24 @@
* \ingroup ConvDiscr
* \author C. Pederson, A. Bueno., and A. Campos.
*/
template <class FlowIndices>
class CUpwScalar : public CNumerics {
protected:
su2double a0 = 0.0; /*!< \brief The maximum of the face-normal velocity and 0 */
su2double a1 = 0.0; /*!< \brief The minimum of the face-normal velocity and 0 */
su2double* Flux = nullptr; /*!< \brief Final result, diffusive flux/residual. */
su2double** Jacobian_i = nullptr; /*!< \brief Flux Jacobian w.r.t. node i. */
su2double** Jacobian_j = nullptr; /*!< \brief Flux Jacobian w.r.t. node j. */
enum : unsigned short {MAXNVAR = 8};

const FlowIndices idx; /*!< \brief Object to manage the access to the flow primitives. */
su2double a0 = 0.0; /*!< \brief The maximum of the face-normal velocity and 0. */
su2double a1 = 0.0; /*!< \brief The minimum of the face-normal velocity and 0. */
su2double Flux[MAXNVAR]; /*!< \brief Final result, diffusive flux/residual. */
su2double* Jacobian_i[MAXNVAR]; /*!< \brief Flux Jacobian w.r.t. node i. */
su2double* Jacobian_j[MAXNVAR]; /*!< \brief Flux Jacobian w.r.t. node j. */
su2double JacobianBuffer[2*MAXNVAR*MAXNVAR]; /*!< \brief Static storage for the two Jacobians. */

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

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

Expand All @@ -69,21 +75,65 @@ class CUpwScalar : public CNumerics {
public:
/*!
* \brief Constructor of the class.
* \param[in] val_nDim - Number of dimensions of the problem.
* \param[in] val_nVar - Number of variables of the problem.
* \param[in] ndim - Number of dimensions of the problem.
* \param[in] nvar - Number of variables of the problem.
* \param[in] config - Definition of the particular problem.
*/
CUpwScalar(unsigned short val_nDim, unsigned short val_nVar, const CConfig* config);

/*!
* \brief Destructor of the class.
*/
~CUpwScalar(void) override;
CUpwScalar(unsigned short ndim, unsigned short nvar, const CConfig* config)
: CNumerics(ndim, nvar, config),
idx(ndim, config->GetnSpecies()),
implicit(config->GetKind_TimeIntScheme_Turb() == EULER_IMPLICIT),
incompressible(config->GetKind_Regime() == ENUM_REGIME::INCOMPRESSIBLE),
dynamic_grid(config->GetDynamic_Grid()) {
if (nVar > MAXNVAR) {
SU2_MPI::Error("Static arrays are too small.", CURRENT_FUNCTION);
}
for (unsigned short iVar = 0; iVar < nVar; iVar++) {
Jacobian_i[iVar] = &JacobianBuffer[iVar * nVar];
Jacobian_j[iVar] = &JacobianBuffer[iVar * nVar + MAXNVAR * MAXNVAR];
}
}

/*!
* \brief Compute the scalar upwind flux between two nodes i and j.
* \param[in] config - Definition of the particular problem.
* \return A lightweight const-view (read-only) of the residual/flux and Jacobians.
*/
ResidualType<> ComputeResidual(const CConfig* config) override;
CNumerics::ResidualType<> ComputeResidual(const CConfig* config) final {
AD::StartPreacc();
AD::SetPreaccIn(Normal, nDim);
AD::SetPreaccIn(ScalarVar_i, nVar);
AD::SetPreaccIn(ScalarVar_j, nVar);
if (dynamic_grid) {
AD::SetPreaccIn(GridVel_i, nDim);
AD::SetPreaccIn(GridVel_j, nDim);
}
AD::SetPreaccIn(&V_i[idx.Velocity()], nDim);
AD::SetPreaccIn(&V_j[idx.Velocity()], nDim);

ExtraADPreaccIn();

su2double q_ij = 0.0;
if (dynamic_grid) {
for (unsigned short iDim = 0; iDim < nDim; iDim++) {
su2double Velocity_i = V_i[iDim + idx.Velocity()] - GridVel_i[iDim];
su2double Velocity_j = V_j[iDim + idx.Velocity()] - GridVel_j[iDim];
q_ij += 0.5 * (Velocity_i + Velocity_j) * Normal[iDim];
}
} else {
for (unsigned short iDim = 0; iDim < nDim; iDim++) {
q_ij += 0.5 * (V_i[iDim + idx.Velocity()] + V_j[iDim + idx.Velocity()]) * Normal[iDim];
}
}

a0 = 0.5 * (q_ij + fabs(q_ij));
a1 = 0.5 * (q_ij - fabs(q_ij));

FinishResidualCalc(config);

AD::SetPreaccOut(Flux, nVar);
AD::EndPreacc();

return ResidualType<>(Flux, Jacobian_i, Jacobian_j);
}
};
72 changes: 59 additions & 13 deletions SU2_CFD/include/numerics/scalar/scalar_diffusion.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -44,14 +44,18 @@
* \ingroup ViscDiscr
* \author C. Pederson, A. Bueno, and F. Palacios
*/
template <class FlowIndices>
class CAvgGrad_Scalar : public CNumerics {
protected:
su2double* Proj_Mean_GradScalarVar_Normal = nullptr; /*!< \brief Mean_gradScalarVar DOT normal. */
su2double* Proj_Mean_GradScalarVar = nullptr; /*!< \brief Mean_gradScalarVar DOT normal, corrected if required. */
su2double proj_vector_ij = 0.0; /*!< \brief (Edge_Vector DOT normal)/|Edge_Vector|^2 */
su2double* Flux = nullptr; /*!< \brief Final result, diffusive flux/residual. */
su2double** Jacobian_i = nullptr; /*!< \brief Flux Jacobian w.r.t. node i. */
su2double** Jacobian_j = nullptr; /*!< \brief Flux Jacobian w.r.t. node j. */
enum : unsigned short {MAXNVAR = 8};

const FlowIndices idx; /*!< \brief Object to manage the access to the flow primitives. */
su2double Proj_Mean_GradScalarVar[MAXNVAR]; /*!< \brief Mean_gradScalarVar DOT normal, corrected if required. */
su2double proj_vector_ij = 0.0; /*!< \brief (Edge_Vector DOT normal)/|Edge_Vector|^2 */
su2double Flux[MAXNVAR]; /*!< \brief Final result, diffusive flux/residual. */
su2double* Jacobian_i[MAXNVAR]; /*!< \brief Flux Jacobian w.r.t. node i. */
su2double* Jacobian_j[MAXNVAR]; /*!< \brief Flux Jacobian w.r.t. node j. */
su2double JacobianBuffer[2*MAXNVAR*MAXNVAR];/*!< \brief Static storage for the two Jacobians. */

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

Expand All @@ -75,17 +79,59 @@ class CAvgGrad_Scalar : public CNumerics {
* \param[in] correct_gradient - Whether to correct gradient for skewness.
* \param[in] config - Definition of the particular problem.
*/
CAvgGrad_Scalar(unsigned short val_nDim, unsigned short val_nVar, bool correct_gradient, const CConfig* config);

/*!
* \brief Destructor of the class.
*/
~CAvgGrad_Scalar(void) override;
CAvgGrad_Scalar(unsigned short val_nDim, unsigned short val_nVar, bool correct_grad,
const CConfig* config)
: CNumerics(val_nDim, val_nVar, config),
idx(val_nDim, config->GetnSpecies()),
correct_gradient(correct_grad),
implicit(config->GetKind_TimeIntScheme_Turb() == EULER_IMPLICIT),
incompressible(config->GetKind_Regime() == ENUM_REGIME::INCOMPRESSIBLE) {
if (nVar > MAXNVAR) {
SU2_MPI::Error("Static arrays are too small.", CURRENT_FUNCTION);
}
for (unsigned short iVar = 0; iVar < nVar; iVar++) {
Jacobian_i[iVar] = &JacobianBuffer[iVar * nVar];
Jacobian_j[iVar] = &JacobianBuffer[iVar * nVar + MAXNVAR * MAXNVAR];
}
}

/*!
* \brief Compute the viscous residual using an average of gradients without correction.
* \param[in] config - Definition of the particular problem.
* \return A lightweight const-view (read-only) of the residual/flux and Jacobians.
*/
ResidualType<> ComputeResidual(const CConfig* config) override;
ResidualType<> ComputeResidual(const CConfig* config) final {
AD::StartPreacc();
AD::SetPreaccIn(Coord_i, nDim);
AD::SetPreaccIn(Coord_j, nDim);
AD::SetPreaccIn(Normal, nDim);
AD::SetPreaccIn(ScalarVar_Grad_i, nVar, nDim);
AD::SetPreaccIn(ScalarVar_Grad_j, nVar, nDim);
if (correct_gradient) {
AD::SetPreaccIn(ScalarVar_i, nVar);
AD::SetPreaccIn(ScalarVar_j, nVar);
}
AD::SetPreaccIn(V_i[idx.Density()], V_i[idx.LaminarViscosity()], V_i[idx.EddyViscosity()]);
AD::SetPreaccIn(V_j[idx.Density()], V_j[idx.LaminarViscosity()], V_j[idx.EddyViscosity()]);

Density_i = V_i[idx.Density()];
Density_j = V_j[idx.Density()];
Laminar_Viscosity_i = V_i[idx.LaminarViscosity()];
Laminar_Viscosity_j = V_j[idx.LaminarViscosity()];
Eddy_Viscosity_i = V_i[idx.EddyViscosity()];
Eddy_Viscosity_j = V_j[idx.EddyViscosity()];

ExtraADPreaccIn();

su2double ProjGradScalarVarNoCorr[MAXNVAR];
proj_vector_ij = ComputeProjectedGradient(nDim, nVar, Normal, Coord_i, Coord_j, ScalarVar_Grad_i, ScalarVar_Grad_j,
correct_gradient, ScalarVar_i, ScalarVar_j, ProjGradScalarVarNoCorr,
Proj_Mean_GradScalarVar);
FinishResidualCalc(config);

AD::SetPreaccOut(Flux, nVar);
AD::EndPreacc();

return ResidualType<>(Flux, Jacobian_i, Jacobian_j);
}
};
69 changes: 58 additions & 11 deletions SU2_CFD/include/numerics/turbulent/turb_convection.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,18 +36,36 @@
* \ingroup ConvDiscr
* \author A. Bueno.
*/
class CUpwSca_TurbSA final : public CUpwScalar {
template <class FlowIndices>
class CUpwSca_TurbSA final : public CUpwScalar<FlowIndices> {
private:
using Base = CUpwScalar<FlowIndices>;
using Base::a0;
using Base::a1;
using Base::Flux;
using Base::Jacobian_i;
using Base::Jacobian_j;
using Base::ScalarVar_i;
using Base::ScalarVar_j;
using Base::implicit;

/*!
* \brief Adds any extra variables to AD
* \brief Adds any extra variables to AD.
*/
void ExtraADPreaccIn() override;
void ExtraADPreaccIn() override {}

/*!
* \brief SA specific steps in the ComputeResidual method
* \param[in] config - Definition of the particular problem.
*/
void FinishResidualCalc(const CConfig* config) override;
void FinishResidualCalc(const CConfig* config) override {
Flux[0] = a0*ScalarVar_i[0] + a1*ScalarVar_j[0];

if (implicit) {
Jacobian_i[0][0] = a0;
Jacobian_j[0][0] = a1;
}
}

public:
/*!
Expand All @@ -56,8 +74,8 @@ class CUpwSca_TurbSA final : public CUpwScalar {
* \param[in] val_nVar - Number of variables of the problem.
* \param[in] config - Definition of the particular problem.
*/
CUpwSca_TurbSA(unsigned short val_nDim, unsigned short val_nVar, const CConfig* config);

CUpwSca_TurbSA(unsigned short val_nDim, unsigned short val_nVar, const CConfig* config)
: CUpwScalar<FlowIndices>(val_nDim, val_nVar, config) {}
};

/*!
Expand All @@ -66,18 +84,47 @@ class CUpwSca_TurbSA final : public CUpwScalar {
* \ingroup ConvDiscr
* \author A. Campos.
*/
class CUpwSca_TurbSST final : public CUpwScalar {
template <class FlowIndices>
class CUpwSca_TurbSST final : public CUpwScalar<FlowIndices> {
private:
using Base = CUpwScalar<FlowIndices>;
using Base::nDim;
using Base::V_i;
using Base::V_j;
using Base::a0;
using Base::a1;
using Base::Flux;
using Base::Jacobian_i;
using Base::Jacobian_j;
using Base::ScalarVar_i;
using Base::ScalarVar_j;
using Base::implicit;
using Base::idx;

/*!
* \brief Adds any extra variables to AD
*/
void ExtraADPreaccIn() override;
void ExtraADPreaccIn() override {
AD::SetPreaccIn(V_i[idx.Density()]);
AD::SetPreaccIn(V_j[idx.Density()]);
}

/*!
* \brief SST specific steps in the ComputeResidual method
* \param[in] config - Definition of the particular problem.
*/
void FinishResidualCalc(const CConfig* config) override;
void FinishResidualCalc(const CConfig* config) override {
Flux[0] = a0*V_i[idx.Density()]*ScalarVar_i[0] + a1*V_j[idx.Density()]*ScalarVar_j[0];
Flux[1] = a0*V_i[idx.Density()]*ScalarVar_i[1] + a1*V_j[idx.Density()]*ScalarVar_j[1];

if (implicit) {
Jacobian_i[0][0] = a0; Jacobian_i[0][1] = 0.0;
Jacobian_i[1][0] = 0.0; Jacobian_i[1][1] = a0;

Jacobian_j[0][0] = a1; Jacobian_j[0][1] = 0.0;
Jacobian_j[1][0] = 0.0; Jacobian_j[1][1] = a1;
}
}

public:
/*!
Expand All @@ -86,6 +133,6 @@ class CUpwSca_TurbSST final : public CUpwScalar {
* \param[in] val_nVar - Number of variables of the problem.
* \param[in] config - Definition of the particular problem.
*/
CUpwSca_TurbSST(unsigned short val_nDim, unsigned short val_nVar, const CConfig* config);

CUpwSca_TurbSST(unsigned short val_nDim, unsigned short val_nVar, const CConfig* config)
: CUpwScalar<FlowIndices>(val_nDim, val_nVar, config) {}
};
Loading

0 comments on commit 641d254

Please sign in to comment.