Multizone py wrapper example (FSI) (#2111)

* checkpoint naming and update some functions

* add FSI example with additional custom load via the python wrapper

* remove mesh files

* fix segfault

* add some arm64 values

* disable polar naca on arm64

Author: pcarruscag

SU2_CFD/include/drivers/CDriverBase.hpp

@@ -54,7 +54,7 @@ class CDriverBase {
       UsedTime;        /*!< \brief Elapsed time between Start and Stop point of the timer. */
 
   unsigned long TimeIter;
-  unsigned short selected_iZone = ZONE_0; /*!< \brief Selected zone for the driver. Defaults to ZONE_0 */
+  unsigned short selected_zone = ZONE_0; /*!< \brief Selected zone for the driver. Defaults to ZONE_0 */
   unsigned short iMesh,  /*!< \brief Iterator on mesh levels. */
       iZone,             /*!< \brief Iterator on zones. */
       nZone,             /*!< \brief Total number of zones in the problem. */
@@ -227,7 +227,7 @@ class CDriverBase {
       SU2_MPI::Error("Initial coordinates are only available with DEFORM_MESH= YES", CURRENT_FUNCTION);
     }
     auto* coords =
-        const_cast<su2activematrix*>(solver_container[selected_iZone][INST_0][MESH_0][MESH_SOL]->GetNodes()->GetMesh_Coord());
+        const_cast<su2activematrix*>(solver_container[selected_zone][INST_0][MESH_0][MESH_SOL]->GetNodes()->GetMesh_Coord());
     return CPyWrapperMatrixView(*coords, "InitialCoordinates", true);
   }
 
@@ -241,7 +241,7 @@ class CDriverBase {
     if (iMarker >= GetNumberMarkers()) SU2_MPI::Error("Marker index exceeds size.", CURRENT_FUNCTION);
 
     auto* coords =
-        const_cast<su2activematrix*>(solver_container[selected_iZone][INST_0][MESH_0][MESH_SOL]->GetNodes()->GetMesh_Coord());
+        const_cast<su2activematrix*>(solver_container[selected_zone][INST_0][MESH_0][MESH_SOL]->GetNodes()->GetMesh_Coord());
     return CPyWrapperMarkerMatrixView(*coords, main_geometry->vertex[iMarker], main_geometry->GetnVertex(iMarker),
                                       "MarkerInitialCoordinates", true);
   }
@@ -550,16 +550,21 @@ class CDriverBase {
   }
 
   /*!
-   * \brief Selects zone to be used for Driver operation
+   * \brief Selects zone to be used for python driver operations.
    * \param[in] iZone - Zone identifier.
    */
   inline void SelectZone(unsigned short iZone) {
     if (iZone >= nZone) SU2_MPI::Error("Zone index out of range", CURRENT_FUNCTION);
-    selected_iZone = iZone;
-    main_geometry = geometry_container[selected_iZone][INST_0][MESH_0];
-    main_config = config_container[selected_iZone];
+    selected_zone = iZone;
+    main_geometry = geometry_container[selected_zone][INST_0][MESH_0];
+    main_config = config_container[selected_zone];
   }
 
+  /*!
+   * \brief Returns the index of the zone selected for python driver operations.
+   */
+  inline unsigned short SelectedZone() const { return selected_zone; }
+
   /*!
    * \brief Get the wall normal heat flux at a vertex on a specified marker of the flow or heat solver.
    * \note This can be the output of a heat or flow solver in a CHT setting.
@@ -707,7 +712,7 @@ class CDriverBase {
     if (iMarker < std::numeric_limits<unsigned short>::max() && iMarker > GetNumberMarkers()) {
       SU2_MPI::Error("Marker index exceeds size.", CURRENT_FUNCTION);
     }
-    auto* solver = solver_container[selected_iZone][INST_0][MESH_0][iSolver];
+    auto* solver = solver_container[selected_zone][INST_0][MESH_0][iSolver];
     if (solver == nullptr) SU2_MPI::Error("The selected solver does not exist.", CURRENT_FUNCTION);
     return solver;
   }

SU2_CFD/include/drivers/CMultizoneDriver.hpp

@@ -122,18 +122,34 @@ class CMultizoneDriver : public CDriver {
   /*!
    * \brief Destructor of the class.
    */
-  ~CMultizoneDriver(void) override;
+  ~CMultizoneDriver() override;
 
   /*!
    * \brief [Overload] Launch the computation for multizone problems.
    */
   void StartSolver() override;
 
   /*!
-   * \brief Preprocess the multizone iteration
+   * \brief Preprocess the multizone iteration.
    */
   void Preprocess(unsigned long TimeIter) override;
 
+  /*!
+   * \brief Solves one time iteration.
+   */
+  void Run() override {
+    switch (driver_config->GetKind_MZSolver()){
+      case ENUM_MULTIZONE::MZ_BLOCK_GAUSS_SEIDEL: RunGaussSeidel(); break;  // Block Gauss-Seidel iteration
+      case ENUM_MULTIZONE::MZ_BLOCK_JACOBI: RunJacobi(); break;             // Block-Jacobi iteration
+    }
+  }
+
+  /*!
+   * \brief Placeholder for post processing operations to make the interface
+   * of this driver identical to CSinglezoneDriver.
+   */
+  void Postprocess() {}
+
   /*!
    * \brief Update the dual-time solution within multiple zones.
    */

SU2_CFD/include/interfaces/fsi/CFlowTractionInterface.hpp

@@ -43,8 +43,11 @@ class CFlowTractionInterface : public CInterface {
   /*!
    * \brief Sets the dimensional factor for pressure and the consistent_interpolation flag.
    * \param[in] flow_config - Definition of the fluid (donor) problem.
+   * \param[in] struct_config - Definition of the structural (target) problem.
+   * \param[in] geometry - FEA geometry.
+   * \param[in] solution - FEA solver.
    */
-  void Preprocess(const CConfig *flow_config);
+  void Preprocess(const CConfig *flow_config, const CConfig *struct_config, CGeometry *geometry, CSolver *solution);
 
   /*!
    * \brief Computes vertex areas (FEA side) for when tractions need to be integrated.

SU2_CFD/include/variables/CFEABoundVariable.hpp

@@ -149,11 +149,6 @@ class CFEABoundVariable final : public CFEAVariable {
     return FlowTraction_n(iPoint,iVar);
   }
 
-  /*!
-   * \brief Clear the flow traction residual
-   */
-  void Clear_FlowTraction() override;
-
   /*!
    * \brief Register the flow tractions as input variable.
    */

SU2_CFD/include/variables/CVariable.hpp

@@ -1508,11 +1508,6 @@ class CVariable {
    */
   inline virtual su2double Get_FlowTraction_n(unsigned long iPoint, unsigned long iVar) const { return 0.0; }
 
-  /*!
-   * \brief A virtual member.
-   */
-  inline virtual void Clear_FlowTraction() {}
-
   /*!
    * \brief A virtual member.
    */

SU2_CFD/src/drivers/CDriverBase.cpp

@@ -391,14 +391,14 @@ vector<passivedouble> CDriverBase::GetMarkerVertexNormals(unsigned short iMarker
 }
 
 void CDriverBase::CommunicateMeshDisplacements() {
-  solver_container[selected_iZone][INST_0][MESH_0][MESH_SOL]->InitiateComms(main_geometry, main_config, MESH_DISPLACEMENTS);
-  solver_container[selected_iZone][INST_0][MESH_0][MESH_SOL]->CompleteComms(main_geometry, main_config, MESH_DISPLACEMENTS);
+  solver_container[selected_zone][INST_0][MESH_0][MESH_SOL]->InitiateComms(main_geometry, main_config, MESH_DISPLACEMENTS);
+  solver_container[selected_zone][INST_0][MESH_0][MESH_SOL]->CompleteComms(main_geometry, main_config, MESH_DISPLACEMENTS);
 }
 
 map<string, unsigned short> CDriverBase::GetSolverIndices() const {
   map<string, unsigned short> indexMap;
   for (auto iSol = 0u; iSol < MAX_SOLS; iSol++) {
-    const auto* solver = solver_container[selected_iZone][INST_0][MESH_0][iSol];
+    const auto* solver = solver_container[selected_zone][INST_0][MESH_0][iSol];
     if (solver != nullptr) {
       if (solver->GetSolverName().empty()) SU2_MPI::Error("Solver name was not defined.", CURRENT_FUNCTION);
       indexMap[solver->GetSolverName()] = iSol;
@@ -408,7 +408,7 @@ map<string, unsigned short> CDriverBase::GetSolverIndices() const {
 }
 
 std::map<string, unsigned short> CDriverBase::GetFEASolutionIndices() const {
-  if (solver_container[selected_iZone][INST_0][MESH_0][FEA_SOL] == nullptr) {
+  if (solver_container[selected_zone][INST_0][MESH_0][FEA_SOL] == nullptr) {
     SU2_MPI::Error("The FEA solver does not exist.", CURRENT_FUNCTION);
   }
   const auto nDim = main_geometry->GetnDim();
@@ -429,7 +429,7 @@ std::map<string, unsigned short> CDriverBase::GetFEASolutionIndices() const {
 }
 
 map<string, unsigned short> CDriverBase::GetPrimitiveIndices() const {
-  if (solver_container[selected_iZone][INST_0][MESH_0][FLOW_SOL] == nullptr) {
+  if (solver_container[selected_zone][INST_0][MESH_0][FLOW_SOL] == nullptr) {
     SU2_MPI::Error("The flow solver does not exist.", CURRENT_FUNCTION);
   }
   return PrimitiveNameToIndexMap(CPrimitiveIndices<unsigned short>(

SU2_CFD/src/drivers/CMultizoneDriver.cpp

@@ -182,10 +182,7 @@ void CMultizoneDriver::StartSolver() {
 
     /*--- Run a block iteration of the multizone problem. ---*/
 
-    switch (driver_config->GetKind_MZSolver()){
-      case ENUM_MULTIZONE::MZ_BLOCK_GAUSS_SEIDEL: RunGaussSeidel(); break;  // Block Gauss-Seidel iteration
-      case ENUM_MULTIZONE::MZ_BLOCK_JACOBI: RunJacobi(); break;             // Block-Jacobi iteration
-    }
+    Run();
 
     /*--- Update the solution for dual time stepping strategy ---*/
 
@@ -646,11 +643,11 @@ void CMultizoneDriver::SetTurboPerformance() {
   }
 }
 
-bool CMultizoneDriver::Monitor(unsigned long TimeIter){
+bool CMultizoneDriver::Monitor(unsigned long TimeIter) {
 
   /*--- Check whether the inner solver has converged --- */
 
-  if (driver_config->GetTime_Domain() == NO){
+  if (driver_config->GetTime_Domain() == NO) {
 
     const auto OuterIter  = driver_config->GetOuterIter();
     const auto nOuterIter = driver_config->GetnOuter_Iter();
@@ -670,29 +667,27 @@ bool CMultizoneDriver::Monitor(unsigned long TimeIter){
   }
   // i.e. unsteady simulation
 
-    /*--- Check whether the outer time integration has reached the final time ---*/
-    const auto TimeConvergence = GetTimeConvergence();
+  /*--- Check whether the outer time integration has reached the final time. ---*/
+  const auto TimeConvergence = GetTimeConvergence();
 
-    const auto nTimeIter = driver_config->GetnTime_Iter();
-    const auto MaxTime   = driver_config->GetMax_Time();
-    const auto CurTime   = driver_output->GetHistoryFieldValue("CUR_TIME");
+  const auto nTimeIter = driver_config->GetnTime_Iter();
+  const auto MaxTime   = driver_config->GetMax_Time();
+  const auto CurTime   = driver_output->GetHistoryFieldValue("CUR_TIME");
 
-    const bool FinalTimeReached = (CurTime >= MaxTime);
-    const bool MaxIterationsReached = (TimeIter+1 >= nTimeIter);
+  const bool FinalTimeReached = (CurTime >= MaxTime);
+  const bool MaxIterationsReached = (TimeIter+1 >= nTimeIter);
 
-    if ((TimeConvergence || FinalTimeReached || MaxIterationsReached) && (rank == MASTER_NODE)){
-      cout << "\n----------------------------- Solver Exit -------------------------------";
-      if (TimeConvergence)  cout << "\nAll windowed time-averaged convergence criteria are fullfilled." << endl;
-      if (FinalTimeReached) cout << "\nMaximum time reached (MAX_TIME = " << MaxTime << "s)." << endl;
-      else cout << "\nMaximum number of time iterations reached (TIME_ITER = " << nTimeIter << ")." << endl;
-      cout << "-------------------------------------------------------------------------" << endl;
-    }
-
-    return (FinalTimeReached || MaxIterationsReached);
- 
+  if ((TimeConvergence || FinalTimeReached || MaxIterationsReached) && (rank == MASTER_NODE)){
+    cout << "\n----------------------------- Solver Exit -------------------------------";
+    if (TimeConvergence)  cout << "\nAll windowed time-averaged convergence criteria are fullfilled." << endl;
+    if (FinalTimeReached) cout << "\nMaximum time reached (MAX_TIME = " << MaxTime << "s)." << endl;
+    else cout << "\nMaximum number of time iterations reached (TIME_ITER = " << nTimeIter << ")." << endl;
+    cout << "-------------------------------------------------------------------------" << endl;
+  }
 
-  if (rank == MASTER_NODE) SetTurboPerformance();
+  if (rank == MASTER_NODE && driver_config->GetBoolTurbomachinery()) SetTurboPerformance();
 
+  return (FinalTimeReached || MaxIterationsReached);
 }
 
 bool CMultizoneDriver::GetTimeConvergence() const{

SU2_CFD/src/drivers/CSinglezoneDriver.cpp

@@ -92,12 +92,6 @@ void CSinglezoneDriver::StartSolver() {
 
     Output(TimeIter);
 
-    /*--- Save iteration solution for libROM ---*/
-    if (config_container[MESH_0]->GetSave_libROM()) {
-      solver_container[ZONE_0][INST_0][MESH_0][FLOW_SOL]->SavelibROM(geometry_container[ZONE_0][INST_0][MESH_0],
-                                                                     config_container[ZONE_0], StopCalc);
-    }
-
     /*--- If the convergence criteria has been met, terminate the simulation. ---*/
 
     if (StopCalc) break;
@@ -208,7 +202,14 @@ void CSinglezoneDriver::Output(unsigned long TimeIter) {
                                                                solver_container[ZONE_0][INST_0][MESH_0],
                                                                TimeIter, StopCalc);
 
-  if (wrote_files){
+  /*--- Save iteration solution for libROM ---*/
+  if (config_container[MESH_0]->GetSave_libROM()) {
+    solver_container[ZONE_0][INST_0][MESH_0][FLOW_SOL]->SavelibROM(geometry_container[ZONE_0][INST_0][MESH_0],
+                                                                   config_container[ZONE_0], StopCalc);
+    wrote_files = true;
+  }
+
+  if (wrote_files) {
 
     StopTime = SU2_MPI::Wtime();
 

SU2_CFD/src/interfaces/fsi/CDiscAdjFlowTractionInterface.cpp

@@ -40,11 +40,7 @@ void CDiscAdjFlowTractionInterface::GetPhysical_Constants(CSolver *flow_solution
                                                           CGeometry *flow_geometry, CGeometry *struct_geometry,
                                                           const CConfig *flow_config, const CConfig *struct_config){
 
-  /*--- We have to clear the traction before applying it, because we are "adding" to node and not "setting" ---*/
-
-  struct_solution->GetNodes()->Clear_FlowTraction();
-
-  Preprocess(flow_config);
+  Preprocess(flow_config, struct_config, struct_geometry, struct_solution);
 
   if (!conservative) ComputeVertexAreas(struct_config, struct_geometry, struct_solution);
 

SU2_CFD/src/interfaces/fsi/CFlowTractionInterface.cpp

@@ -31,6 +31,7 @@
 #include "../../../../Common/include/geometry/CGeometry.hpp"
 #include "../../../include/solvers/CSolver.hpp"
 #include "../../../../Common/include/toolboxes/geometry_toolbox.hpp"
+#include "../../../Common/include/interface_interpolation/CInterpolator.hpp"
 #include <unordered_set>
 
 CFlowTractionInterface::CFlowTractionInterface(unsigned short val_nVar, unsigned short val_nConst,
@@ -39,7 +40,26 @@ CFlowTractionInterface::CFlowTractionInterface(unsigned short val_nVar, unsigned
   conservative(conservative_) {
 }
 
-void CFlowTractionInterface::Preprocess(const CConfig *flow_config) {
+void CFlowTractionInterface::Preprocess(const CConfig *flow_config, const CConfig *struct_config,
+                                        CGeometry *struct_geometry, CSolver *struct_solution) {
+
+  /*--- Clear the tractions only on the markers involved in interface, fluid tractions
+   * on other markers can be specified via e.g. the python wrapper. ---*/
+
+  for (auto iMarkerInt = 0u; iMarkerInt < struct_config->GetMarker_n_ZoneInterface()/2; iMarkerInt++) {
+    const auto markDonor = flow_config->FindInterfaceMarker(iMarkerInt);
+    const auto markTarget = struct_config->FindInterfaceMarker(iMarkerInt);
+
+    if(!CInterpolator::CheckInterfaceBoundary(markDonor, markTarget)) continue;
+    if (markTarget < 0) continue;
+
+    for (auto iVertex = 0ul; iVertex < struct_geometry->GetnVertex(markTarget); iVertex++) {
+      const auto iPoint = struct_geometry->vertex[markTarget][iVertex]->GetNode();
+      if (!struct_geometry->nodes->GetDomain(iPoint)) continue;
+      su2double zeros[3] = {};
+      struct_solution->GetNodes()->Set_FlowTraction(iPoint, zeros);
+    }
+  }
 
   /*--- Compute the constant factor to dimensionalize pressure and shear stress. ---*/
   const su2double *Velocity_ND, *Velocity_Real;
@@ -121,11 +141,7 @@ void CFlowTractionInterface::GetPhysical_Constants(CSolver *flow_solution, CSolv
                                                    CGeometry *flow_geometry, CGeometry *struct_geometry,
                                                    const CConfig *flow_config, const CConfig *struct_config) {
 
-  /*--- We have to clear the traction before applying it, because we are "adding" to node and not "setting" ---*/
-
-  struct_solution->GetNodes()->Clear_FlowTraction();
-
-  Preprocess(flow_config);
+  Preprocess(flow_config, struct_config, struct_geometry, struct_solution);
 
   if (!conservative) ComputeVertexAreas(struct_config, struct_geometry, struct_solution);