Fix curved symmetries (#2460)

* fix curved symmetries

Co-authored-by: Copilot Autofix powered by AI <62310815+github-advanced-security[bot]@users.noreply.github.com>
Co-authored-by: Pedro Gomes <[email protected]>
Co-authored-by: Pedro Gomes <[email protected]>

Author: 3 people

Common/include/geometry/CGeometry.hpp

@@ -243,6 +243,8 @@ class CGeometry {
   /*!< \brief Corrected normals on nodes with shared symmetry markers. */
   vector<std::unordered_map<unsigned long, std::array<su2double, MAXNDIM>>> symmetryNormals;
 
+  /*!< \brief Bool if boundary-marker is straight(2D)/plane(3D) for each local marker. */
+  vector<bool> boundIsStraight;
   vector<su2double> SurfaceAreaCfgFile; /*!< \brief Total Surface area for all markers. */
 
   /*--- Partitioning-specific variables ---*/
@@ -1005,6 +1007,15 @@ class CGeometry {
    */
   su2double GetSurfaceArea(const CConfig* config, unsigned short val_marker) const;
 
+  /*!
+   * \brief Check if a boundary is straight(2D) / plane(3D) for EULER_WALL and SYMMETRY_PLANE
+   *        only and store the information in boundIsStraight. For all other boundary types
+   *        this will return false and could therfore be wrong. Used ultimately for BC_Slip_Wall.
+   * \param[in] config - Definition of the particular problem.
+   * \param[in] print_on_screen - Boolean whether to print result on screen.
+   */
+  void ComputeSurfStraightness(const CConfig* config, bool print_on_screen);
+
   /*!
    * \brief Find and store all vertices on a sharp corner in the geometry.
    * \param[in] config - Definition of the particular problem.

Common/src/geometry/CGeometry.cpp

@@ -2456,16 +2456,66 @@ su2double CGeometry::GetSurfaceArea(const CConfig* config, unsigned short val_ma
 void CGeometry::ComputeModifiedSymmetryNormals(const CConfig* config) {
   /* Check how many symmetry planes there are and use the first (lowest ID) as the basis to orthogonalize against.
    * All nodes that are shared by multiple symmetries have to get a corrected normal. */
+
+  /*--- Compute if markers are straight lines or planes. ---*/
+  ComputeSurfStraightness(config, false);
+
+  symmetryNormals.clear();
   symmetryNormals.resize(nMarker);
-  std::vector<unsigned short> symMarkers;
+  std::vector<unsigned short> symMarkers, curvedSymMarkers;
 
+  /*--- Check which markers are symmetry or Euler wall and store in a list. ---*/
   for (auto iMarker = 0u; iMarker < nMarker; ++iMarker) {
     if ((config->GetMarker_All_KindBC(iMarker) == SYMMETRY_PLANE) ||
         (config->GetMarker_All_KindBC(iMarker) == EULER_WALL)) {
       symMarkers.push_back(iMarker);
+      if (!boundIsStraight[iMarker]) curvedSymMarkers.push_back(iMarker);
+    }
+  }
+
+  /*--- Merge the normals of curved markers to be independent of how surfaces are divided. ---*/
+
+  std::unordered_map<unsigned long, std::array<su2double, MAXNDIM>> mergedNormals;
+
+  for (const auto iMarker : curvedSymMarkers) {
+    for (auto iVertex = 0ul; iVertex < nVertex[iMarker]; iVertex++) {
+      const auto iPoint = vertex[iMarker][iVertex]->GetNode();
+
+      /*--- Determine if this point is shared with other curved symmetries. ---*/
+      int count = 0;
+      for (const auto jMarker : curvedSymMarkers) {
+        count += static_cast<int>(nodes->GetVertex(iPoint, jMarker) >= 0);
+      }
+      if (count < 2) continue;
+
+      std::array<su2double, MAXNDIM> normal = {};
+      vertex[iMarker][iVertex]->GetNormal(normal.data());
+
+      auto result = mergedNormals.emplace(iPoint, normal);
+      const auto inserted = result.second;
+      auto it = result.first;
+      if (!inserted) {
+        for (auto iDim = 0ul; iDim < MAXNDIM; ++iDim) {
+          it->second[iDim] += normal[iDim];
+        }
+      }
     }
   }
 
+  for (const auto& item : mergedNormals) {
+    const auto iPoint = item.first;
+    auto normal = item.second;
+    const su2double area = GeometryToolbox::Norm(int(MAXNDIM), normal.data());
+    for (auto iDim = 0ul; iDim < MAXNDIM; ++iDim) normal[iDim] /= area;
+
+    for (const auto iMarker : curvedSymMarkers) {
+      const auto iVertex = nodes->GetVertex(iPoint, iMarker);
+      if (iVertex >= 0) symmetryNormals[iMarker][iVertex] = normal;
+    }
+  }
+
+  /*--- Now do Gramm-Schmidt Process for symmetries that are not both curved lines or planes. ---*/
+
   /*--- Loop over all markers and find nodes on symmetry planes that are shared with other symmetries. ---*/
   /*--- The first symmetry does not need a corrected normal vector, hence start at 1. ---*/
   for (size_t i = 1; i < symMarkers.size(); ++i) {
@@ -2477,49 +2527,166 @@ void CGeometry::ComputeModifiedSymmetryNormals(const CConfig* config) {
       /*--- Halo points do not need to be considered. ---*/
       if (!nodes->GetDomain(iPoint)) continue;
 
-      /*--- Get the vertex normal on the current symmetry. ---*/
+      /*--- Get the vertex normal on the current symmetry, which may be a merged normal. ---*/
+      auto GetNormal = [&](unsigned long iMarker, unsigned long iVertex, std::array<su2double, MAXNDIM>& normal) {
+        const auto it = symmetryNormals[iMarker].find(iVertex);
+        if (it != symmetryNormals[iMarker].end()) {
+          normal = it->second;
+        } else {
+          vertex[iMarker][iVertex]->GetNormal(normal.data());
+          const su2double jarea = GeometryToolbox::Norm(int(MAXNDIM), normal.data());
+          for (auto iDim = 0ul; iDim < MAXNDIM; ++iDim) normal[iDim] /= jarea;
+        }
+      };
       std::array<su2double, MAXNDIM> iNormal = {};
-      vertex[iMarker][iVertex]->GetNormal(iNormal.data());
+      GetNormal(iMarker, iVertex, iNormal);
 
-      /*--- Loop over previous symmetries and if this point shares them, make this normal orthogonal to them. ---*/
-      bool isShared = false;
+      /*--- Loop over previous symmetries and if this point shares them, make this normal orthogonal to them.
+       * It's ok if we normalize merged normals against themselves, we get 0 area and this becomes a no-op. ---*/
 
       for (size_t j = 0; j < i; ++j) {
         const auto jMarker = symMarkers[j];
         const auto jVertex = nodes->GetVertex(iPoint, jMarker);
         if (jVertex < 0) continue;
 
-        isShared = true;
-
         std::array<su2double, MAXNDIM> jNormal = {};
-        const auto it = symmetryNormals[jMarker].find(jVertex);
-
-        if (it != symmetryNormals[jMarker].end()) {
-          jNormal = it->second;
-        } else {
-          vertex[jMarker][jVertex]->GetNormal(jNormal.data());
-          const su2double area = GeometryToolbox::Norm(nDim, jNormal.data());
-          for (auto iDim = 0u; iDim < nDim; iDim++) jNormal[iDim] /= area;
-        }
+        GetNormal(jMarker, jVertex, jNormal);
 
-        const auto proj = GeometryToolbox::DotProduct(nDim, jNormal.data(), iNormal.data());
-        for (auto iDim = 0u; iDim < nDim; iDim++) iNormal[iDim] -= proj * jNormal[iDim];
+        const su2double proj = GeometryToolbox::DotProduct(int(MAXNDIM), jNormal.data(), iNormal.data());
+        for (auto iDim = 0ul; iDim < MAXNDIM; ++iDim) iNormal[iDim] -= proj * jNormal[iDim];
       }
 
-      if (!isShared) continue;
-
       /*--- Normalize. If the norm is close to zero it means the normal is a linear combination of previous
        * normals, in this case we don't need to store the corrected normal, using the original in the gradient
        * correction will have no effect since previous corrections will remove components in this direction). ---*/
-      const su2double area = GeometryToolbox::Norm(nDim, iNormal.data());
-      if (area > EPS) {
-        for (auto iDim = 0u; iDim < nDim; iDim++) iNormal[iDim] /= area;
+      const su2double area = GeometryToolbox::Norm(int(MAXNDIM), iNormal.data());
+      if (area > 1e-12) {
+        for (auto iDim = 0ul; iDim < MAXNDIM; ++iDim) iNormal[iDim] /= area;
         symmetryNormals[iMarker][iVertex] = iNormal;
       }
     }
   }
 }
 
+void CGeometry::ComputeSurfStraightness(const CConfig* config, bool print_on_screen) {
+  bool RefUnitNormal_defined;
+  unsigned short iDim, iMarker, iMarker_Global, nMarker_Global = config->GetnMarker_CfgFile();
+  unsigned long iVertex;
+  constexpr passivedouble epsilon = 1.0e-6;
+  su2double Area;
+  string Local_TagBound, Global_TagBound;
+
+  vector<su2double> Normal(nDim), UnitNormal(nDim), RefUnitNormal(nDim);
+
+  /*--- Assume now that this boundary marker is straight. As soon as one
+        AreaElement is found that is not aligend with a Reference then it is
+        certain that the boundary marker is not straight and one can stop
+        searching. Another possibility is that this process doesn't own
+        any nodes of that boundary, in that case we also have to assume the
+        boundary is straight.
+        Any boundary type other than SYMMETRY_PLANE or EULER_WALL gets
+        the value false (or see cases specified in the conditional below)
+        which could be wrong. ---*/
+  boundIsStraight.resize(nMarker);
+  fill(boundIsStraight.begin(), boundIsStraight.end(), true);
+
+  /*--- Loop over all local markers ---*/
+  for (iMarker = 0; iMarker < nMarker; iMarker++) {
+    Local_TagBound = config->GetMarker_All_TagBound(iMarker);
+
+    /*--- Marker has to be Symmetry or Euler. Additionally marker can't be a
+          moving surface and Grid Movement Elasticity is forbidden as well. All
+          other GridMovements are rigid. ---*/
+    if ((config->GetMarker_All_KindBC(iMarker) == SYMMETRY_PLANE ||
+         config->GetMarker_All_KindBC(iMarker) == EULER_WALL) &&
+        !config->GetMarker_Moving_Bool(Local_TagBound) && !config->GetMarker_Deform_Mesh_Bool(Local_TagBound)) {
+      /*--- Loop over all global markers, and find the local-global pair via
+            matching unique string tags. ---*/
+      for (iMarker_Global = 0; iMarker_Global < nMarker_Global; iMarker_Global++) {
+        Global_TagBound = config->GetMarker_CfgFile_TagBound(iMarker_Global);
+        if (Local_TagBound == Global_TagBound) {
+          RefUnitNormal_defined = false;
+          iVertex = 0;
+
+          while (boundIsStraight[iMarker] && iVertex < nVertex[iMarker]) {
+            vertex[iMarker][iVertex]->GetNormal(Normal.data());
+            UnitNormal = Normal;
+
+            /*--- Compute unit normal. ---*/
+            Area = 0.0;
+            for (iDim = 0; iDim < nDim; iDim++) Area += Normal[iDim] * Normal[iDim];
+            Area = sqrt(Area);
+
+            /*--- Negate for outward convention. ---*/
+            for (iDim = 0; iDim < nDim; iDim++) UnitNormal[iDim] /= -Area;
+
+            /*--- Check if unit normal is within tolerance of the Reference unit normal.
+                  Reference unit normal = first unit normal found. ---*/
+            if (RefUnitNormal_defined) {
+              for (iDim = 0; iDim < nDim; iDim++) {
+                if (abs(RefUnitNormal[iDim] - UnitNormal[iDim]) > epsilon) {
+                  boundIsStraight[iMarker] = false;
+                  break;
+                }
+              }
+            } else {
+              RefUnitNormal = UnitNormal;  // deep copy of values
+              RefUnitNormal_defined = true;
+            }
+
+            iVertex++;
+          }  // while iVertex
+        }    // if Local == Global
+      }      // for iMarker_Global
+    } else {
+      /*--- Enforce default value: false ---*/
+      boundIsStraight[iMarker] = false;
+    }  // if sym or euler ...
+  }    // for iMarker
+
+  /*--- Communicate results and print on screen. ---*/
+  if (print_on_screen) {
+    /*--- Additional vector which can later be MPI::Allreduce(d) to pring the results
+          on screen as nMarker (local) can vary across ranks. Default 'true' as it can
+          happen that a local rank does not contain an element of each surface marker.  ---*/
+    vector<bool> boundIsStraightGlobal(nMarker_Global, true);
+    /*--- Match local with global tag bound and fill a Global Marker vector. ---*/
+    for (iMarker = 0; iMarker < nMarker; iMarker++) {
+      Local_TagBound = config->GetMarker_All_TagBound(iMarker);
+      for (iMarker_Global = 0; iMarker_Global < nMarker_Global; iMarker_Global++) {
+        Global_TagBound = config->GetMarker_CfgFile_TagBound(iMarker_Global);
+
+        if (Local_TagBound == Global_TagBound) boundIsStraightGlobal[iMarker_Global] = boundIsStraight[iMarker];
+
+      }  // for iMarker_Global
+    }    // for iMarker
+
+    vector<int> Buff_Send_isStraight(nMarker_Global), Buff_Recv_isStraight(nMarker_Global);
+
+    /*--- Cast to int as std::vector<boolean> can be a special construct. MPI handling using <int>
+          is more straight-forward. ---*/
+    for (iMarker_Global = 0; iMarker_Global < nMarker_Global; iMarker_Global++)
+      Buff_Send_isStraight[iMarker_Global] = static_cast<int>(boundIsStraightGlobal[iMarker_Global]);
+
+    /*--- Product of type <int>(bool) is equivalnt to a 'logical and' ---*/
+    SU2_MPI::Allreduce(Buff_Send_isStraight.data(), Buff_Recv_isStraight.data(), nMarker_Global, MPI_INT, MPI_PROD,
+                       SU2_MPI::GetComm());
+
+    /*--- Print results on screen. ---*/
+    if (rank == MASTER_NODE) {
+      for (iMarker_Global = 0; iMarker_Global < nMarker_Global; iMarker_Global++) {
+        if (config->GetMarker_CfgFile_KindBC(config->GetMarker_CfgFile_TagBound(iMarker_Global)) == SYMMETRY_PLANE ||
+            config->GetMarker_CfgFile_KindBC(config->GetMarker_CfgFile_TagBound(iMarker_Global)) == EULER_WALL) {
+          cout << "Boundary marker " << config->GetMarker_CfgFile_TagBound(iMarker_Global) << " is";
+          if (!static_cast<bool>(Buff_Recv_isStraight[iMarker_Global])) cout << " NOT";
+          if (nDim == 2) cout << " a single straight." << endl;
+          if (nDim == 3) cout << " a single plane." << endl;
+        }  // if sym or euler
+      }    // for iMarker_Global
+    }      // if rank==MASTER
+  }        // if print_on_scren
+}
+
 void CGeometry::ComputeSurf_Curvature(CConfig* config) {
   unsigned short iMarker, iNeigh_Point, iDim, iNode, iNeighbor_Nodes, Neighbor_Node;
   unsigned long Neighbor_Point, iVertex, iPoint, jPoint, iElem_Bound, iEdge, nLocalVertex, MaxLocalVertex,

TestCases/disc_adj_fsi/dyn_fsi/grad_dv.opt.ref

@@ -1,9 +1,9 @@
 INDEX	GRAD
-0	-5.239274142145693e-03
-1	-2.822866261832595e-03
-2	-1.251166631047082e-03
-3	-4.746312021502222e-04
-4	-4.726421390156573e-04
-5	-1.242307153964931e-03
-6	-2.797799332331688e-03
-7	-5.181304648315885e-03
+0	-3.533163525398154e-03
+1	-1.879941699367452e-03
+2	-8.078436962531523e-04
+3	-2.794845540524826e-04
+4	-2.790760674146328e-04
+5	-8.045490756184532e-04
+6	-1.867166054373919e-03
+7	-3.498473613679959e-03

TestCases/hybrid_regression.py

@@ -471,7 +471,7 @@ def main():
     cosine_gust.cfg_dir   = "gust"
     cosine_gust.cfg_file  = "cosine_gust_zdir.cfg"
     cosine_gust.test_iter = 79
-    cosine_gust.test_vals = [-2.418805, 0.002013, -0.001504, 0.000445, -0.000558]
+    cosine_gust.test_vals = [-2.418806, 0.001963, -0.001264, 0.000415, -0.000592]
     cosine_gust.unsteady  = True
     cosine_gust.enabled_with_tsan = False
     test_list.append(cosine_gust)
@@ -543,7 +543,7 @@ def main():
     edge_PPR.cfg_dir   = "nicf/edge"
     edge_PPR.cfg_file  = "edge_PPR.cfg"
     edge_PPR.test_iter = 40
-    edge_PPR.test_vals         = [-6.922907, -0.757793, -0.000034, 0.000000]
+    edge_PPR.test_vals         = [-6.923001, -0.757902, -0.000034, 0.000000]
     edge_PPR.test_vals_aarch64 = [-8.573595, -2.391849, -0.000034, 0.000000]
     test_list.append(edge_PPR)
 

TestCases/hybrid_regression_AD.py

@@ -58,7 +58,7 @@ def main():
     discadj_cylinder3D.cfg_dir   = "disc_adj_euler/cylinder3D"
     discadj_cylinder3D.cfg_file  = "inv_cylinder3D.cfg"
     discadj_cylinder3D.test_iter = 5
-    discadj_cylinder3D.test_vals = [-3.764562, -3.719982, -0.000000, 0.000000]
+    discadj_cylinder3D.test_vals = [-3.773755, -3.737672, -0.000000, 0.000000]
     test_list.append(discadj_cylinder3D)
 
     # Arina nozzle 2D