Fix bug in SA-neg diffusion term (and generalize indices of flow variables for use by scalar solvers and numerics)

SU2_CFD/include/drivers/CDriver.hpp

@@ -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).

SU2_CFD/include/numerics/scalar/scalar_convection.hpp

@@ -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;
 
@@ -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);
+  }
 };

SU2_CFD/include/numerics/scalar/scalar_diffusion.hpp

@@ -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;
 
@@ -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);
+  }
 };

SU2_CFD/include/numerics/turbulent/turb_convection.hpp

@@ -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:
   /*!
@@ -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) {}
 };
 
 /*!
@@ -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:
   /*!
@@ -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) {}
 };