[WIP] Multizone adjoints for turbomachinery

Common/include/CConfig.hpp

@@ -449,6 +449,11 @@ class CConfig {
   TURBOMACHINERY_TYPE *Kind_TurboMachinery;
   su2vector<TURBO_INTERFACE_KIND> Kind_TurboInterface;
 
+    /* Turbomachinery objective functions */
+  su2double *EntropyGeneration;
+  su2double *TotalPressureLoss;
+  su2double *KineticEnergyLoss;
+
   /* Gradient smoothing options */
   su2double SmoothingEps1;          /*!< \brief Parameter for the identity part in gradient smoothing. */
   su2double SmoothingEps2;          /*!< \brief Parameter for the Laplace part in gradient smoothing. */
@@ -8147,6 +8152,27 @@ class CConfig {
    * \param[in] val_surface_species_variance - Value of the species variance.
    */
   void SetSurface_Species_Variance(unsigned short val_marker, su2double val_surface_species_variance) { Surface_Species_Variance[val_marker] = val_surface_species_variance; }
+
+  /*!
+   * \brief Set entropy generation for a turbomachinery zone
+   * \param[in] val_iZone - zone index
+   * \param[in] val_entropy_generation - value of entropy generation
+  */
+  void SetEntropyGeneration(unsigned short val_iZone, su2double val_entropy_generation) { EntropyGeneration[val_iZone] = val_entropy_generation; }
+
+  /*!
+   * \brief Get total pressure loss for a turbomachinery zone
+   * \param[in] val_iZone - zone index
+   * \param[in] val_total_pressure_loss - value of total pressure loss
+   */
+  void SetTotalPressureLoss(unsigned short val_iZone, su2double val_total_pressure_loss) { TotalPressureLoss[val_iZone] = val_total_pressure_loss; }
+
+  /*!
+   * \brief Get kinetic energy loss for a turbomachinery zone
+   * \param[in] val_iZone - zone index
+   * \param[in] val_kinetic_energy_loss - value of kinetic energy loss
+   */
+  void SetKineticEnergyLoss(unsigned short val_iZone, su2double val_kinetic_energy_loss) { KineticEnergyLoss[val_iZone] = val_kinetic_energy_loss; }
 
   /*!
    * \brief Get the back pressure (static) at an outlet boundary.
@@ -10009,6 +10035,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) const;
+
   /*!
    * \brief Get whether or not to save solution data to libROM.
    * \return True if specified in config file.

Common/include/geometry/CGeometry.hpp

@@ -776,7 +776,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

Common/include/interface_interpolation/CInterpolator.hpp

@@ -101,6 +101,12 @@ class CInterpolator {
   };
   vector<vector<CDonorInfo> > targetVertices; /*! \brief Donor information per marker per vertex of the target. */
 
+  struct CSpanDonorInfo {
+      unsigned long donorSpan; // Refers to donor span
+      su2double coefficient; // Refers to coefficient
+  };
+  vector<vector<CSpanDonorInfo>> targetSpans; // <iMarkerInt<TargetSpan<CSpanDonorInfo>>
+
   /*!
    * \brief Constructor of the class.
    * \param[in] geometry_container - Geometrical definition of the problem.

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/CMixingPlane.hpp

@@ -0,0 +1,43 @@
+/*!
+ * \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"
+
+/*!
+ * \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(const CConfig* const* config) override;
+};

Common/include/option_structure.hpp

@@ -2067,6 +2067,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)
@@ -2109,6 +2112,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)
 };
 
 /*!

Common/src/CConfig.cpp

@@ -1049,6 +1049,11 @@ void CConfig::SetPointersNull() {
   nBlades                      = nullptr;
   FreeStreamTurboNormal        = nullptr;
 
+  /*--- Turbomachinery Objective Functions ---*/
+  EntropyGeneration = nullptr;
+  TotalPressureLoss = nullptr;
+  KineticEnergyLoss = nullptr;
+
   top_optim_kernels       = nullptr;
   top_optim_kernel_params = nullptr;
   top_optim_filter_radius = nullptr;
@@ -3479,6 +3484,10 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i
         VolumeOutputFrequencies[iVolumeFreq] = Time_Domain ? 1 : 250;
       }
     }
+  } else if (Multizone_Problem && DiscreteAdjoint) {
+      SU2_MPI::Error(string("OUTPUT_WRT_FREQ cannot be specified for this solver "
+                            "(writing of restart and sensitivity files not possible for multizone discrete adjoint during runtime yet).\n"
+                            "Please remove this option from the config file, output files will be written when solver finalizes.\n"), CURRENT_FUNCTION);
   } else if (nVolumeOutputFrequencies < nVolumeOutputFiles) {
     /*--- If there are fewer frequencies than files, repeat the last frequency.
      *    This is useful to define 1 frequency for the restart file and 1 frequency for all the visualization files.  ---*/
@@ -6172,6 +6181,11 @@ void CConfig::SetMarkers(SU2_COMPONENT val_software) {
         Marker_CfgFile_ZoneInterface[iMarker_CfgFile] = YES;
   }
 
+  /*--- Allocate memory for turbomachinery objective functions ---*/
+  EntropyGeneration = new su2double[nZone] ();
+  TotalPressureLoss = new su2double[nZone] ();
+  KineticEnergyLoss = new su2double[nZone] ();
+
   /*--- Identification of Turbomachinery markers and flag them---*/
 
   for (iMarker_CfgFile = 0; iMarker_CfgFile < nMarker_CfgFile; iMarker_CfgFile++) {
@@ -6956,6 +6970,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 {
@@ -8521,6 +8536,10 @@ CConfig::~CConfig() {
 
   delete [] nBlades;
   delete [] FreeStreamTurboNormal;
+
+  delete [] EntropyGeneration;
+  delete [] TotalPressureLoss;
+  delete [] KineticEnergyLoss;
 }
 
 /*--- Input is the filename base, output is the completed filename. ---*/
@@ -8664,6 +8683,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{
@@ -9942,6 +9964,23 @@ short CConfig::FindInterfaceMarker(unsigned short iInterface) const {
   return -1;
 }
 
+short CConfig::FindMixingPlaneInterfaceMarker(unsigned short nMarker) 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)){
+          /*--- 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/CPhysicalGeometry.cpp

@@ -50,6 +50,8 @@
 #include "../../include/geometry/primal_grid/CPrism.hpp"
 #include "../../include/geometry/primal_grid/CVertexMPI.hpp"
 
+#include "../../../Common/include/tracy_structure.hpp"
+
 #include <sys/types.h>
 #include <sys/stat.h>
 #include <iterator>
@@ -6326,6 +6328,7 @@ void CPhysicalGeometry::GatherInOutAverageValues(CConfig* config, bool allocate)
 
 void CPhysicalGeometry::SetAvgTurboGeoValues(const CConfig* donor_config, CGeometry* donor_geometry,
                                              unsigned short donorZone) {
+
   unsigned short iSpan;
   unsigned short nSpanMaxAllZones = donor_config->GetnSpanMaxAllZones();
 

Common/src/interface_interpolation/CInterpolatorFactory.cpp

@@ -32,11 +32,12 @@
 #include "../../include/interface_interpolation/CNearestNeighbor.hpp"
 #include "../../include/interface_interpolation/CRadialBasisFunction.hpp"
 #include "../../include/interface_interpolation/CSlidingMesh.hpp"
+#include "../../include/interface_interpolation/CMixingPlane.hpp"
 
 namespace CInterpolatorFactory {
 CInterpolator* CreateInterpolator(CGeometry**** geometry_container, const CConfig* const* config,
                                   const CInterpolator* transpInterpolator, unsigned iZone, unsigned jZone,
-                                  bool verbose) {
+                                  bool mixing_plane, bool verbose) {
   CInterpolator* interpolator = nullptr;
 
   /*--- Only print information on master node. ---*/
@@ -47,36 +48,41 @@ CInterpolator* CreateInterpolator(CGeometry**** geometry_container, const CConfi
 
   if (verbose) cout << " Setting coupling ";
 
-  /*--- Conservative interpolation is not applicable to the sliding
-   *    mesh approach so that case is handled first. Then we either
-   *    return a CMirror if the target requires conservative inter-
-   *    polation, or the type of interpolator defined by "type". ---*/
+  if (mixing_plane) {
+    if (verbose) cout << "using a mixing plane interpolation." << endl;
+    interpolator = new CMixingPlane(geometry_container, config, iZone, jZone);
+  } else { // Really awful thing to do
+    /*--- Conservative interpolation is not applicable to the sliding
+    *    mesh approach so that case is handled first. Then we either
+    *    return a CMirror if the target requires conservative inter-
+    *    polation, or the type of interpolator defined by "type". ---*/
 
-  if (type == INTERFACE_INTERPOLATOR::WEIGHTED_AVERAGE) {
-    if (verbose) cout << "using a sliding mesh approach." << endl;
-    interpolator = new CSlidingMesh(geometry_container, config, iZone, jZone);
-  } else if (config[jZone]->GetConservativeInterpolation()) {
-    if (verbose) cout << "using the mirror approach, \"transposing\" coefficients from opposite mesh." << endl;
-    interpolator = new CMirror(geometry_container, config, transpInterpolator, iZone, jZone);
-  } else {
-    switch (type) {
-      case INTERFACE_INTERPOLATOR::ISOPARAMETRIC:
-        if (verbose) cout << "using the isoparametric approach." << endl;
-        interpolator = new CIsoparametric(geometry_container, config, iZone, jZone);
-        break;
+    if (type == INTERFACE_INTERPOLATOR::WEIGHTED_AVERAGE) {
+      if (verbose) cout << "using a sliding mesh approach." << endl;
+      interpolator = new CSlidingMesh(geometry_container, config, iZone, jZone);
+    } else if (config[jZone]->GetConservativeInterpolation()) {
+      if (verbose) cout << "using the mirror approach, \"transposing\" coefficients from opposite mesh." << endl;
+      interpolator = new CMirror(geometry_container, config, transpInterpolator, iZone, jZone);
+    } else {
+      switch (type) {
+        case INTERFACE_INTERPOLATOR::ISOPARAMETRIC:
+          if (verbose) cout << "using the isoparametric approach." << endl;
+          interpolator = new CIsoparametric(geometry_container, config, iZone, jZone);
+          break;
 
-      case INTERFACE_INTERPOLATOR::NEAREST_NEIGHBOR:
-        if (verbose) cout << "using a nearest neighbor approach." << endl;
-        interpolator = new CNearestNeighbor(geometry_container, config, iZone, jZone);
-        break;
+        case INTERFACE_INTERPOLATOR::NEAREST_NEIGHBOR:
+          if (verbose) cout << "using a nearest neighbor approach." << endl;
+          interpolator = new CNearestNeighbor(geometry_container, config, iZone, jZone);
+          break;
 
-      case INTERFACE_INTERPOLATOR::RADIAL_BASIS_FUNCTION:
-        if (verbose) cout << "using a radial basis function approach." << endl;
-        interpolator = new CRadialBasisFunction(geometry_container, config, iZone, jZone);
-        break;
+        case INTERFACE_INTERPOLATOR::RADIAL_BASIS_FUNCTION:
+          if (verbose) cout << "using a radial basis function approach." << endl;
+          interpolator = new CRadialBasisFunction(geometry_container, config, iZone, jZone);
+          break;
 
-      default:
-        SU2_MPI::Error("Unknown type of interpolation.", CURRENT_FUNCTION);
+        default:
+          SU2_MPI::Error("Unknown type of interpolation.", CURRENT_FUNCTION);
+      }
     }
   }