Efficient Parallel-in-Time Solution of Time-Periodic Problems Using a Multi-Harmonic Coarse Grid Correction

TL;DR

This paper introduces a parallel-in-time algorithm for nonlinear time-periodic problems that leverages multi-harmonic correction and a simplified Newton method to enable efficient, highly-parallelized solutions in both frequency and time domains.

Contribution

It proposes a novel multi-harmonic coarse grid correction method that enhances parallelization in solving nonlinear time-periodic problems, improving efficiency over existing approaches.

Findings

The method converges effectively for a 1D model problem.

It outperforms existing methods in solving eddy current problems.

Demonstrates scalability on a 3D transformer model.

Abstract

This paper presents a highly-parallelizable parallel-in-time algorithm for efficient solution of nonlinear time-periodic problems. It is based on the time-periodic extension of the Parareal method, known to accelerate sequential computations via parallelization on the fine grid. The proposed approach reduces the complexity of the periodic Parareal solution by introducing a simplified Newton algorithm, which allows an additional parallelization on the coarse grid. In particular, at each Newton iteration a multi-harmonic correction is performed, which converts the block-cyclic periodic system in the time domain into a block-diagonal system in the frequency domain, thereby solving for each frequency component in parallel. The convergence analysis of the method is discussed for a one-dimensional model problem. The introduced algorithm and several existing solution approaches are compared via their application to the eddy current problem for both linear and nonlinear models of a coaxial cable. Performance of the considered methods is also illustrated for a three-dimensional transformer model.

Tables

Table 1. Table 1: Percentage of the wall clock time of PP-PC MH and TP MH with respect to the sequential time stepping.

	$N=5$	$N=10$	$N=20$	$N=50$
Sequential	$100\%$
PP-PC MH	$6.6\%$	$2.5\%$	$0.9\%$	$0.4\%$
TP MH	$5.3\%$	$2.5\%$	$1.3\%$	$0.6\%$

Equations

$$\displaystyle\mathbf{M}\mbox{\boldmath$u$}^{\prime}(t)+\mathbf{K}\bigl{(}\mbox{\boldmath$u$}(t)\bigr{)}\mbox{\boldmath$u$}(t)$$

$$\displaystyle\qquad\mbox{\boldmath$u$}(0)$$

$$\displaystyle\mathbf{M}\mbox{\boldmath$u$}_{n}^{\prime}(t)+\mathbf{K}\bigl{(}\mbox{\boldmath$u$}_{n}(t)\bigr{)}\mbox{\boldmath$u$}_{n}(t)$$

$$\displaystyle\mbox{\boldmath$u$}_{n}(T_{n-1})$$

$${\mathbf{U}}_0^{(k+1)}$$

$${\mathbf{U}}_n^{(k+1)}$$

$${\mathbf{U}}_0^{(k+1)}$$

$${\mathbf{U}}_n^{(k+1)}$$

$$\mathbf{b}_{n}^{(k)}=\mathcal{F}\bigl{(}T_{n},T_{n-1},\mathbf{U}^{(k)}_{n-1}\bigr{)}-\mathcal{G}\bigl{(}T_{n},T_{n-1},\mathbf{U}^{(k)}_{n-1}\bigr{)},\quad n=1,\dots,N$$

$$\left[\begin{array}{cccc}\mathbf{I}& 0& \dots & -\mathcal{G}\left(T_N,T_{N-1},\cdot \right)\\ -\mathcal{G}\left(T_1,T_0,\cdot \right)& \mathbf{I}& & 0\\ ⋮& \ddots & \ddots & ⋮\\ 0& \dots & -\mathcal{G}\left(T_{N-1},T_{N-2},\cdot \right)& \mathbf{I}\end{array}\right]\left[\begin{array}{c}{\mathbf{U}}_0^{(k+1)}\\ {\mathbf{U}}_1^{(k+1)}\\ ⋮\\ {\mathbf{U}}_{N-1}^{(k+1)}\end{array}\right]=\left[\begin{array}{c}{\mathbf{b}}_N^{(k)}\\ {\mathbf{b}}_1^{(k)}\\ ⋮\\ {\mathbf{b}}_{N-1}^{(k)}\end{array}\right],$$

$$\left[\begin{array}{c}{\mathbf{U}}_0^{(k+1,s+1)}\\ {\mathbf{U}}_1^{(k+1,s+1)}\\ ⋮\\ {\mathbf{U}}_{N-1}^{(k+1,s+1)}\end{array}\right]=\left[\begin{array}{cccc}0& 0& \dots & \mathcal{G}\left(T_N,T_{N-1},\cdot \right)\\ \mathcal{G}\left(T_1,T_0,\cdot \right)& 0& & 0\\ ⋮& \ddots & \ddots & ⋮\\ 0& \dots & \mathcal{G}\left(T_{N-1},T_{N-2},\cdot \right)& 0\end{array}\right]\left[\begin{array}{c}{\mathbf{U}}_0^{(k+1,s)}\\ {\mathbf{U}}_1^{(k+1,s)}\\ ⋮\\ {\mathbf{U}}_{N-1}^{(k+1,s)}\end{array}\right]+\left[\begin{array}{c}{\mathbf{b}}_N^{(k)}\\ {\mathbf{b}}_1^{(k)}\\ ⋮\\ {\mathbf{b}}_{N-1}^{(k)}\end{array}\right]$$

$$u(t)=\sum _{j=-\infty }^{\infty }{\hat{u}}_j{\psi }_j(t),t\in [0,T]$$

$${\hat{u}}_j:=\frac{1}{T}\underset{0}{\overset{T}{\int }}u(t){\psi }_j^{\ast }(t)dt=\frac{1}{T}\underset{0}{\overset{T}{\int }}u(t)e^{-ı2\pi jt/T}dt,j\in \mathbb{Z}.$$

$$F(u^{\prime },u,t)=0,t\in (0,T),$$

$$0=(F(u^{\prime },u,t),\psi {)}_{L^2([0,T],\mathbb{C})}:=\frac{1}{T}\underset{0}{\overset{T}{\int }}F(u^{\prime },u,t){\psi }^{\ast }(t)dt\forall \psi \in L^2([0,T],\mathbb{C}),$$

$$u(t)\approx \sum _{j=1}^N{\hat{u}}_j{e}^{ı{\omega }_{j}t},t\in [0,T]$$

$${\omega }_j=2\pi p_j/T,j=1,\dots ,N$$

$$\vec{p}=\big{[}-\left\lfloor{N/2}\right\rfloor+\delta,\dots,\left\lfloor{N/2}\right\rfloor\big{]}^{\!\top}\in\mathbb{R}^{N},$$

$$\delta =\{\begin{array}{ll}1,& \mathrm{if}N\mathrm{is}\mathrm{even}\\ 0,& \mathrm{if}N\mathrm{is}\mathrm{odd}\end{array}$$

$$0=\frac{1}{T}\underset{0}{\overset{T}{\int }}\hat{F}(\hat{\mathbf{u}},t)e^{-ı{\omega }_{j}t}dt,j=1,\dots ,N$$

$${\hat{u}}_j=\frac{1}{N}\sum _{n=0}^{N-1}u(T_n)e^{-ı{\omega }_j{T}_n},j=1,\dots ,N$$

$$\displaystyle\mathbf{M}\mbox{\boldmath$u$}^{\prime}(t)+\mathbf{K}\mbox{\boldmath$u$}(t)$$

$$\displaystyle\mbox{\boldmath$u$}(0)$$

$$\underbrace{\left[\frac{1}{\Delta T}\mathbf{M}+\mathbf{K}\right]}_{=:\mathbf{Q}}\mathcal{G}\left(T_{n},T_{n-1},\mathbf{U}_{n-1}\right)=\underbrace{\frac{1}{\Delta T}\mathbf{M}}_{=:\mathbf{C}}\mathbf{U}_{n-1}+\mbox{\boldmath$j$}(T_{n}),\quad n=1,\dots,N$$

$$\underset{=:\mathbf{G}}{\underbrace{\left[\begin{array}{cccc}\mathbf{Q}& & & -\mathbf{C}\\ -\mathbf{C}& \mathbf{Q}& & \\ & \ddots & \ddots & \\ & & -\mathbf{C}& \mathbf{Q}\end{array}\right]}}\left[\begin{array}{c}{\mathbf{U}}_0^{(k+1)}\\ {\mathbf{U}}_1^{(k+1)}\\ ⋮\\ {\mathbf{U}}_{N-1}^{(k+1)}\end{array}\right]=\underset{=:{\mathbf{r}}^{(k)}}{\underbrace{\left[\begin{array}{c}{\mathbf{r}}_N^{(k)}\\ {\mathbf{r}}_1^{(k)}\\ ⋮\\ {\mathbf{r}}_{N-1}^{(k)}\end{array}\right]}},$$

$$\mathbf{r}_{n}^{(k)}:=\mathbf{Q}\mathbf{b}_{n}^{(k)}+\mbox{\boldmath$j$}(T_{n}),$$

$${\mathbf{U}}_n^{(k+1)}=\sum _{j=1}^N{\hat{\mathbf{U}}}_j^{(k+1)}{e}^{ı{\omega }_j{T}_n},n=0,\dots ,N-1,$$

$$\underset{=:\hat{\mathbf{G}}}{\underbrace{\left(\tilde{\mathbf{F}}\mathbf{G}{\tilde{\mathbf{F}}}^{\mathrm{H}}\right)}}{\hat{\mathbf{U}}}^{(k+1)}=\underset{=:{\hat{\mathbf{r}}}^{(k)}}{\underbrace{\tilde{\mathbf{F}}{\mathbf{r}}^{(k)}}}$$

$${\mathbf{F}}_{jq}=e^{-ı{\omega }_j{T}_{q-1}}/\sqrt{N},j,q=1,\dots ,N,$$

$${\hat{\mathbf{G}}}_{jj}=\mathbf{Q}-\mathbf{C}{e}^{-ı\Delta T{\omega }_j},j=1,\dots ,N.$$

$${\hat{\mathbf{G}}}_{jj}{\hat{\mathbf{U}}}_j^{(k+1)}={\hat{\mathbf{r}}}_j^{(k)},j=1,\dots ,N.$$

$${\mathbf{U}}^{(k+1)}={\tilde{\mathbf{F}}}^{\mathrm{H}}{\hat{\mathbf{U}}}^{(k+1)}.$$

$$\left[\frac{1}{\delta T}\mathbf{M}+\mathbf{K}\right]\mbox{\boldmath$u$}_{n}=\frac{1}{\delta T}\mathbf{M}\mbox{\boldmath$u$}_{n-1}+\mbox{\boldmath$j$}(t_{n}),\quad n=1,\dots,N_{\mathrm{f}},$$

Full text

{AMS}

34A34, 34A36, 34A37, 65L20, 78M10

1 Introduction

The computation of the time-periodic solution of a stiff evolution problem becomes particularly time consuming, when one is interested in the steady-state behavior e.g., of an electrical machine. The numerical treatment of time-periodic problems is often challenging since it requires the simultaneous solution of the underlying system on the whole period. Indeed, a classical methodology involves a discretization in space, e.g., using the finite element method (FEM), followed by a finite difference discretization in time, which couples all the discrete values on one period through the periodicity condition. This space-time discretization approach of time-periodic problems is known in engineering as the time-periodic (TP) FEM [17]. It leads to a large system of algebraic equations whose solution would usually be very computationally expensive or even prohibitive, especially due to its troublesome block-cyclic structure. These obstacles on the way to solve time-periodic problems cause an urgent need to develop efficient parallel algorithms, able to simplify and accelerate the time-domain computations.

Based on the idea of the multiple shooting method [30], the Parareal algorithm was originally introduced to speed up sequential solution of initial value problems (IVPs) using parallelization in the time domain [22]. It involves computations on two grids (fine and coarse). While the coarse solver gives rough information about the solution and is applied sequentially, the acceleration is obtained via parallel calculation of the accurate fine solution. Convergence properties of the method were investigated in [15], [9]. Application of Parareal to the simulation of an electrical machine in [29] illustrated its significant acceleration capabilities. The original Parareal approach for IVPs was subsequently extended to the class of time-periodic problems and analyzed in [12].

This work introduces a novel Parareal-based iterative algorithm able to solve nonlinear time-periodic problems efficiently. In contrast to the fixed point iteration introduced in [12] a special linearization via the simplified Newton method is proposed here. A specific choice of the initial approximation for Newton’s iteration gives a block-cyclic system matrix, which can be diagonalized using the frequency domain approach from [5]. Following the terminology of [2] we call it the multi-harmonic (MH) solver in this paper. This technique was recently used for parallel-in-time solution of linear time-periodic problems in [20]. It was applied on the coarse level of the time-periodic Parareal system, calculating each frequency component separately and in parallel. Diagonalization using the fast Fourier transform (FFT) [31] makes the MH approach particularly attractive due to its low complexity. The Fourier-based diagonalization was exploited to find a time-periodic solution of a fractional diffusion equation in [33], where each separate system of equations was solved using the multigrid method. Besides, a similar idea of diagonalization and parallelization was proposed for solution of IVPs in [10], [11], where non-equidistant partition of the considered time interval was required. Finally, the MH approach was recently applied to the parallel-in-time solution of IVPs in [14]. There an additional restructuring to a block-cyclic form was performed.

Our paper is organized as follows. At first, a couple of existing parallel-in-time approaches for time-periodic problems [12] are described in Section 2. Section 3 recalls the main ideas of the frequency domain solution of linear time-periodic problems. In Section 4 the MH approach presented in [5] as well as its parallel-in-time extension [20] are described. A special Newton-type solver is then introduced in Section 5 for solution of nonlinear time-periodic problems with the MH approach. Its convergence properties are analyzed using the theory from [8] and illustrated for a scalar model problem. Application of the proposed method to numerical treatment of the eddy current problem for a coaxial cable model, as well as its performance compared to the conventional sequential time stepping, to the time-parallel algorithms from [12], and to the MH solution of [5] are presented in Section 6. Section 7 describes performance of the considered approaches applied to a three-dimensional model of a transformer. The paper is finally concluded with Section 8.

2 Overview of the periodic Parareal-based parallel-in-time algorithms

Consider a time-periodic problem for a system of ordinary differential equations (ODEs) on time interval $(0,T)$

[TABLE]

with unknown \mbox{\boldmathu}:[0,T]\to\mathbb{R}^{d}, matrices $\mathbf{M}$ and \mathbf{K}(\mbox{\boldmathu}(t)) such that the matrix pencil \bigl{(}\mathbf{M},\mathbf{K}(\mathbf{v})\bigr{)}=\mathbf{K}(\mathbf{v})+\lambda\mathbf{M}, with $\lambda\in\mathbb{R}$ is regular for each $\mathbf{v}\in\mathbb{R}^{d},$ and a $T$-periodic right-hand side (RHS) \mbox{\boldmathj}:[0,T]\to\mathbb{R}^{d}. System (1) could stem, e.g., from a spatial discretization of a parabolic partial differential equation (PDE) discretized by FEM [7].

As in the classical Parareal method [22], we initially split the time interval $[0,T]$ into $N$ subintervals $[T_{n-1},T_{n}],$ $n=1,\dots,N$ using partition $0=T_{0}<T_{1}<\dots<T_{N}=T.$ An IVP for unknown \mbox{\boldmathu}_{n}:[T_{n-1},T_{n}]\to\mathbb{R}^{d} given by

[TABLE]

is considered on the $n$th subinterval, $n=1,\dots,N$. Within the Parareal setting, two propagators are applied to the $n$th IVP: the fine \mathcal{F}\bigl{(}t,T_{n-1},{\mathbf{U}_{n-1}}\bigr{)} and the coarse \mathcal{G}\bigl{(}t,T_{n-1},{\mathbf{U}_{n-1}}\bigr{)}, [15]. Both operators calculate a solution to (2) at $t\in(T_{n-1},T_{n}],$ starting from the initial value \mbox{\boldmathU}_{n-1} at $T_{n-1},$ $n=1,\dots,N.$ However, accuracies and therefore computational costs of the two solvers differ. In particular, $\mathcal{F}$ solves the IVP with a very high precision, e.g., via time stepping over a very fine grid, while $\mathcal{G}$ provides only a rough approximation of the solution, e.g., by using a very coarse discretization or a lower-order method, compared to the fine propagator.

Based on the idea of Parareal, Gander et al. [12] introduce two parallel-in-time algorithms for efficient treatment of time-periodic problems. The first one, called periodic Parareal algorithm with initial-value coarse problem (PP-IC), calculates the periodic solution using the iteration

[TABLE]

for $k=0,1,\dots$ and $n=1,\dots,N.$ While in the classical Parareal method the initial value (at $T_{0}$) is fixed, (3)-(4) updates the solution at $T_{0}$ with the solution at $T_{N},$ obtained from the previous iteration. In this way PP-IC weakens the periodic coupling among the discrete values at the synchronization points $T_{n},$ $n=0,\dots,N-1$ and, like Parareal, it involves the solution of a sequence of IVPs (and not of a single time-periodic problem) on the coarse level. The method has been recently applied to simulation of an induction motor for an electric vehicle drive by the authors in [3] and already delivered a significant speedup. However, since PP-IC imposes a relaxed periodicity constraint one may expect a rather slow convergence to the periodic solution, especially when the underlying dynamical system possesses a very long settling time. This is observed in the original paper [12] as well as within our numerical experiments in Section 6.

The second approach presented in [12] is PP-PC: periodic Parareal algorithm with periodic coarse problem. In contrast to PP-IC, it maintains the prescribed periodic coupling between the first and the last values of the period $[T_{0},T_{N}]$. The PP-PC iterations are written for $k=0,1,\dots$ as follows:

[TABLE]

with $n=1,\dots,N-1.$ It can be seen that the fine solutions only involve values from iteration $k$ and can therefore be computed in parallel, in the same way as within Parareal and PP-IC. On the other hand, the coarse propagation can neither be performed sequentially nor in parallel, since the periodicity condition yields interconnection of the solution at all synchronization points. Hence, PP-PC requires the solution of a periodic coarse problem, which implies a joint computation of the coarse grid values over the whole period.

Denoting the difference between the fine and the coarse solutions of IVP (2) at $T_{n},$ calculated starting from the initial value $\mathbf{U}^{(k)}_{n-1}$ at $T_{n-1}$, by

[TABLE]

we write PP-PC the iteration (5)-(6) in a matrix-vector (operator) form

[TABLE]

where $\mathbf{I}\in\mathbb{R}^{d\times d}$ denotes the identity matrix. Note that the system matrix of (8) might have large dimensions, since each of the unknown $\mathbf{U}_{n}^{(k+1)},$ $n=0,\dots,N-1$ possibly consists of many degrees of freedom, e.g., coming from a spatial discretization. The matrix also possesses a block-cyclic-like structure, which is typical for time-periodic systems but inconvenient for many linear solvers.

System (8) is written in an implicit operator form and cannot be solved directly in the nonlinear case. A Jacobi-type fixed point iteration was applied in [12] to iteratively solve the PP-PC system. It is written for $s=0,1,\dots$ as

[TABLE]

at each PP-PC iteration $k+1$. The fixed point iteration (9) decouples the calculation of the values $\mathbf{U}_{n}^{(k+1,s+1)},$ $n=0,\dots,N-1$ completely, thereby allowing parallel solution even on the coarse level. Another advantage of this iterative method is that the is fully non-intrusive, i.e., the propagators $\mathcal{F}$ and $\mathcal{G}$ can be considered as black-box solvers. On the other hand, such a lineariazation relaxes the periodicity constraint similarly to PP-IC and thus a loss of efficiency may be expected. This becomes obvious in the linear case where (9) is equivalent to the block Jacobi method for linear systems.

Our aim is to construct an iterative method, which benefits from the present block-cyclic structure and provides fast convergence.

3 MH frequency domain solution of time-periodic problems

The frequency domain representation of (1) is often exploited in (electrical) engineering since it transforms a differential equation into an algebraic one [24]. In contrast to the possibly lengthy time stepping in the time domain, MH frequency domain computations allow to obtain the steady-state solution directly via calculating its harmonic components. However, it may require many basis functions if the solution exhibits local features, e.g., due to pulsed excitations [4, 23]. This frequency domain solution approach is also called the harmonic balance method [16], [1] and can be interpreted as the Fourier collocation [8] or the Fourier spectral method [31].

Let $u$ be a periodic function from the Lebesgue space $L^{2}([0,T],\mathbb{C}),$ i.e., $u:[0,T]\to\mathbb{C}$ is a measurable function with a square integrable norm $\|u(t)\|_{\mathbb{C}}$ on $[0,T]$ (see [34, Chapter 23.2]) and $u(0)=u(T)$. The main idea of the MH solution [5] is to represent $u$ with the infinite sum

[TABLE]

of the spectral (orthonormal) basis functions $\psi_{j}(t)=e^{\imath{{2\pi}jt/{T}}}$ on time interval $[0,T],$ with (Fourier) coefficients

[TABLE]

Searching for the $T$-periodic solution $u$ of an ODE

[TABLE]

we write its variational formulation

[TABLE]

where it is assumed that $u^{\prime}\in L^{2}([0,T],\mathbb{C}^{*}),$ with $\mathbb{C}^{*}$ denoting the dual space of the space of complex numbers $\mathbb{C}.$ The unknown solution $u$ is then approximated with the finite Fourier series

[TABLE]

using frequencies

[TABLE]

from the double-sided spectrum given by $N$-dimensional vector

[TABLE]

where

[TABLE]

and $\lfloor\cdot\rfloor$ denotes the floor function, which returns the greatest integer less than or equal to its argument. This yields a finite-dimensional system of equations

[TABLE]

with respect to the unknown coefficients $\hat{\mathbf{u}}=\left[\hat{u}_{{1}},\dots,\hat{u}_{{N}}\right]^{\!\top}$ in the frequency domain. Here $\hat{F}(\hat{\mathbf{u}},t)$ denotes the restriction of $F(u^{\prime},u,t)$ in (11) after substitution of $u$ from (13) therein. Using the calculated frequency components $\hat{u}_{j},$ expansion (13) gives the solution in the time domain. Note that the approximation of (10) using the left rectangle quadrature rule on a partition $0=T_{0}<T_{1}<\dots<T_{{N-1}}<T$

[TABLE]

is the discrete Fourier transform (DFT) of vector $\mathbf{u}=\left[u(T_{0}),\dots,u(T_{{N-1}})\right]^{\!\top},$ while (13) evaluated at $t=T_{n},$ $n=0,\dots,{N-1}$ gives the inverse DFT of $\hat{\mathbf{u}}.$ In the following section we exploit the MH discretization idea to solve PP-PC system (8) for a linear time-periodic problem.

4 MH solver for linear time-periodic systems

In this section we describe a combination of the time-domain solution approaches such as PP-PC (5)-(6) and direct TP discretization with the frequency-domain representation from Section 3 for linear time-periodic problems, presented in [20] and [5], respectively.

4.1 Linear PP-PC MH approach

Based on [20] we recall the application of the MH approach [5] to parallel-in-time solution of a linear non-autonomous problem of the form

[TABLE]

with matrices $\mathbf{M}$ and $\mathbf{K}$ having a regular matrix pencil $(\mathbf{M},\mathbf{K})$ and RHS \mbox{\boldmathj}(t). Discretization by the implicit Euler method on the equidistant coarse grid $0=T_{0}<T_{1}<\dots<T_{N}=T$ defines the coarse solution $\mathcal{G}\left(T_{n},T_{n-1},\mathbf{U}_{n-1}\right)$ at $T_{n}$ by

[TABLE]

with step size $\Delta T=T/N.$ This yields PP-PC system (8) in the explicit matrix-vector form

[TABLE]

where the RHS is defined by

[TABLE]

and vector $\mathbf{b}_{n}^{(k)},$ $n=1,\dots,N$ is given by (7). Originating from (8), system matrix $\mathbf{G}$ has the inconvenient block-cyclic structure. The MH solution approach [5] overcomes this difficulty by decoupling the solution on the coarse level without introducing an additional iteration, in contrast to (9).

Since $\mathbf{U}_{0}^{(k+1)},\dots,\mathbf{U}_{N-1}^{(k+1)}$ are values of a periodic function at $T_{0},\dots,T_{N-1}$ one could express them with the finite Fourier series expansion (13) discretized at the synchronization points, i.e.,

[TABLE]

with frequencies $\omega_{j}$ given by (14). Plugging (21) into the PP-PC system (20) and applying DFT (17) gives an equivalent system in frequency domain

[TABLE]

in terms of the unknown Fourier coefficients $\hat{\mathbf{U}}_{j}^{(k+1)},$ $j=1,\dots,N$ contained in the joint vector $\hat{\mathbf{U}}^{(k+1)}\in{\mathbb{C}}^{Nd}$. Matrix $\tilde{\mathbf{F}}$ is defined by $\tilde{\mathbf{F}}=\mathbf{F}\otimes\mathbf{I},$ where $\mathbf{F}$ is the DFT matrix with elements

[TABLE]

$\mathbf{I}$ is a $(d\times d)$-dimensional identity matrix, and ‘$\otimes$’ denotes the Kronecker product of two matrices. $\tilde{\mathbf{F}}^{\mathrm{H}}$ is the Hermite conjugate matrix of $\tilde{\mathbf{F}}$.

We emphasize that the transformation from (20) to (22) comes along with an important property: it transforms the block-cyclic system matrix $\mathbf{G}$ into the block-diagonal system matrix $\hat{\textbf{G}}$. Matrices on the diagonal of $\hat{\textbf{G}}$ can be explicitly calculated as

[TABLE]

Subsequently, the harmonic components can be found solving

[TABLE]

Therefore, additional parallelization can be obtained on the coarse level, i.e., solution of $N$ systems of $d$ linear equations could be performed in parallel, without introduction of an additional inner loop. Note that matrices $\tilde{\mathbf{F}}$ and $\tilde{\mathbf{F}}^{\mathrm{H}}$ do not have to be explicitly constructed, since the system matrices in (25) are determined by (24) for each frequency $\omega_{j}$.

Finally, the solution in the time domain can be obtained by application of the inverse DFT to the calculated solution vector $\hat{\mathbf{U}}^{(k+1)},$ i.e.,

[TABLE]

Note that the DFT and its inverse can be efficiently applied using the FFT algorithm [31], thereby further reducing the complexity of the transformation.

4.2 Linear TP MH approach

Clearly, the transformation into frequency domain of the form (22) can be also applied directly within the discrete TP setting as it was performed in [5]. In this case application of implicit Euler’s method to (18) on the fine grid $0=T_{0}<t_{1}<\dots<t_{N_{\mathrm{f}}}=T$

[TABLE]

together with the periodicity condition \mbox{\boldmathu}_{0}=\mbox{\boldmathu}_{N_{\mathrm{f}}} leads to the block-cyclic system

[TABLE]

where $\mathbf{C}_{\mathrm{f}}:={1}/{\delta T}\cdot\mathbf{M}$ and $\delta T=T/N_{\mathrm{f}}.$ Here \mbox{\boldmathu}_{n}\approx\mbox{\boldmathu}(t_{n}) denotes the discrete solution at the time instant $t_{n},$ $n=1,\dots,N_{\mathrm{f}}.$ Solution of system (28) arises in the TP FEM framework, when using the terminology of [17]. We note that the fine Euler time step in (28) is much smaller than the coarse one in the PP-PC system (20), i.e., $\delta T\ll\Delta T,$ which means that system (28) has a much bigger size compared to that of (20).

Similarly to (22) we can apply the MH solver to (28) and compute each harmonic coefficient \hat{\mbox{\boldmathu}}_{j}, $j=1,\dots,N_{\mathrm{f}}$ separately. This way to solve the time-periodic problem we call TP MH. In this (linear) case the TP MH approach effectively solves only one $d$-dimensional linear system, provided sufficient parallelization is possible, i.e., when (at least) $N_{\mathrm{f}}$ cores are available. Besides, this approach does not include any iteration, in contrast to PP-PC (20), but allows to parallelize the computations only through the decoupling of the frequency components.

5 A novel iterative solver for nonlinear time-periodic systems

We now introduce an iterative algorithm which allows to treat also nonlinear PP-PC and discrete TP formulations in an efficient manner using the MH diagonalization technique.

5.1 Nonlinear PP-PC MH approach

We now consider the following time-periodic problem for a system of nonlinear ODEs

[TABLE]

with nonlinearity present in matrix \mathbf{K}(\mbox{\boldmathu}(t)), e.g., due to nonlinear material characteristics. Similarly to (20), coarse discretization with implicit Euler’s method on an equidistant grid leads to the PP-PC system

[TABLE]

with the following definitions

[TABLE]

where $\mathbf{b}_{n}^{(k)},$ $n=1,\dots,N$ is given by (7) and $\Delta T=T/N$. We point out the nonlinear dependence of $\mathbf{K}$ on \mathcal{G}\bigl{(}T_{n},T_{n-1},\mathbf{U}_{n-1}\bigr{)} in \mathbf{Q}_{n}\bigl{(}\mathbf{U}_{n-1}\bigr{)} in contrast to the implicit Euler discretization, where the nonlinear dependence on the variable \mbox{\boldmathu}_{n} directly is present (see Section 5.2). This is due to the application of the coarse propagator $\mathcal{G}$ within PP-PC, seen in the operator equation (8). Substituting \mathcal{G}\bigl{(}T_{n},T_{n-1},\mathbf{U}_{n-1}^{(k+1)}\bigr{)}, $n=1,\dots,N$ from PP-PC iteration (5)-(6) into (30) and omitting superscript $k+1$ we obtain the following nonlinear system for $\mathbf{U}=\left[\mathbf{U}_{0}^{\!\top},\dots,\mathbf{U}_{N-1}^{\!\top}\right]^{\!\top}$

[TABLE]

where

[TABLE]

with $\mathbf{b}_{n}^{(k)},$ $n=1,\dots,N$ from (7). To solve system (31) at PP-PC iteration $k+1$ we search for the root of mapping \mathbf{R}^{{\color[rgb]{0,0,0}k}}:\mathbb{R}^{Nd}\rightarrow\mathbb{R}^{Nd} defined by

[TABLE]

using the simplified Newton iteration [8]: for $s=0,1,\dots$ compute $\mathbf{U}^{(s+1)}=\mathbf{U}^{(k+1,s+1)}$ from

[TABLE]

with given initial approximation $\mathbf{U}^{(0)}=\mathbf{U}^{(k+1,0)}.$ Jacobian matrix \mathbf{J}_{\mathbf{R}}^{{\color[rgb]{0,0,0}k}}=\displaystyle\frac{\mathrm{d}}{\mathrm{d}\mathbf{U}}\mathbf{R}^{{\color[rgb]{0,0,0}k}} is defined by

[TABLE]

where

[TABLE]

We note that applying the simplified Newton method we have to construct the system matrix only once at each PP-PC iteration. Now, in order to apply the MH coarse grid correction to the iteration (33) as in Section 4.1 the block-matrices $\mathbf{Q}_{\mathrm{d}}(\cdot)$ on the diagonal have to be all equal. We therefore choose the initial approximation for the simplified Newton iteration (33) at PP-PC iteration $k+1$ as

[TABLE]

for a given vector $\mathbf{Z}\in\mathbb{R}^{d}.$ The simplified Newton iteration (33) then reads

[TABLE]

with system matrix

[TABLE]

and RHS

[TABLE]

The matrix in (38) remains constant over the Newton iterations and has the same block-cyclic structure as that of the PP-PC system (20) for a linear problem. Hence, a MH solver can be also applied to the Newton system (37). Analogously to (22) we obtain

[TABLE]

where matrix $\hat{\mathbf{G}}_{\mathrm{d}}$ is block-diagonal with each $(d\times d)$-dimensional block given by

[TABLE]

and the DFT matrix $\tilde{\mathbf{F}}$ is as previously defined by (23). The frequency domain solution $\hat{\mathbf{U}}^{(s+1)}$ at Newton iteration $s+1$ is then transformed into the time domain by the inverse DFT, i.e.,

[TABLE]

Remark 5.1**.**

We note that the simplified Newton iteration (33) could be also interpreted as a linear iterative method based on an additive splitting of the system matrix in (31). Indeed, using any constant matrix $\mathbf{H}^{k}\in\mathbb{R}^{Nd\times Nd}$ we can introduce at PP-PC iteration $k+1$ a linear iteration for $s=0,1,\dots$

[TABLE]

which in case of \mathbf{H}^{k}\!:=\mathbf{J}_{\mathbf{R}}^{{\color[rgb]{0,0,0}k}}\big{(}\mathbf{U}^{(0)}\big{)} becomes exactly the proposed simplified Newton iteration (33). We note that that the choice

[TABLE]

*yields Gander’s Jacobi-like fixed-point iteration (9). The idea of choosing a modified fixed point iteration operator to ensure a convenient block structure is well known and is for example applied by Biro [5] and Gander [10]. However, both do not discuss convergence. *

Remark 5.2**.**

We point out that the MH approach is valid not only when applying the implicit Euler method but also other time intergation algorithms. For instance, using the implicit trapezoidal rule (also called the Crank-Nicolson method) in the linear case, the system (20) maintains the same block-cyclic structure with the following definitions of the block matrices and the RHS:

[TABLE]

5.2 Nonlinear TP MH approach

Similarly as in Section 4.2, one could exploit the proposed methodology to solve a time-periodic problem, without applying a parallel-in-time method. More specifically, the fine discretization of the periodic system (29) on $0=T_{0}<t_{1}<\dots<t_{N_{\mathrm{f}}}=T$ using implicit Euler’s method with time step size $\delta T=T/N_{\mathrm{f}}$ leads to the system of nonlinear algebraic equations

[TABLE]

where $\mathbf{C}_{\mathrm{f}}:={1}/{\delta T}\cdot\mathbf{M}.$ Choosing vector $\mathbf{z}\in\mathbb{R}^{d}$ in the initial approximation

[TABLE]

we obtain the simplified Newton iteration of form (37) applied to (44)

[TABLE]

with \mathbf{Q}_{\mathrm{f}}\left(\mbox{\boldmathz}\right):=\mathbf{C}_{\mathrm{f}}+\mathbf{K}_{\mathrm{d}}\left(\mbox{\boldmathz}\right) and $\mathbf{K}_{\mathrm{d}}$ defined in (35). One could apply the MH solver (40) at Newton iteration $s+1$ and solve a linear algebraic system for each harmonic coefficient \hat{\mbox{\boldmathu}}_{j}^{(s+1)}, $j=1,\dots,N_{\mathrm{f}}$ separately. This is to our best knowledge essentially the same idea as the fixed point method proposed in [5].

When many central processing units (CPUs) are available one may expect the parallelization of the MH frequency domain solution at each Newton iteration (46) to outperform the speed up provided by the parallel-in-time solution on the fine grid in PP-PC (30). However, we note that the application of the MH solver to the PP-PC system and not to (44) is especially beneficial when the number of CPUs is limited ($<N_{\mathrm{f}}$), that is, when one could not calculate each harmonic coefficient \hat{\mbox{\boldmathu}}_{j}^{(s+1)}, $j=1,\dots,N_{\mathrm{f}}$ in parallel.

5.3 Pseudocode and implementation details

The error used in the termination criterion for the periodic Parareal-based algorithms (PP-IC and PP-PC) calculates the maximum among the jumps at the synchronization points $T_{n},$ with $n=1,\dots,N-1,$ as well as the periodicity jump of the solution at $T_{0}$ and $T_{N}$, i.e., the error at the periodic Parareal (PP) iteration $k$ is computed as

[TABLE]

where for two $d$-dimensional vectors $\mathbf{u}$ and $\mathbf{v}$ the $\|\cdot\|_{*}$-norm is defined by

[TABLE]

and $\text{aTol}=10^{-6}$ and $\text{rTol}=10^{-3}$ are chosen as the absolute and the relative tolerances, respectively. Here $\|\cdot\|_{2}$ denotes the Euclidean norm in the $d$-dimensinal vector space.

Furthermore, at PP iteration $k+1$ the errors of the inner fixed point (9) and the simplified Newton (33) iterations are obtained from the difference of the coarse solution vectors at the two subsequent iterations $s+1$ and $s$ for $s=0,1,\dots$ through

[TABLE]

This is similar tothe error of the simplified Newton iteration (46), applied directly to the discrete TP formulation, where the error is defined from the (fine) solutions by

[TABLE]

Finally, for $k\geq 1$ the classical sequential time-stepping considers the periodicity error

[TABLE]

after each period $T$ until the periodic steady state is reached. The numerical calculations are terminated when the values of the corresponding errors from (47), (49)-(51) become smaller than $1$.

The considered methods are implemented in GNU Octave version 3.8.2 using the ‘parallel’ package111http://octave.sourceforge.io/parallel for parallelization. We use Octave’s Fourier analysis implemented by the parallel version of FFTW3222http://www.fftw.org/parallel using $80$ threads. The code is executed on an Intel Xeon cluster with $80\times 2.00\text{\,}\mathrm{GHz}$ cores, i.e., $8\times$E7-8850 and 1TB DDR3 memory.

In Algorithm 1 we present the pseudocode of the nonlinear PP-PC MH method, proposed in Section 5.1.

5.4 Convergence analysis of the proposed algorithm

Convergence of the simplified Newton method (33) is determined by the result presented in [8, Theorem 2.5], which originates from [25]. We recall the theorem here in terms of the notations introduced above.

Theorem 5.3** (Convergence of the simplified Newton method [8]).**

Let $\mathbf{R}:D\rightarrow\mathbb{R}^{N}$ be a continuously differentiable mapping with $D\subset\mathbb{R}^{N}$ open and convex. Let $\mathbf{U}^{(0)}\in D$ denote a given starting point so that \mathbf{J}_{\mathbf{R}}^{{\color[rgb]{0,0,0}k}}\left(\mathbf{U}^{(0)}\right) is invertible. Assume the affine covariant Lipschitz condition

[TABLE]

holds for all $\mathbf{U}\in D.$ Let

[TABLE]

and assume that the closure $\bar{S}\left(\mathbf{U}^{(0)},\rho\right)$ of a ball with center in $\mathbf{U}^{(0)}$ and radius $\rho$ is a subset of $D,$ with

[TABLE]

*Then the simplified Newton iterates $\mathbf{U}^{(s+1)},$ $s=0,1,\dots$ generated by (33) remain in $\bar{S}\left(\mathbf{U}^{(0)},\rho\right)$ and converge to some $\mathbf{U}^{*}$ with \mathbf{R}^{{\color[rgb]{0,0,0}k}}\left(\mathbf{U}^{*}\right)=\mathbf{0}. *

5.4.1 Numerical analysis of a model problem

We discuss the convergence analysis of the simplified Newton iteration (33) for the PP-PC algorithm applied to a time-periodic problem for a single ($d=1$) ODE. Motivated by the properties of the nonlinear material curves proposed in [18] we consider the following model problem

[TABLE]

with $m\in\mathbb{R}_{0}^{+},$ $\kappa:\mathbb{R}_{0}^{+}\rightarrow\mathbb{R}_{0}^{+},$ unknown function $u:[0,T]\rightarrow\mathbb{R},$ and $T$-periodic input $j,$ i.e., $j(0)=j(T).$ We assume that the nonlinear function $\kappa$ satisfies the following conditions.

Assumption 5.4**.**

Let for $\kappa:\mathbb{R}_{0}^{+}\rightarrow\mathbb{R}_{0}^{+}$ it holds:

• $\bullet$

$\kappa(t)\geq c_{1}>0,$* $t\in\mathbb{R}_{0}^{+},$*

• $\bullet$

function {\color[rgb]{0,0,0}f}:t\rightarrow\kappa(t)t is strongly monotone with monotonicity constant $c_{1},$ i.e.,

[TABLE]

• $\bullet$

$\kappa\in C^{1}(\mathbb{R}_{0}^{+})$* and $\lim\limits_{t\to\infty}\kappa^{\prime}(t)=0,$*

• $\bullet$

function {\color[rgb]{0,0,0}g}:t\rightarrow\kappa^{\prime}(t)t is Lipschitz continuous, i.e.,

[TABLE]

We now apply the simplified Newton iteration (33) to (55). Analogously to (34) the Jacobian matrix for the one-dimensional (1D) problem (55) is determined at PP-PC iteration $k$ by

[TABLE]

where $\mathbf{U}=\left[u_{0},\dots,u_{N-1}\right]^{\!\top},$ $c=m/\Delta T$ with $\Delta T=T/N,$ and function $\kappa_{\mathrm{d}}:\mathbb{R}\mapsto\mathbb{R}_{0}^{+}$ given by

[TABLE]

As in (36) choosing a fixed value $z\in\mathbb{R}$ we define

[TABLE]

to be the initial approximation for the Newton algorithm at PP-PC iteration $k+1$. We denote the Jacobian matrix \mathbf{J}_{\mathbf{R}}^{{\color[rgb]{0,0,0}k}} at the chosen $\mathbf{U}^{(0)}$ from (60) by $\mathbf{G}_{\mathrm{d}}$ which is given by

[TABLE]

Before demonstrating the convergence results based on the Theorem 5.3 we derive some properties of the function $\kappa_{\mathrm{d}}$ in the following lemma.

Lemma 5.5**.**

Under the Assumption 5.4 on $\kappa$ we have that function $\kappa_{\mathrm{d}}:\mathbb{R}\mapsto\mathbb{R}_{0}^{+}$ defined in (59) is

bounded from below by $c_{1}$, i.e.,

[TABLE]

• 2)

Lipschitz continuous with Lipschitz constant $L_{2}=\kappa_{2}+L_{1},$ where $\kappa_{2}:=\sup\limits_{t\in\mathbb{R}_{0}^{+}}|\kappa^{\prime}(t)|,$ i.e.,

[TABLE]

Proof.

The proof is based on [18], [26].

Based on the proof of Lemma 2.8 in [26] strong monotonicity of $f:t\rightarrow\kappa(t)t,$ $t\in\mathbb{R}_{0}^{+}$ with constant $c_{1}$ implies strong monotonicity of $f_{\text{abs}}:t\rightarrow\kappa(|t|)t,$ $t\in\mathbb{R}$ with the same constant $c_{1},$ which yields $\kappa_{\mathrm{d}}(t)=f_{\text{abs}}^{\prime}(t)\geq c_{1},$ $t\in\mathbb{R}.$

• 2)

First, continuity of $\kappa^{\prime}$ on every closed interval $[0,q],$ $0<q<\infty$ and its vanishing on infinity implies its boundedness, i.e., $|\kappa^{\prime}(t)|\leq\kappa_{2}<\infty,$ $t\in\mathbb{R}_{0}^{+}.$ Then using the mean value theorem we have that $\kappa$ is Lipschitz continuous with Lipschitz constant $\kappa_{2},$ i.e.,

[TABLE]

Using this and the assumption on $g$, we see that $\kappa_{\mathrm{d}}$ is Lipschitz continuous with Lipschitz constant $L_{2}=\kappa_{2}+L_{1}.$

Theorem 5.6**.**

*If the nonlinear function $\kappa$ in (55) satisfies the Assumption 5.4, then the affine covariant Lipschitz condition (52) holds with constant $\delta_{0}={L_{2}}/{c_{1}},$ where $L_{2}$ is the Lipschitz constant of $\kappa_{\mathrm{d}}.$ *

Proof.

First, we observe that for a chosen $\mathbf{U}^{(0)}$ in (60) the Jacobian \mathbf{J}_{\mathbf{R}}^{{\color[rgb]{0,0,0}k}}\bigl{(}\mathbf{U}^{(0)}\bigr{)}=\mathbf{G}_{\mathrm{d}} given in (61) can be decomposed as

[TABLE]

with unitary DFT matrix $\mathbf{F}$ defined in (23) and diagonal matrix $\hat{\mathbf{G}}_{\mathrm{d}}$ defined in (41) containing eigenvalues of $\mathbf{G}_{\mathrm{d}}$. This implies

[TABLE]

where $\|\cdot\|_{2}=\sigma_{\mathrm{max}}(\cdot)$ denotes the largest singular value of a matrix, which for a normal matrix such as $\mathbf{G}_{\mathrm{d}}^{-1}$ is equal to its spectral radius. Using (66) we have for the Lipschitz condition (52) that

[TABLE]

Now, to show convergence of the simplified Newton method (33) one has to prove that expression (67) can be estimated by $\delta_{0}\left\|\mathbf{U}-\mathbf{U}^{(0)}\right\|_{2},$ $\delta_{0}>0$ from above. Indeed, since due to the Lemma 5.5 boundedness condition (62) holds for $\kappa_{\mathrm{d}}$ and since $c\in\mathbb{R}_{0}^{+},$ we have for $1\leq j\leq N$

[TABLE]

Due to the Lipschitz condition (63) and using (68) in expression (67) we directly obtain the estimate (52) with $\delta_{0}={L_{2}}/{c_{1}}$.

Having the constant $\delta_{0}$ from the Lipschitz condition (52) by the Theorem 5.3 the convergence of the simplified Newton method (33) is then implied, provided that the bound (53) holds for the selected initial approximation $\mathbf{U}^{(0)}$.

From the Theorem 5.6 we have seen that the convergence result of the simplified Newton iteration (33) is derived directly from the properties of the nonlinearity of the problem. However, one should be aware of the fact that it may not always be possible to find an initial guess $\mathbf{U}^{(0)}$ which is at the same time close enough to the solution for convergence of the Newton iteration and has the form (60) for applicability of the MH solver to system (37) with constant matrix $\mathbf{G}_{\mathrm{d}}$.

Remark 5.7**.**

We note that an extension of our convergence results to the $d$-dimensional problem (29), obtained from a spatial discretization of a PDE, could be based on the operator theory provided in [26]. In this case one would need to consider the properties not only of the function $\kappa$ but also of a nonlinear operator $A\!:V\mapsto V^{*},$ which, e.g., for a diffusion equation with nonlinear diffusion coefficient $\kappa(|\nabla u|)$ is defined as

[TABLE]

*with a reflexive Banach space $V$ and a Hilbert space $H,$ e.g., $V=H_{0}^{1}(\Omega)$ and $H=L^{2}(\Omega)$ on domain $\Omega.$ Nevertheless, the necessary conditions for the operator $A$ can be derived directly from the properties of $\kappa$, which are exactly those imposed in the Assumption 5.4. We refer to [26, Corollary 2.13, Theorem 4.2 and Theorem 4.3] for the details. *

5.4.2 Numerical study of a model problem

We now discuss convergence of the proposed simplified Newton method applied to a model 1D problem (55) with $m=10^{-1},$ RHS $j(t)=10^{-3}\sin(2\pi t/T),$ $T=0.02,$ and nonlinearity $\kappa$ defined by a piecewise polynomial

[TABLE]

Such a problem can be interpreted as a nonlinear RL-circuit model with unknown magnetic flux $u$ as in [13]. The nonlinear functions $\kappa$ and $\kappa_{\mathrm{d}}$ are depicted in Figure 1 (left). This shows that the lower bound for both $\kappa(|x|)$ and $\kappa_{\mathrm{d}}(x)$, $x\in\mathbb{R}$ is $c_{1}=1$. Figure 1 (right) depicts the corresponding derivatives $\kappa^{\prime}(|x|)$ and $|\kappa_{\mathrm{d}}^{\prime}(x)|,$ $x\in\mathbb{R}.$ Based on the observed characteristics of $\kappa$ one can clearly deduce that the Assumption 5.4 holds. Then by Lemma 5.5 the Lipschitz condition (63) is satisfied for $\kappa_{\mathrm{d}}$. We estimate the Lipschitz constant $L_{2}$ as the upper bound for $|\kappa_{\mathrm{d}}^{\prime}(x)|\leq L_{2}=0.6.$ Based on Theorem (5.6) we have the estimate for $\delta_{0}=L_{2}/c_{1}=6.$

The periodic steady-state solution of (55) with $\kappa$ from (70) is depicted in Figure 2 (left). It is obtained after sequential calculation with step size $\delta T=10^{-5}$ starting from the initial value $u_{\text{init}}=0$ at $t=0$ over $k^{*}=10$ periods when the tolerance defined in (51) is reached. We now describe performance of the simplified Newton algorithm (within TP MH setting) for the 1D model problem depending on the choice of the initial approximation \mbox{\boldmathu}^{(0)}=[z,z,\dots,z]^{\!\top}{\in D=(-0.25,0.25)^{N}\subset\mathbb{R}^{N}} with $N=10.$ In Figure 2 on the right we show the number of iterations required until convergence with respect to the norm (50). We observed that in the neighborhood of $z=0$ as well as for $|z|\geq 0.2$ the algorithm converged in $2$ iterations, while for $0.1<|z|<0.2$ the number of iterations was between $3$ and $5$. This shows us that for the considered example the simplified Newton method is convergent for all the considered initial approximations \mbox{\boldmathu}^{(0)}. The calculated periodic solution is \mbox{\boldmathu}^{*}\approx 10^{-5}\cdot[-3.02,-1.8,0.1,1.9,3.07,3.01,1.8,-0.1,-1.96,-3.07]^{\!\top} in all cases.

For each \mbox{\boldmathu}^{(0)} we have calculated the values of $h_{0}=\delta_{0}\rho_{1}=0.6\rho_{1}$ from (53) with \rho_{1}=\|\mbox{\boldmathu}^{(1)}-\mbox{\boldmathu}^{(0)}\| which we depicted in Figure 3. We see that condition (53) holds, since for all cases $h_{0}<0.5.$ Then by the Theorem 5.3 the iterates \mbox{\boldmathu}^{(s)} remain in the ball \bar{S}\left(\mbox{\boldmathu}^{(0)},\rho\right) of radius $\rho$ defined in (54), which we have also observed in practice.

6 Numerical experiments for a two-dimensional coaxial cable model

In this section we present the results obtained from application of the proposed algorithm with MH coarse correction to solve the time-periodic problem for a coaxial cable model. Figure 3 (left) illustrates a sketch of a two-dimensional (2D) domain

[TABLE]

It represents the cross-section of the cable, which consists of steel tube $\Omega_{\mathrm{Fe}},$ conducting wire $\Omega_{\mathrm{Cu}},$ and air gap $\Omega_{\mathrm{Air}}.$

The electromagnetic phenomena in the coaxial cable when the inner wire is supplied by a sinusoidal current source are mathematically described by the eddy current problem [19]. In the following subsection we discuss the mathematical modeling of the underlying physical characteristics and describe the properties of the nonlinearity.

6.1 Eddy current problem

A time-periodic eddy current formulation in terms of the magnetic vector potential $\vec{A}$ [19] searches for unknown $\vec{A}(\vec{x},t),$ with $(\vec{x},t)\in\Omega\times[0,T],$ which solves

[TABLE]

where $\partial\Omega$ denotes the boundary of $\Omega$ and $\vec{n}$ is the outward normal vector to $\partial\Omega$. The function $\sigma(\vec{x})\geq 0,$ $\vec{x}\in\Omega$ describes the electric conductivity of the materials. It is positive only for $\vec{x}\in\Omega_{\mathrm{Fe}}$ and is equal to zero in $\Omega\setminus{\Omega}_{\mathrm{Fe}}.$ This gives a parabolic-elliptic character to (72). The sinusoidal current excitation defines the RHS

[TABLE]

where $\mathbbm{1}_{\Omega_{\mathrm{Cu}}}$ denotes the indicator function of subdomain $\Omega_{\mathrm{Cu}}.$ The input $\vec{j}$ is described using winding functions [28]. The magnetic reluctivity \nu\bigl{(}\vec{x},|\nabla\times\vec{A}|\bigr{)}>0 is defined by the ferromagnetic material in $\Omega_{\mathrm{Fe}}$ and equals the reluctivity of vacuum $\nu_{0}=$10^{7}/(4\pi)\text{,}\mathrm{m}\text{,}{\mathrm{H}}^{-1}$$ in $\Omega\setminus{\Omega}_{\mathrm{Fe}}.$

The reluctivity function in $\Omega_{\mathrm{Fe}}$ is determined by the magnetization curve $b:\mathbb{R}^{+}_{0}\to\mathbb{R}^{+}_{0}$, shown in the Figure 3 (right). It is commonly obtained from measured tuples $(|\vec{B}|,|\vec{H}|)$ so that |\vec{B}|=b\bigl{(}|\vec{H}|\bigr{)}, where $|\vec{B}|$ denotes the magnitude of the magnetic flux density $\vec{B}$ and $|\vec{H}|$ is the magnitude of the magnetic field intensity $\vec{H}$. To obtain a continuously differentiable curve interpolation of the measurement data is performed, e.g., using classical splines or specific approximation techniques as in [27] such that the monotonicity is preserved. We can then define the reluctivity

[TABLE]

For each $\vec{x}\in{\Omega}_{\mathrm{Fe}}$ we can construct a continuously differentiable reluctivity function $\nu:\mathbb{R}^{+}_{0}\to\mathbb{R}^{+}_{0}$ from the magnetization curve as [18]

[TABLE]

Based on the physical nature of magnetization curve $b$ we can naturally impose the following assumptions on the reluctivity $\nu$ [18].

Assumption 6.1** (Conditions on the magnetic reluctivity).**

Let for $\nu:\mathbb{R}^{+}_{0}\to\mathbb{R}^{+}_{0}$ it holds:

• $\bullet$

$\nu(t)\geq c_{1}>0,$* $t\in\mathbb{R}_{0}^{+},$*

• $\bullet$

function {\color[rgb]{0,0,0}f}:t\rightarrow\nu(t)t is strongly monotone with monotonicity constant $c_{1},$ see (56),

• $\bullet$

$\nu\in C^{1}(\mathbb{R}_{0}^{+})$* and $\lim\limits_{t\to\infty}\nu^{\prime}(t)=0,$*

• $\bullet$

function {\color[rgb]{0,0,0}g}:t\rightarrow\nu^{\prime}(t)t satisfies Lipschitz condition (57).

We note that the conditions on $\nu$ are the same as those imposed on the function $\kappa$ in Assumption 5.4. Analagously to the Remark 5.7 these properties of the nonlinear reluctivity function $\nu$ define conditions on the nonlinear operator in (72), necessary for convergence of the simplified Newton algorithm (37).

For the numerical solution of (72)-(74) the spatial discretization is performed, e.g., using FEM with linear shape functions [7]. This leads to a time-dependent system of differential-algebraic equations (DAEs) [21] (due to the presence of the non-conducting domain $\Omega\setminus{\Omega}_{\mathrm{Fe}}$), which together with the periodicity constraint has the form of (29), with a (singular) mass matrix $\mathbf{M}\in\mathbb{R}^{d\times d},$ nonlinear stiffness matrix $\mathbf{K}(\cdot):\mathbb{R}^{d}\to\mathbb{R}^{d\times d}$ and unknown \mbox{\boldmathu}:[0,T]\to\mathbb{R}^{d}.

We consider $T=0.02$ s to be the period in the eddy current formulation (72)-(74). In our numerical simulations a FEM discretization of the domain $\Omega$ with $d=2269$ degrees of freedom is exploited. The time integration is performed with the implicit Euler method. The DAE is of index-1, therefore in this case it does not need any extra care, e.g., with respect to consistent initial values [29]. The nonlinear algebraic system at each time step is solved using the Newton method. The time step size $\delta T=10^{-5}$ s is used for the sequential calculation, for TP discretization (44), and also as a step size on the fine level of the considered parallel-in-time methods. The coarse propagator solves only one time step per window $[T_{n-1},T_{n}],$ $n=1,\dots,N,$ thereby using step size $\Delta T=T/N.$

We now illustrate the performance of the MH approach for linear and nonlinear models of the coaxial cable and compare it with existing standard and parallel-in-time approaches. We discuss the performance firstly in terms of effective linear system solved. This number corresponds to the solver calls executed sequentially which, in case of the parallel computations, is the maximum number of linear system solutions performed by a single core. This is a convenient and easily reproducible measure because of its independence with respect to the software implementation. However, it disregards aspects like communication costs and thus also wallclock times are compared.

6.2 Linear model

Let us assume that the reluctivity function $\nu$ in (72) is linear, i.e., it does not depend on the solution $\vec{A}.$ We therefore deal with problem (18), where the matrix $\mathbf{K}$ is constant. One could then construct a linear (PP-PC) system of form (20) and directly transform it into the frequency domain (22) [20], thereby solving $N$ separate $d$-dimensional systems in parallel.

The computational costs of several solution approaches are compared for the linear model in terms of the effective number of required linear system solutions. Figure 5 on the left illustrates how the costs of the parallel algorithms applied to the linear model scale when the number of cores $N$ increases, while the precision of the fine solution remains unchanged (strong scaling). The classical time-stepping approach, when the solution of the initial value problem is performed sequentially starting from zero initial condition, required time integration over {\color[rgb]{0,0,0}k^{*}=42} periods (which corresponds to $84\,000$ time steps and solution of linear systems) until the error (51) of the obtained periodic steady-state solution reached the tolerances specified in (48). We consider this sequential time-stepping solution as a benchmark to compare performance of the other applied approaches.

The application of PP-IC (3)-(4) reduced the number of effective linear system solutions by a factor of $14$ compared to the sequential calculation when using $N=50$ cores and by a factor of $3$ with $N=5$ CPUs. The PP-PC method with the Jacobi-like iteration scheme (9) calculates the periodic solution effectively solving $15$ times less linear systems than the classical time stepping already when exploiting only $N=5$ processing units. However, the disadvantage of this approach is that the computational costs stagnate and there is no further acceleration when applying additional computational power.

In contrast to this, the MH approach applied to the PP-PC system (20) as in [20] delivers the steady-state solution effectively solving $500$ times less linear systems when using $N=50$ and $27$ times less linear systems when using $N=5$ cores, compared to the standard time integration. Finally, the direct MH solution of the TP system (44) is performed on the fine grid (with step size $\delta T=10^{-5}$ s on $[0,0.02]$ s). It requires $N_{\mathrm{f}}=2\,000$ separate linear system solves. Clearly, in the linear case the performance gain obtained by this approach scales perfectly linear as $N$ grows. In particular, when distributing the workload among $N=5$ CPUs the computational efforts of the steady-state calculations in terms of the effective linear system solves are reduced by a factor of $210$ and for $N=50$ by a factor of $2\,100.$ Therefore, as expected the TP solution with the MH approach [5] yields optimal scaling for the solution of the linear time-periodic eddy current problem.

6.3 Nonlinear model

We now present the performance of the iterative method introduced in Section 5 for parallel-in-time solution of a nonlinear time-periodic problem (72)-(74). In this case the nonlinear time-periodic system (30) is solved with the simplified Newton method (33), accompanied by the frequency domain solution at each Newton iteration. The initial guess in (36) for the simplified Newton iteration uses $\mathbf{Z}=\mathbf{0}$, i.e.,

[TABLE]

is chosen at PP-PC iteration $k+1$. Almost identical results have been obtained for the choice

[TABLE]

i.e., the average among the $N$ time steps at PP-PC iteration $k$. Since such choices are typically problem-specific, there is no further discussion of ‘sufficiently good’ initial guesses here. The interested reader is referred to, e.g., [8]

As in Section 6.2, Figure 5 (on the right) illustrates the computational efforts of all the previously described methods to obtain the periodic steady-state solution for the nonlinear coaxial cable model. Sequential time stepping starting from zero initial value reached the steady-state solution, periodic up to the tolerances given in (48) using the norm (51), after calculation over {\color[rgb]{0,0,0}k^{*}=31} periods, which corresponded to solution of $124\,000$ linear systems. PP-IC computed the periodic solution effectively solving $19$ times less linear systems when using $N=50$ cores and $5$ times less systems when $N=5.$ The reduction of the effective number of solved linear systems using PP-PC (9) amounts to a factor of $14$ with $N=5$ CPUs and to a factor of $19$ with $N=50,$ as within PP-IC. As previously seen, the costs of the PP-PC solution do not significantly decrease with the growth of $N$ and even increase slightly, when comparing the results for $N=20$ and $N=50.$

On the other hand, the simplified Newton iteration (33) for PP-PC (30) with MH correction on the coarse grid effectively solves about $16$ times less linear systems for $N=5$ and $420$ times less for $N=50.$ Now, in contrast to the linear case, the TP MH solution (with \mbox{\boldmathu}^{(0)}=\mathbf{0} in (45)) performs worse than PP-PC MH for almost all the considered values of $N$ except $N=5.$ More specifically, for $10$, $20$, and $50$ cores the number of effective linear solutions within TP MH is bigger than that of PP-PC MH, e.g., reduction factor is equal to $206$ (versus factor $420$ of PP-PC MH) when using $N=50.$ In case of $N=5,$ $20$ times less systems were solved effectively with TP MH which is a better result than that of PP-PC MH (where factor $16$ is observed). In Figure 6 on the left we see that the percentages of the corresponding calculations with respect to the sequential solution are considerably small (less than $1\%$ for $N=50$ and only $5$-$6\%$ for $N=5$).

The Figure 6 (right) illustrates the total number of linear systems solved within PP-PC MH and TP MH, i.e., the numbers of linear systems solved by each CPU are added. We see that TP MH always solves the same number of linear systems, since increasing the number of cores does not influence the convergence. The method converged in $15$ iterations, each calculating $2\,000$ steps, which in total gives $15\times 2\,000=30\,000$ systems solves. In contrast to this, the reduction of the computational costs within PP-PC MH comes not only from distributing the workload among CPUs but also from its faster convergence with bigger $N$ due to the increasing accuracy of the coarse solver. Therefore, even when the code is not truly parallelized, both approaches solve much less linear systems compared to the classical sequential solution ($4$ times less for TP MH and from $3$ up to $8$ times less systems for PP-PC MH).

Finally, we compare the actual computational time of the simplified Newton approach, applied to the PP-PC and the TP systems, followed by the MH correction. As before, the wall clock time measurements illustrate better performance of PP-PC MH in case of $N=10$, $20$ and $50$, and of TP MH when $N=5$ (see Figure 7 on the left). For $N=5$ PP-PC MH took $4.4$ hours, while TP MH required $3.5$ hours. In case of $N=50$ calculation with PP-PC MH lasted $14.7$ minutes, whereas TP MH took $23.5$ minutes. In comparison, the sequential time stepping over one period was executed in $2.15$ hours, which would extend to a simulation of about $66.7$ hours (almost $3$ days) calculating over $31$ periods until the steady state is reached. Therefore, for PP-PC MH the speed up factor is $15$ for $N=5$, $40$ for $N=10$, $108$ for $N=20$, and $271.5$ for $N=50.$ In its turn, TP MH accelerates the sequential time stepping by a factor of: $19$ for $N=5$, $40$ for $N=10$, $79$ for $N=20$, and $170$ for $N=50.$ The corresponding percentages with respect to the sequential execution time are presented in Table 1.

On the right in Figure 7 the time needed to calculate the Fourier transform together with its inverse is illustrated. We see that within TP MH the DFT occupies approximately a third part of the total computational time for all $N$ (from $33\%$ to $37\%$), while the DFT cost of PP-PC MH is extremely small (less than $1\%$ of the total execution time.) Besides, we see that for TP MH the time for DFT decreases with the increase of $N$, as expected due to the better parallelization capability. The opposite happens within PP-PC MH: although the cost of DFT still remains tiny it grows together with $N$, since the size of the DFT matrix $\mathbf{F}$ (23) increases. Finally, we would like to comment that the communicatioin costs of the both approaches are within the range of $14$ s and $48$ s for the considered values of $N.$

On the left in Figure 8 we show the convergence of the PP-PC MH approach for the strong scaling. As expected for a Parareal-based method, the more subintervals $N$ are used, the less iterations are required until convergence of PP-PC due to the increasing precision of the coarse solver. This is also represented on the right in Figure 8, where we illustrate the percentage of the principal components of PP-PC MH regarding the total calculation time. In particular, each iteration of PP-PC includes the parallel fine and coarse solutions and the simplified Newton iterations. As a result of the decreasing number of PP-PC iterations with a bigger $N$ the timings of each part of the algorithm decreases. The most time-consuming part is the parallelized fine solution, whose contribution amounts to $95\%$ of the total time when $N=5$ and to $69\%$ for $N=50$. The calculation timings of the coarse solution and of simplified Newton also decrease, while their parts of the total time are respectively $4.3\%$ and $0.7\%$ for $N=5$; and $28\%$ and $3\%$ when $N=50.$

We note that in both linear and nonlinear cases the errors (47), (49)-(51) were respectively calculated for the corresponding solution approaches. We have additionally performed comparison of the obtained results to the benchmark sequential solution using the norm (48) at each fine time step and calculating the maximum among the values. For all the considered methods the calculated solutions were close to the benchmark solution in the norm (48) up to the tolerances of $\text{rTol}=2.5\cdot 10^{-2}$ and $\text{aTol}=2.5\cdot 10^{-5}.$

7 Application to a three-dimensional transformer model

We now apply the proposed algorithm with the MH correction to a three-dimensional (3D) model of a transformer, whose cross-section and dimensions are illustrated in Figure 9 on the left. The computational domain $\Omega$ consists of a steel core and two copper coils, surrounded by air with prescribed homogeneous Dirichlet conditions on the outer boundaries. The total depth of the coils in the dimension $x_{2}$ is $10$ cm, while the width of the core is $3$ cm. The two coils have $358$ and $206$ windings wound around the core, respectively.

The nonlinear magnetic reluctivity $\nu$ is given in the ferromagnetic material of the core by the Brauer’s curve [6]

[TABLE]

with parameters $k_{1}=0.3774,$ $k_{2}=2.970,$ and $k_{3}=388.33.$ The reluctivities in the air and the copper regions are given by the reluctivity of vacuum $\nu_{0}.$

The electric conductivity $\sigma(\vec{x})$ is nonzero only in the steel part and is equal to $5\cdot 10^{5}\text{\,}\mathrm{S}\text{\,}{\mathrm{m}}^{-1}$. The coils provide the sinusoidal current excitation

[TABLE]

where $\vec{\chi}_{i}(\vec{x}),$ $i=1,2$ denote the winding functions of the two coils [28] and $T=0.02$ s.

Using these input data we consider the eddy current problem of the form (72)-(74) here. The spatial discretization is performed using the finite integration technique (FIT) [32] with $d=48\,417$ degrees of freedom, which leads to a time-periodic problem of the form (29). Although the system size is considerably bigger than that of a coaxial cable model from Section 6, we note that the considered grid illustrated in Figure 9 on the right is visually relatively coarse.

The time integration is achieved with the implicit Euler method. The nonlinear system of equations at each time step is solved using the successive substitution method, i.e., for fixed $n\geq 1$ and for $s=0,1,\dots$ the solution \mbox{\boldmathu}_{n}^{(s+1)} of the linear system

[TABLE]

is calculated. The fine step size $\delta T=10^{-5}$ s and the coarse step size $\Delta T=10^{-3}$ s are used within the (parallel-in-time) time-domain simulations. Parallelization among $N=20$ CPUs is exploited within the considered parallelization approaches. The solutions are calculated up to the tolerances specified in (48). The resulting magnetic flux density distribution at $t=0.02$ s is depicted in the Figure 9 on the right. Figure 10 (left) illustrates the induced time-periodic voltages in the two coils, calculated with the time step $\delta T$ over the period $[0,0.02]$ s.

Linearization of the periodic system (31) is performed using the fixed point iteration (43) with constant matrix

[TABLE]

where $\mathbf{C}={1}/{\Delta T}\cdot\mathbf{M}$ and $\bar{\mathbf{K}}$ is obtained for the fixed reluctivity in the core $\bar{\nu}=$\nu_{0}\cdot 10^{-3}\text{,}\mathrm{m}\text{,}{\mathrm{H}}^{-1}$$. An analogous linearization is applied also to the TP system (44), using the matrix $\mathbf{H}$ from (82) with $\mathbf{C}={1}/{\delta T}\cdot\mathbf{M}.$ We now describe the performance of PP-PC MH and TP MH in terms of the effective and total numbers of the solved linear systems, as well as in terms of the wall clock time and compare it to the sequential time stepping.

The standard sequential time integration with step size $\delta T=10^{-5}$ s required solution over $k^{*}=21$ periods until the steady state was reached in terms of the error (51). This corresponds to solution of $84\,000$ linear algebraic systems. It its turn, within the linearization of the form (43) applied to nonlinear TP MH $600$ effective linear solutions were calculated during $6$ iterations. It therefore delivered the steady-state solution solving $140$ times less effective linear systems than the classical sequential approach. In turn, PP-PC MH converged in two iterations, each requiring $6$ inner fixed point iterations (43). It effectively needed $424$ linear systems solutions, which is $198$ less than with the sequential time stepping. The total number of solved linear systems is equal to $6\times 2\,000=12\,000$ for TP MH and $8\,468$ for PP-PC MH. The corresponding percentages of the number of total and effective system solutions calculated by PP-PC MH and TP MH with regard to the sequential time stepping are illustrated Figure 10 on the right.

Finally, measurement of the wall clock time showed that the TP MH calculation lasted slightly more than $9$ hours, while PP-PC MH only needed $2.6$ hours. The sequential computation over one period took almost $22$ hours, which over $21$ periods would extend to about $19$ days. The estimated speed up factors of TP MH and PP-PC MH are therefore about $51$ and $176$, respectively. We point out the considerable difference in the factors when comparing the effective number of solved linear systems and when comparing the actual execution time of the TP MH ($140$ vs. $51$) to the sequential calculation. This is due to the fact that the solution of the linear systems took about $60\%$ of the total computational time (while the DFT required $30\%$). In contrast to this, the measure based on the number of solved effective linear systems worked much better for estimating the computational effort of PP-PC MH, since the linear system solutions took the major part of the execution time.

8 Conclusions

This paper presents an iterative parallel-in-time solution approach for nonlinear time-periodic problems combined with the MH coarse grid correction. The transformation into frequency domain introduces an additional parallelizability on the coarse grid, since the frequency components become decoupled. The proposed simplified Newton method with a special choice of the initial approximation keeps the nonlinearity constant over the Newton iterations and among the time instants, thereby letting the MH frequency domain approach be applicable at each iteration. Convergence of the iterative algorithm has been analyzed in detail for a 1D model problem and derived from the assumptions imposed on the nonlinearity. Performance of the method is illustrated for the time-periodic eddy current problem for both linear and nonlinear 2D models of a coaxial cable. Superiority of the MH solver over several existing (parallel-in-time) time-domain approaches is shown based on the computational costs calculated in terms of the number of required effective and total linear system solutions, as well as when measuring the wall clock time. Numerical results demonstrate that the frequency domain transformation applied directly to the fine time-periodic system (within TP MH) deliver optimal scaling for the linear problem, while the MH correction performed at the simplified Newton iteration for the coarse PP-PC system performed better than TP MH in the nonlinear case for the majority of the considered settings (except $N=5$ CPUs). Finally, the considered solution approaches with MH correction were applied to the steady-state analysis of a 3D transformer model. It was observed that compared to the classical sequential time integration TP MH had a speedup up to factor $51$ when using $N=20$ cores. A better result was obtained with the PP-PC MH algorithm, which was able to calculate the steady-state solution using $176$ times less computational time than the sequential approach.

Acknowledgements

The authors would like to thank Herbert De Gersem for the multiple fruitful discussions on the multi-harmonic solution approach.

Funding: This research was supported by the Excellence Initiative of the German Federal and State Governments and the Graduate School of Computational Engineering at Technische Universität Darmstadt, as well as by DFG grant SCHO1562/1-2 and BMBF grant 05M2018RDA (PASIROM).