Multizone with species transport (#1821)

Co-authored-by: Pedro Gomes <38071223+pcarruscag@users.noreply.github.com>
* move slidingstates to scalarsolver
* regression testcase for multizone
* added bc_fluid_interface to speciessolver
Co-authored-by: EvertBunschoten <evertbun@xs4all.nl>

Author: 3 people

Common/src/CConfig.cpp

@@ -5354,7 +5354,8 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i
         Kind_Solver != MAIN_SOLVER::NAVIER_STOKES &&
         Kind_Solver != MAIN_SOLVER::RANS &&
         Kind_Solver != MAIN_SOLVER::DISC_ADJ_NAVIER_STOKES &&
-        Kind_Solver != MAIN_SOLVER::DISC_ADJ_RANS)
+        Kind_Solver != MAIN_SOLVER::DISC_ADJ_RANS &&
+        Kind_Solver != MAIN_SOLVER::MULTIPHYSICS)
       SU2_MPI::Error("Species transport currently only available for compressible and incompressible flow.", CURRENT_FUNCTION);
 
     /*--- Species specific OF currently can only handle one entry in Marker_Analyze. ---*/

SU2_CFD/include/solvers/CScalarSolver.hpp

@@ -29,6 +29,7 @@
 #include <vector>
 
 #include "../../../Common/include/parallelization/omp_structure.hpp"
+#include "../../../Common/include/toolboxes/geometry_toolbox.hpp"
 #include "../variables/CScalarVariable.hpp"
 #include "CSolver.hpp"
 
@@ -57,6 +58,9 @@ class CScalarSolver : public CSolver {
 
   const bool Conservative; /*!< \brief Transported Variable is conservative. Solution has to be multiplied with rho. */
 
+  vector<su2matrix<su2double*> > SlidingState; // vector of matrix of pointers... inner dim alloc'd elsewhere (welcome, to the twilight zone)
+  vector<vector<int> > SlidingStateNodes;
+
   /*--- Shallow copy of grid coloring for OpenMP parallelization. ---*/
 
 #ifdef HAVE_OMP
@@ -135,6 +139,122 @@ class CScalarSolver : public CSolver {
     }
   }
 
+  /*!
+   * \brief Generic implementation of the fluid interface boundary condition for scalar solvers.
+   * \tparam SolverSpecificNumericsFunc - lambda that implements solver specific contributions to viscous numerics.
+   * \note The functor has to implement (iPoint)
+   * \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.
+   */
+  template <class SolverSpecificNumericsFunc>
+  void BC_Fluid_Interface_impl(const SolverSpecificNumericsFunc& SolverSpecificNumerics, CGeometry *geometry,
+                               CSolver **solver_container, CNumerics *conv_numerics, CNumerics *visc_numerics,
+                               CConfig *config) {
+    const auto nPrimVar = solver_container[FLOW_SOL]->GetnPrimVar();
+    su2activevector PrimVar_j(nPrimVar);
+    su2double solution_j[MAXNVAR] = {0.0};
+
+    for (auto iMarker = 0u; iMarker < config->GetnMarker_All(); iMarker++) {
+
+      if (config->GetMarker_All_KindBC(iMarker) != FLUID_INTERFACE) continue;
+
+      SU2_OMP_FOR_STAT(OMP_MIN_SIZE)
+      for (auto iVertex = 0u; iVertex < geometry->nVertex[iMarker]; iVertex++) {
+
+        const auto iPoint = geometry->vertex[iMarker][iVertex]->GetNode();
+
+        if (!geometry->nodes->GetDomain(iPoint)) continue;
+
+        const auto Point_Normal = geometry->vertex[iMarker][iVertex]->GetNormal_Neighbor();
+        const auto nDonorVertex = GetnSlidingStates(iMarker,iVertex);
+
+        su2double Normal[MAXNDIM] = {0.0};
+        for (auto iDim = 0u; iDim < nDim; iDim++)
+          Normal[iDim] = -geometry->vertex[iMarker][iVertex]->GetNormal()[iDim];
+
+        su2double* PrimVar_i = solver_container[FLOW_SOL]->GetNodes()->GetPrimitive(iPoint);
+
+        auto Jacobian_i = Jacobian.GetBlock(iPoint,iPoint);
+
+        /*--- Loop over the nDonorVertexes and compute the averaged flux ---*/
+
+        for (auto jVertex = 0; jVertex < nDonorVertex; jVertex++) {
+
+          for (auto iVar = 0u; iVar < nPrimVar; iVar++)
+            PrimVar_j[iVar] = solver_container[FLOW_SOL]->GetSlidingState(iMarker, iVertex, iVar, jVertex);
+
+          /*--- Get the weight computed in the interpolator class for the j-th donor vertex ---*/
+
+          const su2double weight = solver_container[FLOW_SOL]->GetSlidingState(iMarker, iVertex, nPrimVar, jVertex);
+
+          /*--- Set primitive variables ---*/
+
+          conv_numerics->SetPrimitive( PrimVar_i, PrimVar_j.data() );
+
+          /*--- Set the scalar variable states ---*/
+
+          for (auto iVar = 0u; iVar < nVar; ++iVar)
+            solution_j[iVar] = GetSlidingState(iMarker, iVertex, iVar, jVertex);
+
+          conv_numerics->SetScalarVar(nodes->GetSolution(iPoint), solution_j);
+
+          /*--- Set the normal vector ---*/
+
+          conv_numerics->SetNormal(Normal);
+
+          if (dynamic_grid)
+            conv_numerics->SetGridVel(geometry->nodes->GetGridVel(iPoint), geometry->nodes->GetGridVel(iPoint));
+
+          auto residual = conv_numerics->ComputeResidual(config);
+
+          /*--- Accumulate the residuals to compute the average ---*/
+
+          for (auto iVar = 0u; iVar < nVar; iVar++) {
+            LinSysRes(iPoint,iVar) += weight*residual[iVar];
+            for (auto jVar = 0u; jVar < nVar; jVar++)
+              Jacobian_i[iVar*nVar+jVar] += SU2_TYPE::GetValue(weight*residual.jacobian_i[iVar][jVar]);
+          }
+        }
+
+        /*--- Set the normal vector and the coordinates ---*/
+
+        visc_numerics->SetNormal(Normal);
+        su2double Coord_Reflected[MAXNDIM];
+        GeometryToolbox::PointPointReflect(nDim, geometry->nodes->GetCoord(Point_Normal),
+                                                geometry->nodes->GetCoord(iPoint), Coord_Reflected);
+        visc_numerics->SetCoord(geometry->nodes->GetCoord(iPoint), Coord_Reflected);
+
+        /*--- Primitive variables ---*/
+
+        visc_numerics->SetPrimitive(PrimVar_i, PrimVar_j.data());
+
+        /*--- Scalar variables and their gradients ---*/
+
+        visc_numerics->SetScalarVar(nodes->GetSolution(iPoint), solution_j);
+        visc_numerics->SetScalarVarGradient(nodes->GetGradient(iPoint), nodes->GetGradient(iPoint));
+
+        /*--- Allow derived solvers to set more variables in numerics. ---*/
+
+        SolverSpecificNumerics(iPoint);
+
+        /*--- Compute and update residual ---*/
+
+        auto residual = visc_numerics->ComputeResidual(config);
+
+        LinSysRes.SubtractBlock(iPoint, residual);
+
+        /*--- Jacobian contribution for implicit integration ---*/
+
+        Jacobian.SubtractBlock2Diag(iPoint, residual.jacobian_i);
+
+      }
+      END_SU2_OMP_FOR
+    }
+  }
+
   /*!
    * \brief Gradient and Limiter computation.
    * \param[in] geometry - Geometrical definition of the problem.
@@ -255,6 +375,20 @@ class CScalarSolver : public CSolver {
    */
   void BC_Periodic(CGeometry* geometry, CSolver** solver_container, CNumerics* numerics, CConfig* config) final;
 
+  /*!
+   * \brief Impose the fluid interface boundary condition using transfer data.
+   * \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.
+   */
+  void BC_Fluid_Interface(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics,
+                          CNumerics *visc_numerics, CConfig *config) override {
+    /*--- By default instantiate the generic implementation w/o extra variables, derived solvers can override. ---*/
+    BC_Fluid_Interface_impl([](unsigned long){}, geometry, solver_container, conv_numerics, visc_numerics, config);
+  }
+
   /*!
    * \brief Set the solution using the Freestream values.
    * \param[in] config - Definition of the particular problem.
@@ -326,4 +460,71 @@ class CScalarSolver : public CSolver {
    * \brief Scalar solvers support OpenMP+MPI.
    */
   inline bool GetHasHybridParallel() const override { return true; }
+
+  /*!
+  * \brief Get the outer state for fluid interface nodes.
+  * \param[in] val_marker - marker index
+  * \param[in] val_vertex - vertex index
+  * \param[in] val_state  - requested state component
+  * \param[in] donor_index- index of the donor node to get
+  */
+  inline su2double GetSlidingState(unsigned short val_marker,
+                                   unsigned long val_vertex,
+                                   unsigned short val_state,
+                                   unsigned long donor_index) const final {
+    return SlidingState[val_marker][val_vertex][val_state][donor_index];
+  }
+
+  /*!
+   * \brief Allocates the final pointer of SlidingState depending on how many donor vertex donate to it. That number is stored in SlidingStateNodes[val_marker][val_vertex].
+   * \param[in] val_marker   - marker index
+   * \param[in] val_vertex   - vertex index
+   */
+  inline void SetSlidingStateStructure(unsigned short val_marker, unsigned long val_vertex) final {
+    int iVar;
+
+    for( iVar = 0; iVar < nVar+1; iVar++){
+      if( SlidingState[val_marker][val_vertex][iVar] != nullptr )
+        delete [] SlidingState[val_marker][val_vertex][iVar];
+    }
+
+    for( iVar = 0; iVar < nVar+1; iVar++)
+      SlidingState[val_marker][val_vertex][iVar] = new su2double[ GetnSlidingStates(val_marker, val_vertex) ];
+  }
+
+  /*!
+   * \brief Set the outer state for fluid interface nodes.
+   * \param[in] val_marker   - marker index
+   * \param[in] val_vertex   - vertex index
+   * \param[in] val_state    - requested state component
+   * \param[in] donor_index  - index of the donor node to set
+   * \param[in] component    - set value
+   */
+  inline void SetSlidingState(unsigned short val_marker,
+                              unsigned long val_vertex,
+                              unsigned short val_state,
+                              unsigned long donor_index,
+                              su2double component) final {
+    SlidingState[val_marker][val_vertex][val_state][donor_index] = component;
+  }
+
+  /*!
+   * \brief Set the number of outer state for fluid interface nodes.
+   * \param[in] val_marker - marker index
+   * \param[in] val_vertex - vertex index
+   * \param[in] value - number of outer states
+   */
+  inline void SetnSlidingStates(unsigned short val_marker,
+                                unsigned long val_vertex,
+                                int value) final { SlidingStateNodes[val_marker][val_vertex] = value; }
+
+  /*!
+   * \brief Get the number of outer state for fluid interface nodes.
+   * \param[in] val_marker - marker index
+   * \param[in] val_vertex - vertex index
+   */
+  inline int GetnSlidingStates(unsigned short val_marker, unsigned long val_vertex) const final {
+    return SlidingStateNodes[val_marker][val_vertex];
+  }
+
 };

SU2_CFD/include/solvers/CSpeciesSolver.hpp

@@ -161,4 +161,25 @@ class CSpeciesSolver : public CScalarSolver<CSpeciesVariable> {
    */
   void Source_Residual(CGeometry* geometry, CSolver** solver_container, CNumerics** numerics_container, CConfig* config,
                        unsigned short iMesh) override;
+
+
+/*!
+   * \brief Impose the fluid interface boundary condition using tranfer data.
+   * \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.
+   */
+  void BC_Fluid_Interface(CGeometry *geometry,
+                          CSolver **solver_container,
+                          CNumerics *conv_numerics,
+                          CNumerics *visc_numerics,
+                          CConfig *config) final {
+    BC_Fluid_Interface_impl(
+      [&](unsigned long iPoint) {
+        visc_numerics->SetDiffusionCoeff(nodes->GetDiffusivity(iPoint), nodes->GetDiffusivity(iPoint));
+      },
+      geometry, solver_container, conv_numerics, visc_numerics, config);
+  }
 };

SU2_CFD/include/solvers/CTurbSSTSolver.hpp

@@ -248,6 +248,26 @@ class CTurbSSTSolver final : public CTurbSolver {
                  CConfig *config,
                  unsigned short val_marker) override;
 
+  /*!
+   * \brief Impose the fluid interface boundary condition using tranfer data.
+   * \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.
+   */
+  void BC_Fluid_Interface(CGeometry *geometry,
+                          CSolver **solver_container,
+                          CNumerics *conv_numerics,
+                          CNumerics *visc_numerics,
+                          CConfig *config) final {
+    BC_Fluid_Interface_impl(
+      [&](unsigned long iPoint) {
+        visc_numerics->SetF1blending(nodes->GetF1blending(iPoint), nodes->GetF1blending(iPoint));
+      },
+      geometry, solver_container, conv_numerics, visc_numerics, config);
+  }
+
   /*!
    * \brief Get the constants for the SST model.
    * \return A pointer to an array containing a set of constants