A time adaptive multirate Dirichlet-Neumann waveform relaxation method for heterogeneous coupled heat equations

We consider partitioned time integration for heterogeneous coupled heat equations. First and second order multirate, as well as time-adaptive Dirichlet-Neumann Waveform relaxation (DNWR) methods are derived. In 1D and for implicit Euler time integration, we analytically determine optimal relaxation parameters for the fully discrete scheme. We test the robustness of the relaxation parameters on the second order multirate method in 2D. DNWR is shown to be very robust and consistently yielding fast convergence rates, whereas the closely related Neumann-Neumann Waveform relaxtion (NNWR) method is slower or even diverges. The waveform approach naturally allows for different timesteps in the subproblems. In a performance comparison for DNWR, the time-adaptive method dominates the multirate method due to automatically finding suitable stepsize ratios. Overall, we obtain a fast, robust, multirate and time adaptive partitioned solver for unsteady conjugate heat transfer.


Introduction
We consider efficient numerical methods for the partitioned time integration of coupled multiphysics problems.In a partitioned approach different codes for the sub-problems are reused and the coupling is done by a coupling code which calls interface functions of the segregated codes [9,10].These algorithms are currently an active research topic driven by certain multiphysics applications where multiple physical models or multiple simultaneous physical phenomena involve solving coupled systems of partial differential equations (PDEs).
An example of this is fluid structure interaction (FSI) [39,7].Our prime motivation is thermal interaction between fluids and structures, also called conjugate heat transfer.There are two domains with jumps in the material coefficients across the connecting interface.Conjugate heat transfer plays an important role in many applications and its simulation has proved essential [2].Examples for thermal fluid structure interaction are cooling of gas-turbine blades, thermal antiicing systems of airplanes [8], supersonic reentry of vehicles from space [27,19], gas quenching, which is an industrial heat treatment of metal workpieces [18,37] or the cooling of rocket nozzles [21,22].
The most common form of coupling is a Dirichlet-Neumann (DN) approach, in which one problem has a Dirichlet boundary condition on the shared interface, while the other one uses a Neumann boundary condition.In the iteration, they provide each other with the suitable boundary data, i.e. a flux or the interface value.Thus, there is a connection to Domain Decomposition methods.
From the partitioned time integration, we require that it allows for variable and adaptive time steps, preserves the order of the time integration in the subsolvers, and it should be robust and fast.A technique that promises to deliver these properties is the so called Waveform relaxation (WR).An iteration requires solving the subproblems on a time window.Thereby, continuous interface functions, obtained via suitable interpolation, are provided from the respective other problem.WR methods were originally introduced by [24] for ordinary differential equation (ODE) systems, and used for the first time to solve time dependent PDEs in [14,15].They allow the use of different spatial and time discretizations for each subdomain.This is especially useful in problems with strong jumps in the material coefficients [12] or the coupling of different models for the subdomains [11].
A key problem is to make the Waveform iteration fast.A black box approach is to make use of Quasi-Newton methods, leading to Quasi-Newton Waveform iterations [35].Here, we follow instead the idea to tailor very fast methods to a specific problem.In particular, we consider the Neumann-Neumann Waveform relaxation (NNWR) and Dirichlet-Neumann Waveform relaxation (DNWR) method of Gander et al. [23,13], which are WR methods based on the classical Neumann-Neumann and Dirichlet-Neumann iterations.The DNWR method is serial, whereas with NNWR, one can solve the subproblems in parallel.Using an optimal relaxation parameter, convergence in two iterations is obtained for the continuous iteration in 1D.In [31], a fully discrete multirate NNWR method for heterogeneous coupled heat equations is presented.Optimal relaxation parameters are determined for the 1D case.The method was extended to the time adaptive case in [33].
However, the NNWR method is extremely sensitive to the choice of the relaxation parameter, leading to a lack of robustness.In this paper, we therefore focus on the DNWR method, see also [32].The standard DN method was known to be a very fast method for thermal interaction between air and steel [6,5,3].This was thoroughly analyzed for the fully discrete case for two coupled heat equations with different material properties in [30] and for coupled Laplace equations in [16].Thus, we can expect the waveform variant to be a fast solver not only when using an optimal relaxation parameter.The technique employed here to determine optimal relaxation parameters follows [31], which considers the fully discrete iteration for a 1D model problem.Then, using the Toeplitz structure of the arising matrices, a formula for the spectral radius of the iteration matrix can be found and the optimal relaxation parameter can be analytically determined.
We present first and second order multirate WR methods, as well as a second order time adaptive method.The time integration methods we use as a base are implicit Euler and a second order singly diagonally implicit Runge-Kutta (SDIRK2) method.The optimal relaxation parameter Θ opt from the analysis of 1D and implicit Euler yields good results for 2D and SDIRK2.We show how to adapt Θ opt for use in the multirate and time-adaptive setting to get good convergence rates.Additionally, we experimentally show that the convergence results also extend to non-square geometries.Overall, our DNWR method yields a fast and robust solver for unsteady conjugate heat transfer.

Model problem
The unsteady transmission problem reads as follows, where we consider a domain Ω ⊂ R d which is cut into two subdomains Ω = Ω 1 ∪Ω 2 with transmission conditions at the interface Γ = ∂Ω 1 ∩∂Ω 2 : where n m is the outward normal to Ω m for m = 1, 2.
The constants λ 1 and λ 2 describe the thermal conductivities of the materials on Ω 1 and Ω 2 respectively.D 1 and D 2 represent the thermal diffusivities of the materials and are defined by where ρ m is the density and c pm the specific heat capacity of the material placed in Ω m , m = 1, 2.

The Dirichlet-Neumann Waveform Relaxation algorithm
The Dirichlet-Neumann Waveform relaxation (DNWR) method is inspired by substructuring methods from Domain Decomposition.The PDEs are solved sequentially using a Dirichletrespectively Neumann boundary condition with data given from the solution of the other problem, c.f. [25,26].
Given an interface solution g (k) (x, t), (x, t) ∈ Γ × [0, T f ], it consists of the following threestep iteration.Imposing continuity of the solution across the interface, one first finds the local solution u (k+1) 1 (x, t) on (x, t) ∈ Ω 1 × [0, T f ] by solving the Dirichlet problem: The typical initial guess is g (0) (x, t) = u(x, 0) Γ , i.e., extrapolation of the interface initial value.Then, imposing continuity of the heat fluxes across the interface, one finds the local solution u (k+1) 2 (x, t) on Ω 2 by solving the Neumann problem: Finally, the interface values are updated with where Θ ∈ (0, 1] is the relaxation parameter.Note that choosing an appropriate relaxation parameter is crucial to get a fast convergence rate.In [13], the optimal relaxation parameter in 1D has been proven to be Θ = 1/2 for λ 1 = λ 2 = α 1 = α 2 = 1 and subdomains of equal size.
If one uses the optimal relaxation parameter, two iterations are enough for subdomains of equal size.

Semidiscrete method
We now describe a rather general space discretization of (2)-( 4).The core property we need is that the meshes of Ω 1 and Ω 2 share the same nodes on Γ as shown in Figure 1.Furthermore, we assume that there is a specific set of unknowns associated with the interface nodes.Otherwise, we allow at this point for arbitrary meshes on both sides.Then, letting u (m) I : [0, T f ] → R Rm where R m is the number of grid points on Ω m , m = 1, 2, and u Γ : [0, T f ] → R s , where s is the number of grid points at the interface Γ, we can write a general discretization of the first two equations in (2) and (3), respectively, in a compact form as: (6) with initial conditions u (m) To close the system, we need an approximation of the normal derivatives on Γ. Letting φ j be a nodal FE basis function on Ω m for a node on Γ we observe that the normal derivative of u m with respect to the interface can be written as a linear functional using Green's formula [38, pp. 3].Thus, the approximation of the normal derivative is given by Consequently, the equation ΓI u (2),( k+1) ΓI u (1),( k+1) is a semi-discrete version of the third equation in (3) and completes the system ( 5)-( 6).
Omitting the iteration indices, the system of IVPs defined by ( 5), ( 6) and ( 7) is a semidiscretization of (1).We refer to it as the (semidiscrete) monolithic system and its solution as the monolithic solution.
Then, the semidiscrete DNWR algorithm is as follows: In each iteration k, one first solves the Dirichlet problem (5), obtaining u , one solves the following Neumann problem that corresponds to equations ( 6) and ( 7): where with the heat flux 1),(k+1) ΓI u (1),( k+1) Finally, the interface solutions are updated by The iteration starts with a function initial guess u Γ , we use u Γ ≡ u Γ (0).Since the iteration is done on functions, one would like to terminate when u is a user defined tolerance.Since we expect the time integration error to grow with t, we only compare the update at T f .Our termination criterion is i.e., the relative update w.r.t. the initial value at the interface.We use the discrete L 2 interface norm, given by for a mesh uniform at the interface.Here, d is the spatial dimension of (1).

Multirate time discretization
In the case of non-matching material parameters α m , λ m , the Dirichlet and the Neumann problems ( 5) and ( 8) have different needs for time-discretization.Consequently, we want the possibility to use different time integration methods and different step-sizes on each subdomain.
To this end, we first define interpolation functions in time at the interface.We then present two fully discrete DNWR methods, which use the implicit Euler, resp. the second order singly diagonally implicit Runge-Kutta method (SDIRK2) for time integration.We formulate these for general step-sizes ∆t n .One obtains the multirate resp.adaptive algorithm by choosing different step-sizes on the different subdomains.

Interpolation
We denote the fully discrete interface solutions resp.heat fluxes by where m is the number of timesteps on Ω m in the k-th iteration.We define the interpolants as follows: omitting the dependencies of the interpolants on the time-grids for readability.Now, when deriving the fully discrete DNWR methods, we replace all evaluations of the interface temperatures in the Dirichlet problem (5) by evaluations of the interpolant I(u . Analogously, we replace heat flux evaluations for the Neumann problem ( 8) by evaluations of I(q (k+1) ).
The interpolation is done by a coupling code, such that a solver for a subdomain need not know about the other time grid.Here, we consider linear polynomial interpolation.
For readability, we omit dependencies of the time-points on (k) in the following derivations.In case of time-grids varying with k, all time-point evaluations correspond to the time-grid of the current iteration.The time-grid of the previous iteration is only used to define the underlying interpolation data.

Implicit Euler
The implicit Euler method is defined by approximating the derivative term via a standard backward differences.We will use this approximation for derivative terms, i.e., 1),(k+1) In each iteration k of the WR algorithm, given the initial value u (1),(k+1),0 I (0) and the function I(u (k) Γ )(t), we apply the implicit Euler scheme to the Dirichlet problem (5), which yields We compute the discrete heat fluxes as output of the Dirichlet solver by approximating the derivatives in (9) using standard backward differences, this yields The initial flux q (k+1),0 , which is required in the interpolation, is analogously computed using standard forward differences to approximate derivative terms.
Next, we rewrite the Neumann problem (8) in terms of the vector of unknowns .
Algorithm 1 shows a pseudocode of this method.
Algorithm 1 Pseudocode of the DNWR IE method.On domain Ω m we do

SDIRK2
We now introduce a higher order version of the same multirate algorithm.Specifically, we consider the second order singly diagonally implicit Runge-Kutta (SDIRK2) method as a basis to discretize the systems ( 5), ( 8) and (10) in time.For a general IVP u(t) = f (t, u(t)), u(0) = u 0 , t ∈ [0, T f ], SDIRK2 is defined as follows: Here, the so-called stage derivatives are In the following we use j , j = 1, 2, to denote the stage solutions resp.derivatives on Ω m , m = 1, 2.
Using the SDIRK2 scheme to solve the Dirichlet problem (5), with initial value u (1),(k+1),0 IΓ I(u We then compute the stage heat flux ΓI k (1) ΓI U (1) ΓI I(u In both (15) and ( 16), we approximate the derivative terms by For the second SDIRK2 stage, one solves II u (1),(k+1),n+1 IΓ I(u and computes the second heat flux ΓI k (1) ΓI u (1),(k+1),n+1 I + M (1) ΓΓ ΓI I(u In the computation of the second stage, we use the derivative approximation This procedure is repeated for all n = 0, . . ., N − 1.The time-point evaluations in our derivative approximations coincide with the time-points for which the SDIRK2 scheme provides discrete solutions, in case of matching time-grids on both subdomains. We now construct the two interpolants I(q (k+1) 1 ) and I(q ), where the former corresponds to the time-points 0, a∆t 0 , t 1 + a∆t 1 , . ... We compute the initial flux q (k+1),0 , which we use in both flux interpolants, using the second order forward differences: Using the SDIRK2 scheme to solve the Neumann problem (8), replacing the heat flux by I(q (k+1) 1 ) resp.I(q ), consists of first solving followed by − 1.The relaxation step and termination check are identical to the implicit Euler method.Here, we use a total of only 3 interpolants: The solution, the heat flux and the stage heat flux.This was achieved by the approximation of U ).As can be seen in the numerical results in Section 8.3, this does not lead to a loss of order.Other options for constructing partitioned time-integration methods based on SDIRK2 are discussed in [28,Chap. 5].
The differences to Algorithm 1 are in SolveDirichlet and SolveNeumann, which require solving two linear systems each.Additionally, SolveDirichlet now returns two heat fluxes, which are both interpolated and passed into SolveNeumann.
Including the relaxation step yields the following iteration (for a single timestep): In the 1D case S (m) and Σ are scalars, thus the optimal relaxation parameter is .
In the following we specifically consider a 1D model problem on Ω = [−1, 1], split at x Γ = 0, and an equidistant discretization using linear finite elements.The matrices M II .Through lengthy but straight forward calculations (see [30,29]), one obtains the following expressions: Using c = ∆t/∆x 2 , Θ opt has the following temporal and spatial limits [30]: These are consistent with the one-dimensional continuous analysis performed in [13,25].There, a convergence analysis using Laplace transforms for the DNWR method ( 2)-( 4) on two identical subdomains Ω 1 and Ω 2 with constant coefficients shows that Θ opt = 1/2.Their result is recovered when approaching the continuous case in the limit ∆t/∆x 2 → ∞ for constant coefficients.Figure 2 shows Θ opt for a few material combinations, see Table 1.One can observe that Θ opt is continuous and bounded by its spatial and temporal limits (20).

Multirate relaxation parameter
For ∆t 1 = ∆t 2 the analysis in the previous section does not apply anymore.Instead, we determine Θ opt based on numerical experiments in Section 8.1.These show that, on average, the optimal choice is to use Θ opt based on the maximum of ∆t 1 and ∆t 2 .This result coincides with the experiments to determine Θ opt for the multirate NNWR method in [33].

Time adaptive method
The goal of adaptivity is to use step-sizes as large as possible and as small as necessary to reduce computational costs, while ensuring a target accuracy.In particular, using adaptive time stepping for both sub-domains separately, one can attain comparable time-integration errors automatically, bypassing the need to determine a suitable step-size ratio for the multirate case.
The basic idea is to control timestepsizes to keep an error estimate at a given tolerance T OL.We use a local error estimate obtained by an embedded technique [17,chap. IV.8].With the SDIRK2 method, we obtain a lower order embedded solution ûn+1 via The local error estimate is n = u n − ûn .We then control timesteps using the proportionalintegral (PI) controller , cf. [1], PI3333.We use this procedure on each subdomain independently.As the initial stepsize we use ∆t , m = 1, 2, c.f. [33].We choose the tolerances T OL (m) = T OL W R /5, m = 1, 2. This choice is motivated by [36] and already used in a similar context in [4,33].We use the discrete L 2 norm where M is the corresponding mass matrix and |Ω u | the area on which u is defined.Using this adaptive method, we get independent time-grids for both sub-domains, that are suitable for the given material parameters.A pseudocode of the adaptive SDIRK2 DNWR method is shown in Algorithm 2.

Relaxation parameter
In the adaptive method the Dirichlet-Neumann operator changes every WR step, since the timestepsizes change.Consequently, we recompute Θ in every iteration.We use Θ as in the multirate case with the average stepsizes from each subdomain.This approach improves upon [33], where Θ was based on the timegrids of the previous WR iteration.
One can solve on the subdomains in parallel.The NNWR algorithm is based on the exact same Dirichlet and Neumann problems as the DNWR algorithm, but with different input data, namely fluxes and initial value for the Neumann problem (22).Hence one can directly use the time-discretizations as described in Section 5, including time-adaptivity, c.f. [33].
Under the same restrictions as in Section 5.4, i.e., Ω = [−1, 1], split at x Γ = 0 and linear finite elements on an equidistant grid, one can analogously calculate Θ opt : with S (m) given by (19), see [31] for details.The spatial and temporal limits based on c = ∆t/∆x 2 are lim c→0 8 Numerical results We now present numerical experiments to illustrate the validity of the theoretical results and to test the robustness of the relaxation parameters in 2D, SDIRK2 and with multirate resp.adaptive time-grids.The methods and algorithms described have been implemented in Python 3.6, the code is available at [34].
As space discretization we use linear finite elements on equidistant grids (1D: ∆x = 1/200, 2D: ∆x = 1/100) as shown in Figure 1, see [4] for more details.The resulting linear equation systems are solved with direct solvers.
We experimentally determine the convergence rate via i.e., the reduction rate in the update.Here, we perform up to k max = 6 iterations and take the mean of the update reductions, but never the last iteration, which could be near machine   precision.This experiment is done using IE for the 1D test case and N = 1, as we aim to determine the asymptotic convergence rates.The results in Figures 3, 4 and 5 show that all options yield comparable results, with "Max" being most consistent, making it our choice for DNWR in the multirate setting.
Numerical experiments in [31] yielded the same conclusion for the NNWR method.As such we use (23) based on the larger step-size.

Optimality of relaxation parameter
We now verify the optimality of the relaxation parameters (18) for DNWR and (23) for NNWR in 1D with implicit Euler and also test the convergence rate and robustness.To this end, we determine the experimental convergences rates as in Section 8.1, for varying Θ in 1D and 2D for both implicit Euler and SDIRK2.Time-integration is done in multirate and non-multirate, up to T f using N = 100 base timesteps.
We expect little variation for implicit Euler and SDIRK2, since SDIRK2 consists of two successive implicit Euler steps.We anticipate more notable differences between 1D and 2D.Lastly, convergence rates might deviate due to transitive effects of WR, since the iteration matrices are non-normal.
In the plots, the blue highlighted range on the x-axis marks the spatial and temporal limits of Θ, see (20) resp.(24).Results are seen in Figure 6.In both 1D and 2D, SDIRK2 rates match those of implicit Euler.In all cases, 1D results closely align with the theoretical result marked by Σ(Θ).2D results are slightly off, but still yield good error reduction rates using the 1D Θ opt .For air-water, the error reduction rate is ≈ 10 −2 , i.e., the coupling residual gains two decimals in accuracy per iteration.The air-steel coupling yields very fast convergence with an error reduction rate of ≈ 10 −4 and water-steel rates are between 0.1 and 0.01 for Θ opt .

DNWR
Additionally, we see that DNWR is convergent for all shown Θ.Thus, in the worst case DNWR convergence is slow, yet not divergent.
We test only the non-multirate setting.Results are shown in Figure 7 and strongly resemble the results for square, identical domains in Figure 6, except for slower convergence rates for the 2D water-steel case.
This similarity can be explained by looking at the Schur-complements (17).M II does affect its inverse, but after multiplication with M (1)

IΓ and M
(1) ΓI /∆t + A (1) ΓI , which are mostly zero, the effect on S (1) is expected to be minor.Results are shown in Figure 8.We additionally mark the divergence limit at 1, showing the range of viable Θ is very small for NNWR and unlike DNWR, relaxation is non-optional for convergence.In particular, one may get divergence for Θ within the range marked by the temporal and spatial limits.Convergence rates for implicit Euler and SDIRK2 results are almost identical.1D convergence rates align well with the theoretical results in all cases.In 2D, the airwater and air-steel results match with the 1D results, yielding rates of about 0.01−0.1.However, water-steel shows divergence in 2D, when using Θ opt .In the convergent cases, the observed error reduction rates are slower than for DNWR, this is particularly pronounced in the air-steel case, with a difference of about 3 orders of magnitude.Overall, NNWR shows a lack of robustness.
One might achieve better convergence rates using macrostepping, i.e., successively performing the algorithm on smaller time-windows.This may speed up convergence on each time-window, but the coupling residual propagates through erroneous initial values.On the other hand, DNWR performs well on the given time-windows.

Multirate -convergence order of time-integration
We show convergence of the error, on the whole domain in the discrete L 2 norm (21) and using T f = 1, for ∆t → 0. Our reference solution is the monolithic solution for sufficiently small step-sizes, thus measuring both the time-integration error and coupling residual.
Results for T OL W R = 10 −13 can be seen in Figure 9 for DNWR and in Figure 10 for NNWR.We attain the expected first and second order convergence rates for ∆t → 0.

Time adaptive results
We consider the time-adaptive DNWR method described in Section 6.The reference for error computation is the solution using T OL W R = 10 −8 in 1D and T OL W R = 10 −7 in 2D.We expect  the errors to be proportional to the tolerance for T OL W R → 0, which is observed in Figure 11.Due to its lacking robustness, we do not consider time-adaptive NNWR.

Error over work comparison
We now compare efficiency of the adaptive and multirate method for the 2D test case with ∆x = 1/200.For this we compare error over work, which we measure as the total number of timesteps.
We choose the stepsize ratios in the multirate setting such that both domains use comparable CFL numbers, which is achieved by c  2 for the resulting stepsize ratios for our material configurations.
To compare the multirate method with the time-adaptive method, we parametrize the former by the number of base timesteps N .Given ∆t m = T f /(c m • N ), m = 1, 2, we compute the associated time integration error e ∆t 1 ,∆t 2 using T OL W R = 10 −12 and a monolithic reference solution with ∆t = min(∆t 1 , ∆t 2 )/2.We then use T OL W R = e ∆t 1 ,∆t 2 /5 in the termination criterion for the multirate method, for its error over work comparison.Finally, our references for the error computations in the error over work comparison are adaptive solutions with T OL W R = 10 −6 .
Results are shown in Figure 12 with the resulting stepsizes ratios for the adaptive case in Table 2.The adaptive method is 4 times more efficient in the water-steel case, of similar efficiency in the air-water case and less efficient in the air-steel case.This can be explained by the stepsize ratios in Table 2.The closer the multirate stepsize ratios correspond to the adaptive ones, the better the performance of multirate in comparison with the adaptive method.As second test case we consider the initial condition u(x, y) = 800 sin((x + 1)π) 2 sin(yπ).( 27) Here, we have u Γ (0) = 0 and thus skip the relative norm for the termination check (11).Results are shown in Figure 13 with stepsize ratios in Table 2.In the air-steel case performance is approximately equal, whereas adaptive performance is about 4 resp.25 times better in the air-water resp.water-steel case.
Overall we see that performance depends on the stepsize ratios.This makes the adaptive method a more robust choice, since it automatically determines suitable stepsize ratios, which vary for e.g., different initial conditions.

Summary and conclusions
We derived first and second order, multirate resp.time-adaptive DNWR methods for heterogeneous coupled heat equations.The optimal relaxation parameter Θ opt for WR is shown to be identical to the one for the basic DN iteration.We experimentally show how to adapt Θ in the multirate case.The observed convergence rates using an analytical Θ opt for 1D implicit Euler  (27).are shown to be very robust, yielding fast convergence rates for a second order method and 2D, for various material combinations and multirate settings on long time intervals.
The same tests for the related NNWR methods employing identical Dirichlet and Neumann subsolvers, using an analytical Θ opt for 1D implicit Euler, show a lack of robustness, possibly resulting in divergence.
The time-adaptive DNWR method is experimentally shown to be favorable over mutirate, due ease of use and superior performance.The latter is due to the resulting stepsizes being more suitably chosen than those of the multirate solver.Overall, we obtain a fast, robust, time adaptive (on each domain), partitioned solver for unsteady conjugate heat transfer.

( 1 )
,(k+1) I .Then, using the function of unknowns u(k+1) Toeplitz structure.Thus, one can write down an exact expression of (17) using an Eigendecomposition to calculate the inverse of M (m) II /∆t + A (m)
II increase in size, but values and structure persist.The matrices M are padded with additional zeros.The increased size of M(1)

Figure 11 :
Figure 11: Left: 1D.Right: 2D.Error over T OL W R for the time-adaptive DNWR method.

Table 2 :
Air-water Air-steel Water-steel multirate (c 1 : c 2 ) Timestep ratios for the multirate and adaptive method (final grid) by materials.u 1 0 is for the initial condition (25) and u 2 0 is for(27).