Source term for axisymmetry in species transport (#1785)

* Source term for axisymmetry in species transport added (need to add test case)

Co-authored-by: Pedro Gomes <38071223+pcarruscag@users.noreply.github.com>
Co-authored-by: Nijso <bigfootedrockmidget@hotmail.com>

Author: 3 people

Common/src/CConfig.cpp

@@ -5354,9 +5354,6 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i
       }
     }
 
-    // For now, do not allow axisymmetric simulations
-    if (Axisymmetric) SU2_MPI::Error("Species transport currently not possible with axissymmetric flow.", CURRENT_FUNCTION);
-
     if(Kind_TimeIntScheme_Species != EULER_IMPLICIT &&
        Kind_TimeIntScheme_Species != EULER_EXPLICIT){
       SU2_MPI::Error("Only TIME_DISCRE_TURB = EULER_IMPLICIT, EULER_EXPLICIT have been implemented in the scalar solver.", CURRENT_FUNCTION);

SU2_CFD/include/numerics/species/species_sources.hpp

@@ -68,7 +68,11 @@ template <class FlowIndices>
 class CSourceAxisymmetric_Species : public CSourceBase_Species {
  protected:
   const FlowIndices idx;  /*!< \brief Object to manage the access to the flow primitives. */
-  bool implicit;
+  const bool implicit;
+  const bool viscous;
+  const bool turbulence;
+  const bool incompressible;
+  const su2double Sc_t;
 
  public:
   /*!

SU2_CFD/include/solvers/CSpeciesSolver.hpp

@@ -38,7 +38,7 @@
  */
 class CSpeciesSolver : public CScalarSolver<CSpeciesVariable> {
  protected:
-  unsigned short Inlet_Position;          /*!< \brief Column index for scalar variables in inlet files. */
+  unsigned short Inlet_Position;             /*!< \brief Column index for scalar variables in inlet files. */
   vector<su2activematrix> Inlet_SpeciesVars; /*!< \brief Species variables at inlet profiles. */
 
  public:
@@ -145,9 +145,20 @@ class CSpeciesSolver : public CScalarSolver<CSpeciesVariable> {
   void BC_Outlet(CGeometry* geometry, CSolver** solver_container, CNumerics* conv_numerics, CNumerics* visc_numerics,
                  CConfig* config, unsigned short val_marker) final;
 
-   /*--- Note that BC_Sym_Plane, BC_HeatFlux_Wall, BC_Isothermal_Wall are all treated as zero-flux BC for the
-    * mass-factions, which can be fulfilled by no additional residual contribution on these nodes.
-    * If a specified mass-fractions flux (like BC_HeatFlux_Wall) or a constant mass-fraction on the boundary
-    * (like BC_Isothermal_Wall) are implemented the respective CHeatSolver implementations can act as a  starting
-    * point ---*/
+  /*--- Note that BC_Sym_Plane, BC_HeatFlux_Wall, BC_Isothermal_Wall are all treated as zero-flux BC for the
+   * mass-factions, which can be fulfilled by no additional residual contribution on these nodes.
+   * If a specified mass-fractions flux (like BC_HeatFlux_Wall) or a constant mass-fraction on the boundary
+   * (like BC_Isothermal_Wall) are implemented the respective CHeatSolver implementations can act as a  starting
+   * point ---*/
+
+  /*!
+   * \brief Source term computation for axisymmetric flow.
+   * \param[in] geometry - Geometrical definition of the problem.
+   * \param[in] solver_container - Container vector with all the solutions.
+   * \param[in] numerics_container - Container for description of the numerical method.
+   * \param[in] config - Definition of the particular problem.
+   * \param[in] iMesh - Index of the mesh in multigrid computations.
+   */
+  void Source_Residual(CGeometry* geometry, CSolver** solver_container, CNumerics** numerics_container, CConfig* config,
+                       unsigned short iMesh) override;
 };

SU2_CFD/src/numerics/species/species_sources.cpp

@@ -57,19 +57,79 @@ template <class T>
 CSourceAxisymmetric_Species<T>::CSourceAxisymmetric_Species(unsigned short val_nDim, unsigned short val_nVar,
                                                          const CConfig* config)
     : CSourceBase_Species(val_nDim, val_nVar, config),
-    idx(val_nDim, config->GetnSpecies()) {
-  implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT);
+      idx(val_nDim, config->GetnSpecies()),
+    implicit(config->GetKind_TimeIntScheme_Species() == EULER_IMPLICIT),
+    viscous(config->GetViscous()),
+    turbulence(config->GetKind_Turb_Model() != TURB_MODEL::NONE),
+    incompressible(config->GetKind_Regime() == ENUM_REGIME::INCOMPRESSIBLE),
+    Sc_t(config->GetSchmidt_Number_Turbulent()) {
 }
 
 template <class T>
 CNumerics::ResidualType<> CSourceAxisymmetric_Species<T>::ComputeResidual(const CConfig* config) {
-  for (unsigned short iVar = 0; iVar < nVar; iVar++) residual[iVar] = 0.0;
+  
 
-  if (implicit) {
-    for (unsigned short iVar = 0; iVar < nVar; iVar++) {
-      for (unsigned short jVar = 0; jVar < nVar; jVar++) jacobian[iVar][jVar] = 0.0;
+  /*--- Preaccumulation ---*/
+  AD::StartPreacc();
+  AD::SetPreaccIn(ScalarVar_i, nVar);
+  AD::SetPreaccIn(Volume); 
+  
+  if (incompressible) {
+    AD::SetPreaccIn(V_i, nDim+6);
+  }
+  else {
+    AD::SetPreaccIn(V_i, nDim+7);
+  }
+
+  /*--- Initialization. ---*/
+  for (auto iVar = 0u; iVar < nVar; iVar++) {
+    residual[iVar] = 0.0;
+    for (auto jVar = 0; jVar < nVar; jVar++) {
+      jacobian[iVar][jVar] = 0.0;
+    }
+  }
+
+  /*--- Contribution due to 2D axisymmetric formulation ---*/
+  if (Coord_i[1] > EPS) {
+
+    AD::SetPreaccIn(Coord_i[1]);
+    AD::SetPreaccIn(Diffusion_Coeff_i, nVar);
+    AD::SetPreaccIn(ScalarVar_Grad_i, nVar, nDim);      
+
+    const su2double yinv = 1.0 / Coord_i[1];
+
+    const su2double Density_i = V_i[idx.Density()];
+    
+    su2double Velocity_i[3];
+    for (auto iDim = 0u; iDim < nDim; iDim++)
+      Velocity_i[iDim] = V_i[idx.Velocity() + iDim];
+
+    /*--- Inviscid component of the source term. ---*/
+    
+    for (auto iVar = 0u; iVar < nVar; iVar++)
+      residual[iVar] -= yinv * Volume * Density_i * ScalarVar_i[iVar] * Velocity_i[1];
+
+    if (implicit) {
+      for (auto iVar = 0u; iVar < nVar; iVar++) {
+        jacobian[iVar][iVar] -= yinv * Volume * Velocity_i[1];
+      }
+    }
+
+    /*--- Add the viscous terms if necessary. ---*/
+    
+    if (config->GetViscous()) {
+      su2double Mass_Diffusivity_Tur = 0.0;
+      if (turbulence)
+        Mass_Diffusivity_Tur = V_i[idx.EddyViscosity()] / Sc_t;
+
+      for (auto iVar=0u; iVar < nVar; iVar++){
+        residual[iVar] += yinv * Volume * (Density_i * Diffusion_Coeff_i[iVar] + Mass_Diffusivity_Tur) * ScalarVar_Grad_i[iVar][1];
+      } 
     }
   }
+  
+  AD::SetPreaccOut(residual, nVar);
+  AD::EndPreacc();
 
   return ResidualType<>(residual, jacobian, nullptr);
 }

SU2_CFD/src/solvers/CSpeciesSolver.cpp

@@ -542,3 +542,54 @@ void CSpeciesSolver::BC_Outlet(CGeometry* geometry, CSolver** solver_container,
   }
   END_SU2_OMP_FOR
 }
+
+void CSpeciesSolver::Source_Residual(CGeometry *geometry, CSolver **solver_container, CNumerics **numerics_container,
+                                      CConfig *config, unsigned short iMesh) {
+
+  const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT);
+  const bool axisymmetric = config->GetAxisymmetric();
+
+  if (axisymmetric) {
+    CNumerics *numerics  = numerics_container[SOURCE_FIRST_TERM  + omp_get_thread_num()*MAX_TERMS];
+  
+    SU2_OMP_FOR_DYN(omp_chunk_size)
+    for (auto iPoint = 0u; iPoint < nPointDomain; iPoint++) {
+      /*--- Set primitive variables w/o reconstruction ---*/
+
+      numerics->SetPrimitive(solver_container[FLOW_SOL]->GetNodes()->GetPrimitive(iPoint), nullptr);
+
+      /*--- Set scalar variables w/o reconstruction ---*/
+
+      numerics->SetScalarVar(nodes->GetSolution(iPoint), nullptr);
+
+      numerics->SetDiffusionCoeff(nodes->GetDiffusivity(iPoint), 0);
+
+      /*--- Set volume of the dual cell. ---*/
+
+      numerics->SetVolume(geometry->nodes->GetVolume(iPoint));
+
+      /*--- Update scalar sources in the fluidmodel ---*/
+
+      /*--- Axisymmetry source term for the scalar equation. ---*/
+      /*--- Set y coordinate ---*/
+      
+      numerics->SetCoord(geometry->nodes->GetCoord(iPoint), nullptr);
+      
+      /*--- Set gradients ---*/
+      
+      numerics->SetScalarVarGradient(nodes->GetGradient(iPoint), nullptr);
+
+      auto residual = numerics->ComputeResidual(config);
+
+      /*--- Add Residual ---*/
+    
+      LinSysRes.SubtractBlock(iPoint, residual);
+    
+      /*--- Implicit part ---*/
+    
+      if (implicit) Jacobian.SubtractBlock2Diag(iPoint, residual.jacobian_i);
+    
+    }
+    END_SU2_OMP_FOR
+  }
+}