Support non-turbo (solid/heat) zones in multizone turbomachinery

Common/include/CConfig.hpp

@@ -480,6 +480,7 @@ class CConfig {
   unsigned short nFFDBox;              /*!< \brief Number of ffd boxes. */
   unsigned short nTurboMachineryKind;  /*!< \brief Number turbomachinery types specified. */
   unsigned short nTurboInterfaces;     /*!< \brief Number of turbomachiery interfaces */
+  unsigned short nTurboZones = 0;      /*!< \brief Number of turbomachinery (blade row) zones, these are the first zones. */
   unsigned short nParamDV;             /*!< \brief Number of parameters of the design variable. */
   string DV_Filename;                  /*!< \brief Filename for providing surface positions from an external parameterization. */
   string DV_Unordered_Sens_Filename;   /*!< \brief Filename of volume sensitivities in an unordered ASCII format. */
@@ -5517,6 +5518,14 @@ class CConfig {
    */
   bool GetBoolTurbomachinery(void) const { return (nMarker_Turbomachinery !=0);}
 
+  /*!
+   * \brief Number of turbomachinery (blade row) zones in the problem, i.e. one per MARKER_TURBOMACHINERY pair.
+   * \note Blade rows are zones 0 to nTurboZones-1, any further zones (e.g. solids coupled via CHT) are not turbo zones.
+   *       Unlike GetBoolTurbomachinery, this is also known in the non-turbo zones.
+   * \return Number of turbomachinery zones.
+   */
+  unsigned short GetnTurboZones(void) const { return nTurboZones;}
+
   /*!
    * \brief Verify if a sliding plane for relative frame is specified from config file.
    * \return boolean.

Common/src/CConfig.cpp

@@ -3707,6 +3707,46 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i
     Multizone_Problem = YES;
   }
 
+  /*--- Turbomachinery options are usually set in the master config and inherited by every zone.
+   *    Only compressible flow zones can be blade rows, other zones (e.g. blade solids coupled
+   *    via CHT) drop the inherited turbomachinery, Giles and mixing-plane markers. The blade
+   *    rows must be the first zones, as the turbo markers and performance are indexed by zone. ---*/
+
+  nTurboZones = nMarker_Turbomachinery;
+
+  if (!base_config && nZone > 1 && nTurboZones > 0) {
+    bool turboCapableZone = false;
+    switch (Kind_Solver) {
+      case MAIN_SOLVER::EULER: case MAIN_SOLVER::NAVIER_STOKES: case MAIN_SOLVER::RANS:
+      case MAIN_SOLVER::DISC_ADJ_EULER: case MAIN_SOLVER::DISC_ADJ_NAVIER_STOKES: case MAIN_SOLVER::DISC_ADJ_RANS:
+        turboCapableZone = true;
+        break;
+      default:
+        break;
+    }
+
+    if (turboCapableZone && iZone >= nTurboZones) {
+      SU2_MPI::Error("Zone " + to_string(iZone) + " is a compressible flow zone but there are only " +
+                     to_string(nTurboZones) + " MARKER_TURBOMACHINERY pairs.\n"
+                     "The blade rows must be the first zones in CONFIG_LIST.", CURRENT_FUNCTION);
+    }
+    if (!turboCapableZone) {
+      if (iZone < nTurboZones) {
+        SU2_MPI::Error("Zone " + to_string(iZone) + " is not a compressible flow zone, but it is one of the first " +
+                       to_string(nTurboZones) + " zones (one per MARKER_TURBOMACHINERY pair).\n"
+                       "The blade rows must be the first zones in CONFIG_LIST.", CURRENT_FUNCTION);
+      }
+      if (Giles_FlowDir != nullptr) {
+        for (iMarker = 0; iMarker < nMarker_Giles; iMarker++) delete [] Giles_FlowDir[iMarker];
+        delete [] Giles_FlowDir;
+        Giles_FlowDir = nullptr;
+      }
+      nMarker_Turbomachinery = 0;
+      nMarker_Giles = 0;
+      nMarker_MixingPlaneInterface = 0;
+    }
+  }
+
   /*--- The solver vectors stay on the device but the halo exchange is host-side, so more than
    * one rank would use stale halos. Use OpenMP for the host parts instead. ---*/
   if (Enable_Cuda && size > 1) {
@@ -4597,7 +4637,7 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i
     }
 
     /*--- Check if turbomachinery performance kind is specified with turbo markers ---*/
-    if (GetBoolTurbomachinery() && !(nTurboMachineryKind/nZone == 1)){
+    if (GetBoolTurbomachinery() && nTurboMachineryKind < nTurboZones){
       SU2_MPI::Error("Insufficient TURBO_PERF_KIND options specified with turbomachinery markers", CURRENT_FUNCTION);
     }
 

SU2_CFD/src/drivers/CDriver.cpp

@@ -2435,7 +2435,10 @@ void CDriver::InitializeInterface(CConfig **config, CSolver***** solver, CGeomet
 
         /*--- Setup the interpolation. ---*/
 
-        if (!config[donor]->GetBoolTurbomachinery()) {
+        /*--- Blade row to blade row interfaces build their own (mixing-plane or sliding) interpolation below. ---*/
+        const bool turbo_interface = config[donor]->GetBoolTurbomachinery() && config[target]->GetBoolTurbomachinery();
+
+        if (!turbo_interface) {
           interpolation[donor][target] = unique_ptr<CInterpolator>(CInterpolatorFactory::CreateInterpolator(
                                        geometry, config, interpolation[target][donor].get(), donor, target, false));
           if (rank == MASTER_NODE) cout << " Transferring ";
@@ -2511,7 +2514,7 @@ void CDriver::InitializeInterface(CConfig **config, CSolver***** solver, CGeomet
         }
         else if (fluid_donor && fluid_target) {
           /*--- Interface handling for turbomachinery applications. ---*/
-          if (config[donor]->GetBoolTurbomachinery()) {
+          if (turbo_interface) {
             auto interfaceIndex = donor+target; // Here we assume that the interfaces at each side are the same kind
             switch (config[donor]->GetKind_TurboInterface(interfaceIndex)) {
               case TURBO_INTERFACE_KIND::MIXING_PLANE: {
@@ -2655,6 +2658,7 @@ void CDriver::PreprocessTurboVertex(CConfig** config, CGeometry**** geometry, CS
 
   unsigned short nSpanMax = 0;
   mixingplane = config[ZONE_0]->GetBoolMixingPlaneInterface();
+  const auto nTurboZones = config[ZONE_0]->GetnTurboZones();
 
   /*--- Create turbovertex structure ---*/
   if (rank == MASTER_NODE) cout<<endl<<"Initialize Turbo Vertex Structure." << endl;
@@ -2667,7 +2671,7 @@ void CDriver::PreprocessTurboVertex(CConfig** config, CGeometry**** geometry, CS
         nSpanMax = config[iZone]->GetnSpanWiseSections();
       }
 
-      config[nZone-1]->SetnSpan_iZones(config[iZone]->GetnSpanWiseSections(), iZone);
+      config[nTurboZones-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);
@@ -2691,25 +2695,28 @@ void CDriver::PreprocessTurbomachinery(CConfig** config, CGeometry**** geometry,
   bool discrete_adjoint = config[ZONE_0]->GetDiscrete_Adjoint();
   su2double areaIn, areaOut, nBlades, flowAngleIn, flowAngleOut;
 
+  /*--- Only the blade rows (the first zones) are turbo zones, other zones (e.g. solids) are skipped. ---*/
+  const auto nTurboZones = config[ZONE_0]->GetnTurboZones();
+
   // TODO(turbo): make it general for turbo HB
   if (rank == MASTER_NODE) cout<<"Compute inflow and outflow average geometric quantities." << endl;
-  for (iZone = 0; iZone < nZone; iZone++) {
+  for (iZone = 0; iZone < nTurboZones; 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 (rank == MASTER_NODE) cout<<"Initialize solver containers for average quantities." << endl;
   if (!dummy){
-    for (iZone = 0; iZone < nZone; iZone++) {
+    for (iZone = 0; iZone < nTurboZones; iZone++) {
       solver[iZone][INST_0][MESH_0][FLOW_SOL]->InitTurboContainers(geometry[iZone][INST_0][MESH_0],config, iZone);
     }
   }
 
   if(mixingplane){
     if (rank == MASTER_NODE) cout << "Set span-wise sections between zones on Mixing-Plane interface." << endl;
-    for (donorZone = 0; donorZone < nZone; donorZone++) {
-      for (targetZone = 0; targetZone < nZone; targetZone++) {
+    for (donorZone = 0; donorZone < nTurboZones; donorZone++) {
+      for (targetZone = 0; targetZone < nTurboZones; targetZone++) {
         if (interface_container[donorZone][targetZone] != nullptr){
           interface[donorZone][targetZone]->SetSpanWiseLevels(config[donorZone], config[targetZone]);
         }
@@ -2718,13 +2725,13 @@ void CDriver::PreprocessTurbomachinery(CConfig** config, CGeometry**** geometry,
   }
 
   /*--- Transfer number of blade to ZONE_0 to correctly compute turbo performance---*/
-  for (iZone = 1; iZone < nZone; iZone++) {
+  for (iZone = 1; iZone < nTurboZones; iZone++) {
     nBlades = config[iZone]->GetnBlades(iZone);
     config[ZONE_0]->SetnBlades(iZone, nBlades);
   }
 
   if (rank == MASTER_NODE && !dummy){
-    for (iZone = 0; iZone < nZone; iZone++) {
+    for (iZone = 0; iZone < nTurboZones; iZone++) {
     areaIn  = geometry[iZone][INST_0][MESH_0]->GetSpanAreaIn(iZone, config[iZone]->GetnSpanWiseSections());
     areaOut = geometry[iZone][INST_0][MESH_0]->GetSpanAreaOut(iZone, config[iZone]->GetnSpanWiseSections());
     nBlades = config[iZone]->GetnBlades(iZone);
@@ -2736,14 +2743,14 @@ void CDriver::PreprocessTurbomachinery(CConfig** config, CGeometry**** geometry,
 
   if(!restart && !discrete_adjoint){
     if (rank == MASTER_NODE) cout<<"Initialize turbomachinery solution quantities." << endl;
-    for(iZone = 0; iZone < nZone; iZone++) {
+    for(iZone = 0; iZone < nTurboZones; iZone++) {
       solver[iZone][INST_0][MESH_0][FLOW_SOL]->SetFreeStream_TurboSolution(config[iZone]);
     }
   }
 
   if (dummy) return; // No need to go further for a dummy run
   if (rank == MASTER_NODE) cout<<"Initialize inflow and outflow average solution quantities." << endl;
-  for(iZone = 0; iZone < nZone; iZone++) {
+  for(iZone = 0; iZone < nTurboZones; iZone++) {
     solver[iZone][INST_0][MESH_0][FLOW_SOL]->PreprocessAverage(solver[iZone][INST_0][MESH_0], geometry[iZone][INST_0][MESH_0],config[iZone],INFLOW);
     solver[iZone][INST_0][MESH_0][FLOW_SOL]->PreprocessAverage(solver[iZone][INST_0][MESH_0], geometry[iZone][INST_0][MESH_0],config[iZone],OUTFLOW);
     solver[iZone][INST_0][MESH_0][FLOW_SOL]->TurboAverageProcess(solver[iZone][INST_0][MESH_0], geometry[iZone][INST_0][MESH_0],config[iZone],INFLOW);

SU2_CFD/src/iteration/CFluidIteration.cpp

@@ -239,9 +239,10 @@ bool CFluidIteration::Monitor(COutput* output, CIntegration**** integration, CGe
 
   /*--- Turbomachinery Specific Montior ---*/
   if (config[ZONE_0]->GetBoolTurbomachinery()){
-    if (val_iZone == config[ZONE_0]->GetnZone()-1) {
+    /*--- Stage performance is computed once all blade rows are updated, i.e. by the last turbo zone. ---*/
+    if (val_iZone == config[ZONE_0]->GetnTurboZones()-1) {
       ComputeTurboPerformance(solver, geometry, config);
-      auto TurbomachineryBladePerformances = GetBladesPerformanceVector(solver, config[val_iZone]->GetnZone());
+      auto TurbomachineryBladePerformances = GetBladesPerformanceVector(solver, config[val_iZone]->GetnTurboZones());
 
       output->SetHistoryOutput(geometry, solver, config, TurbomachineryStagePerformance, TurbomachineryBladePerformances,
             val_iZone, config[val_iZone]->GetTimeIter(), config[val_iZone]->GetOuterIter(),

SU2_CFD/src/iteration/CIteration.cpp

@@ -220,16 +220,18 @@ void CIteration::InitTurboPerformance(CGeometry* geometry, CConfig** config, CFl
 
 void CIteration::ComputeTurboPerformance(CSolver***** solver, CGeometry**** geometry_container, CConfig** config_container) {
   // Computes the turboperformance per blade in zone iBlade
-  const auto nZone = config_container[ZONE_0]->GetnZone();
+  /*--- The stage inlet is the first blade row and the outlet the last, non-turbo zones come after them. ---*/
+  const auto nTurboZones = config_container[ZONE_0]->GetnTurboZones();
+  const auto lastTurboZone = nTurboZones - 1;
 
   if (rank == MASTER_NODE) {
-    auto TurbomachineryBladePerformances = GetBladesPerformanceVector(solver, nZone);
+    auto TurbomachineryBladePerformances = GetBladesPerformanceVector(solver, nTurboZones);
 
     auto nSpan = config_container[ZONE_0]->GetnSpanWiseSections();
     auto InState = TurbomachineryBladePerformances[ZONE_0]->GetBladesPerformances().at(nSpan)->GetInletState();
-    nSpan = config_container[nZone-1]->GetnSpanWiseSections();
-    auto OutState =  TurbomachineryBladePerformances[nZone-1]->GetBladesPerformances().at(nSpan)->GetOutletState();
+    nSpan = config_container[lastTurboZone]->GetnSpanWiseSections();
+    auto OutState =  TurbomachineryBladePerformances[lastTurboZone]->GetBladesPerformances().at(nSpan)->GetOutletState();
 
-    TurbomachineryStagePerformance->ComputePerformanceStage(InState, OutState, config_container[nZone-1]);
+    TurbomachineryStagePerformance->ComputePerformanceStage(InState, OutState, config_container[lastTurboZone]);
   }
 }

SU2_CFD/src/output/CFlowCompOutput.cpp

@@ -206,7 +206,7 @@ void CFlowCompOutput::SetHistoryOutputFields(CConfig *config){
 
   AddNearfieldInverseDesignOutput();
 
-  if (config->GetBoolTurbomachinery()) AddTurboOutput(config->GetnZone());
+  if (config->GetBoolTurbomachinery()) AddTurboOutput(config->GetnTurboZones());
 
 }
 
@@ -535,7 +535,7 @@ void CFlowCompOutput::SetTurboPerformance_Output(su2vector<std::shared_ptr<CTurb
   TurboInOut.SetAlign(PrintingToolbox::CTablePrinter::RIGHT);
   TurboInOut.PrintHeader();
 
-  for (unsigned short iZone = 0; iZone <= config->GetnZone()-1; iZone++) {
+  for (unsigned short iZone = 0; iZone < config->GetnTurboZones(); iZone++) {
     auto nSpan = config->GetnSpan_iZones(iZone);
     const auto& BladePerf = TurboBladePerfs[iZone]->GetBladesPerformances().at(nSpan);
 
@@ -589,7 +589,7 @@ void CFlowCompOutput::SetTurboMultiZonePerformance_Output(std::shared_ptr<CTurbo
 }
 
 void CFlowCompOutput::LoadTurboHistoryData(std::shared_ptr<CTurbomachineryStagePerformance> TurboStagePerf, su2vector<std::shared_ptr<CTurboOutput>> TurboBladePerfs, CConfig *config) {
-  for (unsigned short iZone = 0; iZone <= config->GetnZone()-1; iZone++) {
+  for (unsigned short iZone = 0; iZone < config->GetnTurboZones(); iZone++) {
     auto nSpan = config->GetnSpan_iZones(iZone);
     const auto& BladePerf = TurboBladePerfs[iZone]->GetBladesPerformances().at(nSpan);
 

SU2_CFD/src/output/COutput.cpp

@@ -249,7 +249,7 @@ void COutput::SetHistoryOutput(CGeometry ****geometry, CSolver *****solver, CCon
     SetTurboMultiZonePerformance_Output(TurboStagePerf, TurboBladePerfs, config[val_iZone]); //Stage performance screen
   }
 
-  for (int iZone = 0; iZone < config[ZONE_0]->GetnZone(); iZone ++){
+  for (int iZone = 0; iZone < config[ZONE_0]->GetnTurboZones(); iZone ++){
     if (rank == MASTER_NODE) {
       WriteTurboSpanwisePerformance(TurboBladePerfs, geometry[iZone][val_iInst][MESH_0], config, iZone); //Spanwise files
     }