Add periodic boundary conditions to the heat-equation solver

SU2_CFD/include/solvers/CHeatSolver.hpp

@@ -247,6 +247,15 @@ class CHeatSolver final : public CSolver {
                                   CConfig *config,
                                   unsigned short val_marker) override;
 
+  /*!
+   * \brief Impose a periodic boundary condition by summing contributions from the complete control volume.
+   * \param[in] geometry - Geometrical definition of the problem.
+   * \param[in] solver_container - Container vector with all the solutions.
+   * \param[in] numerics - Description of the numerical method.
+   * \param[in] config - Definition of the particular problem.
+   */
+  void BC_Periodic(CGeometry* geometry, CSolver** solver_container, CNumerics* numerics, CConfig* config) final;
+
   /*!
    * \brief Set the conjugate heat variables.
    * \param[in] val_marker        - marker index

SU2_CFD/src/solvers/CHeatSolver.cpp

@@ -159,6 +159,19 @@ CHeatSolver::CHeatSolver(CGeometry *geometry, CConfig *config, unsigned short iM
 
   SetBaseClassPointerToNodes();
 
+  /*--- Communicate and store volume and the number of neighbors for any dual CVs that lie on on periodic markers. ---*/
+  for (unsigned short iPeriodic = 1; iPeriodic <= config->GetnMarker_Periodic() / 2; iPeriodic++) {
+    InitiatePeriodicComms(geometry, config, iPeriodic, PERIODIC_VOLUME);
+    CompletePeriodicComms(geometry, config, iPeriodic, PERIODIC_VOLUME);
+    InitiatePeriodicComms(geometry, config, iPeriodic, PERIODIC_NEIGHBORS);
+    CompletePeriodicComms(geometry, config, iPeriodic, PERIODIC_NEIGHBORS);
+  }
+  /*--- Store if implicit scheme is used. This has implications on the Residual and Jacobian handling for periodic
+   * boundaries  ---*/
+  const bool euler_implicit = (config->GetKind_TimeIntScheme_Heat() == EULER_IMPLICIT);
+  SetImplicitPeriodic(euler_implicit);
+
+
   /*--- MPI solution ---*/
 
   InitiateComms(geometry, config, SOLUTION);
@@ -893,6 +906,20 @@ void CHeatSolver::BC_ConjugateHeat_Interface(CGeometry *geometry, CSolver **solv
   }
 }
 
+void CHeatSolver::BC_Periodic(CGeometry* geometry, CSolver** solver_container, CNumerics* numerics,
+                                           CConfig* config) {
+  /*--- Complete residuals for periodic boundary conditions. We loop over
+   the periodic BCs in matching pairs so that, in the event that there are
+   adjacent periodic markers, the repeated points will have their residuals
+   accumulated corectly during the communications. For implicit calculations
+   the Jacobians and linear system are also correctly adjusted here. ---*/
+
+  for (unsigned short iPeriodic = 1; iPeriodic <= config->GetnMarker_Periodic() / 2; iPeriodic++) {
+    InitiatePeriodicComms(geometry, config, iPeriodic, PERIODIC_RESIDUAL);
+    CompletePeriodicComms(geometry, config, iPeriodic, PERIODIC_RESIDUAL);
+  }
+}
+
 void CHeatSolver::Heat_Fluxes(CGeometry *geometry, CSolver **solver_container, CConfig *config) {
 
   unsigned long iPointNormal;
@@ -1274,11 +1301,12 @@ void CHeatSolver::ImplicitEuler_Iteration(CGeometry *geometry, CSolver **solver_
 
     /*--- Modify matrix diagonal to assure diagonal dominance ---*/
 
-    if (nodes->GetDelta_Time(iPoint) != 0.0) {
-
+    const su2double dt = nodes->GetDelta_Time(iPoint);
+    if (dt != 0.0) {
+      /*--- For nodes on periodic boundaries, add the respective partner volume. ---*/
       // Identical for flow and heat
-      const su2double Delta = geometry->nodes->GetVolume(iPoint) / nodes->GetDelta_Time(iPoint);
-      Jacobian.AddVal2Diag(iPoint, Delta);
+      const su2double Vol = geometry->nodes->GetVolume(iPoint) + geometry->nodes->GetPeriodicVolume(iPoint);
+      Jacobian.AddVal2Diag(iPoint, Vol / dt);
 
     } else {
       Jacobian.SetVal2Diag(iPoint, 1.0);
@@ -1326,6 +1354,12 @@ void CHeatSolver::ImplicitEuler_Iteration(CGeometry *geometry, CSolver **solver_
     }
   }
 
+  /*--- Synchronize the solution between master and passive periodic nodes after the linear solve. ---*/
+  for (unsigned short iPeriodic = 1; iPeriodic <= config->GetnMarker_Periodic() / 2; iPeriodic++) {
+    InitiatePeriodicComms(geometry, config, iPeriodic, PERIODIC_IMPLICIT);
+    CompletePeriodicComms(geometry, config, iPeriodic, PERIODIC_IMPLICIT);
+  }
+
   /*--- MPI solution ---*/
 
   InitiateComms(geometry, config, SOLUTION);

TestCases/parallel_regression.py

@@ -1258,6 +1258,22 @@ def main():
     p1rad.tol       = 0.00001
     test_list.append(p1rad)
 
+
+    # #############################
+    # ### Solid Heat Conduction ###
+    # #############################
+
+    # 2D pins, periodically connected
+    solid_periodic_pins           = TestCase('solid_periodic_pins')
+    solid_periodic_pins.cfg_dir   = "solid_heat_conduction/periodic_pins"
+    solid_periodic_pins.cfg_file  = "configSolid.cfg"
+    solid_periodic_pins.test_iter = 750
+    solid_periodic_pins.test_vals = [-15.739745, -14.448665, 300.900000, 425.320000, 0.000000, 5.000000, -1.448445] #last 7 lines
+    solid_periodic_pins.su2_exec  = "mpirun -n 2 SU2_CFD"
+    solid_periodic_pins.timeout   = 1600
+    solid_periodic_pins.tol       = 0.00001
+    test_list.append(solid_periodic_pins)
+
     # ###############################
     # ### Conjugate heat transfer ###
     # ###############################

TestCases/solid_heat_conduction/periodic_pins/README.md

@@ -0,0 +1,10 @@
+# Periodic Pins
+
+This a solid heat conduction testcase in which the `HEAT_EQUATION` solver runs standalone (i.e. not as CHT).
+The simulation domain is the solid domain of the `incomp_navierstokes/streamwise_periodic/chtPinArray_2d`-Testcase.
+Therefore the provided gmsh `.geo` file contains the full CHT mesh but only writes out the solid zone when called.
+
+Note that using periodic boundary conditions for the solid zone made the solution take ~10x more iterations to converge , compared to the same setup using adiabatic walls.
+This was found for solid only as well as CHT cases.
+
+Expected results are perfectly matched Temperatures at the periodic interface. Compare e.g. using Paraview's `Transform`-Filter with domain length 0.0111544m.

TestCases/solid_heat_conduction/periodic_pins/configSolid.cfg

@@ -0,0 +1,82 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%                                                                              %
+% SU2 configuration file                                                       %
+% Case description: Solid-only Heated Pins, periodically connected             %
+% Author: T. Kattmann                                                          %
+% Institution: Bosch Thermotechniek B.V.                                       %
+% Date: 2021.09.27                                                             %
+% File Version 7.2.0 "Blackbird"                                               %
+%                                                                              %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% ------------- DIRECT, ADJOINT, AND LINEARIZED PROBLEM DEFINITION ------------%
+%
+SOLVER= HEAT_EQUATION
+%
+OBJECTIVE_FUNCTION= AVG_TEMPERATURE
+OBJECTIVE_WEIGHT= 1.0
+%
+OPT_OBJECTIVE= AVG_TOTALTEMP
+%
+% ---------------- (SOLIDS) CONDUCTION CONDITION DEFINITION -------------------%
+%
+INC_NONDIM= DIMENSIONAL
+FREESTREAM_TEMPERATURE= 345.0
+MATERIAL_DENSITY= 2719
+SPECIFIC_HEAT_CP = 871.0
+THERMAL_CONDUCTIVITY_CONSTANT= 200
+%
+% -------------------- BOUNDARY CONDITION DEFINITION --------------------------%
+%
+MARKER_HEATFLUX= (solid_pin1_inner, 5e5, \
+                  solid_pin2_inner, 5e5, \
+                  solid_pin3_inner, 0.0, \
+                  solid_pin1_walls, 0.0, \
+                  solid_pin2_walls, 0.0, \
+                  solid_pin3_walls, 0.0)
+%
+MARKER_PERIODIC= ( solid_pin1_periodic, solid_pin3_periodic, 0.0,0.0,0.0, 0.0,0.0,0.0, 0.0111544,0.0,0.0 )
+%
+MARKER_ISOTHERMAL= ( solid_pin1_interface, 300, solid_pin2_interface, 300, solid_pin3_interface, 300 )
+%
+% ------------------------ SURFACES IDENTIFICATION ----------------------------%
+%
+MARKER_MONITORING = ( solid_pin1_inner, solid_pin3_inner )
+%
+% ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------%
+%
+%NUM_METHOD_GRAD= GREEN_GAUSS
+NUM_METHOD_GRAD= WEIGHTED_LEAST_SQUARES
+CFL_NUMBER= 1e8
+%
+% ------------------------ LINEAR SOLVER DEFINITION ---------------------------%
+%
+LINEAR_SOLVER= FGMRES
+LINEAR_SOLVER_PREC= ILU
+LINEAR_SOLVER_ERROR= 1E-15
+LINEAR_SOLVER_ITER= 5
+%
+% --------------------------- CONVERGENCE PARAMETERS --------------------------%
+%
+ITER= 1000
+%
+CONV_FIELD= RMS_TEMPERATURE
+CONV_RESIDUAL_MINVAL= -16
+CONV_STARTITER= 10
+%
+% -------------------- HEAT NUMERICAL METHOD DEFINITION -----------------------%
+%
+TIME_DISCRE_HEAT= EULER_IMPLICIT
+%
+% ------------------------- INPUT/OUTPUT INFORMATION --------------------------%
+%
+MESH_FILENAME= solid.su2
+%
+SCREEN_OUTPUT= INNER_ITER, RMS_TEMPERATURE, MAX_TEMPERATURE, AVG_TEMPERATURE, TOTAL_HEATFLUX, MAXIMUM_HEATFLUX, LINSOL_ITER, LINSOL_RESIDUAL
+SCREEN_WRT_FREQ_INNER= 50
+%
+HISTORY_OUTPUT= (ITER, RMS_RES, HEAT, LINSOL)
+%
+OUTPUT_FILES= RESTART, PARAVIEW_MULTIBLOCK
+VOLUME_OUTPUT= RESIDUAL
+OUTPUT_WRT_FREQ= 1000