Merge pull request #2392 from su2code/fix_cornernode_normal

Fix cornernode normal

Author: pcarruscag

Common/include/geometry/CMultiGridGeometry.hpp

@@ -123,7 +123,7 @@ class CMultiGridGeometry final : public CGeometry {
   void SetRestricted_GridVelocity(const CGeometry* fine_grid) override;
 
   /*!
-   * \brief Find and store the closest neighbor to a vertex.
+   * \brief Find and store the closest, most normal, neighbor to a vertex.
    * \param[in] config - Definition of the particular problem.
    */
   void FindNormal_Neighbor(const CConfig* config) override;

Common/include/geometry/CPhysicalGeometry.hpp

@@ -607,7 +607,7 @@ class CPhysicalGeometry final : public CGeometry {
   void ComputeMeshQualityStatistics(const CConfig* config) override;
 
   /*!
-   * \brief Find and store the closest neighbor to a vertex.
+   * \brief Find and store the closest, most normal, neighbor to a vertex.
    * \param[in] config - Definition of the particular problem.
    */
   void FindNormal_Neighbor(const CConfig* config) override;

Common/include/geometry/dual_grid/CVertex.hpp

@@ -45,7 +45,8 @@ class CVertex : public CDualGrid {
   long PeriodicPoint[5] = {-1};  /*!< \brief Store the periodic point of a boundary (iProcessor, iPoint) */
   bool ActDisk_Perimeter = false;      /*!< \brief Identify nodes at the perimeter of the actuator disk */
   short Rotation_Type;                 /*!< \brief Type of rotation associated with the vertex (MPI and periodic) */
-  unsigned long Normal_Neighbor;       /*!< \brief Index of the closest neighbor. */
+  unsigned long Normal_Neighbor;       /*!< \brief Index of the closest, most normal, neighbor. */
+  su2double NearestNeighborDist = 0.0; /*!< \brief Distance to nearest internal neighbor that is not a wall. */
   su2double Basis_Function[3] = {0.0}; /*!< \brief Basis function values for interpolation across zones. */
 
  public:
@@ -98,7 +99,7 @@ class CVertex : public CDualGrid {
    * \brief Get the normal to a face of the control volume asociated with a vertex.
    * \return Dimensional normal vector, the modulus is the area of the face.
    */
-  inline su2double* GetNormal(void) override { return Normal; }
+  inline su2double* GetNormal() override { return Normal; }
 
   /*!
    * \brief Get the ith component of the normal.
@@ -108,7 +109,7 @@ class CVertex : public CDualGrid {
   /*!
    * \brief Initialize normal vector.
    */
-  inline void SetZeroValues(void) override {
+  inline void SetZeroValues() override {
     for (unsigned short iDim = 0; iDim < nDim; iDim++) Normal[iDim] = 0.0;
   }
 
@@ -122,7 +123,7 @@ class CVertex : public CDualGrid {
    * \brief Get the value of an auxiliary variable for gradient computation.
    * \return Value of the auxiliar variable.
    */
-  inline su2double GetAuxVar(void) const { return Aux_Var; }
+  inline su2double GetAuxVar() const { return Aux_Var; }
 
   /*!
    * \brief Add the value of an auxiliary variable for gradient computation.
@@ -167,7 +168,7 @@ class CVertex : public CDualGrid {
    * \brief Get the value of the coordinate variation due to a surface modification.
    * \return Variation of the coordinate.
    */
-  inline su2double* GetVarCoord(void) { return VarCoord; }
+  inline su2double* GetVarCoord() { return VarCoord; }
 
   /*!
    * \brief Set the value of the cartesian coordinate for the vertex.
@@ -181,7 +182,7 @@ class CVertex : public CDualGrid {
    * \brief Get the value of the cartesian coordinate for the vertex.
    * \return Value of the cartesian coordinate of the vertex.
    */
-  inline su2double* GetCoord(void) override { return CartCoord; }
+  inline su2double* GetCoord() override { return CartCoord; }
 
   /*!
    * \brief Get the value of the cartesian coordinate for the vertex.
@@ -200,7 +201,7 @@ class CVertex : public CDualGrid {
    * \brief Get the type of rotation associated to the vertex.
    * \return Value of the rotation that must be applied to the solution of the vertex
    */
-  inline short GetRotation_Type(void) const { return Rotation_Type; }
+  inline short GetRotation_Type() const { return Rotation_Type; }
 
   /*!
    * \overload
@@ -250,49 +251,49 @@ class CVertex : public CDualGrid {
    * \brief Get the value of the periodic point of a vertex.
    * \return Value of the periodic point of a vertex.
    */
-  inline long GetDonorPoint(void) const { return PeriodicPoint[0]; }
+  inline long GetDonorPoint() const { return PeriodicPoint[0]; }
 
   /*!
    * \brief Get the value of the periodic point of a vertex.
    * \return Value of the periodic point of a vertex.
    */
-  inline long GetDonorMarker(void) const { return PeriodicPoint[4]; }
+  inline long GetDonorMarker() const { return PeriodicPoint[4]; }
 
   /*!
    * \brief Get the value of the periodic point of a vertex.
    * \return Value of the periodic point of a vertex.
    */
-  inline long GetDonorVertex(void) const { return PeriodicPoint[3]; }
+  inline long GetDonorVertex() const { return PeriodicPoint[3]; }
 
   /*!
    * \brief Get the value of the periodic point of a vertex.
    * \return Value of the periodic point of a vertex.
    */
-  inline long GetDonorGlobalIndex(void) const { return PeriodicPoint[2]; }
+  inline long GetDonorGlobalIndex() const { return PeriodicPoint[2]; }
 
   /*!
    * \brief Get the value of the periodic point of a vertex.
    * \return Value of the periodic point of a vertex.
    */
-  inline long GetGlobalDonorPoint(void) const { return PeriodicPoint[2]; }
+  inline long GetGlobalDonorPoint() const { return PeriodicPoint[2]; }
 
   /*!
    * \brief Get the value of the periodic point of a vertex.
    * \return Value of the periodic point of a vertex.
    */
-  inline long GetDonorProcessor(void) const { return PeriodicPoint[1]; }
+  inline long GetDonorProcessor() const { return PeriodicPoint[1]; }
 
   /*!
    * \brief Get the value of the periodic point of a vertex, and its somain
    * \return Value of the periodic point of a vertex, and the domain.
    */
-  inline long* GetPeriodicPointDomain(void) { return PeriodicPoint; }
+  inline long* GetPeriodicPointDomain() { return PeriodicPoint; }
 
   /*!
    * \brief Get the value of the periodic point of a vertex, and its somain
    * \return Value of the periodic point of a vertex, and the domain.
    */
-  inline bool GetActDisk_Perimeter(void) const { return ActDisk_Perimeter; }
+  inline bool GetActDisk_Perimeter() const { return ActDisk_Perimeter; }
 
   /*!
    * \brief Set the finite element basis functions needed for interpolation.
@@ -318,5 +319,15 @@ class CVertex : public CDualGrid {
    * \brief Get the value of the closest neighbor.
    * \return Index of the closest neighbor.
    */
-  inline unsigned long GetNormal_Neighbor(void) const { return Normal_Neighbor; }
+  inline unsigned long GetNormal_Neighbor() const { return Normal_Neighbor; }
+
+  /*!
+   * \brief Set the distance of the nearest internal neighbor that is not a wall.
+   */
+  inline void SetNearestNeighborDistance(const su2double& val) { NearestNeighborDist = val; }
+
+  /*!
+   * \brief Get the distance to the nearest internal neighbor that is not a wall.
+   */
+  inline const su2double& GetNearestNeighborDistance() const { return NearestNeighborDist; }
 };

Common/src/geometry/CGeometry.cpp

@@ -3902,6 +3902,36 @@ const CGeometry::CLineletInfo& CGeometry::GetLineletInfo(const CConfig* config)
   return li;
 }
 
+namespace {
+su2double NearestNeighborDistance(CGeometry* geometry, const CConfig* config, const unsigned long iPoint) {
+  const su2double max = std::numeric_limits<su2double>::max();
+  su2double distance = max;
+  for (const auto jPoint : geometry->nodes->GetPoints(iPoint)) {
+    const su2double dist = geometry->nodes->GetWall_Distance(jPoint);
+    if (dist > EPS) distance = fmin(distance, dist);
+  }
+  if (distance > 0 && distance < max) return distance;
+
+  /*--- The point only has wall neighbors, which all have 0 wall distance.
+   *    Compute an alternative distance based on volume and wall area. ---*/
+
+  const auto nDim = geometry->GetnDim();
+  su2double Normal[3] = {};
+  for (auto iMarker = 0u; iMarker < config->GetnMarker_All(); iMarker++) {
+    if (!config->GetViscous_Wall(iMarker)) continue;
+
+    const auto iVertex = geometry->nodes->GetVertex(iPoint, iMarker);
+    if (iVertex < 0) continue;
+
+    for (auto iDim = 0u; iDim < nDim; iDim++) {
+      Normal[iDim] += geometry->vertex[iMarker][iVertex]->GetNormal(iDim);
+    }
+  }
+  const su2double Vol = geometry->nodes->GetVolume(iPoint) + geometry->nodes->GetPeriodicVolume(iPoint);
+  return 2 * Vol / GeometryToolbox::Norm(nDim, Normal);
+}
+}  // namespace
+
 void CGeometry::ComputeWallDistance(const CConfig* const* config_container, CGeometry**** geometry_container) {
   int nZone = config_container[ZONE_0]->GetnZone();
   bool allEmpty = true;
@@ -3946,27 +3976,53 @@ void CGeometry::ComputeWallDistance(const CConfig* const* config_container, CGeo
         CGeometry* geometry = geometry_container[iZone][iInst][MESH_0];
         geometry->SetWallDistance(0.0);
       }
+      continue;
     }
-    /*--- Otherwise, set wall roughnesses. ---*/
-    if (!allEmpty) {
-      /*--- Store all wall roughnesses in a common data structure. ---*/
-      // [iZone][iMarker] -> roughness, for this rank
-      auto roughness_f = make_pair(nZone, [config_container, geometry_container, iInst](unsigned long iZone) {
-        const CConfig* config = config_container[iZone];
-        const auto nMarker = geometry_container[iZone][iInst][MESH_0]->GetnMarker();
-
-        return make_pair(nMarker, [config](unsigned long iMarker) {
-          return config->GetWallRoughnessProperties(config->GetMarker_All_TagBound(iMarker)).second;
-        });
+
+    /*--- Otherwise, set wall roughnesses, storing them in a common data structure. ---*/
+    // [iZone][iMarker] -> roughness, for this rank
+    auto roughness_f = make_pair(nZone, [config_container, geometry_container, iInst](unsigned long iZone) {
+      const CConfig* config = config_container[iZone];
+      const auto nMarker = geometry_container[iZone][iInst][MESH_0]->GetnMarker();
+
+      return make_pair(nMarker, [config](unsigned long iMarker) {
+        return config->GetWallRoughnessProperties(config->GetMarker_All_TagBound(iMarker)).second;
       });
-      NdFlattener<2> roughness_local(roughness_f);
-      // [rank][iZone][iMarker] -> roughness
-      NdFlattener<3> roughness_global(Nd_MPI_Environment(), roughness_local);
-      // use it to update roughnesses
-      for (int jZone = 0; jZone < nZone; jZone++) {
-        if (wallDistanceNeeded[jZone] && config_container[jZone]->GetnRoughWall() > 0) {
-          geometry_container[jZone][iInst][MESH_0]->nodes->SetWallRoughness(roughness_global);
+    });
+    NdFlattener<2> roughness_local(roughness_f);
+    // [rank][iZone][iMarker] -> roughness
+    NdFlattener<3> roughness_global(Nd_MPI_Environment(), roughness_local);
+    // use it to update roughnesses
+    for (int jZone = 0; jZone < nZone; jZone++) {
+      if (wallDistanceNeeded[jZone] && config_container[jZone]->GetnRoughWall() > 0) {
+        geometry_container[jZone][iInst][MESH_0]->nodes->SetWallRoughness(roughness_global);
+      }
+    }
+
+    for (int iZone = 0; iZone < nZone; iZone++) {
+      /*--- For the FEM solver, we use a different mesh structure ---*/
+      MAIN_SOLVER kindSolver = config_container[iZone]->GetKind_Solver();
+      if (!wallDistanceNeeded[iZone] || kindSolver == MAIN_SOLVER::FEM_LES || kindSolver == MAIN_SOLVER::FEM_RANS) {
+        continue;
+      }
+      const auto* config = config_container[iZone];
+      auto* geometry = geometry_container[iZone][iInst][MESH_0];
+
+      for (unsigned short iMarker = 0; iMarker < config->GetnMarker_All(); ++iMarker) {
+        const auto viscous = config->GetViscous_Wall(iMarker);
+
+        SU2_OMP_FOR_STAT(OMP_MIN_SIZE)
+        for (auto iVertex = 0u; iVertex < geometry->nVertex[iMarker]; iVertex++) {
+          const auto iPoint = geometry->vertex[iMarker][iVertex]->GetNode();
+          su2double dist = 0;
+          if (viscous && geometry->nodes->GetDomain(iPoint)) {
+            dist = NearestNeighborDistance(geometry, config, iPoint);
+          } else {
+            dist = geometry->nodes->GetWall_Distance(iPoint);
+          }
+          geometry->vertex[iMarker][iVertex]->SetNearestNeighborDistance(dist);
         }
+        END_SU2_OMP_FOR
       }
     }
   }

SU2_CFD/include/solvers/CFVMFlowSolverBase.inl

@@ -2350,7 +2350,7 @@ void CFVMFlowSolverBase<V, FlowRegime>::Friction_Forces(const CGeometry* geometr
 
   unsigned long iVertex, iPoint, iPointNormal;
   unsigned short iMarker, iMarker_Monitoring, iDim, jDim;
-  su2double Viscosity = 0.0, Area, Density = 0.0, WallDistMod, FrictionVel,
+  su2double Viscosity = 0.0, Area, Density = 0.0, FrictionVel,
             UnitNormal[3] = {0.0}, TauElem[3] = {0.0}, Tau[3][3] = {{0.0}}, Cp,
             thermal_conductivity, MaxNorm = 8.0, Grad_Vel[3][3] = {{0.0}}, Grad_Temp[3] = {0.0},
             Grad_Temp_ve[3] = {0.0}, AxiFactor;
@@ -2497,11 +2497,9 @@ void CFVMFlowSolverBase<V, FlowRegime>::Friction_Forces(const CGeometry* geometr
 
       FrictionVel = sqrt(fabs(WallShearStress[iMarker][iVertex]) / Density);
 
-      if (!wallfunctions && (MGLevel == MESH_0 || geometry->nodes->GetDomain(iPoint))) {
-        // for CMultiGridGeometry, the normal neighbor of halo nodes in not set
-        iPointNormal = geometry->vertex[iMarker][iVertex]->GetNormal_Neighbor();
-        Coord_Normal = geometry->nodes->GetCoord(iPointNormal);
-        WallDistMod = GeometryToolbox::Distance(nDim, Coord, Coord_Normal);
+      if (!wallfunctions && MGLevel == MESH_0 && geometry->nodes->GetDomain(iPoint)) {
+        // for CMultiGridGeometry and halos, the nearest neighbor distance is not set
+        const su2double WallDistMod = geometry->vertex[iMarker][iVertex]->GetNearestNeighborDistance();
         YPlus[iMarker][iVertex] = WallDistMod * FrictionVel / (Viscosity / Density);
       }
 

SU2_CFD/include/solvers/CTurbSSTSolver.hpp

@@ -53,7 +53,6 @@ class CTurbSSTSolver final : public CTurbSolver {
                      CSolver **solver_container,
                      const CConfig *config,
                      unsigned short val_marker);
-
 public:
   /*!
    * \brief Constructor.
@@ -302,5 +301,4 @@ class CTurbSSTSolver final : public CTurbSolver {
    */
   inline su2double GetOmega_Inf(void) const override { return Solution_Inf[1]; }
 
-
 };

SU2_CFD/include/variables/CFlowVariable.hpp

@@ -260,7 +260,7 @@ class CFlowVariable : public CVariable {
    * \param[in] iPoint - Point index.
    * \return Value of magnitude.
    */
-  inline su2double GetStrainMag(unsigned long iPoint) const final { return StrainMag(iPoint); }
+  inline su2double GetStrainMag(unsigned long iPoint) const { return StrainMag(iPoint); }
 
   /*!
    * \brief Get the entire vector of the rate of strain magnitude.

SU2_CFD/include/variables/CVariable.hpp

@@ -1147,13 +1147,6 @@ class CVariable {
   inline virtual su2double *GetVorticity(unsigned long iPoint) { return nullptr; }
   inline virtual const su2double *GetVorticity(unsigned long iPoint) const { return nullptr; }
 
-  /*!
-   * \brief A virtual member.
-   * \param[in] iPoint - Point index.
-   * \return Value of the rate of strain magnitude.
-   */
-  inline virtual su2double GetStrainMag(unsigned long iPoint) const { return 0.0; }
-
   /*!
    * \brief A virtual member.
    * \param[in] iPoint - Point index.

SU2_CFD/src/drivers/CDriver.cpp

@@ -833,7 +833,7 @@ void CDriver::InitializeGeometryFVM(CConfig *config, CGeometry **&geometry) {
     geometry[iMGlevel]->SetBoundControlVolume(geometry[iMGlevel-1], config, ALLOCATE);
     geometry[iMGlevel]->SetCoord(geometry[iMGlevel-1]);
 
-    /*--- Find closest neighbor to a surface point ---*/
+    /*--- Find closest, most normal, neighbor to a surface point ---*/
 
     geometry[iMGlevel]->FindNormal_Neighbor(config);
 
@@ -2517,7 +2517,7 @@ void CDriver::InitializeInterface(CConfig **config, CSolver***** solver, CGeomet
             else
               interface_type = NO_TRANSFER;
           }
-          
+
           if (interface_type != NO_TRANSFER) {
             auto nVar = 4;
             interface[donor][target] = new CConjugateHeatInterface(nVar, 0);

SU2_CFD/src/solvers/CIncNSSolver.cpp

@@ -685,11 +685,7 @@ void CIncNSSolver::SetTau_Wall_WF(CGeometry *geometry, CSolver **solver_containe
       const su2double VelTangMod = GeometryToolbox::Norm(int(MAXNDIM), VelTang);
 
       /*--- Compute normal distance of the interior point from the wall ---*/
-
-      su2double WallDist[MAXNDIM] = {0.0};
-      GeometryToolbox::Distance(nDim, Coord, Coord_Normal, WallDist);
-
-      su2double WallDistMod = GeometryToolbox::Norm(int(MAXNDIM), WallDist);
+      const su2double WallDistMod = GeometryToolbox::Distance(nDim, Coord, Coord_Normal);
 
       su2double Density_Wall = nodes->GetDensity(iPoint);