Merge pull request #2736 from su2code/pedro/adjoint_thermoelastic

Adjoint thermo elasticity

Author: pcarruscag

SU2_CFD/include/iteration/CDiscAdjFEAIteration.hpp

@@ -42,6 +42,9 @@ class CDiscAdjFEAIteration final : public CIteration {
  private:
   unsigned short CurrentRecording; /*!< \brief Stores the current status of the recording. */
 
+  /*!< \brief Used for thermo elastic problems to avoid duplicating code. */
+  CIteration* DiscAdjHeatIteration = nullptr;
+
   /*!
    * \brief load solution for dynamic problems
    * \param[in] geometry - Geometrical definition of the problem.
@@ -64,7 +67,7 @@ class CDiscAdjFEAIteration final : public CIteration {
   /*!
    * \brief Destructor of the class.
    */
-  ~CDiscAdjFEAIteration(void) override;
+  ~CDiscAdjFEAIteration() override;
 
   /*!
    * \brief Preprocessing to prepare for an iteration of the physics.

SU2_CFD/include/output/CAdjElasticityOutput.hpp

@@ -36,6 +36,7 @@
  */
 class CAdjElasticityOutput final : public COutput {
 private:
+  bool coupled_heat; //!< Boolean indicating a thermoelastic analysis
   unsigned short nVar_FEM; //!< Number of FEM variables
 
 public:

SU2_CFD/include/output/CAdjHeatOutput.hpp

@@ -46,14 +46,24 @@ class CAdjHeatOutput final: public COutput {
   /*!
    * \brief Destructor of the class.
    */
-   ~CAdjHeatOutput(void) override;
+   ~CAdjHeatOutput() override;
+
+  /*!
+   * \brief Set the history output field values in another output instance.
+   */
+  static void LoadHistoryDataImpl(CConfig *config, CGeometry *geometry, CSolver **solver, 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 available history output fields in another output instance.
+   */
+  static void SetHistoryOutputFieldsImpl(CConfig *config, COutput* output);
+
   /*!
    * \brief Set the available history output fields
    * \param[in] config - Definition of the particular problem.

SU2_CFD/include/output/CElasticityOutput.hpp

@@ -37,7 +37,6 @@
 class CElasticityOutput final: public COutput {
 protected:
 
-  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

SU2_CFD/include/output/COutput.hpp

@@ -56,6 +56,7 @@ class CFileWriter;
 class CParallelDataSorter;
 class CConfig;
 class CHeatOutput;
+class CAdjHeatOutput;
 
 using namespace std;
 
@@ -67,6 +68,7 @@ using namespace std;
 class COutput {
 protected:
   friend class CHeatOutput;
+  friend class CAdjHeatOutput;
 
   /*----------------------------- General ----------------------------*/
 

SU2_CFD/include/solvers/CFEASolver.hpp

@@ -550,7 +550,7 @@ class CFEASolver : public CFEASolverBase {
    * \param[in] config - Definition of the problem.
    * \param[in] solver - Container vector with all the solutions.
    */
-  void Evaluate_ObjFunc(const CConfig *config, CSolver**) final {
+  void Evaluate_ObjFunc(const CConfig *config, CSolver** solver) final {
     Total_ComboObj = 0.0;
     switch (config->GetKind_ObjFunc()) {
       case REFERENCE_GEOMETRY:
@@ -575,6 +575,10 @@ class CFEASolver : public CFEASolverBase {
         Total_ComboObj = Total_Custom_ObjFunc;
         break;
     }
+    if (config->GetWeakly_Coupled_Heat()) {
+      solver[HEAT_SOL]->Evaluate_ObjFunc(config, solver);
+      Total_ComboObj += solver[HEAT_SOL]->GetTotal_ComboObj();
+    }
   }
 
   /*!

SU2_CFD/src/drivers/CDriver.cpp

@@ -2769,9 +2769,10 @@ void CDriver::PrintDirectResidual(RECORDING kind_recording) {
 
   /*--- Helper lambda func to return lenghty [iVar][iZone] string.  ---*/
   auto iVar_iZone2string = [&](unsigned short ivar, unsigned short izone) {
-    if (multizone)
+    if (multizone) {
       return "[" + std::to_string(ivar) + "][" + std::to_string(izone) + "]";
-          return "[" + std::to_string(ivar) + "]";
+    }
+    return "[" + std::to_string(ivar) + "]";
   };
 
   /*--- Print residuals in the first iteration ---*/
@@ -2827,37 +2828,33 @@ void CDriver::PrintDirectResidual(RECORDING kind_recording) {
           if (!addVals) RMSTable.AddColumn("rms_Rad" + iVar_iZone2string(0, iZone), fieldWidth);
           else RMSTable << log10(solvers[RAD_SOL]->GetRes_RMS(0));
         }
-
-      }
-      else if (configs->GetStructuralProblem()) {
-
+      } else if (configs->GetStructuralProblem()) {
         if (configs->GetGeometricConditions() == STRUCT_DEFORMATION::LARGE){
           if (!addVals) {
             RMSTable.AddColumn("UTOL-A", fieldWidth);
             RMSTable.AddColumn("RTOL-A", fieldWidth);
             RMSTable.AddColumn("ETOL-A", fieldWidth);
-          }
-          else {
+          } else {
             RMSTable << log10(solvers[FEA_SOL]->GetRes_FEM(0))
                      << log10(solvers[FEA_SOL]->GetRes_FEM(1))
                      << log10(solvers[FEA_SOL]->GetRes_FEM(2));
           }
-        }
-        else{
+        } else {
           if (!addVals) {
             RMSTable.AddColumn("log10[RMS Ux]", fieldWidth);
             RMSTable.AddColumn("log10[RMS Uy]", fieldWidth);
             if (nDim == 3) RMSTable.AddColumn("log10[RMS Uz]", fieldWidth);
-          }
-          else {
+          } else {
             RMSTable << log10(solvers[FEA_SOL]->GetRes_FEM(0))
                      << log10(solvers[FEA_SOL]->GetRes_FEM(1));
             if (nDim == 3) RMSTable << log10(solvers[FEA_SOL]->GetRes_FEM(2));
           }
         }
-
-      }
-      else if (configs->GetHeatProblem()) {
+        if (configs->GetWeakly_Coupled_Heat()){
+          if (!addVals) RMSTable.AddColumn("rms_Heat", fieldWidth);
+          else RMSTable << log10(solvers[HEAT_SOL]->GetRes_RMS(0));
+        }
+      } else if (configs->GetHeatProblem()) {
 
         if (!addVals) RMSTable.AddColumn("rms_Heat" + iVar_iZone2string(0, iZone), fieldWidth);
         else RMSTable << log10(solvers[HEAT_SOL]->GetRes_RMS(0));

SU2_CFD/src/iteration/CDiscAdjFEAIteration.cpp

@@ -26,13 +26,20 @@
  */
 
 #include "../../include/iteration/CDiscAdjFEAIteration.hpp"
+#include "../../include/iteration/CDiscAdjHeatIteration.hpp"
 #include "../../include/iteration/CFEAIteration.hpp"
 #include "../../include/solvers/CFEASolver.hpp"
 #include "../../include/output/COutput.hpp"
 
-CDiscAdjFEAIteration::CDiscAdjFEAIteration(const CConfig *config) : CIteration(config), CurrentRecording(NONE) {}
+CDiscAdjFEAIteration::CDiscAdjFEAIteration(const CConfig *config) : CIteration(config), CurrentRecording(NONE) {
+  if (config->GetWeakly_Coupled_Heat()) {
+    DiscAdjHeatIteration = new CDiscAdjHeatIteration(config);
+  }
+}
 
-CDiscAdjFEAIteration::~CDiscAdjFEAIteration() {}
+CDiscAdjFEAIteration::~CDiscAdjFEAIteration() {
+  delete DiscAdjHeatIteration;
+}
 
 void CDiscAdjFEAIteration::Preprocess(COutput* output, CIntegration**** integration, CGeometry**** geometry,
                                       CSolver***** solver, CNumerics****** numerics, CConfig** config,
@@ -90,6 +97,10 @@ void CDiscAdjFEAIteration::Preprocess(COutput* output, CIntegration**** integrat
 
   solvers0[ADJFEA_SOL]->Preprocessing(geometry0, solvers0, config[iZone], MESH_0, 0, RUNTIME_ADJFEA_SYS, false);
 
+  if (DiscAdjHeatIteration) {
+    DiscAdjHeatIteration->Preprocess(output, integration, geometry, solver, numerics, config, surface_movement,
+                                     grid_movement, FFDBox, iZone, iInst);
+  }
 }
 
 void CDiscAdjFEAIteration::LoadDynamic_Solution(CGeometry**** geometry, CSolver***** solver, CConfig** config,
@@ -115,23 +126,23 @@ void CDiscAdjFEAIteration::IterateDiscAdj(CGeometry**** geometry, CSolver***** s
                                           unsigned short iZone, unsigned short iInst, bool CrossTerm) {
 
   /*--- Extract the adjoints of the conservative input variables and store them for the next iteration ---*/
-
-  solver[iZone][iInst][MESH_0][ADJFEA_SOL]->ExtractAdjoint_Solution(geometry[iZone][iInst][MESH_0], config[iZone],
-                                                                    CrossTerm);
-
-  solver[iZone][iInst][MESH_0][ADJFEA_SOL]->ExtractAdjoint_Variables(geometry[iZone][iInst][MESH_0], config[iZone]);
+  for (const auto iSol : {ADJFEA_SOL, ADJHEAT_SOL}) {
+    if (auto* sol = solver[iZone][iInst][MESH_0][iSol]; sol != nullptr) {
+      sol->ExtractAdjoint_Solution(geometry[iZone][iInst][MESH_0], config[iZone], CrossTerm);
+      sol->ExtractAdjoint_Variables(geometry[iZone][iInst][MESH_0], config[iZone]);
+    }
+  }
 }
 
 void CDiscAdjFEAIteration::RegisterInput(CSolver***** solver, CGeometry**** geometry, CConfig** config,
                                          unsigned short iZone, unsigned short iInst, RECORDING kind_recording) {
   if (kind_recording != RECORDING::MESH_COORDS) {
-    /*--- Register structural displacements as input ---*/
-
-    solver[iZone][iInst][MESH_0][ADJFEA_SOL]->RegisterSolution(geometry[iZone][iInst][MESH_0], config[iZone]);
-
-    /*--- Register variables as input ---*/
-
-    solver[iZone][iInst][MESH_0][ADJFEA_SOL]->RegisterVariables(geometry[iZone][iInst][MESH_0], config[iZone]);
+    for (const auto iSol : {ADJFEA_SOL, ADJHEAT_SOL}) {
+      if (auto* sol = solver[iZone][iInst][MESH_0][iSol]; sol != nullptr) {
+        sol->RegisterSolution(geometry[iZone][iInst][MESH_0], config[iZone]);
+        sol->RegisterVariables(geometry[iZone][iInst][MESH_0], config[iZone]);
+      }
+    }
   } else {
     /*--- Register topology optimization densities (note direct solver) ---*/
 
@@ -146,6 +157,10 @@ void CDiscAdjFEAIteration::RegisterInput(CSolver***** solver, CGeometry**** geom
 void CDiscAdjFEAIteration::SetDependencies(CSolver***** solver, CGeometry**** geometry, CNumerics****** numerics,
                                            CConfig** config, unsigned short iZone, unsigned short iInst,
                                            RECORDING kind_recording) {
+  if (DiscAdjHeatIteration) {
+    DiscAdjHeatIteration->SetDependencies(solver, geometry, numerics, config, iZone, iInst, kind_recording);
+  }
+
   auto dir_solver = solver[iZone][iInst][MESH_0][FEA_SOL];
   auto adj_solver = solver[iZone][iInst][MESH_0][ADJFEA_SOL];
   auto structural_geometry = geometry[iZone][iInst][MESH_0];
@@ -263,18 +278,25 @@ void CDiscAdjFEAIteration::SetDependencies(CSolver***** solver, CGeometry**** ge
 
 void CDiscAdjFEAIteration::RegisterOutput(CSolver***** solver, CGeometry**** geometry, CConfig** config,
                                           unsigned short iZone, unsigned short iInst) {
-  /*--- Register conservative variables as output of the iteration ---*/
-
-  solver[iZone][iInst][MESH_0][ADJFEA_SOL]->RegisterOutput(geometry[iZone][iInst][MESH_0], config[iZone]);
+  /*--- Register solution variables as output of the iteration. ---*/
+  for (const auto iSol : {ADJFEA_SOL, ADJHEAT_SOL}) {
+    if (auto* sol = solver[iZone][iInst][MESH_0][iSol]; sol != nullptr) {
+      sol->RegisterOutput(geometry[iZone][iInst][MESH_0], config[iZone]);
+    }
+  }
 }
 
 void CDiscAdjFEAIteration::InitializeAdjoint(CSolver***** solver, CGeometry**** geometry, CConfig** config,
                                              unsigned short iZone, unsigned short iInst) {
-  /*--- Initialize the adjoints the conservative variables ---*/
+  /*--- Initialize the adjoints of the solution variables. ---*/
 
   AD::ResizeAdjoints();
   AD::BeginUseAdjoints();
-  solver[iZone][iInst][MESH_0][ADJFEA_SOL]->SetAdjoint_Output(geometry[iZone][iInst][MESH_0], config[iZone]);
+  for (const auto iSol : {ADJFEA_SOL, ADJHEAT_SOL}) {
+    if (auto* sol = solver[iZone][iInst][MESH_0][iSol]; sol != nullptr) {
+      sol->SetAdjoint_Output(geometry[iZone][iInst][MESH_0], config[iZone]);
+    }
+  }
   AD::EndUseAdjoints();
 }
 
@@ -285,8 +307,8 @@ bool CDiscAdjFEAIteration::Monitor(COutput* output, CIntegration**** integration
   /*--- Write the convergence history (only screen output) ---*/
 
   output->SetHistoryOutput(geometry[iZone][INST_0][MESH_0], solver[iZone][INST_0][MESH_0], config[iZone],
-                            config[iZone]->GetTimeIter(), config[iZone]->GetOuterIter(),
-                            config[iZone]->GetInnerIter());
+                           config[iZone]->GetTimeIter(), config[iZone]->GetOuterIter(),
+                           config[iZone]->GetInnerIter());
 
   return output->GetConvergence();
 }

SU2_CFD/src/output/CAdjElasticityOutput.cpp

@@ -27,13 +27,16 @@
 
 
 #include "../../include/output/CAdjElasticityOutput.hpp"
+#include "../../include/output/CAdjHeatOutput.hpp"
 #include <string>
 
 #include "../../../Common/include/geometry/CGeometry.hpp"
 #include "../../include/solvers/CSolver.hpp"
 
 CAdjElasticityOutput::CAdjElasticityOutput(CConfig *config, unsigned short nDim) : COutput(config, nDim, false) {
 
+  coupled_heat = config->GetWeakly_Coupled_Heat();
+
   /*--- Initialize number of variables ---*/
   nVar_FEM = nDim;
 
@@ -51,8 +54,10 @@ CAdjElasticityOutput::CAdjElasticityOutput(CConfig *config, unsigned short nDim)
     requestedScreenFields.emplace_back("INNER_ITER");
     requestedScreenFields.emplace_back("ADJOINT_DISP_X");
     requestedScreenFields.emplace_back("ADJOINT_DISP_Y");
+    if (coupled_heat) requestedScreenFields.emplace_back("RMS_ADJ_TEMPERATURE");
     requestedScreenFields.emplace_back("SENS_E_0");
     requestedScreenFields.emplace_back("SENS_NU_0");
+    if (coupled_heat) requestedScreenFields.emplace_back("SENS_GEO");
     nRequestedScreenFields = requestedScreenFields.size();
   }
 
@@ -134,6 +139,12 @@ void CAdjElasticityOutput::SetHistoryOutputFields(CConfig *config){
       AddHistoryOutput("BGS_ADJ_DISP_Z", "bgs[A_Uz]", ScreenOutputFormat::FIXED, "BGS_RES", "BGS residual of the adjoint Z displacement.", HistoryFieldType::RESIDUAL);
     }
   }
+
+  if (coupled_heat) {
+    CAdjHeatOutput::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.");
+  }
 }
 
 inline void CAdjElasticityOutput::LoadHistoryData(CConfig *config, CGeometry *geometry, CSolver **solver) {
@@ -173,6 +184,13 @@ inline void CAdjElasticityOutput::LoadHistoryData(CConfig *config, CGeometry *ge
     }
   }
 
+  /*--- Add heat solver data if available. ---*/
+  if (coupled_heat) {
+    CAdjHeatOutput::LoadHistoryDataImpl(config, geometry, solver, this);
+    SetHistoryOutputValue("LINSOL_ITER_HEAT", solver[ADJHEAT_SOL]->GetIterLinSolver());
+    SetHistoryOutputValue("LINSOL_RESIDUAL_HEAT", log10(solver[ADJHEAT_SOL]->GetResLinSolver()));
+  }
+
   ComputeSimpleCustomOutputs(config);
 }
 
@@ -191,6 +209,11 @@ void CAdjElasticityOutput::LoadVolumeData(CConfig *config, CGeometry *geometry,
   if (nVar_FEM == 3)
     SetVolumeOutputValue("ADJOINT-Z", iPoint, Node_Struc->GetSolution(iPoint, 2));
 
+  if (coupled_heat) {
+    const auto* Node_Heat = solver[ADJHEAT_SOL]->GetNodes();
+    SetVolumeOutputValue("ADJ_TEMPERATURE", iPoint, Node_Heat->GetSolution(iPoint, 0));
+  }
+
   SetVolumeOutputValue("SENSITIVITY-X", iPoint, Node_Struc->GetSensitivity(iPoint, 0));
   SetVolumeOutputValue("SENSITIVITY-Y", iPoint, Node_Struc->GetSensitivity(iPoint, 1));
   if (nDim == 3)
@@ -212,6 +235,14 @@ void CAdjElasticityOutput::LoadVolumeData(CConfig *config, CGeometry *geometry,
   SetVolumeOutputValue("SENS_ACCEL-Y", iPoint, Node_Struc->GetSolution_time_n(iPoint, 2 * nDim + 1));
   if (nDim == 3)
     SetVolumeOutputValue("SENS_ACCEL-Z", iPoint, Node_Struc->GetSolution_time_n(iPoint, 8));
+
+  if (coupled_heat) {
+    const auto* Node_Heat = solver[ADJHEAT_SOL]->GetNodes();
+    SetVolumeOutputValue("SENS_TEMP_N", iPoint, Node_Heat->GetSolution_time_n(iPoint, 0));
+    if (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_2ND) {
+      SetVolumeOutputValue("SENS_TEMP_N1", iPoint, Node_Heat->GetSolution_time_n1(iPoint, 0));
+    }
+  }
 }
 
 void CAdjElasticityOutput::SetVolumeOutputFields(CConfig *config){
@@ -232,6 +263,10 @@ void CAdjElasticityOutput::SetVolumeOutputFields(CConfig *config){
     AddVolumeOutput("ADJOINT-Z", "Adjoint_z", "SOLUTION", "adjoint of displacement in the z direction");
   /// END_GROUP
 
+  if (coupled_heat) {
+    AddVolumeOutput("ADJ_TEMPERATURE", "Adjoint_Temperature", "SOLUTION", "Adjoint Temperature");
+  }
+
   /// BEGIN_GROUP: SENSITIVITY, DESCRIPTION: Geometrical sensitivities of the current objective function.
   /// DESCRIPTION: Sensitivity x-component.
   AddVolumeOutput("SENSITIVITY-X", "Sensitivity_x", "SENSITIVITY", "geometric sensitivity in the x direction");
@@ -261,4 +296,10 @@ void CAdjElasticityOutput::SetVolumeOutputFields(CConfig *config){
   if (nDim == 3)
     AddVolumeOutput("SENS_ACCEL-Z", "SensitivityAccelN_z", "SENSITIVITY_N", "sensitivity to the previous z acceleration");
 
+  if (coupled_heat) {
+    AddVolumeOutput("SENS_TEMP_N", "SensitivityTempN", "SENSITIVITY_N", "sensitivity to the previous temperature");
+    if (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_2ND) {
+      AddVolumeOutput("SENS_TEMP_N1", "SensitivityTempN1", "SENSITIVITY_N", "sensitivity to the previous-1 temperature");
+    }
+  }
 }