Allow FSI+CHT (aerothermoelasticity) (#2807)

* allow fsi+cht

* Apply suggestions from code review

Co-authored-by: Pedro Gomes <38071223+pcarruscag@users.noreply.github.com>

* set regression values

* back to develop

Author: pcarruscag

Common/src/CConfig.cpp

@@ -1311,9 +1311,9 @@ void CConfig::SetConfig_Options() {
   /* DESCRIPTION: Definition of the turbulent thermal conductivity model (CONSTANT_PRANDTL_TURB (default), NONE). */
   addEnumOption("TURBULENT_CONDUCTIVITY_MODEL", Kind_ConductivityModel_Turb, TurbConductivityModel_Map, CONDUCTIVITYMODEL_TURB::CONSTANT_PRANDTL);
 
- /*--- Options related to Constant Thermal Conductivity Model ---*/
+  /*--- Options related to Constant Thermal Conductivity Model ---*/
 
- /* DESCRIPTION: default value for AIR */
+  /* DESCRIPTION: default value for AIR */
   addDoubleListOption("THERMAL_CONDUCTIVITY_CONSTANT", nThermal_Conductivity_Constant , Thermal_Conductivity_Constant);
 
   /*--- Options related to temperature polynomial coefficients for fluid models. ---*/

SU2_CFD/include/interfaces/CInterface.hpp

@@ -29,6 +29,7 @@
 #pragma once
 
 #include "../../../Common/include/parallelization/mpi_structure.hpp"
+#include "../../../Common/include/option_structure.hpp"
 
 #include <cmath>
 #include <string>
@@ -77,7 +78,7 @@ class CInterface {
   /*!
    * \brief Constructor of the class.
    */
-  CInterface(void);
+  CInterface();
 
   /*!
    * \overload
@@ -89,7 +90,7 @@ class CInterface {
   /*!
    * \brief Destructor of the class.
    */
-  virtual ~CInterface(void);
+  virtual ~CInterface();
 
   /*!
    * \brief Interpolate data and broadcast it into all processors, for nonmatching meshes.
@@ -224,4 +225,11 @@ class CInterface {
    * \param[in] val_contact_resistance - Contact resistance value in m^2/W
    */
   inline virtual void SetContactResistance(su2double val_contact_resistance) {};
+
+  /*!
+   * \brief These can be used to chain interfaces between the same zones but for other variables,
+   * without having to mix physics in the interface classes. Currently this is used for FSI+CHT.
+   */
+  ENUM_TRANSFER NextInterfaceType = ENUM_TRANSFER::NO_TRANSFER;
+  CInterface* NextInterface = nullptr;
 };

SU2_CFD/src/drivers/CDriver.cpp

@@ -2443,6 +2443,31 @@ void CDriver::InitializeInterface(CConfig **config, CSolver***** solver, CGeomet
         interpolation[donor][target] = unique_ptr<CInterpolator>(CInterpolatorFactory::CreateInterpolator(
                                        geometry, config, interpolation[target][donor].get(), donor, target));
 
+        /*--- Helpers with logic to create CHT interfaces. ---*/
+
+        auto GetChtInterfaceType = [donor, target, config](bool heat_donor, bool heat_target) {
+          if (heat_donor && heat_target) return CONJUGATE_HEAT_SS;
+
+          const auto fluidZone = heat_target ? donor : target;
+          if (config[fluidZone]->GetEnergy_Equation() ||
+              config[fluidZone]->GetKind_Regime() == ENUM_REGIME::COMPRESSIBLE ||
+              config[fluidZone]->GetKind_FluidModel() == ENUM_FLUIDMODEL::FLUID_FLAMELET) {
+            return heat_target ? CONJUGATE_HEAT_FS : CONJUGATE_HEAT_SF;
+          } else if (config[fluidZone]->GetWeakly_Coupled_Heat()) {
+            return heat_target ? CONJUGATE_HEAT_WEAKLY_FS : CONJUGATE_HEAT_WEAKLY_SF;
+          }
+          return NO_TRANSFER;
+        };
+
+        auto MakeChtInterface = [&](const auto type) {
+          if (type != NO_TRANSFER) {
+            if (rank == MASTER_NODE) cout << " Conjugate heat variables." << endl;
+            return new CConjugateHeatInterface(4, 0);
+          }
+          if (rank == MASTER_NODE) cout << " NO heat variables." << endl;
+          return static_cast<CConjugateHeatInterface*>(nullptr);
+        };
+
         /*--- The type of variables transferred depends on the donor/target physics. ---*/
 
         const bool heat_target = config[target]->GetHeatProblem();
@@ -2467,16 +2492,26 @@ void CDriver::InitializeInterface(CConfig **config, CSolver***** solver, CGeomet
             interface[donor][target] = new CDiscAdjFlowTractionInterface(nDim, nConst, config[donor], conservative);
           }
           if (rank == MASTER_NODE) cout << "fluid " << (conservative? "forces." : "tractions.") << endl;
+
+          if (config[target]->GetWeakly_Coupled_Heat()) {
+            interface[donor][target]->NextInterfaceType = GetChtInterfaceType(false, true);
+            interface[donor][target]->NextInterface = MakeChtInterface(interface[donor][target]->NextInterfaceType);
+          }
         }
         else if (structural_donor && (fluid_target || heat_target)) {
           if (solver_container[target][INST_0][MESH_0][MESH_SOL] == nullptr) {
             SU2_MPI::Error("Mesh deformation was not correctly specified for the fluid/heat zone.\n"
                            "Use DEFORM_MESH=YES, and setup MARKER_DEFORM_MESH=(...)", CURRENT_FUNCTION);
           }
           interface_type = BOUNDARY_DISPLACEMENTS;
-          if (!config[donor]->GetTime_Domain()) interface[donor][target] = new CDisplacementsInterface(nDim, 0);
-          else interface[donor][target] = new CDisplacementsInterface(2*nDim, 0);
+          const auto nVar = config[donor]->GetTime_Domain() ? 2 * nDim : nDim;
+          interface[donor][target] = new CDisplacementsInterface(nVar, 0);
           if (rank == MASTER_NODE) cout << "boundary displacements from the structural solver." << endl;
+
+          if (fluid_target && config[donor]->GetWeakly_Coupled_Heat()) {
+            interface[donor][target]->NextInterfaceType = GetChtInterfaceType(true, false);
+            interface[donor][target]->NextInterface = MakeChtInterface(interface[donor][target]->NextInterfaceType);
+          }
         }
         else if (fluid_donor && fluid_target) {
           /*--- Interface handling for turbomachinery applications. ---*/
@@ -2487,48 +2522,27 @@ void CDriver::InitializeInterface(CConfig **config, CSolver***** solver, CGeomet
                 interface_type = MIXING_PLANE;
                 auto nVar = solver[donor][INST_0][MESH_0][FLOW_SOL]->GetnVar();
                 interface[donor][target] = new CMixingPlaneInterface(nVar, 0);
-                if (rank == MASTER_NODE) cout << "using a mixing-plane interface from donor zone " << donor << " to target zone " << target << "." << endl;
+                if (rank == MASTER_NODE) cout << " Using a mixing-plane interface from donor zone " << donor << " to target zone " << target << "." << endl;
                 break;
               }
               case TURBO_INTERFACE_KIND::FROZEN_ROTOR: {
                 auto nVar = solver[donor][INST_0][MESH_0][FLOW_SOL]->GetnPrimVar();
                 interface_type = SLIDING_INTERFACE;
                 interface[donor][target] = new CSlidingInterface(nVar, 0);
-                if (rank == MASTER_NODE) cout << "using a fluid interface interface from donor zone " << donor << " to target zone " << target << "." << endl;
+                if (rank == MASTER_NODE) cout << " Using a fluid interface interface from donor zone " << donor << " to target zone " << target << "." << endl;
               }
             }
           }
           else{
             auto nVar = solver[donor][INST_0][MESH_0][FLOW_SOL]->GetnPrimVar();
               interface_type = SLIDING_INTERFACE;
               interface[donor][target] = new CSlidingInterface(nVar, 0);
-              if (rank == MASTER_NODE) cout << "sliding interface." << endl;
+              if (rank == MASTER_NODE) cout << " Sliding interface." << endl;
           }
         }
         else if (heat_donor || heat_target) {
-          if (heat_donor && heat_target){
-            interface_type = CONJUGATE_HEAT_SS;
-
-          } else {
-
-            const auto fluidZone = heat_target? donor : target;
-            if (config[fluidZone]->GetEnergy_Equation() || (config[fluidZone]->GetKind_Regime() == ENUM_REGIME::COMPRESSIBLE)
-                || (config[fluidZone]->GetKind_FluidModel() == ENUM_FLUIDMODEL::FLUID_FLAMELET))
-              interface_type = heat_target? CONJUGATE_HEAT_FS : CONJUGATE_HEAT_SF;
-            else if (config[fluidZone]->GetWeakly_Coupled_Heat())
-              interface_type = heat_target? CONJUGATE_HEAT_WEAKLY_FS : CONJUGATE_HEAT_WEAKLY_SF;
-            else
-              interface_type = NO_TRANSFER;
-          }
-
-          if (interface_type != NO_TRANSFER) {
-            auto nVar = 4;
-            interface[donor][target] = new CConjugateHeatInterface(nVar, 0);
-            if (rank == MASTER_NODE) cout << "conjugate heat variables." << endl;
-          }
-          else {
-            if (rank == MASTER_NODE) cout << "NO heat variables." << endl;
-          }
+          interface_type = GetChtInterfaceType(heat_donor, heat_target);
+          interface[donor][target] = MakeChtInterface(interface_type);
         }
         else {
           if (solver[donor][INST_0][MESH_0][FLOW_SOL] == nullptr)
@@ -2537,7 +2551,7 @@ void CDriver::InitializeInterface(CConfig **config, CSolver***** solver, CGeomet
           auto nVar = solver[donor][INST_0][MESH_0][FLOW_SOL]->GetnVar();
           interface_type = CONSERVATIVE_VARIABLES;
           interface[donor][target] = new CConservativeVarsInterface(nVar, 0);
-          if (rank == MASTER_NODE) cout << "generic conservative variables." << endl;
+          if (rank == MASTER_NODE) cout << " Generic conservative variables." << endl;
         }
       }
 

SU2_CFD/src/drivers/CMultizoneDriver.cpp

@@ -549,19 +549,19 @@ bool CMultizoneDriver::TransferData(unsigned short donorZone, unsigned short tar
 
   /*--- Select the transfer method according to the magnitudes being transferred ---*/
 
-  auto BroadcastData = [&](int donorSol, int targetSol) {
-    interface_container[donorZone][targetZone]->BroadcastData(
-      *interpolator_container[donorZone][targetZone].get(),
-      solver_container[donorZone][INST_0][MESH_0][donorSol],
-      solver_container[targetZone][INST_0][MESH_0][targetSol],
-      geometry_container[donorZone][INST_0][MESH_0],
-      geometry_container[targetZone][INST_0][MESH_0],
-      config_container[donorZone],
-      config_container[targetZone]);
-  };
-
-  switch (interface_types[donorZone][targetZone]) {
-
+  auto HandleInterfaceType = [&] (const auto interface_type, auto* interface) {
+    auto BroadcastData = [&](int donorSol, int targetSol) {
+      interface->BroadcastData(
+        *interpolator_container[donorZone][targetZone],
+        solver_container[donorZone][INST_0][MESH_0][donorSol],
+        solver_container[targetZone][INST_0][MESH_0][targetSol],
+        geometry_container[donorZone][INST_0][MESH_0],
+        geometry_container[targetZone][INST_0][MESH_0],
+        config_container[donorZone],
+        config_container[targetZone]);
+    };
+
+    switch (interface_type) {
     case SLIDING_INTERFACE:
       BroadcastData(FLOW_SOL, FLOW_SOL);
 
@@ -598,35 +598,49 @@ bool CMultizoneDriver::TransferData(unsigned short donorZone, unsigned short tar
     case FLOW_TRACTION:
       BroadcastData(FLOW_SOL, FEA_SOL);
       break;
-    case MIXING_PLANE:
-    {
+    case MIXING_PLANE: {
       const auto nMarkerInt = config_container[donorZone]->GetnMarker_MixingPlaneInterface() / 2;
 
-      /*--- Transfer the average value from the donorZone to the targetZone ---*/
-      /*--- Loops over the mixing planes defined in the config file to find the correct mixing plane for the donor-target combination ---*/
+      /*--- Transfer the average value from the donorZone to the targetZone
+       * Loops over the mixing planes defined in the config file to find the
+       * correct mixing plane for the donor-target combination ---*/
       for (auto iMarkerInt = 1; iMarkerInt <= nMarkerInt; iMarkerInt++) {
-            interface_container[donorZone][targetZone]->AllgatherAverage(solver_container[donorZone][INST_0][MESH_0][FLOW_SOL],solver_container[targetZone][INST_0][MESH_0][FLOW_SOL],
-                geometry_container[donorZone][INST_0][MESH_0],geometry_container[targetZone][INST_0][MESH_0],
-                config_container[donorZone], config_container[targetZone], iMarkerInt );
+        interface->AllgatherAverage(
+          solver_container[donorZone][INST_0][MESH_0][FLOW_SOL],
+          solver_container[targetZone][INST_0][MESH_0][FLOW_SOL],
+          geometry_container[donorZone][INST_0][MESH_0],
+          geometry_container[targetZone][INST_0][MESH_0],
+          config_container[donorZone], config_container[targetZone], iMarkerInt);
       }
 
       /*--- Set average value donorZone->targetZone ---*/
-      interface_container[donorZone][targetZone]->SetAverageValues(solver_container[donorZone][INST_0][MESH_0][FLOW_SOL],solver_container[targetZone][INST_0][MESH_0][FLOW_SOL], donorZone);
+      interface->SetAverageValues(solver_container[donorZone][INST_0][MESH_0][FLOW_SOL],
+                                  solver_container[targetZone][INST_0][MESH_0][FLOW_SOL], donorZone);
 
       /*--- Set average geometrical properties FROM donorZone IN targetZone ---*/
-      geometry_container[targetZone][INST_0][MESH_0]->SetAvgTurboGeoValues(config_container[iZone],geometry_container[iZone][INST_0][MESH_0], iZone);
-
-      break;
-    }
+      geometry_container[targetZone][INST_0][MESH_0]->SetAvgTurboGeoValues(
+          config_container[iZone], geometry_container[iZone][INST_0][MESH_0], iZone);
+     } break;
     case NO_TRANSFER:
     case ZONES_ARE_EQUAL:
     case NO_COMMON_INTERFACE:
       break;
     default:
-      if(rank == MASTER_NODE)
+      if (rank == MASTER_NODE) {
         cout << "WARNING: One of the intended interface transfer routines is not "
              << "known to the chosen driver and has not been executed." << endl;
+      }
       break;
+    }
+  };
+
+  auto type = interface_types[donorZone][targetZone];
+  auto* interface = interface_container[donorZone][targetZone];
+
+  while (interface != nullptr) {
+    HandleInterfaceType(type, interface);
+    type = interface->NextInterfaceType;
+    interface = interface->NextInterface;
   }
 
   return UpdateMesh;

SU2_CFD/src/interfaces/CInterface.cpp

@@ -57,6 +57,8 @@ CInterface::~CInterface() {
 
   delete[] SpanValueCoeffTarget;
   delete[] SpanLevelDonor;
+
+  delete NextInterface;
 }
 
 void CInterface::BroadcastData(const CInterpolator& interpolator,

SU2_CFD/src/interfaces/fsi/CFlowTractionInterface.cpp

@@ -158,9 +158,9 @@ void CFlowTractionInterface::GetPhysical_Constants(CSolver *flow_solution, CSolv
   Physical_Constants[1] = ModAmpl;
 
   /*--- For static FSI, we cannot apply the ramp like this ---*/
-  if ((!flow_config->GetTime_Domain())){
+  if (!flow_config->GetTime_Domain()) {
     Physical_Constants[1] = 1.0;
-    if (Ramp_Load){
+    if (Ramp_Load) {
       CurrentTime = static_cast<su2double>(struct_config->GetOuterIter());
       Ramp_Time = static_cast<su2double>(struct_config->GetnIterFSI_Ramp() - 1);
 

SU2_CFD/src/solvers/CHeatSolver.cpp

@@ -257,6 +257,7 @@ void CHeatSolver::LoadRestart(CGeometry **geometry, CSolver ***solver, CConfig *
   solver[MESH_0][HEAT_SOL]->InitiateComms(geometry[MESH_0], config, MPI_QUANTITIES::SOLUTION);
   solver[MESH_0][HEAT_SOL]->CompleteComms(geometry[MESH_0], config, MPI_QUANTITIES::SOLUTION);
 
+  SU2_OMP_SAFE_GLOBAL_ACCESS(config->SetGlobalParam(MAIN_SOLVER::HEAT_EQUATION, RUNTIME_HEAT_SYS);)
   solver[MESH_0][HEAT_SOL]->Preprocessing(geometry[MESH_0], solver[MESH_0], config, MESH_0, NO_RK_ITER, RUNTIME_HEAT_SYS, false);
 
   /*--- Interpolate the solution down to the coarse multigrid levels ---*/

TestCases/fea_fsi/stat_fsi/config.cfg

@@ -1,24 +1,32 @@
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%                                                                              %
 % SU2 configuration file                                                       %
-% Case description: Fluid structure interaction - Beam in channel - 2D - FEM   %
-% Author: R.Sanchez                                                            %
-% Institution: Imperial College London                                         %
-%                                                                              %
+% Case description: Aero-thermo-elasticity                                     %
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
+% FSI settings
 SOLVER= MULTIPHYSICS
+MULTIZONE_SOLVER= BLOCK_GAUSS_SEIDEL
+RAMP_LOADING= YES
+RAMP_FSI_ITER= 5
+BGS_RELAXATION= FIXED_PARAMETER
+STAT_RELAX_PARAMETER= 0.9
 
-CONFIG_LIST = (configFlow.cfg, configFEA.cfg)
+CONFIG_LIST= ( configFlow.cfg, configFEA.cfg )
 
-MULTIZONE_SOLVER = BLOCK_GAUSS_SEIDEL
+MARKER_ZONE_INTERFACE= ( wallF, wallS )
 
-MARKER_ZONE_INTERFACE = (wallF, wallS)
+SCREEN_OUTPUT= \
+  OUTER_ITER, INNER_ITER[0], AVG_BGS_RES[0], AVG_BGS_RES[1],\
+  RMS_DENSITY[0], RMS_UTOL[1], RMS_TEMPERATURE[1], VMS[1],\
+  TOTAL_HEATFLUX[1], DEFORM_MIN_VOLUME[0], DEFORM_ITER[0]
 
-MULTIZONE_MESH = NO
-SCREEN_OUTPUT=(OUTER_ITER, BGS_DENSITY[0], AVG_BGS_RES[1], DEFORM_MIN_VOLUME[0], DEFORM_ITER[0])
 RESTART_SOL= NO
-RESTART_ITER = 0
+RESTART_ITER= 0
+
+OUTER_ITER= 30
+CONV_RESIDUAL_MINVAL= -5
+
+MULTIZONE_MESH= NO
+WRT_ZONE_CONV= NO
+OUTPUT_WRT_FREQ= 100
 
-TIME_DOMAIN = NO
-OUTER_ITER = 8