Consistent mass fluxes between bulk flow and scalar solvers and boundedness correction

SU2_CFD/include/solvers/CScalarSolver.hpp

@@ -31,6 +31,8 @@
 #include "../../../Common/include/parallelization/omp_structure.hpp"
 #include "../../../Common/include/toolboxes/geometry_toolbox.hpp"
 #include "../variables/CScalarVariable.hpp"
+#include "../variables/CFlowVariable.hpp"
+#include "../variables/CPrimitiveIndices.hpp"
 #include "CSolver.hpp"
 
 /*!
@@ -50,14 +52,15 @@ class CScalarSolver : public CSolver {
 
   unsigned long omp_chunk_size; /*!< \brief Chunk size used in light point loops. */
 
-  su2double lowerlimit[MAXNVAR]; /*!< \brief contains lower limits for turbulence variables. Note that ::min()
-                                             returns the smallest positive value for floats. */
-  su2double upperlimit[MAXNVAR]; /*!< \brief contains upper limits for turbulence variables. */
+  su2double lowerlimit[MAXNVAR]; /*!< \brief contains lower limits for scalar variables. */
+  su2double upperlimit[MAXNVAR]; /*!< \brief contains upper limits for scalar variables. */
 
   su2double Solution_Inf[MAXNVAR]; /*!< \brief Far-field solution. */
 
   const bool Conservative; /*!< \brief Transported Variable is conservative. Solution has to be multiplied with rho. */
 
+  const CPrimitiveIndices<unsigned short> prim_idx; /*!< \brief Indices of the primitive flow variables. */
+
   vector<su2matrix<su2double*> > SlidingState; // vector of matrix of pointers... inner dim alloc'd elsewhere (welcome, to the twilight zone)
   vector<vector<int> > SlidingStateNodes;
 
@@ -86,7 +89,7 @@ class CScalarSolver : public CSolver {
   inline CVariable* GetBaseClassPointerToNodes() final { return nodes; }
 
   /*!
-   * \brief Compute the viscous flux for the turbulent equation at a particular edge.
+   * \brief Compute the viscous flux for the scalar equation at a particular edge.
    * \tparam SolverSpecificNumericsFunc - lambda-function, that implements solver specific contributions to numerics.
    * \note The functor has to implement (iPoint, jPoint)
    * \param[in] iEdge - Edge for which we want to compute the flux
@@ -100,7 +103,7 @@ class CScalarSolver : public CSolver {
                                          CGeometry* geometry, CSolver** solver_container, CNumerics* numerics,
                                          CConfig* config) {
     const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT);
-    CVariable* flowNodes = solver_container[FLOW_SOL]->GetNodes();
+    auto* flowNodes = su2staticcast_p<CFlowVariable*>(solver_container[FLOW_SOL]->GetNodes());
 
     /*--- Points in edge ---*/
 
@@ -177,7 +180,7 @@ class CScalarSolver : public CSolver {
 
         su2double* PrimVar_i = solver_container[FLOW_SOL]->GetNodes()->GetPrimitive(iPoint);
 
-        auto Jacobian_i = Jacobian.GetBlock(iPoint,iPoint);
+        auto* Jacobian_i = Jacobian.GetBlock(iPoint, iPoint);
 
         /*--- Loop over the nDonorVertexes and compute the averaged flux ---*/
 
@@ -208,6 +211,12 @@ class CScalarSolver : public CSolver {
           if (dynamic_grid)
             conv_numerics->SetGridVel(geometry->nodes->GetGridVel(iPoint), geometry->nodes->GetGridVel(iPoint));
 
+          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));
+          }
+
           auto residual = conv_numerics->ComputeResidual(config);
 
           /*--- Accumulate the residuals to compute the average ---*/
@@ -255,6 +264,20 @@ 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-poiting normal (-1 * normal of vertex).
+   * \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);
+    return edgeMassFlux;
+  }
+
   /*!
    * \brief Gradient and Limiter computation.
    * \param[in] geometry - Geometrical definition of the problem.
@@ -265,7 +288,7 @@ class CScalarSolver : public CSolver {
 
  private:
   /*!
-   * \brief Compute the viscous flux for the turbulent equation at a particular edge.
+   * \brief Compute the viscous flux for the scalar equation at a particular edge.
    * \param[in] iEdge - Edge for which we want to compute the flux
    * \param[in] geometry - Geometrical definition of the problem.
    * \param[in] solver_container - Container vector with all the solutions.
@@ -319,8 +342,20 @@ class CScalarSolver : public CSolver {
    * \param[in] config - Definition of the particular problem.
    * \param[in] iMesh - Index of the mesh in multigrid computations.
    */
-  void Upwind_Residual(CGeometry* geometry, CSolver** solver_container, CNumerics** numerics_container, CConfig* config,
-                       unsigned short iMesh) override;
+  void Upwind_Residual(CGeometry* geometry, CSolver** solver_container, CNumerics** numerics_container,
+                       CConfig* config, unsigned short iMesh) override;
+
+  /*!
+   * \brief Impose the Far Field boundary condition.
+   * \param[in] geometry - Geometrical definition of the problem.
+   * \param[in] solver_container - Container vector with all the solutions.
+   * \param[in] conv_numerics - Description of the numerical method.
+   * \param[in] visc_numerics - Description of the numerical method.
+   * \param[in] config - Definition of the particular problem.
+   * \param[in] val_marker - Surface marker where the boundary condition is applied.
+   */
+  void BC_Far_Field(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics,
+                    CNumerics *visc_numerics, CConfig *config, unsigned short val_marker) final;
 
   /*!
    * \brief Impose the Symmetry Plane boundary condition.

SU2_CFD/include/solvers/CScalarSolver.inl

@@ -31,7 +31,9 @@
 
 template <class VariableType>
 CScalarSolver<VariableType>::CScalarSolver(CGeometry* geometry, CConfig* config, bool conservative)
-    : CSolver(), Conservative(conservative) {
+    : CSolver(), Conservative(conservative),
+      prim_idx(config->GetKind_Regime() == ENUM_REGIME::INCOMPRESSIBLE,
+               config->GetNEMOProblem(), geometry->GetnDim(), config->GetnSpecies()) {
   nMarker = config->GetnMarker_All();
 
   /*--- Store the number of vertices on each marker for deallocation later ---*/
@@ -281,13 +283,9 @@ void CScalarSolver<VariableType>::Upwind_Residual(CGeometry* geometry, CSolver**
        * to avoid race conditions in accessing nodes shared by edges handled by different threads. ---*/
 
       if (bounded_scalar && !ReducerStrategy) {
-        for (iVar = 0; iVar < nVar; iVar++) {
-          const su2double FluxCorrection_i = nodes->GetSolution(iPoint, iVar) * EdgeMassFlux;
-          const su2double FluxCorrection_j = nodes->GetSolution(jPoint, iVar) * EdgeMassFlux;
+        LinSysRes.AddBlock(iPoint, nodes->GetSolution(iPoint), -EdgeMassFlux);
+        LinSysRes.AddBlock(jPoint, nodes->GetSolution(jPoint), EdgeMassFlux);
 
-          LinSysRes(iPoint, iVar) -= FluxCorrection_i;
-          LinSysRes(jPoint, iVar) += FluxCorrection_j;
-        }
         if (implicit) {
           Jacobian.AddVal2Diag(iPoint, -EdgeMassFlux);
           Jacobian.AddVal2Diag(jPoint, EdgeMassFlux);
@@ -363,6 +361,67 @@ void CScalarSolver<VariableType>::BC_Periodic(CGeometry* geometry, CSolver** sol
   }
 }
 
+template <class VariableType>
+void CScalarSolver<VariableType>::BC_Far_Field(CGeometry* geometry, CSolver** solver_container,
+                                               CNumerics* conv_numerics, CNumerics*, CConfig *config,
+                                               unsigned short val_marker) {
+
+  const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT);
+
+  SU2_OMP_FOR_STAT(OMP_MIN_SIZE)
+  for (auto iVertex = 0u; iVertex < geometry->nVertex[val_marker]; iVertex++) {
+
+    const auto iPoint = geometry->vertex[val_marker][iVertex]->GetNode();
+
+    /*--- Check if the node belongs to the domain (i.e, not a halo node) ---*/
+
+    if (geometry->nodes->GetDomain(iPoint)) {
+
+      /*--- Allocate the value at the infinity ---*/
+
+      auto V_infty = solver_container[FLOW_SOL]->GetCharacPrimVar(val_marker, iVertex);
+
+      /*--- Retrieve solution at the farfield boundary node ---*/
+
+      auto V_domain = solver_container[FLOW_SOL]->GetNodes()->GetPrimitive(iPoint);
+
+      /*--- Grid Movement ---*/
+
+      if (dynamic_grid)
+        conv_numerics->SetGridVel(geometry->nodes->GetGridVel(iPoint), geometry->nodes->GetGridVel(iPoint));
+
+      conv_numerics->SetPrimitive(V_domain, V_infty);
+
+      /*--- Set turbulent variable at the wall, and at infinity ---*/
+
+      conv_numerics->SetScalarVar(nodes->GetSolution(iPoint), Solution_Inf);
+
+      /*--- Set Normal (it is necessary to change the sign) ---*/
+
+      su2double Normal[MAXNDIM] = {0.0};
+      for (auto iDim = 0u; iDim < nDim; iDim++)
+        Normal[iDim] = -geometry->vertex[val_marker][iVertex]->GetNormal(iDim);
+      conv_numerics->SetNormal(Normal);
+
+      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));
+      }
+
+      /*--- Compute residuals and Jacobians ---*/
+
+      auto residual = conv_numerics->ComputeResidual(config);
+
+      /*--- Add residuals and Jacobians ---*/
+
+      LinSysRes.AddBlock(iPoint, residual);
+      if (implicit) Jacobian.AddBlock2Diag(iPoint, residual.jacobian_i);
+    }
+  }
+  END_SU2_OMP_FOR
+}
+
 template <class VariableType>
 void CScalarSolver<VariableType>::PrepareImplicitIteration(CGeometry* geometry, CSolver** solver_container,
                                                            CConfig* config) {

SU2_CFD/include/solvers/CTransLMSolver.hpp

@@ -153,23 +153,7 @@ class CTransLMSolver final : public CTurbSolver {
                           CNumerics *conv_numerics,
                           CNumerics *visc_numerics,
                           CConfig *config,
-                          unsigned short val_marker) override; 
-
-  /*!
-   * \brief Impose the Far Field boundary condition.
-   * \param[in] geometry - Geometrical definition of the problem.
-   * \param[in] solver_container - Container vector with all the solutions.
-   * \param[in] conv_numerics - Description of the numerical method.
-   * \param[in] visc_numerics - Description of the numerical method.
-   * \param[in] config - Definition of the particular problem.
-   * \param[in] val_marker - Surface marker where the boundary condition is applied.
-   */
-  void BC_Far_Field(CGeometry *geometry,
-                    CSolver **solver_container,
-                    CNumerics *conv_numerics,
-                    CNumerics *visc_numerics,
-                    CConfig *config,
-                    unsigned short val_marker) override;
+                          unsigned short val_marker) override;
 
   /*!
    * \brief Impose the inlet boundary condition.

SU2_CFD/include/solvers/CTurbSASolver.hpp

@@ -188,22 +188,6 @@ class CTurbSASolver final : public CTurbSolver {
                           CConfig *config,
                           unsigned short val_marker) override;
 
-  /*!
-   * \brief Impose the Far Field boundary condition.
-   * \param[in] geometry - Geometrical definition of the problem.
-   * \param[in] solver_container - Container vector with all the solutions.
-   * \param[in] conv_numerics - Description of the numerical method.
-   * \param[in] visc_numerics - Description of the numerical method.
-   * \param[in] config - Definition of the particular problem.
-   * \param[in] val_marker - Surface marker where the boundary condition is applied.
-   */
-  void BC_Far_Field(CGeometry *geometry,
-                    CSolver **solver_container,
-                    CNumerics *conv_numerics,
-                    CNumerics *visc_numerics,
-                    CConfig *config,
-                    unsigned short val_marker) override;
-
   /*!
    * \brief Impose the inlet boundary condition.
    * \param[in] geometry - Geometrical definition of the problem.

SU2_CFD/include/solvers/CTurbSSTSolver.hpp

@@ -169,22 +169,6 @@ class CTurbSSTSolver final : public CTurbSolver {
                           CConfig *config,
                           unsigned short val_marker) override;
 
-  /*!
-   * \brief Impose the Far Field boundary condition.
-   * \param[in] geometry - Geometrical definition of the problem.
-   * \param[in] solver_container - Container vector with all the solutions.
-   * \param[in] conv_numerics - Description of the numerical method.
-   * \param[in] visc_numerics - Description of the numerical method.
-   * \param[in] config - Definition of the particular problem.
-   * \param[in] val_marker - Surface marker where the boundary condition is applied.
-   */
-  void BC_Far_Field(CGeometry *geometry,
-                    CSolver **solver_container,
-                    CNumerics *conv_numerics,
-                    CNumerics *visc_numerics,
-                    CConfig *config,
-                    unsigned short val_marker) override;
-
   /*!
    * \brief Impose the inlet boundary condition.
    * \param[in] geometry - Geometrical definition of the problem.