[WIP] Fix symmetry plane issues in SU2, change SU2-NEMO sym plane to match flow solver

SU2_CFD/include/solvers/CFVMFlowSolverBase.inl

@@ -1062,23 +1062,90 @@ void CFVMFlowSolverBase<V, R>::PushSolutionBackInTime(unsigned long TimeIter, bo
 template <class V, ENUM_REGIME R>
 void CFVMFlowSolverBase<V, R>::BC_Sym_Plane(CGeometry* geometry, CSolver** solver_container, CNumerics* conv_numerics,
                                             CNumerics* visc_numerics, CConfig* config, unsigned short val_marker) {
-  unsigned short iDim, iVar;
-  unsigned long iVertex, iPoint;
+  unsigned short iDim, iVar, iNeigh;
+  unsigned long iVertex, iPoint, jPoint, iEdge, iMarker;
+  su2double Tangent[MAXNDIM]  = {0.0};
+  su2double Normal_Sym[MAXNDIM] = {0.0}, UnitNormal_Sym[MAXNDIM] = {0.0};
+  su2double Normal_Product, Product, tol = 1e-16;
 
-  bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT);
-  bool viscous = config->GetViscous();
+  const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT);
+  const bool viscous = config->GetViscous();
+  const bool nemo = config->GetNEMOProblem();
   bool preprocessed = false;
+  const unsigned short nSpecies = config->GetnSpecies();
 
   /*--- Allocation of variables necessary for convective fluxes. ---*/
   su2double Area, ProjVelocity_i, *V_reflected, *V_domain, Normal[MAXNDIM] = {0.0}, UnitNormal[MAXNDIM] = {0.0};
 
   /*--- Allocation of variables necessary for viscous fluxes. ---*/
   su2double ProjGradient, ProjNormVelGrad, ProjTangVelGrad, TangentialNorm,
-      Tangential[MAXNDIM] = {0.0}, GradNormVel[MAXNDIM] = {0.0}, GradTangVel[MAXNDIM] = {0.0};
+      Tangential[MAXNDIM] = {0.0}, GradNormVel[MAXNDIM] = {0.0}, GradTangVel[MAXNDIM] = {0.0}, Residual[MAXNVAR] = {0.0};
 
   /*--- Allocation of primitive gradient arrays for viscous fluxes. ---*/
   su2activematrix Grad_Reflected(nPrimVarGrad, nDim);
 
+/*--- Count number of symmetry planes where each Vertex is inserted ---*/
+  for (iMarker = 0; iMarker < geometry->GetnMarker(); iMarker++) {
+    if (config->GetMarker_All_KindBC(iMarker) == SYMMETRY_PLANE){
+      for (iVertex = 0; iVertex < geometry->GetnVertex(iMarker); iVertex++) {
+        iPoint = geometry->vertex[iMarker][iVertex]->GetNode();
+        nodes->SetSymmetry(iPoint);
+      }
+    }
+  }
+
+  /*--- Correct normal directions of edges ---*/
+  for (iMarker = 0; iMarker < geometry->GetnMarker(); iMarker++) {
+    if (config->GetMarker_All_KindBC(iMarker) == SYMMETRY_PLANE){
+
+      for (iVertex = 0; iVertex < geometry->GetnVertex(iMarker); iVertex++) {
+
+        iPoint = geometry->vertex[iMarker][iVertex]->GetNode();
+
+        geometry->vertex[iMarker][iVertex]->GetNormal(Normal_Sym);
+
+        Area = GeometryToolbox::Norm(nDim, Normal_Sym);
+
+        for(iDim = 0; iDim<nDim; iDim++){
+          UnitNormal_Sym[iDim] = Normal_Sym[iDim]/Area;
+        }
+
+        for (iNeigh = 0; iNeigh < geometry->nodes->GetnPoint(iPoint); ++iNeigh){
+          Product = 0.0;
+
+          jPoint = geometry->nodes->GetPoint(iPoint,iNeigh);
+
+          /*---Check if neighbour point is on the same plane as the Symmetry_Plane
+               by computing the internal product and of the Normal Vertex vector and
+               the vector connecting iPoint and jPoint. If the product is lower than
+               estabilished tolerance (to account for Numerical errors) both points are
+               in the same plane as SYMMETRY_PLANE---*/
+
+          for(iDim = 0; iDim<nDim; iDim++){
+            Tangent[iDim] = geometry->nodes->GetCoord(jPoint,iDim) - geometry->nodes->GetCoord(iPoint,iDim);
+            Product += Tangent[iDim] * Normal_Sym[iDim];
+          }
+
+          if (abs(Product) < tol) {
+            Product = 0.0;
+
+            iEdge = geometry->nodes->GetEdge(iPoint,iNeigh);
+
+            geometry->edges->GetNormal(iEdge,Normal);
+
+            for(iDim = 0; iDim<nDim; iDim++)
+              Product += Normal[iDim]*UnitNormal_Sym[iDim];
+
+            for(iDim = 0; iDim<nDim; iDim++)
+              Normal[iDim]-=Product*UnitNormal_Sym[iDim];
+
+            geometry->edges->SetNormal(iEdge,Normal);
+          }
+        }
+      }
+    }
+  }
+
   /*--- Loop over all the vertices on this boundary marker. ---*/
 
   SU2_OMP_FOR_DYN(OMP_MIN_SIZE)
@@ -1143,9 +1210,10 @@ void CFVMFlowSolverBase<V, R>::BC_Sym_Plane(CGeometry* geometry, CSolver** solve
     }      // if bound_is_straight
 
     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)) {
+      Normal_Product = 0.0;
       /*-------------------------------------------------------------------------------*/
       /*--- Step 1: For the convective fluxes, create a reflected state of the      ---*/
       /*---         Primitive variables by copying all interior values to the       ---*/
@@ -1172,25 +1240,30 @@ void CFVMFlowSolverBase<V, R>::BC_Sym_Plane(CGeometry* geometry, CSolver** solve
             the velocity is mirrored along the symmetry boundary, i.e. the velocity in
             normal direction is substracted twice. ---*/
       for (iVar = 0; iVar < nPrimVar; iVar++) V_reflected[iVar] = nodes->GetPrimitive(iPoint, iVar);
-
+        
       /*--- Compute velocity in normal direction (ProjVelcity_i=(v*n)) und substract twice from
             velocity in normal direction: v_r = v - 2 (v*n)n ---*/
       ProjVelocity_i = nodes->GetProjVel(iPoint, UnitNormal);
-
+        
       /*--- Adjustment to v.n due to grid movement. ---*/
       if (dynamic_grid) {
         ProjVelocity_i -= GeometryToolbox::DotProduct(nDim, geometry->nodes->GetGridVel(iPoint), UnitNormal);
       }
 
       for (iDim = 0; iDim < nDim; iDim++)
-        V_reflected[iDim + 1] = nodes->GetVelocity(iPoint, iDim) - 2.0 * ProjVelocity_i * UnitNormal[iDim];
+        V_reflected[prim_idx.Velocity() + iDim] = nodes->GetVelocity(iPoint, iDim) - 2.0 * ProjVelocity_i * UnitNormal[iDim];
 
       /*--- Set Primitive and Secondary for numerics class. ---*/
       conv_numerics->SetPrimitive(V_domain, V_reflected);
       conv_numerics->SetSecondary(nodes->GetSecondary(iPoint), nodes->GetSecondary(iPoint));
 
-      /*--- Compute the residual using an upwind scheme. ---*/
+      if (nemo) {
+        conv_numerics->SetEve(nodes->GetEve(iPoint),    nodes->GetEve(iPoint));
+        conv_numerics->SetCvve(nodes->GetCvve(iPoint),   nodes->GetCvve(iPoint));
+        conv_numerics->SetGamma(nodes->GetGamma(iPoint),  nodes->GetGamma(iPoint));
+      }
 
+      /*--- Compute the residual using an upwind scheme. ---*/
       auto residual = conv_numerics->ComputeResidual(config);
 
       /*--- Update residual value ---*/

SU2_CFD/include/solvers/CNEMOEulerSolver.hpp

@@ -238,18 +238,6 @@ class CNEMOEulerSolver : public CFVMFlowSolverBase<CNEMOEulerVariable, ENUM_REGI
   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 symmetry boundary condition using the residual.
-   * \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 for convective terms.
-   * \param[in] visc_numerics - Description of the numerical method for viscous terms.
-   * \param[in] config - Definition of the particular problem.
-   * \param[in] val_marker - Surface marker where the boundary condition is applied.
-   */
-  void BC_Sym_Plane(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics,
-                    CNumerics *visc_numerics, CConfig *config, unsigned short val_marker) final;
-
   /*!
    * \brief Impose a subsonic inlet boundary condition.
    * \param[in] geometry - Geometrical definition of the problem.

SU2_CFD/include/variables/CEulerVariable.hpp

@@ -71,6 +71,8 @@ class CEulerVariable : public CFlowVariable {
   MatrixType WindGust;      /*! < \brief Wind gust value */
   MatrixType WindGustDer;   /*! < \brief Wind gust derivatives value */
 
+  VectorType symmetry;      /*!< \brief Nodes in symmetry planes. */
+
  public:
   /*!
    * \brief Constructor of the class.
@@ -321,4 +323,30 @@ class CEulerVariable : public CFlowVariable {
     for (unsigned long iDim = 0; iDim < nDim; iDim++) Solution(iPoint, iDim+1) = GetDensity(iPoint) * val_vector[iDim];
   }
 
+ /*!
+   * \brief Returns the stored value of Eve at the specified node
+   */
+  su2double *GetEve(unsigned long iPoint) {return nullptr;}
+
+  /*!
+   * \brief Returns the value of Cvve at the specified node
+   */
+  su2double *GetCvve(unsigned long iPoint) {return nullptr;}
+
+  /*!
+   * \brief Returns the stored value of Gamma at the specified node
+   */
+  su2double GetGamma(unsigned long iPoint) {return 0;}
+
+    /*!
+   * \brief Retrieves the number of symmetry planes at the specified node.
+   */  
+  inline su2double GetSymmetry(unsigned long iPoint) { return symmetry[iPoint]; }
+
+  /*!
+   * \brief Increases the number of symmetry planes at the specified node by one.
+   */  
+  inline void SetSymmetry(unsigned long iPoint) {symmetry[iPoint] += 1.0;}
+
+
 };

SU2_CFD/include/variables/CIncEulerVariable.hpp

@@ -68,6 +68,9 @@ class CIncEulerVariable : public CFlowVariable {
 
   VectorType Streamwise_Periodic_RecoveredPressure,    /*!< \brief Recovered/Physical pressure [Pa] for streamwise periodic flow. */
              Streamwise_Periodic_RecoveredTemperature; /*!< \brief Recovered/Physical temperature [K] for streamwise periodic flow. */
+
+  VectorType symmetry;      /*!< \brief Nodes in symmetry planes. */
+
  public:
   /*!
    * \brief Constructor of the class.
@@ -271,5 +274,29 @@ class CIncEulerVariable : public CFlowVariable {
   inline void SetVelSolutionVector(unsigned long iPoint, const su2double *val_vector) final {
     for (unsigned long iDim = 0; iDim < nDim; iDim++) Solution(iPoint, iDim+1) = val_vector[iDim];
   }
+ /*!
+   * \brief Returns the stored value of Eve at the specified node
+   */
+  su2double *GetEve(unsigned long iPoint) {return nullptr;}
+
+  /*!
+   * \brief Returns the value of Cvve at the specified node
+   */
+  su2double *GetCvve(unsigned long iPoint) {return nullptr;}
+
+  /*!
+   * \brief Returns the stored value of Gamma at the specified node
+   */
+  su2double GetGamma(unsigned long iPoint) {return 0;}
+
+    /*!
+   * \brief Retrieves the number of symmetry planes at the specified node.
+   */  
+  inline su2double GetSymmetry(unsigned long iPoint) { return symmetry[iPoint]; }
+
+  /*!
+   * \brief Increases the number of symmetry planes at the specified node by one.
+   */  
+  inline void SetSymmetry(unsigned long iPoint) {symmetry[iPoint] += 1.0;}
 
 };

SU2_CFD/include/variables/CNEMOEulerVariable.hpp

@@ -95,6 +95,8 @@ class CNEMOEulerVariable : public CFlowVariable {
   LAM_VISC_INDEX, EDDY_VISC_INDEX, nSpecies;
 
   su2double Tve_Freestream; /*!< \brief Freestream vib-el temperature. */
+  VectorType symmetry;
+
   const bool implicit;      /*!< \brief Implicit flag. */
 
  public:
@@ -420,4 +422,14 @@ class CNEMOEulerVariable : public CFlowVariable {
     for (unsigned long iDim = 0; iDim < nDim; iDim++) Res_TruncError(iPoint,nSpecies+iDim) = 0.0;
   }
 
+    /*!
+   * \brief Increases the number of symmetry planes at the specified node by one.
+   */  
+  inline void SetSymmetry(unsigned long iPoint) {symmetry[iPoint] += 1.0;}
+
+  /*!
+   * \brief Retrieves the number of symmetry planes at the specified node.
+   */  
+  inline su2double GetSymmetry(unsigned long iPoint) { return symmetry[iPoint]; }
+
 };

SU2_CFD/src/solvers/CNEMOEulerSolver.cpp

@@ -1408,107 +1408,6 @@ void CNEMOEulerSolver::SetReferenceValues(const CConfig& config) {
 
 }
 
-void CNEMOEulerSolver::BC_Sym_Plane(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics,
-                                    CNumerics *visc_numerics, CConfig *config, unsigned short val_marker) {
-
-  unsigned short iDim, jDim, iSpecies, iVar, jVar;
-  unsigned long iPoint, iVertex;
-
-  /*--- Allocate the necessary vector structures ---*/
-  su2double Normal[MAXNDIM] = {0.0}, UnitNormal[MAXNDIM] = {0.0};
-
-  bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT);
-
-  /*--- Loop over all the vertices on this boundary (val_marker) ---*/
-  for (iVertex = 0; iVertex < geometry->nVertex[val_marker]; iVertex++) {
-    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)) {
-
-      /*--- Normal vector for this vertex (negative for outward convention) ---*/
-      geometry->vertex[val_marker][iVertex]->GetNormal(Normal);
-
-      /*--- Calculate parameters from the geometry ---*/
-      su2double Area = GeometryToolbox::Norm(nDim, Normal);
-
-      for (iDim = 0; iDim < nDim; iDim++){
-        UnitNormal[iDim] = -Normal[iDim]/Area;
-      }
-
-      /*--- Retrieve the pressure on the vertex ---*/
-      su2double P = nodes->GetPressure(iPoint);
-
-      /*--- Apply the flow-tangency b.c. to the convective flux ---*/
-      for (iSpecies = 0; iSpecies < nSpecies; iSpecies++)
-        Residual[iSpecies] = 0.0;
-      for (iDim = 0; iDim < nDim; iDim++){
-        Residual[nSpecies+iDim] = P * UnitNormal[iDim] * Area;
-      }
-      Residual[nSpecies+nDim]   = 0.0;
-      Residual[nSpecies+nDim+1] = 0.0;
-
-      /*--- Add value to the residual ---*/
-      LinSysRes.AddBlock(iPoint, Residual);
-
-      /*--- If using implicit time-stepping, calculate b.c. contribution to Jacobian ---*/
-      if (implicit) {
-
-        /*--- Allocate arrays needed for implicit solver ---*/
-        su2double Velocity[MAXNDIM] = {0.0};
-
-        /*--- Get species molar mass ---*/
-        auto& Ms = FluidModel->GetSpeciesMolarMass();
-
-        /*--- Initialize Jacobian ---*/
-        for (iVar = 0; iVar < nVar; iVar++)
-          for (jVar = 0; jVar < nVar; jVar++)
-            Jacobian_i[iVar][jVar] = 0.0;
-
-        /*--- Calculate state i ---*/
-        su2double rho     = nodes->GetDensity(iPoint);
-        su2double rhoE    = nodes->GetSolution(iPoint,nSpecies+nDim);
-        su2double rhoEve  = nodes->GetSolution(iPoint,nSpecies+nDim+1);
-        auto dPdU    = nodes->GetdPdU(iPoint);
-        for (iDim = 0; iDim < nDim; iDim++)
-          Velocity[iDim] = nodes->GetVelocity(iPoint,iDim);
-
-        su2double conc = 0.0;
-        for (iSpecies = 0; iSpecies < nSpecies; iSpecies++) {
-          su2double cs = nodes->GetMassFraction(iPoint,iSpecies);
-          conc += cs * rho/Ms[iSpecies];
-
-          for (iDim = 0; iDim < nDim; iDim++) {
-            Jacobian_i[nSpecies+iDim][iSpecies] = dPdU[iSpecies] * UnitNormal[iDim];
-            Jacobian_i[iSpecies][nSpecies+iDim] = cs * UnitNormal[iDim];
-          }
-        }
-
-        for (iDim = 0; iDim < nDim; iDim++) {
-          for (jDim = 0; jDim < nDim; jDim++) {
-            Jacobian_i[nSpecies+iDim][nSpecies+jDim] = Velocity[iDim]*UnitNormal[jDim]
-                + dPdU[nSpecies+jDim]*UnitNormal[iDim];
-          }
-          Jacobian_i[nSpecies+iDim][nSpecies+nDim]   = dPdU[nSpecies+nDim]  *UnitNormal[iDim];
-          Jacobian_i[nSpecies+iDim][nSpecies+nDim+1] = dPdU[nSpecies+nDim+1]*UnitNormal[iDim];
-
-          Jacobian_i[nSpecies+nDim][nSpecies+iDim]   = (rhoE+P)/rho * UnitNormal[iDim];
-          Jacobian_i[nSpecies+nDim+1][nSpecies+iDim] = rhoEve/rho   * UnitNormal[iDim];
-        }
-
-        /*--- Integrate over the dual-grid area ---*/
-        for (iVar = 0; iVar < nVar; iVar++)
-          for (jVar = 0; jVar < nVar; jVar++)
-            Jacobian_i[iVar][jVar] = Jacobian_i[iVar][jVar] * Area;
-
-        /*--- Apply the contribution to the system ---*/
-        Jacobian.AddBlock2Diag(iPoint, Jacobian_i);
-
-      }
-    }
-  }
-}
-
 void CNEMOEulerSolver::BC_Far_Field(CGeometry *geometry, CSolver **solver_container,
                                     CNumerics *conv_numerics, CNumerics *visc_numerics,
                                     CConfig *config, unsigned short val_marker) {

SU2_CFD/src/variables/CNEMOEulerVariable.cpp

@@ -90,6 +90,9 @@ CNEMOEulerVariable::CNEMOEulerVariable(su2double val_pressure,
   eves.resize(nPoint, nSpecies)  = su2double(0.0);
   Gamma.resize(nPoint)           = su2double(0.0);
 
+  /* Vector to count number of symmetry planes at each node. */
+  symmetry.resize(nPoint) = su2double(0.0);
+
   /*--- Set mixture state ---*/
   fluidmodel->SetTDStatePTTv(val_pressure, val_massfrac, val_temperature, val_temperature_ve);
 

SU2_IDE/Xcode/SU2_CFD.xcodeproj/project.pbxproj

@@ -2790,4 +2790,4 @@
 /* End XCConfigurationList section */
 	};
 	rootObject = 08FB7793FE84155DC02AAC07 /* Project object */;
-}
+}