Reworked consecutive handling for multiple dv markers.

thomasdick · thomasdick · commit 576d52a460dd · 2022-01-11T14:44:38.000+01:00
diff --git a/SU2_CFD/include/solvers/CGradientSmoothingSolver.hpp b/SU2_CFD/include/solvers/CGradientSmoothingSolver.hpp
@@ -66,6 +66,8 @@ class CGradientSmoothingSolver final : public CFEASolverBase {
 
   CSquareMatrixCM hessian;                   /*!< \brief The approximated Hessian with respect to the design variables. */
 
+  std::vector<bool> visited;                 /*! <\brief Stores already visited points for surface applications with multiple markers. */
+
   /*!
    * \brief The highest level in the variable hierarchy all derived solvers can safely use,
    * CVariable is the common denominator between the FEA and Mesh deformationd variables.
@@ -190,9 +192,10 @@ class CGradientSmoothingSolver final : public CFEASolverBase {
 
   /*!
    * \brief Record a tape containing the parameter Jacobian.
-   * \param geometry
-   * \param config
-   * \param surface_movement
+   * \param geometry - Geometrical definition of the problem.
+   * \param surface_movement - Surface movement of the problem.
+   * \param registeredCoord - Indexes of the affected coordinates.
+   * \param config - Definition of the particular problem.
    */
   void RecordParameterizationJacobian(CGeometry *geometry,
                                       CSurfaceMovement *surface_movement,
@@ -201,9 +204,6 @@ class CGradientSmoothingSolver final : public CFEASolverBase {
 
   /*!
    * \brief Forward evaluate parameterization Jacobian.
-   * \param geometry
-   * \param config
-   * \param surface_movement
    */
   void ProjectDVtoMesh(CGeometry *geometry,
                        std::vector<su2double>& seeding,
@@ -213,9 +213,6 @@ class CGradientSmoothingSolver final : public CFEASolverBase {
 
   /*!
    * \brief Reverse evaluate parameterization Jacobian.
-   * \param geometry
-   * \param config
-   * \param surface_movement
    */
   void ProjectMeshToDV(CGeometry *geometry,
                        CSysVector<su2matvecscalar>& sensitivity,
@@ -235,6 +232,14 @@ class CGradientSmoothingSolver final : public CFEASolverBase {
    */
   void ReadSensFromGeometry(const CGeometry *geometry) override;
 
+  /*!
+   * \brief Helper function for derivative treatment.
+   *        to determine if a boundary marker is a DVmarker on any mpi rank.
+   * \param[in] iMarker - the current boundary marker.
+   * \param[in] config - Definition of the particular problem.
+   */
+  bool MarkerIsDVMarker(unsigned short iMarker, const CConfig *config) const override;
+
  private:
   /*!
    * \brief Write the solution of the linear solver into the sensitivities of the nodes
@@ -268,6 +273,12 @@ class CGradientSmoothingSolver final : public CFEASolverBase {
    */
   void Set_VertexEliminationSchedule(CGeometry* geometry, const CConfig* config);
 
+  /*!
+   * \brief Complete the calculation of the surface stiffness matrix
+   * \param[in] geometry - Geometrical definition of the problem.
+   */
+  void Complete_Surface_StiffMatrix(const CGeometry* geometry);
+
   /*!
    * \brief Get the element container index and number of nodes of a given VTK type.
    * \param[in] VTK_Type - Type of element.
diff --git a/SU2_CFD/include/solvers/CSolver.hpp b/SU2_CFD/include/solvers/CSolver.hpp
@@ -4188,6 +4188,14 @@ class CSolver {
    */
   virtual void WriteSensToGeometry(CGeometry* geometry) const {}
 
+  /*!
+   * \brief Helper function for derivative treatment.
+   *        to determine if a boundary marker is a DVmarker on any mpi rank.
+   * \param[in] imMarker - the current boundary marker.
+   * \param[in] config - Definition of the particular problem.
+   */
+  virtual bool MarkerIsDVMarker(unsigned short iMarker, const CConfig *config) const { return false; };
+
   /*!
    * \brief Routine that sets the flag controlling implicit treatment for periodic BCs.
    * \param[in] val_implicit_periodic - Flag controlling implicit treatment for periodic BCs.
diff --git a/SU2_CFD/src/solvers/CGradientSmoothingSolver.cpp b/SU2_CFD/src/solvers/CGradientSmoothingSolver.cpp
@@ -95,6 +95,7 @@ CGradientSmoothingSolver::CGradientSmoothingSolver(CGeometry *geometry, CConfig
       LinSysRes.Initialize(nPoint, nPointDomain, 1, 0.0);
       Jacobian.Initialize(nPoint, nPointDomain, 1, 1, false, geometry, config, false, true);
     }
+    visited.resize(geometry->GetnPoint(), false);
   }
 
   /*--- vectors needed for projection when working on the complete system ---*/
@@ -201,9 +202,11 @@ void CGradientSmoothingSolver::ApplyGradientSmoothingSurface(CGeometry* geometry
   LinSysSol.SetValZero();
   LinSysRes.SetValZero();
   Jacobian.SetValZero();
+  std::fill(visited.begin(), visited.end(), false);
 
   /*--- Compute the stiffness matrix for the smoothing operator. ---*/
   Compute_Surface_StiffMatrix(geometry, numerics, config, val_marker);
+  Complete_Surface_StiffMatrix(geometry);
 
   Compute_Surface_Residual(geometry, config, val_marker);
 
@@ -397,18 +400,12 @@ void CGradientSmoothingSolver::Compute_Surface_StiffMatrix(CGeometry* geometry,
   std::array<unsigned long, MAXNNODE_2D> indexVertex;
   int EL_KIND = 0;
   su2double val_Coord, HiHj = 0.0;
-  su2activematrix DHiDHj, Jacobian_block, mId_Aux;
+  su2activematrix DHiDHj, Jacobian_block;
 
   Jacobian_block.resize(nDim, nDim) = su2double(0.0);
-  mId_Aux.resize(nDim, nDim) = su2double(0.0);
-  for (iDim = 0; iDim < nDim; iDim++) {
-    mId_Aux[iDim][iDim] = su2double(1.0);
-  }
-
-  vector<bool> visited(geometry->GetnPoint(), false);
 
   /*--- Check if the current MPI rank has a part of the marker ---*/
-  if (val_marker!=NOT_AVAILABLE) {
+  if (val_marker!=BC_TYPE::NOT_AVAILABLE) {
 
     /*--- Loops over all the elements ---*/
     for (iElem = 0; iElem < geometry->GetnElem_Bound(val_marker); iElem++) {
@@ -463,13 +460,6 @@ void CGradientSmoothingSolver::Compute_Surface_StiffMatrix(CGeometry* geometry,
     }
   }
 
-  /*--- Assembling the stiffness matrix on the design surface means the Jacobian is the identity for nodes inside the domain. ---*/
-  for (iPoint = 0; iPoint < geometry->GetnPointDomain(); iPoint++){
-    if (visited[iPoint]==false) {
-      Jacobian.AddBlock(iPoint, iPoint, mId_Aux);
-    }
-  }
-
 }
 
 void CGradientSmoothingSolver::Compute_Residual(CGeometry* geometry, const CConfig* config) {
@@ -539,7 +529,7 @@ void CGradientSmoothingSolver::Compute_Surface_Residual(CGeometry* geometry, con
   su2double* normal = NULL;
 
   /*--- Check if the current MPI rank has a part of the marker ---*/
-  if (val_marker!=NOT_AVAILABLE) {
+  if (val_marker!=BC_TYPE::NOT_AVAILABLE) {
 
     for (iElem = 0; iElem < geometry->GetnElem_Bound(val_marker); iElem++) {
 
@@ -637,7 +627,7 @@ void CGradientSmoothingSolver::BC_Surface_Dirichlet(const CGeometry* geometry, c
   const su2double zeros[MAXNDIM] = {0.0};
 
   /*--- Check if the current MPI rank has a part of the marker ---*/
-  if (val_marker!=NOT_AVAILABLE) {
+  if (val_marker!=BC_TYPE::NOT_AVAILABLE) {
     for (iVertex = 0; iVertex < geometry->nVertex[val_marker]; iVertex++) {
 
       /*--- Get node index ---*/
@@ -680,17 +670,15 @@ template <typename scalar_type>
 CSysMatrixVectorProduct<scalar_type> CGradientSmoothingSolver::GetStiffnessMatrixVectorProduct(CGeometry* geometry,
                                                                                                CNumerics* numerics,
                                                                                                const CConfig* config) {
-  bool surf = config->GetSmoothOnSurface();
 
   /*--- Compute the sparse stiffness matrix ---*/
-  if (surf) {
-    unsigned long dvMarker=NOT_AVAILABLE;
-    for (unsigned int iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) {
-      if ( config->GetMarker_All_DV(iMarker) == YES ) {
-        dvMarker=iMarker;
+  if (config->GetSmoothOnSurface()) {
+    for (unsigned int iMarker = 0; iMarker < config->GetnMarker_CfgFile(); iMarker++) {
+      if (config->GetMarker_All_DV(iMarker) == YES) {
+        Compute_Surface_StiffMatrix(geometry, numerics, config, iMarker, nDim);
       }
     }
-    Compute_Surface_StiffMatrix(geometry, numerics, config, dvMarker, nDim);
+    Complete_Surface_StiffMatrix(geometry);
   } else {
     Compute_StiffMatrix(geometry, numerics, config);
   }
@@ -917,7 +905,7 @@ void CGradientSmoothingSolver::WriteSensitivity(CGeometry* geometry, const CConf
 
   /*--- split between surface and volume first, to avoid mpi ranks with no part of the marker to write back nonphysical solutions in the surface case ---*/
   if ( config->GetSmoothOnSurface() ) {
-    if( val_marker!=NOT_AVAILABLE ) {
+    if( val_marker!=BC_TYPE::NOT_AVAILABLE ) {
       for (unsigned long iVertex =0; iVertex<geometry->nVertex[val_marker]; iVertex++)  {
         iPoint = geometry->vertex[val_marker][iVertex]->GetNode();
         normal = geometry->vertex[val_marker][iVertex]->GetNormal();
@@ -980,23 +968,25 @@ void CGradientSmoothingSolver::Set_VertexEliminationSchedule(CGeometry* geometry
   if (config->GetSmoothOnSurface()) {
 
     /*--- Find Marker_DV if this node has any ---*/
-    unsigned long dvMarker=NOT_AVAILABLE;
+    unsigned long dvMarker=BC_TYPE::NOT_AVAILABLE;
     for (unsigned short iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) {
-      if (config->GetMarker_All_DV(iMarker) == YES) {
-        dvMarker = iMarker;
-      }
-    }
-
-    /*--- Surface case:
-     * Fix design boundary border if Dirichlet condition is set and the current rank holds part of it. ---*/
-    if ( config->GetDirichletSurfaceBound() && dvMarker!=NOT_AVAILABLE) {
-      for (iVertex = 0; iVertex < geometry->nVertex[dvMarker]; iVertex++) {
-        /*--- Get node index ---*/
-        iPoint = geometry->vertex[dvMarker][iVertex]->GetNode();
-        if (nodes->GetIsBoundaryPoint(iPoint)) {
-          myPoints.push_back(geometry->nodes->GetGlobalIndex(iPoint));
+      //if (MarkerIsDVMarker(iMarker, config)) {
+        //dvMarker = BC_TYPE::NOT_AVAILABLE;
+        if (config->GetMarker_All_DV(iMarker) == YES) dvMarker = iMarker;
+      //}
+
+      /*--- Surface case:
+       * Fix design boundary border if Dirichlet condition is set and the current rank holds part of it. ---*/
+      if ( config->GetDirichletSurfaceBound() && dvMarker!=BC_TYPE::NOT_AVAILABLE) {
+        for (iVertex = 0; iVertex < geometry->nVertex[dvMarker]; iVertex++) {
+          /*--- Get node index ---*/
+          iPoint = geometry->vertex[dvMarker][iVertex]->GetNode();
+          if (nodes->GetIsBoundaryPoint(iPoint)) {
+            myPoints.push_back(geometry->nodes->GetGlobalIndex(iPoint));
+          }
         }
       }
+
     }
 
   } else {
@@ -1024,3 +1014,36 @@ void CGradientSmoothingSolver::Set_VertexEliminationSchedule(CGeometry* geometry
     Jacobian.EnforceSolutionAtNode(iPoint, LinSysSol.GetBlock(iPoint), LinSysRes);
   }
 }
+
+bool CGradientSmoothingSolver::MarkerIsDVMarker(unsigned short iMarker, const CConfig *config) const {
+
+  int isLocalDVmarker, isGlobalDVmarker;
+  bool isDVmarker=false;
+
+  if ( config->GetMarker_All_DV(iMarker) == YES ) {
+    isLocalDVmarker = 1;
+  } else {
+    isLocalDVmarker = 0;
+  }
+  SU2_MPI::Allreduce(&isLocalDVmarker, &isGlobalDVmarker, 1, MPI_UNSIGNED_LONG, MPI_SUM, SU2_MPI::GetComm());
+  if (isGlobalDVmarker >= 1) isDVmarker=true;
+
+  return isDVmarker;
+}
+
+void CGradientSmoothingSolver::Complete_Surface_StiffMatrix(const CGeometry* geometry) {
+
+  su2activematrix mId_Aux;
+  mId_Aux.resize(nDim, nDim) = su2double(0.0);
+  for (unsigned iDim = 0; iDim < nDim; iDim++) {
+    mId_Aux[iDim][iDim] = su2double(1.0);
+  }
+
+  /*--- Assembling the stiffness matrix on the design surface means the Jacobian is the identity for nodes inside the domain. ---*/
+  for (unsigned long iPoint = 0; iPoint < geometry->GetnPointDomain(); iPoint++){
+    if (visited[iPoint]==false) {
+      Jacobian.AddBlock(iPoint, iPoint, mId_Aux);
+    }
+  }
+
+}
diff --git a/SU2_DOT/src/SU2_DOT.cpp b/SU2_DOT/src/SU2_DOT.cpp
@@ -1003,19 +1003,15 @@ void DerivativeTreatment_MeshSensitivity(CGeometry *geometry, CConfig *config, C
       /*--- Get the sensitivities from the geometry class to work with. ---*/
       solver->ReadSensFromGeometry(geometry);
 
-      /*--- Select DV marker, or NOT_AVAILABLE if none is specified. ---*/
-      unsigned long ndvMarker=0, dvMarker = BC_TYPE::NOT_AVAILABLE;
-      for (unsigned short iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) {
-        if ( config->GetMarker_All_DV(iMarker) == YES ) {
-          dvMarker = iMarker;
-          ++ndvMarker;
+      /*--- Perform the smoothing procedure sequentially for each DV marker. ---*/
+      unsigned long dvMarker;
+      for (unsigned short iMarker = 0; iMarker < config->GetnMarker_CfgFile(); iMarker++) {
+        if (solver->MarkerIsDVMarker(iMarker, config)) {
+          dvMarker = BC_TYPE::NOT_AVAILABLE;
+          if (config->GetMarker_All_DV(iMarker) == YES) dvMarker = iMarker;
+          solver->ApplyGradientSmoothingSurface(geometry, numerics, config, dvMarker);
         }
       }
-      if(ndvMarker>1) {
-        SU2_MPI::Error("Sobolev smoothing in surface mode only works with one marker at a time.",CURRENT_FUNCTION);
-      }
-
-      solver->ApplyGradientSmoothingSurface(geometry, numerics, config, dvMarker);
 
       /*--- After appling the solver write the results back ---*/
       solver->WriteSensToGeometry(geometry);
