Fix the second-order convection of turbulence and species

Common/include/CConfig.hpp

@@ -4859,6 +4859,15 @@ class CConfig {
    */
   bool GetMUSCL_Flow(void) const { return MUSCL_Flow; }
 
+  /*! \brief Whether the turbulence equations reconstruct. */
+  bool GetMUSCL_Turb(void) const { return MUSCL_Turb; }
+
+  /*! \brief Whether the species equations reconstruct. */
+  bool GetMUSCL_Species(void) const { return MUSCL_Species; }
+
+  /*! \brief Whether the heat equation reconstructs. */
+  bool GetMUSCL_Heat(void) const { return MUSCL_Heat; }
+
   /*!
    * \brief Get if the upwind scheme used MUSCL or not.
    * \return MUSCL scheme.

SU2_CFD/include/numerics/scalar/scalar_edge_flux.hpp

@@ -62,6 +62,12 @@ struct ScalarFluxOptions {
   bool viscous = false;   /*!< \brief Whether the diffusion term contributes. */
   bool oneSided = false;  /*!< \brief Whether only the row of i is assembled. */
   bool muscl = false;     /*!< \brief Whether the convective scheme reconstructs. */
+  /*!
+   * \brief Bounds of the transported variables (those that clip the solution). A reconstructed face
+   *        value outside them falls back to the value of its point (first order). Null: no bounds.
+   */
+  const su2double* lowerLimit = nullptr;
+  const su2double* upperLimit = nullptr;
 
   /*!
    * \brief Options of the interior edge loop: both terms, both rows, and reconstruction and
@@ -349,10 +355,9 @@ class CUpwScalarBase : public CUpwScalarFlux<Double_, Derived, FlowIndices, nDim
 
   /*!
    * \brief MUSCL reconstruction parameters, read from CConfig once per construction (i.e. once
-   *        per nonlinear iteration, see CScalarSolver::EdgeFluxResidual) instead of per edge;
-   *        this is also where the scalar limiter's freezing (GetLimiterIter) is resolved, by
-   *        collapsing its type to NONE once frozen. The flow limiter is not frozen this way: once
-   *        the flow solver stops recomputing it, it keeps applying the last values it has.
+   *        per nonlinear iteration, see CScalarSolver::EdgeFluxResidual) instead of per edge.
+   *        Point limiters keep applying their last stored values after LIMITER_ITER.
+   *        VAN_ALBADA_EDGE is computed on every edge and is never frozen.
    */
   const su2double kappa, umusclRamp, kappaFlow;
   const LIMITER limiterType, limiterTypeFlow;
@@ -372,9 +377,8 @@ class CUpwScalarBase : public CUpwScalarFlux<Double_, Derived, FlowIndices, nDim
         kappa(config.GetMUSCL_Kappa()),
         umusclRamp(config.GetMUSCLRampValue()),
         kappaFlow(config.GetMUSCL_Kappa_Flow()),
-        limiterType(config.GetInnerIter() <= config.GetLimiterIter() ? config.GetKind_SlopeLimit() : LIMITER::NONE),
-        limiterTypeFlow(config.GetKind_SlopeLimit_Flow() != LIMITER::VAN_ALBADA_EDGE ? config.GetKind_SlopeLimit_Flow()
-                                                                                     : LIMITER::NONE),
+        limiterType(config.GetKind_SlopeLimit()),
+        limiterTypeFlow(config.GetKind_SlopeLimit_Flow()),
         musclFlow(config.GetMUSCL_Flow() && config.GetKind_ConvNumScheme_Flow() == SPACE_UPWIND) {
     if (nEqn > Size) {
       SU2_MPI::Error("The number of equations exceeds MAX_TRANSPORTED_SPECIES, the size of the static arrays. Increase it\n"
@@ -474,14 +478,29 @@ class CUpwScalarBase : public CUpwScalarFlux<Double_, Derived, FlowIndices, nDim
       }
 
       if (opt.muscl) {
+        const auto phi1st = phi;
         if constexpr (nVar != Dynamic) {
-          reconstruct<nVar>(iPoint, jPoint, vector_ij, side_i.scalarNodes.GetGradient_Reconstruction(),
-                            side_i.scalarNodes.GetLimiter(), limiterType, 0, phi, kappa, umusclRamp);
+          reconstruct<nVar, MusclLimiter::EDGE_SCALAR>(
+              iPoint, jPoint, vector_ij, side_i.scalarNodes.GetGradient_Reconstruction(),
+              side_i.scalarNodes.GetLimiter(), limiterType, 0, phi, kappa, umusclRamp);
         } else {
           /*--- A dynamic model's equation count is only known at runtime, so the reconstructed
            * width is passed as an argument instead of a template parameter. ---*/
-          reconstruct(iPoint, jPoint, vector_ij, side_i.scalarNodes.GetGradient_Reconstruction(),
-                      side_i.scalarNodes.GetLimiter(), limiterType, 0, phi, kappa, umusclRamp, res.nVar);
+          reconstruct<0, MusclLimiter::EDGE_SCALAR>(
+              iPoint, jPoint, vector_ij, side_i.scalarNodes.GetGradient_Reconstruction(),
+              side_i.scalarNodes.GetLimiter(), limiterType, 0, phi, kappa, umusclRamp, res.nVar);
+        }
+
+        /*--- A face value outside the bounds of the variable (e.g. a negative k or omega next to a wall,
+         * where omega varies by orders of magnitude) falls back to the value of its point. ---*/
+        if (opt.lowerLimit) {
+          for (size_t iVar = 0; iVar < res.nVar; ++iVar) {
+            const su2double lo = opt.lowerLimit[iVar], hi = opt.upperLimit[iVar];
+            const Double out_i = fmax(phi.i.all(iVar) < lo, phi.i.all(iVar) > hi);
+            const Double out_j = fmax(phi.j.all(iVar) < lo, phi.j.all(iVar) > hi);
+            phi.i.all(iVar) = out_i * phi1st.i.all(iVar) + (1 - out_i) * phi.i.all(iVar);
+            phi.j.all(iVar) = out_j * phi1st.j.all(iVar) + (1 - out_j) * phi.j.all(iVar);
+          }
         }
       }
 

SU2_CFD/include/numerics/util.hpp

@@ -34,6 +34,8 @@
 #include "../../../Common/include/containers/C2DContainer.hpp"
 #include "../../../Common/include/linear_algebra/CSysVector.hpp"
 #include "../../../Common/include/linear_algebra/CSysMatrix.hpp"
+#include "../../../Common/include/geometry/CGeometry.hpp"
+#include "../limiters/CLimiterDetails.hpp"
 
 /*!
  * \enum UpdateType
@@ -498,8 +500,11 @@ FORCEINLINE auto gradientReader(typename CLaneTraits<Double>::Int iPoint, const
 
 /*!
  * \brief How the reconstructed differences are limited.
+ * \note EDGE is the van Albada edge limiter of the flow solver, which the scalar solvers also use to
+ *       reconstruct the flow velocity so that both see the same face velocity. EDGE_SCALAR is the
+ *       van Albada function of the point limiters, used for the transported scalar variables.
  */
-enum class MusclLimiter { NONE, EDGE, POINT };
+enum class MusclLimiter { NONE, EDGE, EDGE_SCALAR, POINT };
 
 /*!
  * \brief U-MUSCL reconstruction of nVarGrad variables, from the gradient rows starting at iRow.
@@ -534,6 +539,9 @@ FORCEINLINE void muscl(typename CLaneTraits<Double>::Int iPoint, typename CLaneT
       const Double delta_ij_2 = pow(delta_ij, 2) + 1e-6;
       lim_i = (delta_ij_2 + proj_i * delta_ij) / (pow(proj_i, 2) + delta_ij_2);
       lim_j = (delta_ij_2 + proj_j * delta_ij) / (pow(proj_j, 2) + delta_ij_2);
+    } else if constexpr (limiterKind == MusclLimiter::EDGE_SCALAR) {
+      lim_i = LimiterHelpers<Double>::vanAlbadaFunction(proj_i, delta_ij, Double(EPS));
+      lim_j = LimiterHelpers<Double>::vanAlbadaFunction(proj_j, delta_ij, Double(EPS));
     } else if constexpr (limiterKind == MusclLimiter::POINT) {
       lim_i = gatherVariables(iPoint, limiter, iRow + iVar);
       lim_j = gatherVariables(jPoint, limiter, iRow + iVar);
@@ -548,21 +556,25 @@ FORCEINLINE void muscl(typename CLaneTraits<Double>::Int iPoint, typename CLaneT
 /*!
  * \brief Reconstruct a slice of nVarGrad variables starting at column iRow, dispatching on the
  *        limiter type; shared by the flow and the scalar reconstructions so both call one body.
+ * \note edgeLimiter is the form that VAN_ALBADA_EDGE uses (EDGE for flow variables, EDGE_SCALAR for
+ *       transported scalars).
  */
-template <size_t nVarGrad_ = 0, size_t nDim, class Double, class VarType, class Limiter_t, class Gradient_t>
+template <size_t nVarGrad_ = 0, MusclLimiter edgeLimiter = MusclLimiter::EDGE, size_t nDim, class Double, class VarType,
+          class Limiter_t, class Gradient_t>
 FORCEINLINE void reconstruct(typename CLaneTraits<Double>::Int iPoint, typename CLaneTraits<Double>::Int jPoint,
                              const Vector<Double, nDim>& vector_ij, const Gradient_t& gradient,
                              const Limiter_t& limiter, LIMITER limiterType, size_t iRow, CPair<VarType>& V,
                              const CNonDeduced<Double>& kappa, const CNonDeduced<Double>& umusclRamp,
                              size_t nVarGradRuntime = 0) {
+  static_assert(edgeLimiter == MusclLimiter::EDGE || edgeLimiter == MusclLimiter::EDGE_SCALAR);
   switch (limiterType) {
     case LIMITER::NONE:
       muscl<MusclLimiter::NONE, nVarGrad_>(iPoint, jPoint, vector_ij, gradient, limiter, iRow, V, kappa, umusclRamp,
                                            nVarGradRuntime);
       break;
     case LIMITER::VAN_ALBADA_EDGE:
-      muscl<MusclLimiter::EDGE, nVarGrad_>(iPoint, jPoint, vector_ij, gradient, limiter, iRow, V, kappa, umusclRamp,
-                                           nVarGradRuntime);
+      muscl<edgeLimiter, nVarGrad_>(iPoint, jPoint, vector_ij, gradient, limiter, iRow, V, kappa, umusclRamp,
+                                    nVarGradRuntime);
       break;
     default:
       muscl<MusclLimiter::POINT, nVarGrad_>(iPoint, jPoint, vector_ij, gradient, limiter, iRow, V, kappa, umusclRamp,

SU2_CFD/include/solvers/CScalarSolver.hpp

@@ -190,7 +190,8 @@ class CScalarSolver : public CSolver {
           if (conv_numerics->GetBoundedScalar()) {
             const su2double* velocity = &PrimVar_j[prim_idx.Velocity()];
             const su2double density = solver_container[FLOW_SOL]->GetNodes()->GetDensity(iPoint);
-            conv_numerics->SetMassFlux(BoundedScalarBCFlux(iPoint, true, density, velocity, Normal));
+            conv_numerics->SetMassFlux(BoundedScalarBCFlux(iPoint, true, density, velocity, Normal,
+                                                         &PrimVar_j[prim_idx.Density()], weight));
           }
 
           auto residual = conv_numerics->ComputeResidual(config);
@@ -244,16 +245,35 @@ class CScalarSolver : public CSolver {
    * \brief Applies a convective flux correction to negate the effects of flow divergence at a BC node.
    * \note This function should be used for nodes that are part of a boundary marker, it computes a mass flux
    * from density and velocity at the node, and the outward-pointing normal (-1 * normal of vertex).
+   * \param[in] weight - Weight of the boundary state, e.g. of a donor at a sliding interface, which also weights
+   *            the flux (the returned mass flux is not weighted).
    * \return The mass flux.
    */
   inline su2double BoundedScalarBCFlux(unsigned long iPoint, bool implicit, const su2double& density,
-                                       const su2double* velocity, const su2double* normal) {
-    const su2double edgeMassFlux = density * GeometryToolbox::DotProduct(nDim, velocity, normal);
-    LinSysRes.AddBlock(iPoint, nodes->GetSolution(iPoint), -edgeMassFlux);
-    if (implicit) Jacobian.AddVal2Diag(iPoint, -edgeMassFlux);
+                                       const su2double* velocity, const su2double* normal,
+                                       const su2double* densityGhost = nullptr, su2double weight = 1.0) {
+    const su2double rhoGhost = densityGhost ? *densityGhost : density;
+    const su2double volumeFlux = GeometryToolbox::DotProduct(nDim, velocity, normal);
+    const su2double edgeMassFlux = volumeFlux * (volumeFlux > 0.0 ? density : rhoGhost);
+    const su2double q = weight * BoundedScalarDivergenceFlux(edgeMassFlux, density, rhoGhost);
+    LinSysRes.AddBlock(iPoint, nodes->GetSolution(iPoint), -q);
+    if (implicit) Jacobian.AddVal2Diag(iPoint, -q / (Conservative ? density : 1.0));
     return edgeMassFlux;
   }
 
+  /*!
+   * \brief Face flux of the bounded scalar divergence correction, -q * phi_i, such that it cancels the
+   *        part of the upwind flux due to the discrete divergence of the mass flux. It uses the same
+   *        coefficients as the flux: the mass flux for a conservative scalar (whose flux is
+   *        mdot+ phi_i + mdot- phi_j), the upwinded volume flux mdot+/rho_i + mdot-/rho_j otherwise.
+   * \note The Jacobian w.r.t. the conserved variable of a conservative scalar, rho * phi, is -q / rho_i.
+   */
+  inline su2double BoundedScalarDivergenceFlux(const su2double& massFlux, const su2double& rho_i,
+                                               const su2double& rho_j) const {
+    if (Conservative) return massFlux;
+    return fmax(massFlux, 0.0) / rho_i + fmin(massFlux, 0.0) / rho_j;
+  }
+
   /*!
    * \brief Move solution to previous time levels (for restarts).
    */

SU2_CFD/include/solvers/CScalarSolver.inl

@@ -118,7 +118,9 @@ void CScalarSolver<VariableType>::CommonPreprocessing(CGeometry *geometry, const
    * before calling these solver functions. ---*/
   const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT);
   const bool muscl = config->GetMUSCL();
+  /*--- The edge limiter (VAN_ALBADA_EDGE) is applied in the flux kernel and needs no point values. ---*/
   const bool limiter = (config->GetKind_SlopeLimit() != LIMITER::NONE) &&
+                       (config->GetKind_SlopeLimit() != LIMITER::VAN_ALBADA_EDGE) &&
                        (config->GetInnerIter() <= config->GetLimiterIter());
 
   /*--- Clear residual and system matrix, not needed for
@@ -132,9 +134,11 @@ void CScalarSolver<VariableType>::CommonPreprocessing(CGeometry *geometry, const
     }
   }
 
-  /*--- Upwind second order reconstruction and gradients ---*/
+  /*--- Upwind second order reconstruction and gradients. The reconstruction gradient is global to the
+   * case (it is needed as soon as any equation system reconstructs), so it is skipped here when this
+   * system does not reconstruct: only the reconstruction and the limiter read it. ---*/
 
-  if (config->GetReconstructionGradientRequired()) {
+  if (muscl && config->GetReconstructionGradientRequired()) {
     switch(config->GetKind_Gradient_Method_Recon()) {
       case GREEN_GAUSS: SetSolution_Gradient_GG(geometry, config, -1, true); break;
       case LEAST_SQUARES: SetSolution_Gradient_LS(geometry, config, -1, true); break;
@@ -174,9 +178,14 @@ void CScalarSolver<VariableType>::SumEdgeFluxes(const CGeometry* geometry) {
 template <class VariableType>
 template <class Scheme>
 void CScalarSolver<VariableType>::EdgeFluxResidual(const CGeometry* geometry, CSolver** solver_container,
-                                                    const CConfig* config, const ScalarFluxOptions& opt) {
+                                                    const CConfig* config, const ScalarFluxOptions& optIn) {
   SU2_ZONE_SCOPED
 
+  /*--- The reconstructed face values are kept within the bounds that clip the solution. ---*/
+  auto opt = optIn;
+  opt.lowerLimit = lowerlimit;
+  opt.upperLimit = upperlimit;
+
   using Double = typename Scheme::Double;
   constexpr int nDim = Scheme::nDim;
 
@@ -217,11 +226,13 @@ void CScalarSolver<VariableType>::EdgeFluxResidual(const CGeometry* geometry, CS
       /*--- Bounded scalar divergence correction, per edge; the ReducerStrategy equivalent runs
        * in a per-point pass below, where the diagonal is not written from the edge loop. ---*/
       if (opt.boundedScalar && !ReducerStrategy) {
-        LinSysRes.AddBlock(iPoint, nodes->GetSolution(iPoint), -massFlux);
-        LinSysRes.AddBlock(jPoint, nodes->GetSolution(jPoint), massFlux);
+        const su2double rho_i = flowNodes->GetDensity(iPoint), rho_j = flowNodes->GetDensity(jPoint);
+        const su2double q = BoundedScalarDivergenceFlux(massFlux, rho_i, rho_j);
+        LinSysRes.AddBlock(iPoint, nodes->GetSolution(iPoint), -q);
+        LinSysRes.AddBlock(jPoint, nodes->GetSolution(jPoint), q);
         if (opt.implicit) {
-          Jacobian.AddVal2Diag(iPoint, -massFlux);
-          Jacobian.AddVal2Diag(jPoint, massFlux);
+          Jacobian.AddVal2Diag(iPoint, -q / (Conservative ? rho_i : 1.0));
+          Jacobian.AddVal2Diag(jPoint, q / (Conservative ? rho_j : 1.0));
         }
       }
     }
@@ -239,15 +250,19 @@ void CScalarSolver<VariableType>::EdgeFluxResidual(const CGeometry* geometry, CS
       SU2_OMP_FOR_STAT(omp_chunk_size)
       for (unsigned long iPoint = 0; iPoint < nPoint; ++iPoint) {
         const auto* solution = nodes->GetSolution(iPoint);
+        const su2double rho_i = flowNodes->GetDensity(iPoint);
         su2double divergence = 0;
 
-        for (auto iEdge : geometry->nodes->GetEdges(iPoint)) {
+        for (auto iNeigh = 0u; iNeigh < geometry->nodes->GetnPoint(iPoint); ++iNeigh) {
+          const auto iEdge = geometry->nodes->GetEdge(iPoint, iNeigh);
+          const auto jPoint = geometry->nodes->GetPoint(iPoint, iNeigh);
           const auto sign = (iPoint == geometry->edges->GetNode(iEdge, 0)) ? 1 : -1;
-          const su2double edgeMassFlux = sign * (*edgeMassFluxes)[iEdge];
-          divergence += edgeMassFlux;
-          LinSysRes.AddBlock(iPoint, solution, -edgeMassFlux);
+          const su2double q = BoundedScalarDivergenceFlux(sign * (*edgeMassFluxes)[iEdge], rho_i,
+                                                          flowNodes->GetDensity(jPoint));
+          divergence += q;
+          LinSysRes.AddBlock(iPoint, solution, -q);
         }
-        if (opt.implicit) Jacobian.AddVal2Diag(iPoint, -divergence);
+        if (opt.implicit) Jacobian.AddVal2Diag(iPoint, -divergence / (Conservative ? rho_i : 1.0));
       }
       END_SU2_OMP_FOR
     }
@@ -282,8 +297,9 @@ void CScalarSolver<VariableType>::BoundaryFluxResidual(const CGeometry* geometry
 
     Double massFlux = 0.0;
     if (opt.boundedScalar) {
+      const auto* ghostPrim = ghostFlowNodes->GetPrimitive(iVertex);
       massFlux = BoundedScalarBCFlux(iPoint, opt.implicit, flowNodes->GetDensity(iPoint),
-                                     &ghostFlowNodes->GetPrimitive(iVertex)[prim_idx.Velocity()], normal.data());
+                                     &ghostPrim[prim_idx.Velocity()], normal.data(), &ghostPrim[prim_idx.Density()]);
     }
 
     const auto res = flux.ComputeFlux(opt, iPoint, side_i, iVertex, side_j, normal, massFlux);
@@ -347,7 +363,8 @@ void CScalarSolver<VariableType>::FluidInterfaceFluxResidual(const CGeometry* ge
         su2double massFlux = 0.0;
         if (optConv.boundedScalar) {
           massFlux = BoundedScalarBCFlux(iPoint, optConv.implicit, flowNodes->GetDensity(iPoint),
-                                         &PrimVar_j[prim_idx.Velocity()], normal.data());
+                                         &PrimVar_j[prim_idx.Velocity()], normal.data(), &PrimVar_j[prim_idx.Density()],
+                                         weight);
         }
 
         const auto res = flux.ComputeFlux(optConv, iPoint, side_i, iVertex, side_j, normal, massFlux);
@@ -446,7 +463,8 @@ void CScalarSolver<VariableType>::BC_Far_Field(CGeometry* geometry, CSolver** so
       if (conv_numerics->GetBoundedScalar()) {
         const su2double* velocity = &V_infty[prim_idx.Velocity()];
         const su2double density = solver_container[FLOW_SOL]->GetNodes()->GetDensity(iPoint);
-        conv_numerics->SetMassFlux(BoundedScalarBCFlux(iPoint, implicit, density, velocity, Normal));
+        conv_numerics->SetMassFlux(BoundedScalarBCFlux(iPoint, implicit, density, velocity, Normal,
+                                                     &V_infty[prim_idx.Density()]));
       }
 
       /*--- Compute residuals and Jacobians ---*/

SU2_CFD/include/variables/CScalarVariable.hpp

@@ -49,8 +49,9 @@ class CScalarVariable : public CVariable {
    * \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.
+   * \param[in] muscl - Whether this scalar system reconstructs.
    */
-  CScalarVariable(unsigned long npoint, unsigned long ndim, unsigned long nvar, const CConfig* config);
+  CScalarVariable(unsigned long npoint, unsigned long ndim, unsigned long nvar, const CConfig* config, bool muscl);
 
   /*!
    * \brief Get the array of the reconstruction variables gradient at a node.