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

Common/include/option_structure.hpp

@@ -826,7 +826,8 @@ enum ENUM_UPWIND {
   LAX_FRIEDRICH = 15,         /*!< \brief Lax-Friedrich numerical method. */
   AUSMPLUSUP = 16,            /*!< \brief AUSM+ -up numerical method (All Speed) */
   AUSMPLUSUP2 = 17,            /*!< \brief AUSM+ -up2 numerical method (All Speed) */
-  AUSMPWPLUS = 18            /*!< \brief AUSMplus numerical method. (MAYBE for TNE2 ONLY)*/
+  AUSMPWPLUS = 18,            /*!< \brief AUSMplus numerical method. (MAYBE for TNE2 ONLY)*/
+  SCALAR_ADVECTION = 19
 };
 static const MapType<std::string, ENUM_UPWIND> Upwind_Map = {
   MakePair("NONE", NO_UPWIND)
@@ -842,6 +843,7 @@ static const MapType<std::string, ENUM_UPWIND> Upwind_Map = {
   MakePair("MSW", MSW)
   MakePair("CUSP", CUSP)
   MakePair("SCALAR_UPWIND", SCALAR_UPWIND)
+  MakePair("SCALAR_ADVECTION", SCALAR_ADVECTION)
   MakePair("CONVECTIVE_TEMPLATE", CONVECTIVE_TEMPLATE)
   MakePair("L2ROE", L2ROE)
   MakePair("LMROE", LMROE)

SU2_CFD/include/numerics/CNumerics.hpp

@@ -55,6 +55,7 @@ class CNumerics {
   su2double Gas_Constant;     /*!< \brief Gas constant. */
   su2double Prandtl_Lam;      /*!< \brief Laminar Prandtl's number. */
   su2double Prandtl_Turb;     /*!< \brief Turbulent Prandtl's number. */
+  su2double MassFlux;         /*!< \brief Mass flux across edge. */
   su2double
   *Proj_Flux_Tensor;  /*!< \brief Flux tensor projected in a direction. */
   su2double **tau;    /*!< \brief Viscous stress tensor. */
@@ -1598,6 +1599,13 @@ class CNumerics {
    * \param[in] SolverSPvals - Struct holding the values.
    */
   virtual void SetStreamwisePeriodicValues(const StreamwisePeriodicValues SolverSPvals) { }
+
+  /*!
+   * \brief SetMassFlux
+   * \param[in] val_MassFlux: Mass flux across the edge
+   */
+  inline void SetMassFlux(const su2double val_MassFlux) {MassFlux = val_MassFlux;}
+
 };
 
 /*!

SU2_CFD/include/numerics/scalar/scalar_convection.hpp

@@ -118,6 +118,7 @@ class CUpwScalar : public CNumerics {
 
     ExtraADPreaccIn();
 
+    /*    
     su2double q_ij = 0.0;
     if (dynamic_grid) {
       for (unsigned short iDim = 0; iDim < nDim; iDim++) {
@@ -130,9 +131,12 @@ class CUpwScalar : public CNumerics {
         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));
+    */
+
+    //a0 = 0.5 * (q_ij + fabs(q_ij));
+    //a1 = 0.5 * (q_ij - fabs(q_ij));
+    a0 = max(0.0, MassFlux);
+    a1 = min(0.0, MassFlux);
 
     FinishResidualCalc(config);
 

SU2_CFD/include/numerics/species/species_convection.hpp

@@ -66,7 +66,7 @@ class CUpwSca_Species final : public CUpwScalar<FlowIndices> {
    */
   void FinishResidualCalc(const CConfig* config) override {
     for (auto iVar = 0u; iVar < nVar; iVar++) {
-      Flux[iVar] = a0 * V_i[idx.Density()] * ScalarVar_i[iVar] + a1 * V_j[idx.Density()] * ScalarVar_j[iVar];
+      Flux[iVar] = a0 * ScalarVar_i[iVar] + a1 * ScalarVar_j[iVar];
 
       /*--- Jacobians are taken wrt rho*Y not Y alone in the species solver. ---*/
       /*--- Off-diagonal entries are zero. ---*/

SU2_CFD/include/solvers/CScalarSolver.hpp

@@ -112,6 +112,9 @@ class CScalarSolver : public CSolver {
 
     numerics->SetPrimitive(flowNodes->GetPrimitive(iPoint), flowNodes->GetPrimitive(jPoint));
 
+    su2double EdgeMassFlux = solver_container[FLOW_SOL]->GetEdgeMassFlux(iEdge);
+    numerics->SetMassFlux(EdgeMassFlux);
+
     /*--- Turbulent variables w/o reconstruction, and its gradients ---*/
 
     numerics->SetScalarVar(nodes->GetSolution(iPoint), nodes->GetSolution(jPoint));
@@ -122,7 +125,7 @@ class CScalarSolver : public CSolver {
     SolverSpecificNumerics(iPoint, jPoint);
 
     /*--- Compute residual, and Jacobians ---*/
-
+    
     auto residual = numerics->ComputeResidual(config);
 
     if (ReducerStrategy) {

SU2_CFD/include/solvers/CScalarSolver.inl

@@ -119,6 +119,20 @@ template <class VariableType>
 void CScalarSolver<VariableType>::Upwind_Residual(CGeometry* geometry, CSolver** solver_container,
                                                   CNumerics** numerics_container, CConfig* config,
                                                   unsigned short iMesh) {
+
+  su2double EdgeMassFlux, Project_Grad_i, Project_Grad_j, FluxCorrection_i, FluxCorrection_j;
+
+  su2double **Jacobian_i_correction, **Jacobian_j_correction;
+  Jacobian_i_correction = new su2double*[nVar];
+  Jacobian_j_correction = new su2double*[nVar];
+  for (unsigned short iVar=0; iVar<nVar; iVar++){
+    Jacobian_i_correction[iVar] = new su2double[nVar];
+    Jacobian_j_correction[iVar] = new su2double[nVar];
+    for(unsigned short jVar=0; jVar<nVar; jVar++){
+      Jacobian_i_correction[iVar][jVar] = 0.0;
+      Jacobian_j_correction[iVar][jVar] = 0.0;
+    }
+  }
   /*--- Define booleans that are solver specific through CConfig's GlobalParams which have to be set in CFluidIteration
    * before calling these solver functions. ---*/
   const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT);
@@ -149,6 +163,10 @@ void CScalarSolver<VariableType>::Upwind_Residual(CGeometry* geometry, CSolver**
   else
     AD::StartNoSharedReading();
 
+  bool advection{false};
+
+  if(config->GetKind_Upwind_Species() == SCALAR_ADVECTION) advection = true;
+
   /*--- Loop over edge colors. ---*/
   for (auto color : EdgeColoring) {
     /*--- Chunk size is at least OMP_MIN_SIZE and a multiple of the color group size. ---*/
@@ -164,7 +182,7 @@ void CScalarSolver<VariableType>::Upwind_Residual(CGeometry* geometry, CSolver**
       auto jPoint = geometry->edges->GetNode(iEdge, 1);
 
       numerics->SetNormal(geometry->edges->GetNormal(iEdge));
-
+ 
       /*--- Primitive variables w/o reconstruction ---*/
 
       const auto V_i = flowNodes->GetPrimitive(iPoint);
@@ -204,7 +222,8 @@ void CScalarSolver<VariableType>::Upwind_Residual(CGeometry* geometry, CSolver**
           }
 
           for (iVar = 0; iVar < solver_container[FLOW_SOL]->GetnPrimVarGrad(); iVar++) {
-            su2double Project_Grad_i = 0.0, Project_Grad_j = 0.0;
+            Project_Grad_i = 0.0;
+            Project_Grad_j = 0.0;
             for (iDim = 0; iDim < nDim; iDim++) {
               Project_Grad_i += Vector_ij[iDim] * Gradient_i[iVar][iDim];
               Project_Grad_j -= Vector_ij[iDim] * Gradient_j[iVar][iDim];
@@ -249,8 +268,11 @@ void CScalarSolver<VariableType>::Upwind_Residual(CGeometry* geometry, CSolver**
         }
       }
 
-      /*--- Update convective residual value ---*/
+      /*--- Convective flux ---*/
+      EdgeMassFlux = solver_container[FLOW_SOL]->GetEdgeMassFlux(iEdge);
+      numerics->SetMassFlux(EdgeMassFlux);
 
+      /*--- Update convective residual value ---*/
       auto residual = numerics->ComputeResidual(config);
 
       if (ReducerStrategy) {
@@ -262,6 +284,18 @@ void CScalarSolver<VariableType>::Upwind_Residual(CGeometry* geometry, CSolver**
         if (implicit) Jacobian.UpdateBlocks(iEdge, iPoint, jPoint, residual.jacobian_i, residual.jacobian_j);
       }
 
+      /*--- Applying convective flux correction to negate the effects of flow divergence ---*/
+      if(advection){
+        for(iVar=0; iVar<nVar; iVar++){
+          FluxCorrection_i = GetNodes()->GetSolution(iPoint, iVar) * EdgeMassFlux;
+          FluxCorrection_j = GetNodes()->GetSolution(jPoint, iVar) * EdgeMassFlux;
+
+          LinSysRes(iPoint, iVar) -= FluxCorrection_i;
+          LinSysRes(jPoint, iVar) += FluxCorrection_j;
+          Jacobian_i_correction[iVar][iVar] = max(0.0, EdgeMassFlux);;
+          Jacobian_j_correction[iVar][iVar] = min(0.0, EdgeMassFlux);
+        }
+      }if (implicit) Jacobian.UpdateBlocks(iEdge, iPoint, jPoint, Jacobian_i_correction, Jacobian_j_correction);
       /*--- Viscous contribution. ---*/
 
       Viscous_Residual(iEdge, geometry, solver_container,
@@ -278,6 +312,13 @@ void CScalarSolver<VariableType>::Upwind_Residual(CGeometry* geometry, CSolver**
     SumEdgeFluxes(geometry);
     if (implicit) Jacobian.SetDiagonalAsColumnSum();
   }
+
+  for(unsigned short iVar=0; iVar<nVar; iVar++){
+    delete [] Jacobian_i_correction[iVar];
+    delete [] Jacobian_j_correction[iVar];
+  }
+  delete [] Jacobian_i_correction;
+  delete [] Jacobian_j_correction;
 }
 
 template <class VariableType>

SU2_CFD/include/solvers/CSolver.hpp

@@ -127,6 +127,8 @@ class CSolver {
   **Jacobian_ji,            /*!< \brief Auxiliary matrices for storing point to point Jacobians. */
   **Jacobian_jj;            /*!< \brief Auxiliary matrices for storing point to point Jacobians. */
 
+  su2double *EdgeMassFluxes;  /*!< \brief Mass fluxes across each edge, for discretization of passive scalars. */
+
   /*--- End variables that need to go. ---*/
 
   su2activevector iPoint_UndLapl;  /*!< \brief Auxiliary variable for the undivided Laplacians. */
@@ -4311,9 +4313,16 @@ class CSolver {
    * \param[in] geometry - Geometrical definition of the problem.
    * \param[in] config - Definition of the particular problem.
    * \param[in] converged - Whether or not solution has converged.
-  */
+   */
   void SavelibROM(CGeometry *geometry, CConfig *config, bool converged);
 
+  /*!
+   * \brief Get the mass flux across an edge (computed and stored during the discretization of convective fluxes).
+   * \param[in] iEdge - Index of the edge.
+   * \return The mass flux across the edge.
+   */
+  inline su2double GetEdgeMassFlux(const unsigned long iEdge) const {return EdgeMassFluxes[iEdge];}
+
 protected:
   /*!
    * \brief Allocate the memory for the verification solution, if necessary.

SU2_CFD/src/solvers/CIncEulerSolver.cpp

@@ -1316,6 +1316,13 @@ void CIncEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_cont
 
     Viscous_Residual(iEdge, geometry, solver_container,
                      numerics_container[VISC_TERM + omp_get_thread_num()*MAX_TERMS], config);
+
+    //if (rans && (iMesh == MESH_0)) EdgeMassFluxes[iEdge] = Res_Conv[0];
+    /* only when we solve additional scalars? So turbulence and species? */
+    if (iMesh == MESH_0) {
+      EdgeMassFluxes[iEdge] = residual[0];
+    }
+
   }
   END_SU2_OMP_FOR
   } // end color loop

SU2_CFD/src/solvers/CIncNSSolver.cpp

@@ -43,6 +43,16 @@ CIncNSSolver::CIncNSSolver(CGeometry *geometry, CConfig *config, unsigned short
   Viscosity_Inf   = config->GetViscosity_FreeStreamND();
   Tke_Inf         = config->GetTke_FreeStreamND();
 
+  //if (rans && (iMesh == MESH_0)) {
+  // only when we solve additional scalars?    
+  if (iMesh == MESH_0) {
+    EdgeMassFluxes = new su2double [geometry->GetnEdge()];
+
+    for(unsigned long iEdge = 0; iEdge < geometry->GetnEdge(); ++iEdge)
+      EdgeMassFluxes[iEdge] = 0.0;
+  }
+
+
   /*--- Initialize the secondary values for direct derivative approximations ---*/
 
   switch (config->GetDirectDiff()) {

SU2_CFD/src/solvers/CSpeciesSolver.cpp

@@ -385,8 +385,11 @@ void CSpeciesSolver::BC_Inlet(CGeometry* geometry, CSolver** solver_container, C
       if (dynamic_grid)
         conv_numerics->SetGridVel(geometry->nodes->GetGridVel(iPoint), geometry->nodes->GetGridVel(iPoint));
 
-      /*--- Compute the residual using an upwind scheme ---*/
+      /* nijso: TODO we have to compute mass flux here as well */ 
+      //su2double EdgeMassFlux = solver_container[FLOW_SOL]->GetEdgeMassFlux(iEdge);
+      //conv_numerics->SetMassFlux(EdgeMassFlux);
 
+      /*--- Compute the residual using an upwind scheme ---*/
       auto residual = conv_numerics->ComputeResidual(config);
       LinSysRes.AddBlock(iPoint, residual);
 
@@ -530,8 +533,11 @@ void CSpeciesSolver::BC_Outlet(CGeometry* geometry, CSolver** solver_container,
       if (dynamic_grid)
         conv_numerics->SetGridVel(geometry->nodes->GetGridVel(iPoint), geometry->nodes->GetGridVel(iPoint));
 
-      /*--- Compute the residual using an upwind scheme ---*/
+      /* nijso: TODO we have to compute mass flux here as well */ 
+      //su2double EdgeMassFlux = solver_container[FLOW_SOL]->GetEdgeMassFlux(iEdge);
+      //conv_numerics->SetMassFlux(EdgeMassFlux);
 
+      /*--- Compute the residual using an upwind scheme ---*/
       auto residual = conv_numerics->ComputeResidual(config);
       LinSysRes.AddBlock(iPoint, residual);
 

TestCases/species_transport/venturi_primitive_3species/DAspecies3_primitiveVenturi.cfg

@@ -53,7 +53,7 @@ INLET_FILENAME= inlet_venturi.dat
 INC_INLET_TYPE=  VELOCITY_INLET VELOCITY_INLET
 MARKER_INLET= ( gas_inlet,       300, 1.0, 1.0,  0.0, 0.0,\
                 air_axial_inlet, 300, 1.0, 0.0, -1.0, 0.0 )
-SPECIES_USE_STRONG_BC= NO
+SPECIES_USE_STRONG_BC= YES
 MARKER_INLET_SPECIES= (gas_inlet,       0.5, 0.5,\
                        air_axial_inlet, 0.6, 0.0 )
 %
@@ -128,7 +128,7 @@ CONV_FILENAME= history
 MARKER_ANALYZE= outlet
 MARKER_ANALYZE_AVERAGE= MASSFLUX
 %
-OUTPUT_FILES= RESTART_ASCII, PARAVIEW_MULTIBLOCK
+OUTPUT_FILES= RESTART_ASCII, PARAVIEW
 VOLUME_OUTPUT= RESIDUAL, PRIMITIVE
 OUTPUT_WRT_FREQ= 1000
 %

config_template.cfg

@@ -715,7 +715,7 @@ MARKER_INLET_SPECIES= (inlet, 0.5, ..., inlet2, 0.6, ...)
 % Use strong inlet and outlet BC in the species solver
 SPECIES_USE_STRONG_BC= NO
 %
-% Convective numerical method for species transport (SCALAR_UPWIND)
+% Convective numerical method for species transport (SCALAR_UPWIND, SCALAR_ADVECTION)
 CONV_NUM_METHOD_SPECIES= SCALAR_UPWIND
 %
 % Monotonic Upwind Scheme for Conservation Laws (TVD) in the species equations.
@@ -1272,7 +1272,7 @@ KIND_MATRIX_COLORING= GREEDY_COLORING
 
 % -------------------- TURBULENT NUMERICAL METHOD DEFINITION ------------------%
 %
-% Convective numerical method (SCALAR_UPWIND)
+% Convective numerical method (SCALAR_UPWIND, SCALAR_ADVECTION)
 CONV_NUM_METHOD_TURB= SCALAR_UPWIND
 %
 % Time discretization (EULER_IMPLICIT, EULER_EXPLICIT)