[WIP] Multizone adjoints for turbomachinery

Common/include/CConfig.hpp

@@ -445,6 +445,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. */
@@ -7990,6 +7995,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.
@@ -8271,6 +8297,24 @@ class CConfig {
    */
   su2double GetSurface_Species_Variance(unsigned short val_marker) const { return Surface_Species_Variance[val_marker]; }
 
+    /*!
+   * \brief Get entropy generation for a turbomachine at a boundary
+   * \param[in] val_iZone - zone index
+   */
+  su2double GetEntropyGeneration(unsigned short val_iZone) const { return EntropyGeneration[val_iZone]; }
+
+  /*!
+   * \brief Get total pressure loss for a turbomachinery zone
+   * \param[in] val_iZone - zone index
+   */
+  su2double GetTotalPressureLoss(unsigned short val_iZone) const { return TotalPressureLoss[val_iZone]; }
+
+  /*!
+   * \brief Get kinetic energy loss for a turbomachinery zone
+   * \param[in] val_iZone - zone index
+   */
+  su2double GetKineticEnergyLoss(unsigned short val_iZone) const { return KineticEnergyLoss[val_iZone]; }
+
   /*!
    * \brief Get the back pressure (static) at an outlet boundary.
    * \param[in] val_index - Index corresponding to the outlet boundary.

Common/include/geometry/CGeometry.hpp

@@ -62,6 +62,7 @@ extern "C" {
 #include "../CConfig.hpp"
 #include "../toolboxes/graph_toolbox.hpp"
 #include "../adt/CADTElemClass.hpp"
+#include "../../../SU2_CFD/include/interfaces/CInterface.hpp"
 
 using namespace std;
 
@@ -776,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
@@ -1352,6 +1353,10 @@ class CGeometry {
    */
   static void UpdateGeometry(CGeometry** geometry_container, CConfig* config);
 
+  static void UpdateTurboGeometry(CGeometry**** geometry, CInterface*** interface, CConfig** config_container);
+
+  static void InitTurboVertexAdj(CGeometry**** geometry, CConfig** config);
+
   /*!
    * \brief Update the multi-grid structure for the customized boundary conditions
    * \param geometry_container - Geometrical definition.

Common/include/option_structure.hpp

@@ -2010,6 +2010,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)
@@ -2052,6 +2055,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

@@ -1048,6 +1048,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;
@@ -3413,8 +3418,17 @@ 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++){
-      VolumeOutputFrequencies[iVolumeFreq] = Time_Domain ? 1 : 250;
+      if (Multizone_Problem && DiscreteAdjoint) {
+        VolumeOutputFrequencies[iVolumeFreq] = nOuterIter;
+      }
+      else {
+        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.  ---*/
@@ -6056,6 +6070,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++) {
@@ -8364,6 +8383,10 @@ CConfig::~CConfig() {
 
   delete [] nBlades;
   delete [] FreeStreamTurboNormal;
+
+  delete [] EntropyGeneration;
+  delete [] TotalPressureLoss;
+  delete [] KineticEnergyLoss;
 }
 
 string CConfig::GetFilename(string filename, const string& ext, int timeIter) const {
@@ -8541,6 +8564,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{

Common/src/geometry/CGeometry.cpp

@@ -2346,6 +2346,62 @@ void CGeometry::UpdateGeometry(CGeometry** geometry_container, CConfig* config)
   geometry_container[MESH_0]->ComputeSurfaceAreaCfgFile(config);
 }
 
+void CGeometry::InitTurboVertexAdj(CGeometry**** geometry, CConfig** config){
+  auto nSpanMax = 0u; 
+  auto nZone = config[ZONE_0]->GetnZone();
+
+  /*--- Create turbovertex ---*/
+  for (auto iZone = 0u; iZone < nZone; iZone++){
+    geometry[iZone][INST_0][MESH_0]->ComputeNSpan(config[iZone], iZone, INFLOW, true);
+    geometry[iZone][INST_0][MESH_0]->ComputeNSpan(config[iZone], iZone, OUTFLOW, true);
+    if (config[iZone]->GetnSpanWiseSections() > nSpanMax){
+      nSpanMax = config[iZone]->GetnSpanWiseSections();
+    }
+
+    config[nZone-1]->SetnSpan_iZones(config[iZone]->GetnSpanWiseSections(), iZone);
+
+    geometry[iZone][INST_0][MESH_0]->SetTurboVertex(config[iZone], iZone, INFLOW, true);
+    geometry[iZone][INST_0][MESH_0]->SetTurboVertex(config[iZone], iZone, OUTFLOW, true);
+  }
+
+  /*--- Set max span in all zones ---*/
+  for (auto iZone = 0u; iZone < nZone; iZone++){
+    if (config[iZone]->GetBoolTurbomachinery()) {
+      config[iZone]->SetnSpanMaxAllZones(nSpanMax);
+    }
+  }
+}
+
+void CGeometry::UpdateTurboGeometry(CGeometry**** geometry, CInterface*** interface, CConfig** config){
+  auto nZone = config[ZONE_0]->GetnZone();
+
+  for (auto iZone = 0u; iZone < nZone; iZone++) {
+    geometry[iZone][INST_0][MESH_0]->SetAvgTurboValue(config[iZone], iZone, INFLOW, true);
+    geometry[iZone][INST_0][MESH_0]->SetAvgTurboValue(config[iZone],iZone, OUTFLOW, true);
+    geometry[iZone][INST_0][MESH_0]->GatherInOutAverageValues(config[iZone], true);
+  }
+
+  if(config[ZONE_0]->GetBoolMixingPlaneInterface()){
+    for (auto donorZone = 0u; donorZone < nZone; donorZone++) {
+      for (auto targetZone = 0u; targetZone < nZone; targetZone++) {
+        if (interface[donorZone][targetZone] != nullptr){
+          interface[donorZone][targetZone]->SetSpanWiseLevels(config[donorZone], config[targetZone]);
+        }
+      }
+    }
+  }
+
+  for (auto iZone = 0u; iZone < nZone-1; iZone++) {
+    geometry[nZone-1][INST_0][MESH_0]->SetAvgTurboGeoValues(config[iZone],geometry[iZone][INST_0][MESH_0], iZone);
+  }
+
+  /*--- Transfer number of blade to ZONE_0 to correctly compute turbo performance---*/
+  for (auto iZone = 1u; iZone < nZone; iZone++) {
+    auto nBlades = config[iZone]->GetnBlades(iZone);
+    config[ZONE_0]->SetnBlades(iZone, nBlades);
+  }
+}
+
 void CGeometry::SetCustomBoundary(CConfig* config) {
   unsigned short iMarker;
   unsigned long iVertex;

SU2_CFD/include/iteration/CDiscAdjFluidIteration.hpp

@@ -157,7 +157,7 @@ class CDiscAdjFluidIteration final : public CIteration {
    * \param[in] iInst - Index of the zone.
    * \param[in] kind_recording - The kind of recording (geometry or flow).
    */
-  void SetDependencies(CSolver***** solver, CGeometry**** geometry, CNumerics****** numerics, CConfig** config,
+  void SetDependencies(CSolver***** solver, CGeometry**** geometry, CNumerics****** numerics, CConfig** config, CInterface*** interface,
                        unsigned short iZone, unsigned short iInst, RECORDING kind_recording) override;
 
 };

SU2_CFD/include/iteration/CFluidIteration.hpp

@@ -116,11 +116,6 @@ class CFluidIteration : public CIteration {
    */
   void TurboMonitor(CGeometry**** geometry_container, CConfig** config_container, unsigned long ExtIter, unsigned short iZone);
 
-  /*!
-   * \brief Computes turboperformance.
-   */
-  void ComputeTurboPerformance(CSolver***** solver, CGeometry**** geometry_container, CConfig** config_container, unsigned long ExtIter);
-
   /*!
    * \brief Postprocesses the fluid system before heading to another physics system or the next iteration.
    * \param[in] solver - Container vector with all the solutions.

SU2_CFD/include/iteration/CIteration.hpp

@@ -293,6 +293,19 @@ class CIteration {
   virtual void SetDependencies(CSolver***** solver, CGeometry**** geometry, CNumerics****** numerics, CConfig** config,
                                unsigned short iZone, unsigned short iInst, RECORDING kind_recording) {}
 
+  virtual void SetDependencies(CSolver***** solver, CGeometry**** geometry, CNumerics****** numerics, CConfig** config, CInterface*** interface,
+                                unsigned short iZone, unsigned short iInst, RECORDING kind_recording) {} 
+
   virtual void RegisterOutput(CSolver***** solver, CGeometry**** geometry, CConfig** config,
                               unsigned short iZone, unsigned short iInst) {}
+
+  /*!
+   * \brief Computes turboperformance.
+   */
+  virtual void ComputeTurboPerformance(CSolver***** solver, CGeometry**** geometry_container, CConfig** config_container){}
+
+  /*!
+   * \brief Initialises turboperformance classes.
+   */
+  virtual void InitTurboPerformance(CGeometry *geometry, CConfig** config, CFluidModel *fluid){};
 };

SU2_CFD/include/iteration/CTurboIteration.hpp

@@ -67,5 +67,10 @@ class CTurboIteration : public CFluidIteration {
   /*!
    * \brief Initialises turboperformance classes.
    */
-  void InitTurboPerformance(CGeometry *geometry, CConfig** config, CFluidModel *fluid);
+  void InitTurboPerformance(CGeometry *geometry, CConfig** config, CFluidModel *fluid) override;
+
+  /*!
+   * \brief Computes turboperformance.
+   */
+  void ComputeTurboPerformance(CSolver***** solver, CGeometry**** geometry_container, CConfig** config_container) override;
 };

SU2_CFD/include/solvers/CEulerSolver.hpp

@@ -155,6 +155,10 @@ class CEulerSolver : public CFVMFlowSolverBase<CEulerVariable, ENUM_REGIME::COMP
 
   vector<su2matrix<complex<su2double> > > CkInflow, CkOutflow1, CkOutflow2;
 
+  vector<su2double> EntropyGeneration;
+  vector<su2double> TotalPressureLoss;
+  vector<su2double> KineticEnergyLoss;
+
   /*--- End of Turbomachinery Solver Variables ---*/
 
   /*!
@@ -1145,6 +1149,35 @@ class CEulerSolver : public CFVMFlowSolverBase<CEulerVariable, ENUM_REGIME::COMP
     }
   }
 
+  inline su2double GetTurboObjectiveFunction(short unsigned int ObjFunc, int bladeRow) const final { 
+    switch (ObjFunc) {
+      case ENTROPY_GENERATION:
+        return EntropyGeneration[bladeRow];
+      case TOTAL_PRESSURE_LOSS:
+        return TotalPressureLoss[bladeRow];
+      case KINETIC_ENERGY_LOSS:
+        return KineticEnergyLoss[bladeRow];
+      default:
+        return 0.0;
+    }
+  }
+
+  inline void SetTurboObjectiveFunction(short unsigned int ObjFunc, int bladeRow, su2double val) final { 
+    switch (ObjFunc) {
+      case ENTROPY_GENERATION:
+        EntropyGeneration[bladeRow] = val;
+        break;
+      case TOTAL_PRESSURE_LOSS:
+        TotalPressureLoss[bladeRow] = val;
+        break;
+      case KINETIC_ENERGY_LOSS:
+        KineticEnergyLoss[bladeRow] = val;
+        break;
+      default:
+        break;
+    }
+  }
+
   /*!
    * \brief it take a velocity in the cartesian reference of framework and transform into the turbomachinery frame of reference.
    * \param[in] cartesianVelocity - cartesian components of velocity vector.

SU2_CFD/include/solvers/CFVMFlowSolverBase.inl

@@ -672,11 +672,12 @@ void CFVMFlowSolverBase<V, R>::ComputeVorticityAndStrainMag(const CConfig& confi
     StrainMag(iPoint) = sqrt(2.0*StrainMag(iPoint));
     AD::SetPreaccOut(StrainMag(iPoint));
 
+    AD::EndPreacc();
+
     /*--- Max is not differentiable, so we not register them for preacc. ---*/
     strainMax = max(strainMax, StrainMag(iPoint));
     omegaMax = max(omegaMax, GeometryToolbox::Norm(3, Vorticity));
 
-    AD::EndPreacc();
   }
   END_SU2_OMP_FOR
 

SU2_CFD/include/solvers/CSolver.hpp

@@ -1089,6 +1089,10 @@ class CSolver {
                                CConfig *config,
                                unsigned short val_marker) { }
 
+  inline virtual void SetTurboObjectiveFunction(short unsigned int ObjFunc, int bladeRow, su2double val)  { }
+
+  inline virtual su2double GetTurboObjectiveFunction(short unsigned int ObjFunc, int bladeRow) const { return 0.0; }
+
   /*!
    * \brief A virtual member.
    * \param[in] geometry - Geometrical definition of the problem.

SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp

@@ -622,7 +622,7 @@ void CDiscAdjMultizoneDriver::SetRecording(RECORDING kind_recording, Kind_Tape t
 
   for (iZone = 0; iZone < nZone; iZone++) {
     iteration_container[iZone][INST_0]->SetDependencies(solver_container, geometry_container, numerics_container,
-                                                        config_container, iZone, INST_0, kind_recording);
+                                                        config_container, interface_container, iZone, INST_0, kind_recording);
   }
 
   AD::Push_TapePosition(); /// DEPENDENCIES
@@ -636,7 +636,7 @@ void CDiscAdjMultizoneDriver::SetRecording(RECORDING kind_recording, Kind_Tape t
     for (iZone = 0; iZone < nZone; iZone++) {
       if (Has_Deformation(iZone)) {
         iteration_container[iZone][INST_0]->SetDependencies(solver_container, geometry_container, numerics_container,
-                                                            config_container, iZone, INST_0, kind_recording);
+                                                            config_container, interface_container, iZone, INST_0, kind_recording);
       }
     }
     SetObjFunction(kind_recording);
@@ -718,6 +718,13 @@ void CDiscAdjMultizoneDriver::SetObjFunction(RECORDING kind_recording) {
         }
 
         direct_output[iZone]->SetHistoryOutput(geometry, solvers, config);
+        // Need to evaluate turbomachinery objective functions here as they are only calculated in the final zone
+        if (config->GetBoolTurbomachinery() && iZone == config->GetnZone()-1){
+          direct_iteration[iZone][INST_0]->InitTurboPerformance(geometry, config_container, solvers[FLOW_SOL]->GetFluidModel());
+          direct_iteration[iZone][INST_0]->ComputeTurboPerformance(solver_container, geometry_container, config_container);
+          const auto weight = config->GetWeight_ObjFunc(0);
+          ObjFunc += weight * solvers[FLOW_SOL]->GetTurboObjectiveFunction(config->GetKind_ObjFunc(), config_container[ZONE_0]->GetnMarker_Turbomachinery());
+        }
         ObjFunc += solvers[FLOW_SOL]->GetTotal_ComboObj();
         break;