[WIP] Multizone discrete adjoints for turbomachinery

.github/workflows/regression.yml

@@ -88,8 +88,6 @@ jobs:
           #  flags: '--buildtype=debug -Denable-directdiff=true -Denable-normal=false -Dwith-omp=true -Denable-mixedprec=true -Denable-pywrapper=true -Denable-tecio=false --warnlevel=3 --werror'
     runs-on: ${{ inputs.runner || 'ubuntu-latest' }}
     steps:
-      - name: Reduce ASLR entropy for tsan
-        run: sudo sysctl -w vm.mmap_rnd_bits=28
       - name: Cache Object Files
         uses: actions/cache@v4
         with:
@@ -130,8 +128,6 @@ jobs:
             flags: '--buildtype=debugoptimized -Denable-autodiff=true -Denable-normal=false -Dwith-mpi=disabled --warnlevel=3 --werror'
     runs-on: ${{ inputs.runner || 'ubuntu-latest' }}
     steps:
-      - name: Reduce ASLR entropy for asan
-        run: sudo sysctl -w vm.mmap_rnd_bits=28
       - name: Cache Object Files
         uses: actions/cache@v4
         with:
@@ -229,8 +225,6 @@ jobs:
       matrix:
         testscript: ['hybrid_regression.py', 'hybrid_regression_AD.py']
     steps:
-      - name: Reduce ASLR entropy for tsan
-        run: sudo sysctl -w vm.mmap_rnd_bits=28
       - name: Pre Cleanup
         uses: docker://ghcr.io/su2code/su2/test-su2-tsan:240320-1536
         with:
@@ -276,8 +270,6 @@ jobs:
       matrix:
         testscript: ['serial_regression.py', 'serial_regression_AD.py']
     steps:
-      - name: Reduce ASLR entropy for asan
-        run: sudo sysctl -w vm.mmap_rnd_bits=28
       - name: Pre Cleanup
         uses: docker://ghcr.io/su2code/su2/test-su2-asan:240320-1536
         with:

Common/include/CConfig.hpp

@@ -441,6 +441,11 @@ class CConfig {
   TURBO_PERF_KIND *Kind_TurboPerf;           /*!< \brief Kind of turbomachynery architecture.*/
   TURBOMACHINERY_TYPE *Kind_TurboMachinery;
 
+  /* 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. */
@@ -7930,7 +7935,28 @@ class CConfig {
   void SetSurface_Species_Variance(unsigned short val_marker, su2double val_surface_species_variance) { Surface_Species_Variance[val_marker] = val_surface_species_variance; }
 
   /*!
-   * \brief Get the back pressure (static) at an outlet boundary.
+   * \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.ß
    * \param[in] val_index - Index corresponding to the outlet boundary.
    * \return The outlet pressure.
    */
@@ -8209,6 +8235,25 @@ 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/option_structure.hpp

@@ -1993,6 +1993,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)
@@ -2035,6 +2038,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

@@ -1043,6 +1043,11 @@
   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;
@@ -3391,8 +3396,17 @@
 
     /*---  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.  ---*/
@@ -3920,6 +3934,13 @@
         SU2_MPI::Error("The use of FLUID_MIXTURE requires the INC_DENSITY_MODEL option to be VARIABLE",
                        CURRENT_FUNCTION);
       }
+      /*--- Check whether the Kind scalar model used is correct, in the case of FLUID_MIXTURE the kind scalar model must
+       be SPECIES_TRANSPORT. Otherwise, if the scalar model is NONE, the species transport equations will not be solved.
+       --- */
+      if (Kind_Species_Model != SPECIES_MODEL::SPECIES_TRANSPORT) {
+        SU2_MPI::Error("The use of FLUID_MIXTURE requires the KIND_SCALAR_MODEL option to be SPECIES_TRANSPORT",
+                       CURRENT_FUNCTION);
+      }
 
       switch (Kind_ViscosityModel) {
         case VISCOSITYMODEL::CONSTANT:
@@ -5984,6 +6005,11 @@
         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++) {
@@ -8417,6 +8443,9 @@
         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{

SU2_CFD/include/output/CTurboOutput.hpp

@@ -87,7 +87,7 @@ class CTurbomachineryCombinedPrimitiveStates {
 class CTurbomachineryState {
  private:
   su2double Density, Pressure, Entropy, Enthalpy, Temperature, TotalTemperature, TotalPressure, TotalEnthalpy;
-  su2double AbsFlowAngle, FlowAngle, MassFlow, Rothalpy, TotalRelPressure;
+  su2double AbsFlowAngle, FlowAngle, MassFlow, Rothalpy, TangVelocity, TotalRelPressure;
   vector<su2double> Velocity, RelVelocity, Mach, RelMach;
   su2double Area, Radius;
 
@@ -124,6 +124,8 @@ class CTurbomachineryState {
 
   const su2double& GetRothalpy() const { return Rothalpy; }
 
+  const su2double& GetTangVelocity() const { return TangVelocity; }
+
   const vector<su2double>& GetVelocity() const { return Velocity; }
 
   const vector<su2double>& GetMach() const { return Mach; }
@@ -207,7 +209,7 @@ class CPropellorBladePerformance : public CTurbomachineryBladePerformance {
  */
 class CTurbomachineryStagePerformance {
  protected:
-  su2double TotalStaticEfficiency, TotalTotalEfficiency, NormEntropyGen, TotalStaticPressureRatio, TotalTotalPressureRatio, EulerianWork;
+  su2double TotalStaticEfficiency, TotalTotalEfficiency, NormEntropyGen, TotalStaticPressureRatio, TotalTotalPressureRatio, EulerianWork, TotalPressureLoss, KineticEnergyLoss;
   CFluidModel& fluidModel;
 
  public:
@@ -232,6 +234,10 @@ class CTurbomachineryStagePerformance {
   su2double GetTotalStaticPressureRatio() const { return TotalStaticPressureRatio; }
 
   su2double GetTotalTotalPressureRatio() const { return TotalTotalPressureRatio; }
+
+  su2double GetTotalPressureLoss() const { return TotalPressureLoss; }
+
+  su2double GetKineticEnergyLoss() const { return KineticEnergyLoss; }
 };
 
 /*!

SU2_CFD/include/solvers/CFVMFlowSolverBase.inl

@@ -2908,6 +2908,12 @@ su2double CFVMFlowSolverBase<V,R>::EvaluateCommonObjFunc(const CConfig& config)
     case SURFACE_SPECIES_VARIANCE:
       objFun += weight * config.GetSurface_Species_Variance(0);
       break;
+    case ENTROPY_GENERATION:
+      objFun += weight * config.GetEntropyGeneration(config.GetnMarker_Turbomachinery());
+    case TOTAL_PRESSURE_LOSS:
+      objFun += weight * config.GetTotalPressureLoss(config.GetnMarker_Turbomachinery());
+    case KINETIC_ENERGY_LOSS:
+      objFun += weight * config.GetKineticEnergyLoss(config.GetnMarker_Turbomachinery());
     case CUSTOM_OBJFUNC:
       objFun += weight * Total_Custom_ObjFunc;
       break;

SU2_CFD/src/iteration/CFluidIteration.cpp

@@ -325,16 +325,15 @@ void CFluidIteration::TurboMonitor(CGeometry**** geometry_container, CConfig** c
 void CFluidIteration::ComputeTurboPerformance(CSolver***** solver, CGeometry**** geometry_container, CConfig** config_container, unsigned long ExtIter) {
   unsigned short nDim = geometry_container[ZONE_0][INST_0][MESH_0]->GetnDim();
   unsigned short nBladesRow = config_container[ZONE_0]->GetnMarker_Turbomachinery();
-  unsigned short iBlade=0, iSpan;
   vector<su2double> TurboPrimitiveIn, TurboPrimitiveOut;
   std::vector<std::vector<CTurbomachineryCombinedPrimitiveStates>> bladesPrimitives;
 
   if (rank == MASTER_NODE) {
-      for (iBlade = 0; iBlade < nBladesRow; iBlade++){
+      for (auto iBlade = 0u; iBlade < nBladesRow; iBlade++){
       /* Blade Primitive initialized per blade */
       std::vector<CTurbomachineryCombinedPrimitiveStates> bladePrimitives;
       auto nSpan = config_container[iBlade]->GetnSpanWiseSections();
-      for (iSpan = 0; iSpan < nSpan + 1; iSpan++) {
+      for (auto iSpan = 0u; iSpan < nSpan + 1; iSpan++) {
         TurboPrimitiveIn= solver[iBlade][INST_0][MESH_0][FLOW_SOL]->GetTurboPrimitive(iBlade, iSpan, true);
         TurboPrimitiveOut= solver[iBlade][INST_0][MESH_0][FLOW_SOL]->GetTurboPrimitive(iBlade, iSpan, false);
         auto spanInletPrimitive = CTurbomachineryPrimitiveState(TurboPrimitiveIn, nDim, geometry_container[iBlade][INST_0][MESH_0]->GetTangGridVelIn(iBlade, iSpan));
@@ -352,6 +351,21 @@ void CFluidIteration::ComputeTurboPerformance(CSolver***** solver, CGeometry****
     auto OutState =  TurbomachineryPerformance->GetBladesPerformances().at(nZone-1).at(nSpan)->GetOutletState();
 
     TurbomachineryStagePerformance->ComputePerformanceStage(InState, OutState, config_container[nZone-1]);
+
+    /*--- Set turbomachinery objective function value in each zone ---*/
+    for (auto iBlade = 0u; iBlade < nBladesRow; iBlade++) {
+      // Should we set in ZONE_0 or nZone-1?
+      auto iBladePerf = TurbomachineryPerformance->GetBladesPerformances().at(iBlade).at(nSpan);
+      InState = iBladePerf->GetInletState();
+      OutState = iBladePerf->GetOutletState();
+      config_container[nZone-1]->SetEntropyGeneration(iBlade, (OutState.GetEntropy() - InState.GetEntropy())/InState.GetEntropy() * 100);
+      config_container[nZone-1]->SetTotalPressureLoss(iBlade, iBladePerf->GetTotalPressureLoss());
+      config_container[nZone-1]->SetKineticEnergyLoss(iBlade, iBladePerf->GetKineticEnergyLoss());
+    }
+    /*--- Set global turbomachinery objective function ---*/
+    config_container[nZone-1]->SetEntropyGeneration(nBladesRow, TurbomachineryStagePerformance->GetNormEntropyGen()*100);
+    config_container[nZone-1]->SetTotalPressureLoss(nBladesRow, TurbomachineryStagePerformance->GetTotalPressureLoss());
+    config_container[nZone-1]->SetKineticEnergyLoss(nBladesRow, TurbomachineryStagePerformance->GetKineticEnergyLoss());
   }
 }
 

SU2_CFD/src/output/CFlowCompOutput.cpp

@@ -577,7 +577,7 @@ void CFlowCompOutput::SetTurboMultiZonePerformance_Output(std::shared_ptr<CTurbo
 
 void CFlowCompOutput::LoadTurboHistoryData(std::shared_ptr<CTurbomachineryStagePerformance> TurboStagePerf, std::shared_ptr<CTurboOutput> TurboPerf, CConfig *config) {
   auto BladePerformance = TurboPerf->GetBladesPerformances();
-  for (unsigned short iZone = 0; iZone <= config->GetnZone()-1; iZone++) {
+  for (auto iZone = 0u; iZone < config->GetnZone(); iZone++) {
     auto nSpan = config->GetnSpan_iZones(iZone);
     const auto& BladePerf = BladePerformance.at(iZone).at(nSpan);
 

SU2_CFD/src/output/CTurboOutput.cpp

@@ -60,7 +60,7 @@ void CTurbomachineryState::ComputeState(CFluidModel& fluidModel, const CTurbomac
   Pressure = primitiveState.GetPressure();
   std::vector<su2double> velocity = primitiveState.GetVelocity();
   Velocity.assign(velocity.begin(), velocity.end());
-  su2double tangVel = primitiveState.GetTangVelocity();
+  TangVelocity = primitiveState.GetTangVelocity();
 
   /*--- Compute static TD quantities ---*/
   fluidModel.SetTDState_Prho(Pressure, Density);
@@ -81,9 +81,9 @@ void CTurbomachineryState::ComputeState(CFluidModel& fluidModel, const CTurbomac
   std::for_each(Mach.begin(), Mach.end(), [&](su2double& el) { el /= soundSpeed; });
 
   /*--- Compute relative kinematic quantities ---*/
-  su2double tangVel2 = tangVel * tangVel;
+  su2double tangVel2 = TangVelocity * TangVelocity;
   RelVelocity.assign(Velocity.begin(), Velocity.end());
-  RelVelocity[1] -= tangVel;
+  RelVelocity[1] -= TangVelocity;
   su2double relVel2 = GetRelVelocityValue();
   FlowAngle = atan(RelVelocity[1] / RelVelocity[0]);
   RelMach.assign(RelVelocity.begin(), RelVelocity.end());
@@ -245,13 +245,20 @@ void CTurbomachineryStagePerformance::ComputeTurbineStagePerformance(const CTurb
   su2double enthalpyOutIs = fluidModel.GetStaticEnergy() + OutState.GetPressure() / fluidModel.GetDensity();
   su2double totEnthalpyOutIs = enthalpyOutIs + 0.5 * OutState.GetVelocityValue() * OutState.GetVelocityValue();
 
+  su2double tangVel = OutState.GetTangVelocity();
+  su2double relVelOutIs2 = 2 * (OutState.GetRothalpy() - enthalpyOutIs) + tangVel * tangVel;
+
   /*--- Compute turbine stage performance ---*/
   NormEntropyGen = (OutState.GetEntropy() - InState.GetEntropy()) / InState.GetEntropy();
   EulerianWork = InState.GetTotalEnthalpy() - OutState.GetTotalEnthalpy();
   TotalStaticEfficiency = EulerianWork / (InState.GetTotalEnthalpy() - enthalpyOutIs);
   TotalTotalEfficiency = EulerianWork / (InState.GetTotalEnthalpy() - totEnthalpyOutIs);
   TotalStaticPressureRatio = InState.GetTotalPressure() / OutState.GetPressure();
   TotalTotalPressureRatio = InState.GetTotalPressure() / OutState.GetTotalPressure();
+
+  TotalPressureLoss = (InState.GetTotalRelPressure() - OutState.GetTotalRelPressure()) /
+                      (OutState.GetTotalRelPressure() - OutState.GetPressure());
+  KineticEnergyLoss = 2 * (OutState.GetEnthalpy() - enthalpyOutIs) / relVelOutIs2;
 }
 
 void CTurbomachineryStagePerformance::ComputeCompressorStagePerformance(const CTurbomachineryState& InState,
@@ -260,6 +267,8 @@ void CTurbomachineryStagePerformance::ComputeCompressorStagePerformance(const CT
   fluidModel.SetTDState_Ps(OutState.GetPressure(), InState.GetEntropy());
   su2double enthalpyOutIs = fluidModel.GetStaticEnergy() + OutState.GetPressure() / fluidModel.GetDensity();
   su2double totEnthalpyOutIs = enthalpyOutIs + 0.5 * OutState.GetVelocityValue() * OutState.GetVelocityValue();
+  su2double tangVel = OutState.GetTangVelocity();
+  su2double relVelOutIs2 = 2 * (OutState.GetRothalpy() - enthalpyOutIs) + tangVel * tangVel;
 
   /*--- Compute compressor stage performance ---*/
   NormEntropyGen = (OutState.GetEntropy() - InState.GetEntropy()) / InState.GetEntropy();
@@ -268,4 +277,8 @@ void CTurbomachineryStagePerformance::ComputeCompressorStagePerformance(const CT
   TotalTotalEfficiency = (totEnthalpyOutIs - InState.GetTotalEnthalpy()) / EulerianWork;
   TotalStaticPressureRatio = OutState.GetPressure() / InState.GetTotalPressure();
   TotalTotalPressureRatio = OutState.GetTotalPressure() / InState.GetTotalPressure();
+
+  TotalPressureLoss = (InState.GetTotalRelPressure() - OutState.GetTotalRelPressure()) /
+                      (InState.GetTotalRelPressure() - InState.GetPressure());
+  KineticEnergyLoss = 2 * (OutState.GetEnthalpy() - enthalpyOutIs) / relVelOutIs2;
 }