Heat solver fixes for primal and adjoint CHT simulations

Common/src/linear_algebra/CSysSolve.cpp

@@ -246,7 +246,7 @@ unsigned long CSysSolve<ScalarType>::CG_LinSolver(const CSysVector<ScalarType> &
     norm_r = r.norm();
     norm0  = b.norm();
     if ((norm_r < tol*norm0) || (norm_r < eps)) {
-      if (master) cout << "CSysSolve::ConjugateGradient(): system solved by initial guess." << endl;
+      if (master && !mesh_deform) cout << "CSysSolve::ConjugateGradient(): system solved by initial guess." << endl;
       return 0;
     }
 

SU2_CFD/include/solvers/CDiscAdjFEASolver.hpp

@@ -176,10 +176,9 @@ class CDiscAdjFEASolver final : public CSolver {
   /*!
    * \brief Extract and set the geometrical sensitivity.
    * \param[in] geometry - Geometrical definition of the problem.
-   * \param[in] solver - The solver container holding all terms of the solution.
    * \param[in] config - Definition of the particular problem.
    */
-  void SetSensitivity(CGeometry *geometry, CSolver **solver, CConfig *config) override;
+  void SetSensitivity(CGeometry *geometry, CConfig *config, CSolver*) override;
 
   /*!
    * \brief Set the objective function.

SU2_CFD/include/solvers/CDiscAdjMeshSolver.hpp

@@ -97,10 +97,10 @@ class CDiscAdjMeshSolver final : public CSolver {
   /*!
    * \brief Extract and set the geometrical sensitivity.
    * \param[in] geometry - Geometrical definition of the problem.
-   * \param[in] solver - The solver container holding all terms of the solution.
    * \param[in] config - Definition of the particular problem.
+   * \param[in] target_solver - The target solver to store the sensitivities.
    */
-  void SetSensitivity(CGeometry *geometry, CSolver **solver, CConfig *config) override;
+  void SetSensitivity(CGeometry *geometry, CConfig *config, CSolver* target_solver) override;
 
   /*!
    * \brief Prepare the solver for a new recording.

SU2_CFD/include/solvers/CDiscAdjSolver.hpp

@@ -156,10 +156,9 @@ class CDiscAdjSolver final : public CSolver {
   /*!
    * \brief Extract and set the geometrical sensitivity.
    * \param[in] geometry - Geometrical definition of the problem.
-   * \param[in] solver - The solver container holding all terms of the solution.
    * \param[in] config - Definition of the particular problem.
    */
-  void SetSensitivity(CGeometry *geometry, CSolver **solver, CConfig *config) override;
+  void SetSensitivity(CGeometry *geometry, CConfig *config, CSolver*) override;
 
   /*!
    * \brief Set the objective function.

SU2_CFD/include/solvers/CSolver.hpp

@@ -3805,10 +3805,10 @@ class CSolver {
   /*!
    * \brief A virtual member. Extract and set the geometrical sensitivity.
    * \param[in] geometry - Geometrical definition of the problem.
-   * \param[in] solver - The solver container holding all terms of the solution.
    * \param[in] config - Definition of the particular problem.
+   * \param[in] solver - The target solver for the sensitivities, optional, for when the mesh solver is used.
    */
-  inline virtual void SetSensitivity(CGeometry *geometry, CSolver **solver, CConfig *config){ }
+  inline virtual void SetSensitivity(CGeometry *geometry, CConfig *config, CSolver *target_solver = nullptr){ }
 
   /*!
    * \brief A virtual member. Extract and set the derivative of objective function.

SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp

@@ -404,20 +404,24 @@ void CDiscAdjMultizoneDriver::EvaluateSensitivities(unsigned long iOuterIter, bo
       case DISC_ADJ_INC_EULER: case DISC_ADJ_INC_NAVIER_STOKES: case DISC_ADJ_INC_RANS:
 
         if(Has_Deformation(iZone)) {
-          solvers[ADJMESH_SOL]->SetSensitivity(geometry, solvers, config);
+          solvers[ADJMESH_SOL]->SetSensitivity(geometry, config, solvers[ADJFLOW_SOL]);
         } else {
-          solvers[ADJFLOW_SOL]->SetSensitivity(geometry, solvers, config);
+          solvers[ADJFLOW_SOL]->SetSensitivity(geometry, config);
         }
         break;
 
       case DISC_ADJ_HEAT:
 
-        solvers[ADJHEAT_SOL]->SetSensitivity(geometry, solvers, config);
+        if(Has_Deformation(iZone)) {
+          solvers[ADJMESH_SOL]->SetSensitivity(geometry, config, solvers[ADJHEAT_SOL]);
+        } else {
+          solvers[ADJHEAT_SOL]->SetSensitivity(geometry, config);
+        }
         break;
 
       case DISC_ADJ_FEM:
 
-        solvers[ADJFEA_SOL]->SetSensitivity(geometry, solvers, config);
+        solvers[ADJFEA_SOL]->SetSensitivity(geometry, config);
         break;
 
       default:

SU2_CFD/src/drivers/CDiscAdjSinglezoneDriver.cpp

@@ -555,20 +555,15 @@ void CDiscAdjSinglezoneDriver::SecondaryRecording(){
 
   /*--- Extract the computed sensitivity values. ---*/
 
-  int IDX_SOL;
-
   if (config->GetKind_Solver() == DISC_ADJ_FEM) {
-    IDX_SOL = ADJFEA_SOL;
-  } else if(SecondaryVariables == MESH_COORDS) {
-    IDX_SOL = ADJFLOW_SOL;
-  } else if(SecondaryVariables == MESH_DEFORM) {
-    IDX_SOL = ADJMESH_SOL;
-  } else {
-    IDX_SOL = -1;
+    solver[ADJFEA_SOL]->SetSensitivity(geometry, config);
+  }
+  else if(SecondaryVariables == MESH_COORDS) {
+    solver[ADJFLOW_SOL]->SetSensitivity(geometry, config);
+  }
+  else { // MESH_DEFORM
+    solver[ADJMESH_SOL]->SetSensitivity(geometry, config, solver[ADJFLOW_SOL]);
   }
-
-  if(IDX_SOL >= 0)
-    solver[IDX_SOL]->SetSensitivity(geometry, solver, config);
 
   /*--- Clear the stored adjoint information to be ready for a new evaluation. ---*/
 

SU2_CFD/src/drivers/CDriver.cpp

@@ -1063,7 +1063,7 @@ void CDriver::Solver_Preprocessing(CConfig* config, CGeometry** geometry, CSolve
 
   /*--- Set up any necessary inlet profiles ---*/
 
-  Inlet_Preprocessing(solver, geometry, config); 
+  Inlet_Preprocessing(solver, geometry, config);
 
 }
 
@@ -1196,7 +1196,7 @@ void CDriver::Solver_Restart(CSolver ***solver, CGeometry **geometry,
                              CConfig *config, bool update_geo) {
 
   bool euler, ns, turbulent,
-  NEMO_euler, NEMO_ns, 
+  NEMO_euler, NEMO_ns,
   adj_euler, adj_ns, adj_turb,
   heat, fem, fem_euler, fem_ns, fem_dg_flow,
   template_solver, disc_adj, disc_adj_fem, disc_adj_turb, disc_adj_heat;
@@ -1447,7 +1447,7 @@ void CDriver::Numerics_Preprocessing(CConfig *config, CGeometry **geometry, CSol
       NEMO_euler = compressible = true; break;
 
     case NEMO_NAVIER_STOKES:
-      NEMO_ns = compressible = true; break;  
+      NEMO_ns = compressible = true; break;
 
     case RANS:
     case DISC_ADJ_RANS:
@@ -3409,7 +3409,7 @@ bool CTurbomachineryDriver::Monitor(unsigned long ExtIter) {
   case INC_EULER: case INC_NAVIER_STOKES: case INC_RANS:
     StopCalc = integration_container[ZONE_0][INST_0][FLOW_SOL]->GetConvergence(); break;
   case NEMO_EULER: case NEMO_NAVIER_STOKES:
-      StopCalc = integration_container[ZONE_0][INST_0][FLOW_SOL]->GetConvergence(); break;  
+      StopCalc = integration_container[ZONE_0][INST_0][FLOW_SOL]->GetConvergence(); break;
   case DISC_ADJ_EULER: case DISC_ADJ_NAVIER_STOKES: case DISC_ADJ_RANS:
   case DISC_ADJ_INC_EULER: case DISC_ADJ_INC_NAVIER_STOKES: case DISC_ADJ_INC_RANS:
   case DISC_ADJ_FEM_EULER: case DISC_ADJ_FEM_NS: case DISC_ADJ_FEM_RANS:

SU2_CFD/src/iteration/CDiscAdjHeatIteration.cpp

@@ -195,6 +195,14 @@ void CDiscAdjHeatIteration::RegisterInput(CSolver***** solver, CGeometry**** geo
 
     geometry[iZone][iInst][MESH_0]->RegisterCoordinates(config[iZone]);
   }
+  /*--- Register the variables of the mesh deformation ---*/
+  if (kind_recording == MESH_DEFORM) {
+    /*--- Undeformed mesh coordinates ---*/
+    solver[iZone][iInst][MESH_0][ADJMESH_SOL]->RegisterSolution(geometry[iZone][iInst][MESH_0], config[iZone]);
+
+    /*--- Boundary displacements ---*/
+    solver[iZone][iInst][MESH_0][ADJMESH_SOL]->RegisterVariables(geometry[iZone][iInst][MESH_0], config[iZone]);
+  }
 }
 
 void CDiscAdjHeatIteration::SetDependencies(CSolver***** solver, CGeometry**** geometry, CNumerics****** numerics,

SU2_CFD/src/solvers/CDiscAdjFEASolver.cpp

@@ -805,7 +805,7 @@ void CDiscAdjFEASolver::Preprocessing(CGeometry *geometry, CSolver **solver_cont
 
 }
 
-void CDiscAdjFEASolver::SetSensitivity(CGeometry *geometry, CSolver **solver, CConfig *config){
+void CDiscAdjFEASolver::SetSensitivity(CGeometry *geometry, CConfig *config, CSolver*){
 
   unsigned short iVar;
 

SU2_CFD/src/solvers/CDiscAdjMeshSolver.cpp

@@ -198,12 +198,10 @@ void CDiscAdjMeshSolver::ExtractAdjoint_Variables(CGeometry *geometry, CConfig *
 
 }
 
-void CDiscAdjMeshSolver::SetSensitivity(CGeometry *geometry, CSolver **solver, CConfig *config) {
+void CDiscAdjMeshSolver::SetSensitivity(CGeometry *geometry, CConfig *config, CSolver *solver) {
 
-  unsigned long iPoint;
-  unsigned short iDim;
-  su2double Sensitivity, eps;
-  bool time_stepping = (config->GetTime_Marching() != STEADY);
+  const bool time_stepping = (config->GetTime_Marching() != STEADY);
+  const auto eps = config->GetAdjSharp_LimiterCoeff()*config->GetRefElemLength();
 
   /*--- Extract the sensitivities ---*/
   ExtractAdjoint_Solution(geometry, config);
@@ -212,30 +210,28 @@ void CDiscAdjMeshSolver::SetSensitivity(CGeometry *geometry, CSolver **solver, C
   ExtractAdjoint_Variables(geometry, config);
 
   /*--- Store the sensitivities in the flow adjoint container ---*/
-  for (iPoint = 0; iPoint < nPoint; iPoint++) {
+  for (auto iPoint = 0u; iPoint < nPoint; iPoint++) {
 
-    for (iDim = 0; iDim < nDim; iDim++) {
+    /*--- If sharp edge, set the sensitivity to 0 on that region ---*/
+    su2double limiter = 1.0;
+    if (config->GetSens_Remove_Sharp() && (geometry->nodes->GetSharpEdge_Distance(iPoint) < eps)) {
+      limiter = 0.0;
+    }
 
+    for (auto iDim = 0u; iDim < nDim; iDim++) {
       /*--- The sensitivity was extracted using ExtractAdjoint_Solution ---*/
-      Sensitivity = nodes->GetSolution(iPoint,iDim);
-
-      /*--- If sharp edge, set the sensitivity to 0 on that region ---*/
-      if (config->GetSens_Remove_Sharp()) {
-        eps = config->GetVenkat_LimiterCoeff()*config->GetRefElemLength();
-        if ( geometry->nodes->GetSharpEdge_Distance(iPoint) < config->GetAdjSharp_LimiterCoeff()*eps )
-          Sensitivity = 0.0;
-      }
+      const su2double Sensitivity = limiter * nodes->GetSolution(iPoint,iDim);
 
       /*--- Store the sensitivities ---*/
       if (!time_stepping) {
-        solver[ADJFLOW_SOL]->GetNodes()->SetSensitivity(iPoint, iDim, Sensitivity);
+        solver->GetNodes()->SetSensitivity(iPoint, iDim, Sensitivity);
       } else {
-        solver[ADJFLOW_SOL]->GetNodes()->SetSensitivity(iPoint, iDim,
-          solver[ADJFLOW_SOL]->GetNodes()->GetSensitivity(iPoint, iDim) + Sensitivity);
+        solver->GetNodes()->SetSensitivity(iPoint, iDim,
+          solver->GetNodes()->GetSensitivity(iPoint, iDim) + Sensitivity);
       }
     }
   }
-  solver[ADJFLOW_SOL]->SetSurface_Sensitivity(geometry, config);
+  solver->SetSurface_Sensitivity(geometry, config);
 
 }
 

SU2_CFD/src/solvers/CDiscAdjSolver.cpp

@@ -739,7 +739,7 @@ void CDiscAdjSolver::SetAdjoint_OutputMesh(CGeometry *geometry, CConfig *config)
 
 }
 
-void CDiscAdjSolver::SetSensitivity(CGeometry *geometry, CSolver **solver, CConfig *config) {
+void CDiscAdjSolver::SetSensitivity(CGeometry *geometry, CConfig *config, CSolver*) {
 
   unsigned long iPoint;
   unsigned short iDim;

SU2_CFD/src/solvers/CHeatSolver.cpp

@@ -770,65 +770,20 @@ void CHeatSolver::Set_Heatflux_Areas(CGeometry *geometry, CConfig *config) {
 void CHeatSolver::BC_Isothermal_Wall(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, CNumerics *visc_numerics, CConfig *config,
                                        unsigned short val_marker) {
 
-  unsigned long iPoint, iVertex, Point_Normal;
-  unsigned short iDim;
-  su2double *Normal, *Coord_i, *Coord_j, Area, dist_ij, laminar_viscosity, thermal_diffusivity, Twall, dTdn, Prandtl_Lam;
-  //su2double Prandtl_Turb;
-  bool implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT);
-
-  bool flow = ((config->GetKind_Solver() == INC_NAVIER_STOKES)
-               || (config->GetKind_Solver() == INC_RANS)
-               || (config->GetKind_Solver() == DISC_ADJ_INC_NAVIER_STOKES)
-               || (config->GetKind_Solver() == DISC_ADJ_INC_RANS));
+  const bool implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT);
+  const auto Marker_Tag = config->GetMarker_All_TagBound(val_marker);
+  const su2double Twall = config->GetIsothermal_Temperature(Marker_Tag)/config->GetTemperature_Ref();
 
-  Prandtl_Lam = config->GetPrandtl_Lam();
-//  Prandtl_Turb = config->GetPrandtl_Turb();
-  laminar_viscosity = config->GetMu_ConstantND();
-  //Prandtl_Turb = config->GetPrandtl_Turb();
-  //laminar_viscosity = config->GetViscosity_FreeStreamND(); // TDE check for consistency for CHT
+  for (auto iVertex = 0u; iVertex < geometry->nVertex[val_marker]; iVertex++) {
+    const auto iPoint = geometry->vertex[val_marker][iVertex]->GetNode();
 
-  string Marker_Tag = config->GetMarker_All_TagBound(val_marker);
-
-  Twall = config->GetIsothermal_Temperature(Marker_Tag)/config->GetTemperature_Ref();
+    if (!geometry->nodes->GetDomain(iPoint)) continue;
 
-  for (iVertex = 0; iVertex < geometry->nVertex[val_marker]; iVertex++) {
+    nodes->SetSolution_Old(iPoint,&Twall);
+    LinSysRes(iPoint, 0) = 0.0;
+    nodes->SetRes_TruncErrorZero(iPoint);
 
-    iPoint = geometry->vertex[val_marker][iVertex]->GetNode();
-
-    if (geometry->nodes->GetDomain(iPoint)) {
-
-        Point_Normal = geometry->vertex[val_marker][iVertex]->GetNormal_Neighbor();
-
-        Normal = geometry->vertex[val_marker][iVertex]->GetNormal();
-        Area = 0.0;
-        for (iDim = 0; iDim < nDim; iDim++) Area += Normal[iDim]*Normal[iDim];
-        Area = sqrt (Area);
-
-        Coord_i = geometry->nodes->GetCoord(iPoint);
-        Coord_j = geometry->nodes->GetCoord(Point_Normal);
-        dist_ij = 0;
-        for (iDim = 0; iDim < nDim; iDim++)
-          dist_ij += (Coord_j[iDim]-Coord_i[iDim])*(Coord_j[iDim]-Coord_i[iDim]);
-        dist_ij = sqrt(dist_ij);
-
-        dTdn = -(nodes->GetSolution(Point_Normal,0) - Twall)/dist_ij;
-
-        if(flow) {
-          thermal_diffusivity = laminar_viscosity/Prandtl_Lam;
-        }
-        else
-          thermal_diffusivity = config->GetThermalDiffusivity_Solid();
-
-        Res_Visc[0] = thermal_diffusivity*dTdn*Area;
-
-        if(implicit) {
-
-          Jacobian_i[0][0] = -thermal_diffusivity/dist_ij * Area;
-        }
-
-        LinSysRes.SubtractBlock(iPoint, Res_Visc);
-        Jacobian.SubtractBlock2Diag(iPoint, Jacobian_i);
-    }
+    if (implicit) Jacobian.DeleteValsRowi(iPoint);
   }
 }
 
@@ -1143,6 +1098,12 @@ void CHeatSolver::BC_ConjugateHeat_Interface(CGeometry *geometry, CSolver **solv
 
           HeatFluxDensity   = thermal_diffusivity*(Tinterface - Tnormal_Conjugate);
           HeatFlux          = HeatFluxDensity * Area;
+
+          if (implicit) {
+
+            Jacobian_i[0][0] = -thermal_diffusivity*Area;
+            Jacobian.SubtractBlock2Diag(iPoint, Jacobian_i);
+          }
         }
         else {
 
@@ -1152,12 +1113,6 @@ void CHeatSolver::BC_ConjugateHeat_Interface(CGeometry *geometry, CSolver **solv
 
         Res_Visc[0] = -HeatFlux;
         LinSysRes.SubtractBlock(iPoint, Res_Visc);
-
-        if (implicit) {
-
-          Jacobian_i[0][0] = thermal_diffusivity*Area;
-          Jacobian.SubtractBlock2Diag(iPoint, Jacobian_i);
-        }
       }
     }
   }
@@ -1642,7 +1597,9 @@ void CHeatSolver::ImplicitEuler_Iteration(CGeometry *geometry, CSolver **solver_
 
   /*--- Solve or smooth the linear system ---*/
 
-  System.Solve(Jacobian, LinSysRes, LinSysSol, geometry, config);
+  auto iter = System.Solve(Jacobian, LinSysRes, LinSysSol, geometry, config);
+  SetIterLinSolver(iter);
+  SetResLinSolver(System.GetResidual());
 
   for (iPoint = 0; iPoint < nPointDomain; iPoint++) {
     for (iVar = 0; iVar < nVar; iVar++) {