[WIP] Multizone adjoints for turbomachinery

.github/workflows/regression.yml

@@ -211,7 +211,7 @@ jobs:
         uses: docker://ghcr.io/su2code/su2/test-su2:250717-1402
         with:
           # -t <Tutorials-branch> -c <Testcases-branch>
-          args: -b ${{github.ref}} -t develop -c develop -s ${{matrix.testscript}}
+          args: -b ${{github.ref}} -t develop -c feature_mz_turbomachinery_adjoint -s ${{matrix.testscript}}
       - name: Cleanup
         uses: docker://ghcr.io/su2code/su2/test-su2:250717-1402
         with:

Common/include/CConfig.hpp

@@ -1973,12 +1973,6 @@ class CConfig {
    */
   su2double GetPressure_FreeStreamND(void) const { return Pressure_FreeStreamND; }
 
-  /*!
-   * \brief Get a reference to the non-dimensionalized freestream pressure (used for AD tracking).
-   * \return Reference to non-dimensionalized freestream pressure.
-   */
-  su2double& GetPressure_FreeStreamND(void) { return Pressure_FreeStreamND; }
-
   /*!
    * \brief Get the value of the thermodynamic pressure.
    * \return Thermodynamic pressure.
@@ -2004,12 +1998,6 @@ class CConfig {
    */
   su2double GetTemperature_FreeStreamND(void) const { return Temperature_FreeStreamND; }
 
-  /*!
-   * \brief Get a reference to the non-dimensionalized freestream temperature (used for AD tracking).
-   * \return Reference to non-dimensionalized freestream temperature.
-   */
-  su2double& GetTemperature_FreeStreamND(void) { return Temperature_FreeStreamND; }
-
   /*!
    * \brief Get the value of the non-dimensionalized vibrational-electronic freestream temperature.
    * \return Non-dimensionalized vibrational-electronic freestream temperature.
@@ -10122,6 +10110,13 @@ class CConfig {
    */
   short FindInterfaceMarker(unsigned short iInterface) const;
 
+  /*!
+   * \brief Find the marker index (if any) that is part of a mixing plane interface pair.
+   * \param[in] nMarker - Number of the marker in a zone being tested, starting at 0.
+   * \return value > 1 if (on this mpi rank) the zone defined by config is part of the mixing plane.
+   */
+  short FindMixingPlaneInterfaceMarker(unsigned short nMarker, unsigned short iMarkerInt) const;
+
   /*!
    * \brief Get whether or not to save solution data to libROM.
    * \return True if specified in config file.

Common/include/geometry/CGeometry.hpp

@@ -777,7 +777,7 @@ class CGeometry {
   inline virtual void GatherInOutAverageValues(CConfig* config, bool allocate) {}
 
   /*!
-   * \brief Store all the turboperformance in the solver in ZONE_0.
+   * \brief Store all the turboperformance in the solver in final zone.
    * \param[in] donor_geometry  - Solution from the donor mesh.
    * \param[in] target_geometry - Solution from the target mesh.
    * \param[in] donorZone       - counter of the donor solution
@@ -1862,7 +1862,8 @@ class CGeometry {
    * \param[in] config_container - Definition of the particular problem.
    * \param[in] geometry_container - Geometrical definition of the problem.
    */
-  static void ComputeWallDistance(const CConfig* const* config_container, CGeometry**** geometry_container);
+  static void ComputeWallDistance(const CConfig* const* config_container, CGeometry**** geometry_container,
+                                  const int record_zone = -1);
 
   /*!
    * \brief Set the amount of nonconvex elements in the mesh.

Common/include/interface_interpolation/CInterpolator.hpp

@@ -99,7 +99,13 @@ class CInterpolator {
       coefficient.resize(nDonor);
     }
   };
-  vector<vector<CDonorInfo> > targetVertices; /*! \brief Donor information per marker per vertex of the target. */
+  vector<vector<CDonorInfo>> targetVertices; /*! \brief Donor information per marker per vertex of the target. */
+
+  struct CSpanDonorInfo {
+    size_t donorSpan;       // Refers to donor span
+    su2double coefficient;  // Refers to coefficient
+  };
+  vector<vector<CSpanDonorInfo>> targetSpans;  // <iMarkerInt<TargetSpan<CSpanDonorInfo>>
 
   /*!
    * \brief Constructor of the class.
@@ -125,13 +131,18 @@ class CInterpolator {
    * \note Main method that derived classes must implement.
    * \param[in] config - Definition of the particular problem.
    */
-  virtual void SetTransferCoeff(const CConfig* const* config) = 0;
+  virtual void SetTransferCoeff(CGeometry**** geometry, const CConfig* const* config) = 0;
 
   /*!
    * \brief Print information about the interpolation.
    */
   virtual void PrintStatistics(void) const {}
 
+  /*!
+   * \brief Write mixing plane interpolation details to file
+   */
+  inline virtual void WriteInterpolationDetails(const string& filename, const CConfig* const* config){};
+
   /*!
    * \brief Check whether an interface should be processed or not, i.e. if it is part of the zones.
    * \param[in] val_markDonor  - Marker tag from donor zone.

Common/include/interface_interpolation/CInterpolatorFactory.hpp

@@ -43,5 +43,5 @@ namespace CInterpolatorFactory {
  */
 CInterpolator* CreateInterpolator(CGeometry**** geometry_container, const CConfig* const* config,
                                   const CInterpolator* transpInterpolator, unsigned iZone, unsigned jZone,
-                                  bool verbose = true);
+                                  bool mixing_plane, bool verbose = true);
 }  // namespace CInterpolatorFactory

Common/include/interface_interpolation/CIsoparametric.hpp

@@ -64,7 +64,7 @@ class CIsoparametric final : public CInterpolator {
    * \brief Set up transfer matrix defining relation between two meshes
    * \param[in] config - Definition of the particular problem.
    */
-  void SetTransferCoeff(const CConfig* const* config) override;
+  void SetTransferCoeff(CGeometry**** geometry, const CConfig* const* config) override;
 
   /*!
    * \brief Print information about the interpolation.

Common/include/interface_interpolation/CMirror.hpp

@@ -54,5 +54,5 @@ class CMirror final : public CInterpolator {
    * \brief Set up transfer matrix defining relation between two meshes
    * \param[in] config - Definition of the particular problem.
    */
-  void SetTransferCoeff(const CConfig* const* config) override;
+  void SetTransferCoeff(CGeometry**** geometry, const CConfig* const* config) override;
 };

Common/include/interface_interpolation/CMixingPlane.hpp

@@ -0,0 +1,101 @@
+/*!
+ * \file CMixingPlane.hpp
+ * \brief Header of mixing plane interpolation methods.
+ * \author J. Kelly
+ * \version 8.3.0 "Harrier"
+ *
+ * SU2 Project Website: https://su2code.github.io
+ *
+ * The SU2 Project is maintained by the SU2 Foundation
+ * (http://su2foundation.org)
+ *
+ * Copyright 2012-2025, SU2 Contributors (cf. AUTHORS.md)
+ *
+ * SU2 is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public
+ * License as published by the Free Software Foundation; either
+ * version 2.1 of the License, or (at your option) any later version.
+ *
+ * SU2 is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with SU2. If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#pragma once
+#include "CInterpolator.hpp"
+#include "../option_structure.hpp"
+
+/*!
+ * \brief Mixing plane interpolation.
+ * \note This contains several interpolation methods used in the mixing plane interpolation
+ * and enables the mixing state class structure for proper recording in AD mode
+ * \ingroup Interfaces
+ */
+class CMixingPlane final : public CInterpolator {
+ public:
+  CMixingPlane(CGeometry**** geometry_container, const CConfig* const* config, unsigned int iZone, unsigned int jZone);
+
+  void SetTransferCoeff(CGeometry**** geometry, const CConfig* const* config) override;
+
+  inline CSpanDonorInfo MapMatchingSpan(unsigned short iSpanTarget) { return {iSpanTarget, 0.0}; }
+
+  inline CSpanDonorInfo MapNearestSpan(const su2double iSpanTargetValue, const su2double* spanValuesDonor,
+                                       unsigned long nSpanDonor) {
+    unsigned short tSpan = 0;  // Nearest donor span index
+    auto minDist = std::numeric_limits<su2double>::max();
+
+    for (auto iSpanDonor = 0u; iSpanDonor < nSpanDonor - 1; iSpanDonor++) {
+      const auto dist = abs(iSpanTargetValue - spanValuesDonor[iSpanDonor]);
+      if (dist < minDist) {
+        minDist = dist;
+        tSpan = iSpanDonor;
+      }
+    }
+    return {tSpan, 0.0};
+  };
+
+  inline CSpanDonorInfo MapLinearInterpolationSpan(const su2double iSpanTargetValue, const su2double* spanValuesDonor,
+                                                   unsigned long nSpanDonor, int rank) {
+    unsigned short kSpan = 0;  // Lower bound donor span for interpolation
+    auto minDist = std::numeric_limits<su2double>::max();
+    su2double coeff = 0.0;  // Interpolation coefficient
+
+    if (iSpanTargetValue <= spanValuesDonor[0]) {
+      PrintClampingWarning(rank, true);
+      return {0, 0.0};
+    }
+
+    if (iSpanTargetValue >= spanValuesDonor[nSpanDonor - 1]) {
+      PrintClampingWarning(rank, false);
+      return {nSpanDonor - 1, 0.0};
+    }
+
+    for (auto iSpanDonor = 0u; iSpanDonor < nSpanDonor - 1; iSpanDonor++) {
+      const auto dist = abs(iSpanTargetValue - spanValuesDonor[iSpanDonor]);
+      if (dist < minDist && iSpanTargetValue >= spanValuesDonor[iSpanDonor]) {
+        kSpan = iSpanDonor;
+        minDist = dist;
+      }
+    }
+    coeff = (iSpanTargetValue - spanValuesDonor[kSpan]) / (spanValuesDonor[kSpan + 1] - spanValuesDonor[kSpan]);
+    return {kSpan, coeff};
+  };
+
+  inline void PrintClampingWarning(int rank, bool atHub) {
+    if (rank != MASTER_NODE) return;
+    cout << "Warning! Target spans exist outside the bounds of donor spans! Clamping interpolator..." << endl;
+    cout << (atHub ? "This is an issue at the hub." : "This is an issue at the shroud.") << endl;
+    cout << "Setting coeff = 0.0 and transferring endwall value!" << endl;
+  };
+
+  /*!
+   * \brief Write interpolation details to file.
+   * \param[in] filename - Name of output file.
+   * \param[in] config - Configuration for all zones.
+   */
+  void WriteInterpolationDetails(const string& filename, const CConfig* const* config) override;
+};

Common/include/interface_interpolation/CNearestNeighbor.hpp

@@ -63,7 +63,7 @@ class CNearestNeighbor final : public CInterpolator {
    * \brief Set up transfer matrix defining relation between two meshes.
    * \param[in] config - Definition of the particular problem.
    */
-  void SetTransferCoeff(const CConfig* const* config) override;
+  void SetTransferCoeff(CGeometry**** geometry, const CConfig* const* config) override;
 
   /*!
    * \brief Print interpolation statistics.

Common/include/interface_interpolation/CRadialBasisFunction.hpp

@@ -56,7 +56,7 @@ class CRadialBasisFunction final : public CInterpolator {
    * \brief Set up transfer matrix defining relation between two meshes
    * \param[in] config - Definition of the particular problem.
    */
-  void SetTransferCoeff(const CConfig* const* config) override;
+  void SetTransferCoeff(CGeometry**** geometry, const CConfig* const* config) override;
 
   /*!
    * \brief Print information about the interpolation.

Common/include/interface_interpolation/CSlidingMesh.hpp

@@ -49,7 +49,7 @@ class CSlidingMesh final : public CInterpolator {
    * \brief Set up transfer matrix defining relation between two meshes
    * \param[in] config - Definition of the particular problem.
    */
-  void SetTransferCoeff(const CConfig* const* config) override;
+  void SetTransferCoeff(CGeometry**** geometry, const CConfig* const* config) override;
 
  private:
   /*!

Common/include/option_structure.hpp

@@ -2100,6 +2100,9 @@ enum ENUM_OBJECTIVE {
   TOPOL_DISCRETENESS = 63,      /*!< \brief Measure of the discreteness of the current topology. */
   TOPOL_COMPLIANCE = 64,        /*!< \brief Measure of the discreteness of the current topology. */
   STRESS_PENALTY = 65,          /*!< \brief Penalty function of VM stresses above a maximum value. */
+  ENTROPY_GENERATION = 80,      /*!< \brief Entropy generation turbomachinery objective function. */
+  TOTAL_PRESSURE_LOSS = 81,     /*!< \brief Total pressure loss turbomachinery objective function. */
+  KINETIC_ENERGY_LOSS = 82      /*!< \breif Kinetic energy loss coefficient turbomachinery objective function. */
 };
 static const MapType<std::string, ENUM_OBJECTIVE> Objective_Map = {
   MakePair("DRAG", DRAG_COEFFICIENT)
@@ -2142,6 +2145,9 @@ static const MapType<std::string, ENUM_OBJECTIVE> Objective_Map = {
   MakePair("TOPOL_DISCRETENESS", TOPOL_DISCRETENESS)
   MakePair("TOPOL_COMPLIANCE", TOPOL_COMPLIANCE)
   MakePair("STRESS_PENALTY", STRESS_PENALTY)
+  MakePair("ENTROPY_GENERATION", ENTROPY_GENERATION)
+  MakePair("TOTAL_PRESSURE_LOSS", TOTAL_PRESSURE_LOSS)
+  MakePair("KINETIC_ENERGY_LOSS", KINETIC_ENERGY_LOSS)
 };
 
 /*!
@@ -2616,19 +2622,15 @@ enum class CHECK_TAPE_VARIABLES {
 };
 static const MapType<std::string, CHECK_TAPE_VARIABLES> CheckTapeVariables_Map = {
     MakePair("SOLVER_VARIABLES", CHECK_TAPE_VARIABLES::SOLVER_VARIABLES)
-    MakePair("SOLVER_VARIABLES_AND_MESH_COORDINATES", CHECK_TAPE_VARIABLES::MESH_COORDINATES)
+    MakePair("MESH_COORDINATES", CHECK_TAPE_VARIABLES::MESH_COORDINATES)
 };
 
 enum class RECORDING {
   CLEAR_INDICES,
   SOLUTION_VARIABLES,
   MESH_COORDS,
   MESH_DEFORM,
-  SOLUTION_AND_MESH,
-  TAG_INIT_SOLVER_VARIABLES,
-  TAG_CHECK_SOLVER_VARIABLES,
-  TAG_INIT_SOLVER_AND_MESH,
-  TAG_CHECK_SOLVER_AND_MESH
+  SOLUTION_AND_MESH
 };
 
 /*!

Common/src/CConfig.cpp

@@ -3514,12 +3514,7 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i
 
     /*---  Using default frequency of 250 for all files when steady, and 1 for unsteady. ---*/
     for (auto iVolumeFreq = 0; iVolumeFreq < nVolumeOutputFrequencies; iVolumeFreq++){
-      if (Multizone_Problem && DiscreteAdjoint) {
-        VolumeOutputFrequencies[iVolumeFreq] = Time_Domain ? 1 : nOuterIter;
-      }
-      else {
-        VolumeOutputFrequencies[iVolumeFreq] = Time_Domain ? 1 : 250;
-      }
+      VolumeOutputFrequencies[iVolumeFreq] = Time_Domain ? 1 : nOuterIter;
     }
   } else if (nVolumeOutputFrequencies < nVolumeOutputFiles) {
     /*--- If there are fewer frequencies than files, repeat the last frequency.
@@ -7002,6 +6997,7 @@ void CConfig::SetOutput(SU2_COMPONENT val_software, unsigned short val_izone) {
         case TOPOL_DISCRETENESS:         cout << "Topology discreteness objective function." << endl; break;
         case TOPOL_COMPLIANCE:           cout << "Topology compliance objective function." << endl; break;
         case STRESS_PENALTY:             cout << "Stress penalty objective function." << endl; break;
+        case ENTROPY_GENERATION:         cout << "Entropy generation objective function." << endl; break;
       }
     }
     else {
@@ -8595,6 +8591,7 @@ CConfig::~CConfig() {
 
   delete [] nBlades;
   delete [] FreeStreamTurboNormal;
+
 }
 
 /*--- Input is the filename base, output is the completed filename. ---*/
@@ -8738,6 +8735,9 @@ string CConfig::GetObjFunc_Extension(string val_filename) const {
         case TOPOL_DISCRETENESS:          AdjExt = "_topdisc";  break;
         case TOPOL_COMPLIANCE:            AdjExt = "_topcomp";  break;
         case STRESS_PENALTY:              AdjExt = "_stress";   break;
+        case ENTROPY_GENERATION:          AdjExt = "_entg";     break;
+        case TOTAL_PRESSURE_LOSS:         AdjExt = "_tot_press_loss"; break;
+        case KINETIC_ENERGY_LOSS:         AdjExt = "_kin_en_loss"; break;
       }
     }
     else{
@@ -10038,6 +10038,23 @@ short CConfig::FindInterfaceMarker(unsigned short iInterface) const {
   return -1;
 }
 
+short CConfig::FindMixingPlaneInterfaceMarker(unsigned short nMarker, unsigned short iMarkerInt) const {
+  short mark;
+  for (auto iMarker = 0; iMarker < nMarker; iMarker++){
+      /*--- If the tag GetMarker_All_MixingPlaneInterface equals the index we are looping at ---*/
+      if (GetMarker_All_MixingPlaneInterface(iMarker) == iMarkerInt){
+          /*--- We have identified the local index of the marker ---*/
+          /*--- Store the identifier for the marker ---*/
+          mark = iMarker;
+          /*--- Exit the for loop: we have found the local index for Mixing-Plane interface ---*/
+          return mark;
+      }
+      /*--- If the tag hasn't matched any tag within the donor markers ---*/
+      mark = -1;
+  }
+  return mark;
+}
+
 void CConfig::Tick(double *val_start_time) {
 
 #ifdef PROFILE

Common/src/geometry/CGeometry.cpp

@@ -4153,7 +4153,8 @@ su2double NearestNeighborDistance(CGeometry* geometry, const CConfig* config, co
 }
 }  // namespace
 
-void CGeometry::ComputeWallDistance(const CConfig* const* config_container, CGeometry**** geometry_container) {
+void CGeometry::ComputeWallDistance(const CConfig* const* config_container, CGeometry**** geometry_container,
+                                    const int record_zone) {
   int nZone = config_container[ZONE_0]->GetnZone();
   bool allEmpty = true;
   vector<bool> wallDistanceNeeded(nZone, false);
@@ -4221,6 +4222,11 @@ void CGeometry::ComputeWallDistance(const CConfig* const* config_container, CGeo
     }
 
     for (int iZone = 0; iZone < nZone; iZone++) {
+      /*--- When recording for a specific zone, only compute nearest-neighbor distances (which read vertex
+       *    normals) for that zone. Reading normals from other zones at this tape position would create
+       *    cross-zone AD dependencies before those zones have had their geometry updated. ---*/
+      if (record_zone >= 0 && iZone != record_zone) continue;
+
       /*--- 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) {