Merge pull request #2429 from su2code/coupled_thermoelasticity

Coupled thermoelasticity

Author: pcarruscag

SU2_CFD/include/output/CElasticityOutput.hpp

@@ -40,6 +40,7 @@ class CElasticityOutput final: public COutput {
   unsigned short nVar_FEM; //!< Number of FEM variables
   bool linear_analysis,    //!< Boolean indicating a linear analysis
        nonlinear_analysis, //!< Boolean indicating a nonlinear analysis
+       coupled_heat,       //!< Boolean indicating a thermoelastic analysis
        dynamic;            //!< Boolean indicating a dynamic analysis
 
 public:
@@ -82,6 +83,11 @@ class CElasticityOutput final: public COutput {
    * \param[in] config - Definition of the particular problem.
    * \return <TRUE> if the residuals should be initialized.
    */
-  bool SetInitResiduals(const CConfig *config) override ;
+  bool SetInitResiduals(const CConfig *config) override;
 
+  /*!
+   * \brief LoadSurfaceData
+   */
+  void LoadSurfaceData(CConfig *config, CGeometry *geometry, CSolver **solver, unsigned long iPoint,
+                       unsigned short iMarker, unsigned long iVertex) override;
 };

SU2_CFD/include/output/CHeatOutput.hpp

@@ -50,12 +50,22 @@ class CHeatOutput final: public CFVMOutput {
    */
   void SetHistoryOutputFields(CConfig *config) override;
 
+  /*!
+   * \brief Set the available history output fields in another output instance.
+   */
+  static void SetHistoryOutputFieldsImpl(CConfig *config, COutput* output);
+
   /*!
    * \brief Load the history output field values
    * \param[in] config - Definition of the particular problem.
    */
   void LoadHistoryData(CConfig *config, CGeometry *geometry, CSolver **solver) override;
 
+  /*!
+   * \brief Set the history output field values in another output instance.
+   */
+  static void LoadHistoryDataImpl(CConfig *config, CGeometry *geometry, CSolver **solver, COutput* output);
+
   /*!
    * \brief Set the available volume output fields
    * \param[in] config - Definition of the particular problem.

SU2_CFD/include/output/COutput.hpp

@@ -55,6 +55,7 @@ class CSolver;
 class CFileWriter;
 class CParallelDataSorter;
 class CConfig;
+class CHeatOutput;
 
 using namespace std;
 
@@ -65,6 +66,7 @@ using namespace std;
  */
 class COutput {
 protected:
+  friend class CHeatOutput;
 
   /*----------------------------- General ----------------------------*/
 

SU2_CFD/include/solvers/CFEASolver.hpp

@@ -100,6 +100,11 @@ 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 nodes to access temperature for coupled simulations.
+   */
+  const CVariable* heat_nodes = 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/include/solvers/CHeatSolver.hpp

@@ -42,7 +42,6 @@ class CHeatSolver final : public CScalarSolver<CHeatVariable> {
   static constexpr size_t MAXNVAR = 1; /*!< \brief Max number of variables, for static arrays. */
 
   const bool flow; /*!< \brief Use solver as a scalar transport equation of Temperature for the inc solver. */
-  const bool heat_equation; /*!< \brief use solver for heat conduction in solids. */
 
   su2double Global_Delta_Time = 0.0, Global_Delta_UnstTimeND = 0.0;
 

SU2_CFD/include/solvers/CScalarSolver.inl

@@ -503,10 +503,10 @@ void CScalarSolver<VariableType>::CompleteImplicitIteration(CGeometry* geometry,
       SU2_OMP_FOR_STAT(omp_chunk_size)
       for (unsigned long iPoint = 0; iPoint < nPointDomain; iPoint++) {
         /*--- Multiply the Solution var with density to get the conservative transported quantity, if necessary. ---*/
-        /* Note that for consistency with residual and jacobian calaulcations, use of current density for conservative variables 
+        /* Note that for consistency with residual and jacobian calaulcations, use of current density for conservative variables
         * of the old solution is used. see pull request https://github.com/su2code/SU2/pull/2458*/
         const su2double density = flowNodes->GetDensity(iPoint);
-        
+
         for (unsigned short iVar = 0; iVar < nVar; iVar++) {
           nodes->AddClippedSolution(iPoint, iVar, nodes->GetUnderRelaxation(iPoint) * LinSysSol(iPoint, iVar),
                                     lowerlimit[iVar], upperlimit[iVar], density, density);

SU2_CFD/src/drivers/CDriver.cpp

@@ -1527,7 +1527,9 @@ void CDriver::InitializeNumerics(CConfig *config, CGeometry **geometry, CSolver
 
     case MAIN_SOLVER::FEM_ELASTICITY:
     case MAIN_SOLVER::DISC_ADJ_FEM:
-      fem = true; break;
+      fem = true;
+      heat = config->GetWeakly_Coupled_Heat();
+      break;
 
     case MAIN_SOLVER::ADJ_EULER:
       adj_euler = euler = compressible = true; break;

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

@@ -24,18 +24,18 @@
  * You should have received a copy of the GNU Lesser General Public
  * License along with SU2. If not, see <http://www.gnu.org/licenses/>.
  */
-
-
 #include "../../include/output/CElasticityOutput.hpp"
+#include "../../include/output/CHeatOutput.hpp"
 
 #include "../../../Common/include/geometry/CGeometry.hpp"
 #include "../../include/solvers/CSolver.hpp"
 
 CElasticityOutput::CElasticityOutput(CConfig *config, unsigned short nDim) : COutput(config, nDim, false) {
 
-  linear_analysis = (config->GetGeometricConditions() == STRUCT_DEFORMATION::SMALL);
-  nonlinear_analysis = (config->GetGeometricConditions() == STRUCT_DEFORMATION::LARGE);
-  dynamic = (config->GetTime_Domain());
+  linear_analysis = config->GetGeometricConditions() == STRUCT_DEFORMATION::SMALL;
+  nonlinear_analysis = config->GetGeometricConditions() == STRUCT_DEFORMATION::LARGE;
+  coupled_heat = config->GetWeakly_Coupled_Heat();
+  dynamic = config->GetTime_Domain();
 
   /*--- Initialize number of variables ---*/
   if (linear_analysis) nVar_FEM = nDim;
@@ -48,20 +48,21 @@ CElasticityOutput::CElasticityOutput(CConfig *config, unsigned short nDim) : COu
   }
 
   /*--- Default fields for screen output ---*/
-  if (nRequestedScreenFields == 0){
+  if (nRequestedScreenFields == 0) {
     if (dynamic) requestedScreenFields.emplace_back("TIME_ITER");
     if (multiZone) requestedScreenFields.emplace_back("OUTER_ITER");
     requestedScreenFields.emplace_back("INNER_ITER");
-    if(linear_analysis){
+    if (linear_analysis) {
       requestedScreenFields.emplace_back("RMS_DISP_X");
       requestedScreenFields.emplace_back("RMS_DISP_Y");
       requestedScreenFields.emplace_back("RMS_DISP_Z");
     }
-    if(nonlinear_analysis){
+    if (nonlinear_analysis) {
       requestedScreenFields.emplace_back("RMS_UTOL");
       requestedScreenFields.emplace_back("RMS_RTOL");
       requestedScreenFields.emplace_back("RMS_ETOL");
     }
+    if (coupled_heat) requestedScreenFields.emplace_back("RMS_TEMPERATURE");
     requestedScreenFields.emplace_back("VMS");
     nRequestedScreenFields = requestedScreenFields.size();
   }
@@ -71,11 +72,8 @@ CElasticityOutput::CElasticityOutput(CConfig *config, unsigned short nDim) : COu
     requestedVolumeFields.emplace_back("COORDINATES");
     requestedVolumeFields.emplace_back("SOLUTION");
     requestedVolumeFields.emplace_back("STRESS");
-    if (dynamic) {
-      requestedVolumeFields.emplace_back("VELOCITY");
-      requestedVolumeFields.emplace_back("ACCELERATION");
-    }
     if (config->GetTopology_Optimization()) requestedVolumeFields.emplace_back("TOPOLOGY");
+    if (coupled_heat) requestedVolumeFields.emplace_back("PRIMITIVE");
     nRequestedVolumeFields = requestedVolumeFields.size();
   }
 
@@ -106,6 +104,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 ---*/
@@ -147,6 +146,13 @@ void CElasticityOutput::LoadHistoryData(CConfig *config, CGeometry *geometry, CS
     SetHistoryOutputValue("TOPOL_DISCRETENESS", fea_solver->GetTotal_OFDiscreteness());
   }
 
+  /*--- Add heat solver data if available. ---*/
+  if (coupled_heat) {
+    CHeatOutput::LoadHistoryDataImpl(config, geometry, solver, this);
+    SetHistoryOutputValue("LINSOL_ITER_HEAT", heat_solver->GetIterLinSolver());
+    SetHistoryOutputValue("LINSOL_RESIDUAL_HEAT", log10(heat_solver->GetResLinSolver()));
+  }
+
   ComputeSimpleCustomOutputs(config);
 
   /*--- Keep this as last, since it uses the history values that were set. ---*/
@@ -189,6 +195,11 @@ void CElasticityOutput::SetHistoryOutputFields(CConfig *config) {
   }
   AddHistoryOutput("COMBO", "ComboObj", ScreenOutputFormat::SCIENTIFIC, "COMBO", "Combined obj. function value.", HistoryFieldType::COEFFICIENT);
 
+  if (coupled_heat) {
+    CHeatOutput::SetHistoryOutputFieldsImpl(config, this);
+    AddHistoryOutput("LINSOL_ITER_HEAT", "LinSolIterHeat", ScreenOutputFormat::INTEGER, "LINSOL", "Number of iterations of the linear solver.");
+    AddHistoryOutput("LINSOL_RESIDUAL_HEAT", "LinSolResHeat", ScreenOutputFormat::FIXED, "LINSOL", "Residual of the linear solver.");
+  }
 }
 
 void CElasticityOutput::LoadVolumeData(CConfig *config, CGeometry *geometry, CSolver **solver, unsigned long iPoint){
@@ -214,6 +225,10 @@ void CElasticityOutput::LoadVolumeData(CConfig *config, CGeometry *geometry, CSo
     SetVolumeOutputValue("ACCELERATION-Y", iPoint, Node_Struc->GetSolution_Accel(iPoint, 1));
     if (nDim == 3) SetVolumeOutputValue("ACCELERATION-Z", iPoint, Node_Struc->GetSolution_Accel(iPoint, 2));
   }
+  if (coupled_heat) {
+    CVariable* Node_Heat = solver[HEAT_SOL]->GetNodes();
+    SetVolumeOutputValue("TEMPERATURE", iPoint, Node_Heat->GetSolution(iPoint, 0));
+  }
 
   SetVolumeOutputValue("STRESS-XX", iPoint, Node_Struc->GetStress_FEM(iPoint)[0]);
   SetVolumeOutputValue("STRESS-YY", iPoint, Node_Struc->GetStress_FEM(iPoint)[1]);
@@ -228,6 +243,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){
@@ -251,6 +274,10 @@ void CElasticityOutput::SetVolumeOutputFields(CConfig *config){
     if (nDim == 3) AddVolumeOutput("ACCELERATION-Z", "Acceleration_z", "SOLUTION", "z-component of the acceleration vector");
   }
 
+  if (coupled_heat) {
+    AddVolumeOutput("TEMPERATURE", "Temperature", "SOLUTION", "Temperature");
+  }
+
   AddVolumeOutput("STRESS-XX",    "Sxx", "STRESS", "x-component of the normal stress vector");
   AddVolumeOutput("STRESS-YY",    "Syy", "STRESS", "y-component of the normal stress vector");
   AddVolumeOutput("STRESS-XY",    "Sxy", "STRESS", "xy shear stress component");
@@ -266,10 +293,25 @@ void CElasticityOutput::SetVolumeOutputFields(CConfig *config){
   if (config->GetTopology_Optimization()) {
     AddVolumeOutput("TOPOL_DENSITY", "Topology_Density", "TOPOLOGY", "filtered topology density");
   }
+
+  if (coupled_heat) {
+    AddVolumeOutput("HEAT_FLUX", "Heat_Flux", "PRIMITIVE", "Heatflux");
+    AddVolumeOutput("RES_TEMPERATURE", "Residual_Temperature", "RESIDUAL", "Residual of the temperature");
+  }
+
 }
 
 bool CElasticityOutput::SetInitResiduals(const CConfig *config){
 
   return (config->GetTime_Domain() == NO && (curInnerIter  == 0));
 
 }
+
+void CElasticityOutput::LoadSurfaceData(CConfig *config, CGeometry *geometry, CSolver **solver, unsigned long iPoint,
+                                        unsigned short iMarker, unsigned long iVertex) {
+  if (!coupled_heat || !config->GetViscous_Wall(iMarker)) return;
+
+  /* Heat flux value at each surface grid node. */
+  SetVolumeOutputValue("HEAT_FLUX", iPoint, solver[HEAT_SOL]->GetHeatFlux(iMarker, iVertex));
+
+}

SU2_CFD/src/output/CHeatOutput.cpp

@@ -72,42 +72,52 @@ CHeatOutput::CHeatOutput(CConfig *config, unsigned short nDim) : CFVMOutput(conf
 
 }
 
+void CHeatOutput::LoadHistoryDataImpl(CConfig *config, CGeometry *geometry, CSolver **solver, COutput* output) {
+
+  CSolver* heat_solver = solver[HEAT_SOL];
+
+  output->SetHistoryOutputValue("TOTAL_HEATFLUX", heat_solver->GetTotal_HeatFlux());
+  output->SetHistoryOutputValue("AVG_TEMPERATURE", heat_solver->GetTotal_AvgTemperature());
+  output->SetHistoryOutputValue("RMS_TEMPERATURE", log10(heat_solver->GetRes_RMS(0)));
+  output->SetHistoryOutputValue("MAX_TEMPERATURE", log10(heat_solver->GetRes_Max(0)));
+  if (config->GetMultizone_Problem()) {
+    output->SetHistoryOutputValue("BGS_TEMPERATURE", log10(heat_solver->GetRes_BGS(0)));
+  }
+  output->SetHistoryOutputValue("CFL_NUMBER", config->GetCFL(MESH_0));
+}
+
 void CHeatOutput::LoadHistoryData(CConfig *config, CGeometry *geometry, CSolver **solver) {
 
   CSolver* heat_solver = solver[HEAT_SOL];
 
-  SetHistoryOutputValue("TOTAL_HEATFLUX", heat_solver->GetTotal_HeatFlux());
-  SetHistoryOutputValue("MAXIMUM_HEATFLUX", heat_solver->GetTotal_MaxHeatFlux());
-  SetHistoryOutputValue("AVG_TEMPERATURE", heat_solver->GetTotal_AvgTemperature());
-  SetHistoryOutputValue("RMS_TEMPERATURE", log10(heat_solver->GetRes_RMS(0)));
-  SetHistoryOutputValue("MAX_TEMPERATURE", log10(heat_solver->GetRes_Max(0)));
-  if (multiZone)
-    SetHistoryOutputValue("BGS_TEMPERATURE", log10(heat_solver->GetRes_BGS(0)));
+  LoadHistoryDataImpl(config, geometry, solver, this);
 
   SetHistoryOutputValue("LINSOL_ITER", heat_solver->GetIterLinSolver());
   SetHistoryOutputValue("LINSOL_RESIDUAL", log10(heat_solver->GetResLinSolver()));
-  SetHistoryOutputValue("CFL_NUMBER", config->GetCFL(MESH_0));
 
   ComputeSimpleCustomOutputs(config);
 
   /*--- Keep this as last, since it uses the history values that were set. ---*/
   SetCustomAndComboObjectives(HEAT_SOL, config, solver);
 }
 
+void CHeatOutput::SetHistoryOutputFieldsImpl(CConfig *config, COutput* output) {
 
-void CHeatOutput::SetHistoryOutputFields(CConfig *config){
+  output->AddHistoryOutput("RMS_TEMPERATURE", "rms[T]", ScreenOutputFormat::FIXED, "RMS_RES", "Root mean square residual of the temperature", HistoryFieldType::RESIDUAL);
+  output->AddHistoryOutput("MAX_TEMPERATURE", "max[T]", ScreenOutputFormat::FIXED, "MAX_RES", "Maximum residual of the temperature", HistoryFieldType::RESIDUAL);
+  output->AddHistoryOutput("BGS_TEMPERATURE", "bgs[T]", ScreenOutputFormat::FIXED, "BGS_RES", "Block-Gauss-Seidel residual of the temperature", HistoryFieldType::RESIDUAL);
 
-  AddHistoryOutput("LINSOL_ITER", "LinSolIter", ScreenOutputFormat::INTEGER, "LINSOL", "Number of iterations of the linear solver.");
-  AddHistoryOutput("LINSOL_RESIDUAL", "LinSolRes", ScreenOutputFormat::FIXED, "LINSOL", "Residual of the linear solver.");
+  output->AddHistoryOutput("TOTAL_HEATFLUX", "HF", ScreenOutputFormat::SCIENTIFIC, "HEAT", "Total heatflux on all surfaces defined in MARKER_MONITORING", HistoryFieldType::COEFFICIENT);
+  output->AddHistoryOutput("AVG_TEMPERATURE", "AvgTemp", ScreenOutputFormat::SCIENTIFIC, "HEAT", "Average temperature on all surfaces defined in MARKER_MONITORING", HistoryFieldType::COEFFICIENT);
+  output->AddHistoryOutput("CFL_NUMBER", "CFL number", ScreenOutputFormat::SCIENTIFIC, "CFL_NUMBER", "Current value of the CFL number");
+}
 
-  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("BGS_TEMPERATURE", "bgs[T]", ScreenOutputFormat::FIXED, "BGS_RES", "Block-Gauss-Seidel residual of the temperature", HistoryFieldType::RESIDUAL);
+void CHeatOutput::SetHistoryOutputFields(CConfig *config) {
 
-  AddHistoryOutput("TOTAL_HEATFLUX", "HF", ScreenOutputFormat::SCIENTIFIC, "HEAT", "Total heatflux on all surfaces defined in MARKER_MONITORING", HistoryFieldType::COEFFICIENT);
-  AddHistoryOutput("MAXIMUM_HEATFLUX", "MaxHF", ScreenOutputFormat::SCIENTIFIC, "HEAT", "Maximum heatflux on all surfaces defined in MARKER_MONITORING", HistoryFieldType::COEFFICIENT);
-  AddHistoryOutput("AVG_TEMPERATURE", "AvgTemp", ScreenOutputFormat::SCIENTIFIC, "HEAT", "Average temperature on all surfaces defined in MARKER_MONITORING", HistoryFieldType::COEFFICIENT);
-  AddHistoryOutput("CFL_NUMBER", "CFL number", ScreenOutputFormat::SCIENTIFIC, "CFL_NUMBER", "Current value of the CFL number");
+  SetHistoryOutputFieldsImpl(config, this);
+
+  AddHistoryOutput("LINSOL_ITER", "LinSolIter", ScreenOutputFormat::INTEGER, "LINSOL", "Number of iterations of the linear solver.");
+  AddHistoryOutput("LINSOL_RESIDUAL", "LinSolRes", ScreenOutputFormat::FIXED, "LINSOL", "Residual of the linear solver.");
 
   AddHistoryOutput("COMBO", "ComboObj", ScreenOutputFormat::SCIENTIFIC, "COMBO", "Combined obj. function value.", HistoryFieldType::COEFFICIENT);
 }