Merge pull request #2404 from Vaish-W/coupled_thermoelasticity

[WIP] Coupled solver for thermoelasticity

Author: pcarruscag

SU2_CFD/include/solvers/CFEASolver.hpp

@@ -100,6 +100,22 @@ class CFEASolver : public CFEASolverBase {
   bool initial_calc = true;    /*!< \brief Becomes false after first call to Preprocessing. */
   bool body_forces = false;    /*!< \brief Whether any body force is active. */
 
+  /*!
+ * \brief Pointer to the heat solver for coupled simulations.
+ *
+ * This member stores a pointer to the heat solver, which handles
+ * the solution of the heat equation in weakly coupled simulations.
+ * It is initialized during the preprocessing step if the configuration
+ * enables the weak coupling of heat and elasticity solvers. This solver
+ * provides temperature information to the finite element elasticity solver
+ * and contributes to the coupled residuals.
+ *
+ * The `heat_solver` pointer remains `nullptr` when the heat solver is not enabled
+ * in the configuration. Memory for the heat solver is dynamically allocated
+ * during initialization and released in the destructor to avoid memory leaks.
+ */
+  CSolver* heat_solver = nullptr;
+
   /*!
    * \brief The highest level in the variable hierarchy this solver can safely use,
    * CVariable is the common denominator between the FEA and Mesh deformation variables.

SU2_CFD/src/iteration/CFEAIteration.cpp

@@ -48,6 +48,13 @@ void CFEAIteration::Iterate(COutput* output, CIntegration**** integration, CGeom
   CIntegration* feaIntegration = integration[val_iZone][val_iInst][FEA_SOL];
   CSolver* feaSolver = solver[val_iZone][val_iInst][MESH_0][FEA_SOL];
 
+  /*--- Add heat solver integration step ---*/
+    if (config[val_iZone]->GetWeakly_Coupled_Heat()) {
+    config[val_iZone]->SetGlobalParam(MAIN_SOLVER::HEAT_EQUATION, RUNTIME_HEAT_SYS);
+    integration[val_iZone][val_iInst][HEAT_SOL]->SingleGrid_Iteration(
+        geometry, solver, numerics, config, RUNTIME_HEAT_SYS, val_iZone, val_iInst);
+    }
+
   /*--- FEA equations ---*/
   config[val_iZone]->SetGlobalParam(MAIN_SOLVER::FEM_ELASTICITY, RUNTIME_FEA_SYS);
 

SU2_CFD/src/output/CElasticityOutput.cpp

@@ -106,6 +106,7 @@ CElasticityOutput::CElasticityOutput(CConfig *config, unsigned short nDim) : COu
 void CElasticityOutput::LoadHistoryData(CConfig *config, CGeometry *geometry, CSolver **solver)  {
 
   CSolver* fea_solver = solver[FEA_SOL];
+  CSolver* heat_solver = solver[HEAT_SOL];
 
   /*--- Residuals: ---*/
   /*--- Linear analysis: RMS of the displacements in the nDim coordinates ---*/
@@ -152,6 +153,15 @@ void CElasticityOutput::LoadHistoryData(CConfig *config, CGeometry *geometry, CS
   /*--- Keep this as last, since it uses the history values that were set. ---*/
   SetCustomAndComboObjectives(FEA_SOL, config, solver);
 
+  /*--- Add heat solver data if available ---*/
+  if (heat_solver) {
+    SetHistoryOutputValue("RMS_TEMPERATURE", log10(heat_solver->GetRes_RMS(0)));
+    SetHistoryOutputValue("MAX_TEMPERATURE", log10(heat_solver->GetRes_Max(0)));
+    SetHistoryOutputValue("AVG_TEMPERATURE", heat_solver->GetTotal_AvgTemperature());
+    SetHistoryOutputValue("TOTAL_HEATFLUX", heat_solver->GetTotal_HeatFlux());
+    SetHistoryOutputValue("MAXIMUM_HEATFLUX", heat_solver->GetTotal_MaxHeatFlux());
+  }
+
 }
 
 void CElasticityOutput::SetHistoryOutputFields(CConfig *config) {
@@ -189,6 +199,14 @@ void CElasticityOutput::SetHistoryOutputFields(CConfig *config) {
   }
   AddHistoryOutput("COMBO", "ComboObj", ScreenOutputFormat::SCIENTIFIC, "COMBO", "Combined obj. function value.", HistoryFieldType::COEFFICIENT);
 
+  if (config->GetWeakly_Coupled_Heat()) {
+    AddHistoryOutput("RMS_TEMPERATURE", "rms[T]", ScreenOutputFormat::FIXED, "RMS_RES", "Root mean square residual of the temperature", HistoryFieldType::RESIDUAL);
+    AddHistoryOutput("MAX_TEMPERATURE", "max[T]", ScreenOutputFormat::FIXED, "MAX_RES", "Maximum residual of the temperature", HistoryFieldType::RESIDUAL);
+    AddHistoryOutput("AVG_TEMPERATURE", "avg[T]", ScreenOutputFormat::SCIENTIFIC, "HEAT", "Average temperature on all surfaces", HistoryFieldType::COEFFICIENT);
+    AddHistoryOutput("TOTAL_HEATFLUX", "HF", ScreenOutputFormat::SCIENTIFIC, "HEAT", "Total heat flux on all surfaces", HistoryFieldType::COEFFICIENT);
+    AddHistoryOutput("MAXIMUM_HEATFLUX", "MaxHF", ScreenOutputFormat::SCIENTIFIC, "HEAT", "Maximum heat flux on all surfaces", HistoryFieldType::COEFFICIENT);
+  }
+
 }
 
 void CElasticityOutput::LoadVolumeData(CConfig *config, CGeometry *geometry, CSolver **solver, unsigned long iPoint){
@@ -228,6 +246,14 @@ void CElasticityOutput::LoadVolumeData(CConfig *config, CGeometry *geometry, CSo
   if (config->GetTopology_Optimization()) {
     SetVolumeOutputValue("TOPOL_DENSITY", iPoint, Node_Struc->GetAuxVar(iPoint));
   }
+
+  CSolver* heat_solver = solver[HEAT_SOL];
+  if (heat_solver) {
+    const auto Node_Heat = heat_solver->GetNodes();
+    SetVolumeOutputValue("TEMPERATURE", iPoint, Node_Heat->GetSolution(iPoint, 0));
+    SetVolumeOutputValue("RES_TEMPERATURE", iPoint, heat_solver->LinSysRes(iPoint, 0));
+  }
+
 }
 
 void CElasticityOutput::SetVolumeOutputFields(CConfig *config){
@@ -266,6 +292,12 @@ void CElasticityOutput::SetVolumeOutputFields(CConfig *config){
   if (config->GetTopology_Optimization()) {
     AddVolumeOutput("TOPOL_DENSITY", "Topology_Density", "TOPOLOGY", "filtered topology density");
   }
+
+  if (config->GetWeakly_Coupled_Heat()) {
+    AddVolumeOutput("TEMPERATURE", "Temperature", "SOLUTION", "Temperature");
+    AddVolumeOutput("RES_TEMPERATURE", "Residual_Temperature", "RESIDUAL", "Residual of the temperature");
+  }
+
 }
 
 bool CElasticityOutput::SetInitResiduals(const CConfig *config){

SU2_CFD/src/solvers/CFEASolver.cpp

@@ -30,6 +30,7 @@
 #include "../../include/numerics/elasticity/CFEAElasticity.hpp"
 #include "../../../Common/include/toolboxes/printing_toolbox.hpp"
 #include "../../../Common/include/toolboxes/geometry_toolbox.hpp"
+#include "../../include/solvers/CHeatSolver.hpp"
 #include <algorithm>
 
 using namespace GeometryToolbox;
@@ -563,6 +564,10 @@ void CFEASolver::Preprocessing(CGeometry *geometry, CSolver **solver_container,
   const bool disc_adj_fem = (config->GetKind_Solver() == MAIN_SOLVER::DISC_ADJ_FEM);
   const bool topology_mode = config->GetTopology_Optimization();
 
+  if (config->GetWeakly_Coupled_Heat()) {
+        heat_solver = solver_container[HEAT_SOL];
+    }
+
   /*
    * For topology optimization we apply a filter on the design density field to avoid
    * numerical issues (checkerboards), ensure mesh independence, and impose a length scale.
@@ -687,6 +692,10 @@ void CFEASolver::Compute_StiffMatrix(CGeometry *geometry, CNumerics **numerics,
             su2double val_Sol = nodes->GetSolution(indexNode[iNode],iDim) + val_Coord;
             element->SetRef_Coord(iNode, iDim, val_Coord);
             element->SetCurr_Coord(iNode, iDim, val_Sol);
+          }  
+          if (heat_solver) {
+                      const su2double temperature = heat_solver->GetNodes()->GetSolution(indexNode[iNode], 0);
+                      element->SetTemperature(iNode, temperature);
           }
         }
 
@@ -795,6 +804,10 @@ void CFEASolver::Compute_StiffMatrix_NodalStressRes(CGeometry *geometry, CNumeri
               de_elem->SetRef_Coord(iNode, iDim, val_Coord);
             }
           }
+          if (heat_solver) {
+            const su2double temperature = heat_solver->GetNodes()->GetSolution(indexNode[iNode], 0);
+            fea_elem->SetTemperature(iNode, temperature);
+          }
         }
 
         /*--- In topology mode determine the penalty to apply to the stiffness. ---*/
@@ -1093,6 +1106,10 @@ void CFEASolver::Compute_NodalStressRes(CGeometry *geometry, CNumerics **numeric
             element->SetCurr_Coord(iNode, iDim, val_Sol);
             element->SetRef_Coord(iNode, iDim, val_Coord);
           }
+          if (heat_solver) {
+            const su2double temperature = heat_solver->GetNodes()->GetSolution(indexNode[iNode], 0);
+            element->SetTemperature(iNode, temperature);
+          }
         }
 
         /*--- In topology mode determine the penalty to apply to the stiffness ---*/
@@ -1209,6 +1226,10 @@ void CFEASolver::Compute_NodalStress(CGeometry *geometry, CNumerics **numerics,
             element->SetCurr_Coord(iNode, iDim, val_Sol);
             element->SetRef_Coord(iNode, iDim, val_Coord);
           }
+          if (heat_solver) {
+            const su2double temperature = heat_solver->GetNodes()->GetSolution(indexNode[iNode], 0);
+            element->SetTemperature(iNode, temperature);
+          }
         }
 
         /*--- In topology mode determine the penalty to apply to the stiffness ---*/

SU2_CFD/src/solvers/CSolverFactory.cpp

@@ -173,6 +173,9 @@ CSolver** CSolverFactory::CreateSolverContainer(MAIN_SOLVER kindMainSolver, CCon
       break;
     case MAIN_SOLVER::FEM_ELASTICITY:
       solver[FEA_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::FEA, solver, geometry, config, iMGLevel);
+      if (config->GetWeakly_Coupled_Heat()) {
+        solver[HEAT_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::HEAT, solver, geometry, config, iMGLevel);
+      }
       break;
     case MAIN_SOLVER::DISC_ADJ_FEM:
       solver[FEA_SOL]    = CreateSubSolver(SUB_SOLVER_TYPE::FEA, solver, geometry, config, iMGLevel);