[WIP] Quasi-Newton correction for multizone couplings

SU2_CFD/include/interfaces/CInterface.hpp

@@ -29,6 +29,7 @@
 #pragma once
 
 #include "../../../Common/include/parallelization/mpi_structure.hpp"
+#include "../../../Common/include/toolboxes/CQuasiNewtonInvLeastSquares.hpp"
 
 #include <cmath>
 #include <string>
@@ -73,6 +74,9 @@
   unsigned short nVar = 0;
   static constexpr size_t MAXNDIM = 3;  /*!< \brief Max number of space dimensions, used in some static arrays. */
 
+  /*--- Quasi-Newton correction for each donor marker. ---*/
+  vector<CQuasiNewtonInvLeastSquares<passivedouble>> QuasiNewtonCorrection;
+
 public:
   /*!
    * \brief Constructor of the class.
@@ -91,6 +95,33 @@
    */
   virtual ~CInterface(void);
 
+  /*!
+   * \brief Initializes the Quasi-Newton correction object for each donor marker.
+   * \param[in] donor_geometry - Geometry of the donor mesh.
+   * \param[in] donor_config - Definition of the problem at the donor mesh, used to get the number of donor markers.
+   */
+  void InitializeQuasiNewtonCorrection(CGeometry *donor_geometry, const CConfig *donor_config);
+
+  /*!
+    * \brief Get the values to be transferred from the donor solution.
+    * \param[in] DonorData - Data container from which the values will be read.
+    * \param[in] TargetData - Data container to which the values will be written.
+    * \param[in] nVertex - Number of vertices on the donor side.
+    * \param[in] nVar - Number of variables.
+    */
+  template <class DonorData, class TargetData>
+  void SetAllValues(DonorData& donor, TargetData& target, unsigned long nVertex, unsigned short nVar) {
+
+    unsigned long index = 0;
+    for(auto iVertex = 0u; iVertex < nVertex; iVertex++) {
+      for (auto iVar = 0u; iVar < nVar; iVar++) {
+        // if (iVertex * nVar + iVar >= donor.size()) {}
+        target(iVertex, iVar) = donor(iVertex, iVar);
+        index++;
+      }
+    }
+  }
+
   /*!
    * \brief Interpolate data and broadcast it into all processors, for nonmatching meshes.
    * \param[in] interpolator - Object defining the interpolation.

SU2_CFD/src/drivers/CMultizoneDriver.cpp

@@ -291,6 +291,19 @@ void CMultizoneDriver::RunGaussSeidel() {
     config_container[iZone]->SetOuterIter(0ul);
   }
 
+  /*--- Initialize a Quasi-Newton correction for data that is broadcasted at physical interfaces. ---*/
+
+  for (iZone = 0; iZone < nZone; iZone++) {
+    for (auto jZone = 0u; jZone < nZone; jZone++) {
+      if (jZone != iZone) {
+        /*--- Donor zone is jzone ---*/
+        if(interface_container[jZone][iZone] != nullptr) {
+          interface_container[jZone][iZone]->InitializeQuasiNewtonCorrection(geometry_container[jZone][INST_0][MESH_0], config_container[jZone]);
+        }
+      }
+    }
+  }
+
   /*--- Loop over the number of outer iterations ---*/
   for (auto iOuter_Iter = 0ul; iOuter_Iter < driver_config->GetnOuter_Iter(); iOuter_Iter++) {
 

SU2_CFD/src/interfaces/CInterface.cpp

@@ -59,6 +59,37 @@ CInterface::~CInterface() {
   delete[] SpanLevelDonor;
 }
 
+void CInterface::InitializeQuasiNewtonCorrection(CGeometry *donor_geometry, const CConfig *donor_config) {
+
+  /*--- One Quasi-Newton correction object for each interface. ---*/
+  QuasiNewtonCorrection.resize(donor_config->GetMarker_n_ZoneInterface()/2);
+
+  /*--- Determine the global size of the interface. ---*/
+  for (auto iMarkerInt = 0u; iMarkerInt < donor_config->GetMarker_n_ZoneInterface()/2; iMarkerInt++) {
+    const auto markDonor = donor_config->FindInterfaceMarker(iMarkerInt);
+    int nLocalVertexDonor = 0;
+    if (markDonor >= 0) {
+      for (auto iVertex = 0ul; iVertex < donor_geometry->GetnVertex(markDonor); iVertex++) {
+        auto Point_Donor = donor_geometry->vertex[markDonor][iVertex]->GetNode();
+        /*--- Only domain points are donors. ---*/
+        nLocalVertexDonor += donor_geometry->nodes->GetDomain(Point_Donor);
+      }
+    }
+
+    vector<int> nAllVertexDonor(size);
+    SU2_MPI::Allgather(&nLocalVertexDonor, 1, MPI_INT, nAllVertexDonor.data(), 1, MPI_INT, SU2_MPI::GetComm());
+
+    int nGlobalVertexDonor = 0;
+    for (int i = 0; i < size; ++i) nGlobalVertexDonor += nAllVertexDonor[i];
+
+    /*--- Initialize the Quasi-Newton correction object for this interface. ---*/
+    if(donor_config->GetnQuasiNewtonSamples()) {
+      std::cout << "Applying a Quasi-Newton correction for zone " << donor_config->GetiZone() << " and interface " << iMarkerInt << " with " << nGlobalVertexDonor << " donor vertices and " << nVar << " variables." << std::endl;
+      QuasiNewtonCorrection[iMarkerInt].resize(donor_config->GetnQuasiNewtonSamples(), nGlobalVertexDonor, 4);
+    }
+  }
+}
+
 void CInterface::BroadcastData(const CInterpolator& interpolator,
                                CSolver *donor_solution, CSolver *target_solution,
                                CGeometry *donor_geometry, CGeometry *target_geometry,
@@ -160,6 +191,14 @@ void CInterface::BroadcastData(const CInterpolator& interpolator,
         swap(donorVar(i,iVar), donorVar(j,iVar));
     }
 
+    /*--- Each rank performs a Quasi-Newton correction step on the global donor data before applying them on the target. ---*/
+
+    if (QuasiNewtonCorrection[iMarkerInt].size()) {
+      SetAllValues(donorVar, QuasiNewtonCorrection[iMarkerInt].FPresult(), nGlobalVertexDonor, nVar);
+      QuasiNewtonCorrection[iMarkerInt].compute();
+      SetAllValues(donorVar, QuasiNewtonCorrection[iMarkerInt], nGlobalVertexDonor, nVar);
+    }
+
     /*--- Loop over target vertices. ---*/
 
     for (auto iVertex = 0ul; iVertex < target_geometry->GetnVertex(markTarget); iVertex++) {