Parallel-in-time integration of Kinematic Dynamos

TL;DR

This paper explores the application of parallel-in-time algorithms, specifically Parareal, to accelerate numerical simulations of kinematic dynamos, achieving significant speedups beyond traditional spatial parallelization methods.

Contribution

It demonstrates the feasibility and effectiveness of using the Parareal parallel-in-time method for dynamo simulations, a novel approach in this context.

Findings

Speed ups beyond spatial parallelization were achieved.

Parareal efficiency was close to 0.3 for Galloway-Proctor flow.

Parallel in space and time speed ups of ~300 with 1600 cores.

Abstract

The precise mechanisms responsible for the natural dynamos in the Earth and Sun are still not fully understood. Numerical simulations of natural dynamos are extremely computationally intensive, and are carried out in parameter regimes many orders of magnitude away from real conditions. Parallelization in space is a common strategy to speed up simulations on high performance computers, but eventually hits a scaling limit. Additional directions of parallelization are desirable to utilise the high number of processor cores now available. Parallel-in-time methods can deliver speed up in addition to that offered by spatial partitioning but have not yet been applied to dynamo simulations. This paper investigates the feasibility of using the parallel-in-time algorithm Parareal to speed up initial value problem simulations of the kinematic dynamo, using the open source Dedalus spectral solver. Both the time independent Roberts and time dependent Galloway-Proctor 2.5D dynamos are investigated over a range of magnetic Reynolds numbers. Speed ups beyond those possible from spatial parallelization are found in both cases. Results for the Galloway-Proctor flow are promising, with Parareal efficiency found to be close to 0.3. Roberts flow results are less efficient, but Parareal still shows some speed up over spatial parallelization alone. Parallel in space and time speed ups of $\sim300$ were found for 1600 cores for the Galloway-Proctor flow, with total parallel efficiency of $\sim0.16$.

Figures

Tables

Table 1. Table 1 : Parameters for simulations. R m subscript 𝑅 𝑚 R_{m} : magnetic Reynolds number, k z subscript 𝑘 𝑧 k_{z} : wave number in z co-ordinate, k x subscript 𝑘 𝑥 k_{x} : wave number in x co-ordinate, N F subscript 𝑁 𝐹 N_{F} : number of modes in fine propagator, N C subscript 𝑁 𝐶 N_{C} : number of modes in coarse propagator, δ ​ t 𝛿 𝑡 \delta t : time step for fine propagator, Δ ​ t Δ 𝑡 \Delta t : time step for coarse propagator, Growth rate indicates growth rate of the magnetic field.

Flow	$R_m$	$k_z$	$k_x$	$N_F$	$N_C$	$\delta t$	$\Delta t$	Growth
								Rate
Roberts	512	2.87		160	80	${10}^{-2}$	$2\times {10}^{-2}$	0.11
	4096	7.5		512	256	$2.5\times {10}^{-3}$	$5\times {10}^{-3}$	0.097
Galloway-	3		0.57	16	8	$5\times {10}^{-3}$	${10}^{-1}$	0.15
Proctor	300		0.57	128	64	${10}^{-3}$	$2\times {10}^{-2}$	0.3
	3000		0.57	512	256	${10}^{-4}$	$5\times {10}^{-3}$	0.3

Equations

$$\frac{\partial \mathbf{B}}{\partial t}=\nabla \times \left(\mathbf{u}\times \mathbf{B}\right)+\frac{1}{R_m}{\nabla }^2\mathbf{B},$$

$$R_m=\frac{UL}{\eta },$$

$$\mathbf{u}=\left(C\sin (z)+B\cos (y),A\sin (x)+C\cos (z),B\sin (y)+A\cos (x)\right),$$

$$\begin{array}{r}{\partial }_t{b}_x(x,y,t)=-b_y\sin (y)-\cos (y){\partial }_x{b}_x-\sin (x){\partial }_y{b}_x\\ -(\sin (y)+\cos (x))\mathrm{i}{k}_z{b}_x+\frac{1}{R_m}\left({\nabla }^2-k_z^2\right)b_x,\end{array}$$

$$\begin{array}{r}{\partial }_t{b}_y(x,y,t)=b_x\cos (x)-\cos (y){\partial }_x{b}_y-\sin (x){\partial }_y{b}_y\\ -(\sin (y)+\cos (x))\mathrm{i}{k}_z{b}_y+\frac{1}{R_m}\left({\nabla }^2-k_z^2\right)b_y,\end{array}$$

$$\begin{array}{r}\frac{\partial b_x}{\partial x}+\frac{\partial b_y}{\partial y}+\mathrm{i}{k}_z{b}_z=0.\end{array}$$

$$\begin{array}{r}\mathbf{u}=(C\sin (z+\sin \omega t)+B\cos (y+\cos \omega t),\\ C\cos (z+\sin \omega t),B\sin (y+\cos \omega t)),\end{array}$$

$${\partial }_t{b}_y(y,z,t)=b_z{\partial }_{z}v-u\mathrm{i}{k}_x{b}_y-v{\partial }_y{b}_y-w{\partial }_z{b}_y+\frac{1}{R_m}\left({\nabla }^2-{k_z}^2\right)b_y,$$

$${\partial }_t{b}_z(y,z,t)=b_y{\partial }_{y}w-u\mathrm{i}{k}_x{b}_z-v{\partial }_y{b}_z-w{\partial }_z{b}_z+\frac{1}{R_m}\left({\nabla }^2-{k_x}^2\right)b_z,$$

$$\mathrm{i}{k}_x{b}_x+\frac{\partial b_y}{\partial y}+\frac{\partial b_z}{\partial z}=0.$$

$$\frac{\partial U}{\partial t}=f\left(U(t),t\right),U(0)=U_0,0\le t\le T$$

$$U_{n+1}^{k+1}=\mathcal{G}(t_{n+1},t_n,U_n^{k+1})+\mathcal{F}(t_{n+1},t_n,U_n^k)-\mathcal{G}(t_{n+1},t_n,U_n^k).$$

$$\sigma =\frac{{\Vert U_{N_p}^k-U_{N_p}^{k-1}\Vert }_2}{{\Vert U_{N_p}^k\Vert }_2}.$$

$$s={\left[\left(1+\frac{k_{con}}{N_p}\right)\frac{R_c}{R_f}+\frac{k_{con}}{N_p}\right]}^{-1},$$

$$s\le \min \left\{\frac{N_p}{k_{con}},\frac{R_f}{R_c}\right\},$$

$$ϵ=\frac{s}{N_p}=\frac{\text{parallel speed up}}{\text{number of processors}},$$

Full text

Speed ups beyond those possible from spatial parallelisation are found in both cases. Results for the Galloway-Proctor flow are promising, with Parareal efficiency found to be close to 0.3. Roberts flow results are less efficient, but Parareal still shows some speed up over spatial parallelisation alone.

Parallel in space and time speed ups of $\sim 300$ were found for 1600 cores for the Galloway-Proctor flow, with total parallel efficiency of $\sim 0.16$.

††journal: Elsevier

1 Introduction

Dynamo theory seeks to explain the processes that generate magnetic fields in stars, planets, and galaxies. These magnetic fields are sustained by currents generated by the flow of a conducting fluid [1, 2], such as molten iron in the Earth’s core, or plasma in the Sun and stars. The fluid velocity, $\boldsymbol{u}$, twists, stretches, and shears field lines, counteracting the Ohmic diffusion of the magnetic field, $\boldsymbol{B}$, which would otherwise cause the field to decay [3].

The dynamo problem is typically modelled using the induction equation

[TABLE]

coupled with the the momentum equation for the fluid, which determines $\boldsymbol{u}$. The magnetic Reynolds number is defined by

[TABLE]

where $U$ is a characteristic speed, $L$ a characteristic length scale, and $\eta$ is the magnetic diffusivity [4]. A solution of this coupled system is called a dynamo if the magnetic field continues to be generated as time $t\rightarrow\infty$. The smallest scale of magnetic features created by the evolution of the induction equation (1) is the resistive scale which is proportional to ${R_{m}}^{-1/2}$ [3], so that the spatial resolution required to capture these structures increases as ${R_{m}}^{1/2}$. This leads to very high computational requirements that limit the parameter regimes that can be studied. Simulations are currently carried out at parameter regimes far removed from those found in natural dynamos of real stars and planets [5].

Owing to the complexity of the full dynamo problem much attention has focused on the simpler problem of kinematic dynamos where $\boldsymbol{u}$ is prescribed. As the induction equation is linear in $\boldsymbol{B}$, solutions to the kinematic dynamo problem are either exponentially growing or decaying with a well defined growth-rate for a statistically stationary $\boldsymbol{u}$. If $\boldsymbol{u}$ is steady, then the growth-rate can be found by solving an eigenvalue problem, whilst if the flow is periodic in time, the growth-rate can be found by solving a Floquet problem. If in addition the flow is independent of one co-ordinate, say the $z$-co-ordinate (2.5D flow), then one can seek monochromatic solutions of the form $\boldsymbol{B}(x,y,z,t)=\boldsymbol{b}(x,y,t)\mathrm{e}^{\mathrm{i}k_{z}z}$, so that $k_{z}$ becomes a parameter and the problem becomes 2D. In particular, there is much interest in the behaviour of dynamos at different $R_{m}$. Low $R_{m}$ systems are dominated by diffusion, while high $R_{m}$ systems are dominated by advection. Laboratory and engineering type flows tend to have $R_{m}<1$ [6], flows in the Earth’s core are characterised by $R_{m}\sim 10^{2}-10^{3}$ [7], whilst flows in the Sun have $R_{m}\sim 10^{6}-10^{10}$ [8].

Behaviour of dynamos at large $R_{m}$ is of particular interest astrophysically. If the growth-rate stays bounded away from zero as $R_{m}\rightarrow\infty$ then a dynamo is called fast. Otherwise it is called slow. It is a necessary (but not sufficient) condition that a fast dynamo have chaotic Lagrangian particle paths. In 2.5D flows, this can only be the case for unsteady flows. In fully 3D flows, the paths tend to be chaotic even without time dependence. No fast dynamos have yet been proven mathematically, but some have been found to act as fast dynamos in numerical simulations.

A particular choice of 2.5D cellular steady flow was considered by Roberts [9]. The Galloway-Proctor circularly polarised (CP) [10] flow is an extension of this flow to include time dependence. While the Roberts flow must act as a slow dynamo, as it is steady, the Galloway-Proctor dynamo is thought to be fast. The Roberts flow dynamo has been used to investigate the experimental dynamo at Karlsruhe [11], whilst the Galloway-Proctor flow was used to investigate the formation of large scale magnetic fields at high $R_{m}$ [12].

High accuracy solutions are needed for dynamo simulations and so the majority of studies use spectral methods to discretise in space [13], expanding in Fourier series, spherical harmonics, or Chebyshev polynomials as best fits the domain, to make use of their spectral convergence properties. Because of the high resolution requirements, very long compute times are found for the full dynamo system, even when simulations are run on large numbers of cores in high performance computing (HPC) facilities. Matsui et al. [14], for example, tested scaling on up to 16,384 cores. Schaeffer et al. [15] ran simulations that required over 10 million cpu hours. Transforming between spectral and spatial coordinates, usually via Fast Fourier Transform (FFT), acts as a limiting factor on the parallel scalability of pseudospectral codes [16]. Owing to the global communications requirements of spectral methods, simulations are currently unable to scale to the huge number of processors available in modern HPC facilities. Further ways to parallelize computations are therefore required to investigate more realistic parameter regimes.

Parallel in time methods can increase scalability of computer simulations beyond saturation of spatial parallelisation [17] and have been investigated for over 50 years [18]. The most widely studied parallel-in-time algorithm is Parareal [19], which was introduced in 2001 and has spurred a renewed interest in the subject. Further parallel in time methods have been proposed recently, like PFASST [20], PITA [21], and Paraexp [22]. The Parareal algorithm has been utilised in fluid flow problems [23, 17], plasma physics [24], financial simulations [25], and in planetary mantle simulations [26]. The method has been extensively analysed mathematically by Gander and Vandewalle [27], and is thought to perform badly for purely hyperbolic and highly advective systems [28], as it primarily corrects amplitude defects, rather than phase or frequency defects [29]. While some recent works demonstrate efficient application of Parareal to hyperbolic PDEs [30], the applicability of the method for dynamo simulations has not yet been studied.

The present paper studies the ability of Parareal to speed up kinematic dynamo simulations. It presents the first demonstration that parallel-in-time methods can deliver speedup for the induction equation beyond the saturation point of spatial parallelization. We introduce an implementation of the Parareal algorithm in the open source spectral solver Dedalus [31]. and investigate its performance for the stationary Roberts [9] as well as the time-dependent Galloway-Proctor [10] flow. Although these flows are relatively simple, they generate complex dynamics in the magnetic field and are good test problems to demonstrate Parareal’s efficiency for dynamo simulations.

2 The kinematic dynamo problem

In this work, the kinematic dynamo is studied for a subset of the ABC (Arnold, Beltrami and Childress) class of flows

[TABLE]

[32], with one of $A$, $B$, or $C$ equal to 0, making the flow 2.5D.

2.1 Roberts Flow

The Roberts flow is found from the ABC flow by setting $A=B=1$, $C=0$, giving $\boldsymbol{u}=\left(\cos(y),\sin(x),\sin(y)+\cos(x)\right)$. Noting that $\boldsymbol{\nabla}\cdot\boldsymbol{B}=\boldsymbol{\nabla}\cdot\boldsymbol{u}=0$, equation (1) becomes

[TABLE]

This flow is periodic in $x$ and $y$, and equations (4) are solved in a $2\pi$ square plane with periodic boundary conditions.

2.2 Galloway-Proctor Flow

The Galloway-Proctor circularly polarised (CP) flow adds time dependence to a Roberts like flow

[TABLE]

with $\omega=1$ and $C=B=\sqrt{3/2}$, $A=0$. We look for solutions of the form $\boldsymbol{B}=\boldsymbol{b}(y,z,t)\mathrm{e}^{\mathrm{i}k_{x}x}$. When put into (1), we have

[TABLE]

2.3 Morphology of the Roberts and Galloway-Proctor fields

Figure 1 shows the contours of the y-component of the magnetic field for $R_{m}=3$ and $3000$ in the Galloway-Proctor flow, and $R_{m}=4$ and $4096$ in the Roberts flow. Larger scale and more diffuse structures are present in the $R_{m}=3$ and $R_{m}=4$ simulations. Finer structures emerge in the $R_{m}=3000$ and $R_{m}=4096$ cases, showing the effect of the ${R_{m}}^{-1/2}$ scaling on the smallest structures. The Galloway-Proctor magnetic field morphology changes over time, with the periodic change in the velocity field, which leads to the more complicated pattern shown in Figure 1(d). For the Roberts flow, the magnetic field has a constant morphology and simply grows exponentially in magnitude.

3 Parareal Algorithm

We give only a brief description of the Parareal algorithm, for a more detailed description see for example the works by Lions et al. or Gander and Vandewalle [19, 27]. Parareal is an iterative method for solving initial value problems (IVPs) of the form

[TABLE]

where the right hand side $f$ typically comes from spatially discretizing a partial differential equation. Parareal uses two different time stepping methods, the coarse method $\mathcal{G}$ and the fine method $\mathcal{F}$. The time domain is split up into $N_{p}$ time slices, defined by the time-points $0=t_{0}<t_{1}<...<t_{N_{p}-1}<t_{N_{p}}=T$, where $N_{p}$ is the number of processors available for time parallelisation. We assume that all time slice are of equal length $\Delta T=T/N_{p}$. The fine time stepping method has time step $\delta t$, while the coarse method has time step $\Delta t$.

Parareal starts with running the coarse solver from $t=0$ to $t=T$, giving an initial approximation to the solution $U$ at time $t_{n}$, $U^{k=0}_{n}$, for every $n=0,1,...,N_{p}$, where $k$ denotes the Parareal iteration number, and $n$ denotes the time slice. Then, the approximation is refined using the iteration

[TABLE]

Computing the fine method can be parallelized across time slices. However, the correction from the coarse propagator has to be computed in serial, propagating the Parareal correction throughout the time domain.

3.1 Convergence and performance

As $k\rightarrow N_{p}$, Parareal converges to the solution that would have been obtained through the use of the fine method in serial over the time domain $[0,T]$ [27]. Achieving speed up, however, is dependent upon converging to an acceptable tolerance ${tol}$ within a much smaller number of iterations than $N_{p}$. We test for convergence by measuring the relative two-norm defect at the last time slice between successive iterations

[TABLE]

We denote as $k_{con}$ the number of iterations required until $\sigma\leq{tol}$.

Speed up $s$ of Parareal can be estimated as

[TABLE]

where $R_{c}$ is the runtime of the coarse method over one time slice and $R_{f}$ is the runtime of the fine method over one time slice. Speed up is bounded by

[TABLE]

[27]. The bound illustrates the trade off that needs to be optimised to gain optimal performance with Parareal. On the one hand, the ratio between the run times of the fine and coarse solvers should be large, so that the second bound is high. This can be achieved by using a very coarse and cheap $\mathcal{G}$. However, the method also needs to converge in a small number of iterations to ensure that the first bound is high. A less accurate coarse solver will typically require a larger number of iterations to converge. Finding a good compromise between these two competing factors is key to gaining high speed up with Parareal. Parallel efficiency is defined as

[TABLE]

where ideal scaling $s=N_{p}$ would give an efficiency of one. Parareal has an efficiency bound of $1/k_{con}$, which highlights the need to keep the number of iterations low. The need to test for convergence between two successive iterations of Parareal means there is an effective limit of $1/2$ for parallel efficiency when using Parareal.

3.2 Coarse solvers

A number of strategies have been proposed to design coarse solvers for Parareal. The most simple to implement is an increase in time step size for the coarse solver compared with the fine solver, whilst leaving the timestepping method and spatial discretization unchanged [17, 33]. A second strategy is to use an implicit time stepping method with large step size for the coarse solver, with explicit time stepping for the fine solver [34]. Further strategies to find a coarse solver include simplifying the physics so that a less complicated model can be used [35, 36, 37]. Another strategy to be considered is to use a reduced spatial resolution, along with a larger time step [38]. Using this strategy, the method of interpolation from coarse to fine grids has been found to be important to Parareal convergence [39]. In section 5.2 we explore different options in the context of the kinematic dynamo problem.

4 Implementation

This work was carried out using the Dedalus [40] open source, spectral solver. We use a collocation-based pseudo-spectral method. Spectral methods benefit from spectral convergence: for a sufficiently smooth solution, the error of the method is $\mathcal{O}[(1/N_{s})^{N_{s}}]$ for a discretisation with $N_{s}$ collocation points. This compares favourably with a standard finite difference discretisation, which has error $\mathcal{O}[(1/N)^{p}]$ for a $p^{\text{th}}$ order method. Because of the periodic domain, Fourier bases are used. Dedalus uses the FFTW library to perform fast parallel transforms between real and spectral space and parallelizes in space over $n-1$ dimensions of an $n$ dimensional domain using the mpi4py library [41].

Dedalus offers a number of different time stepping methods up to \nth4 order, including multi-step and Runge Kutta IMEX methods [42, 43, 44]. Because multi-step methods in Parareal require restarting in every time slice and every iteration, we focus on Runge-Kutta methods. For optimal serial performance, we rely on implicit/explicit (IMEX) Runge-Kutta methods [42]. Here, we integrate the terms $(\boldsymbol{u}\cdot\boldsymbol{\nabla})\boldsymbol{B}$ and $(\boldsymbol{B}\cdot\boldsymbol{\nabla})\boldsymbol{u}$ explicitly, whilst the diffusion term will be treated implicitly. Below, we compare the IMEX Runge-Kutta methods RK111 (1 implicit stage, 1 explicit stage, \nth1 order), RK222 (2 implicit stages, 2 explicit stages, \nth2 order), and RK443 (4 implicit stages, 4 explicit stages, \nth3 order) to determine the most efficient serial fine solver against which to compare Parareal.

Simulations are carried out in a periodic domain with side length $2\pi$. The spatial resolution is measured in terms of the number of spectral modes in $x$ and $y$, ($N_{x}$, and $N_{y}$ respectively) for the Roberts flow, or number of modes in $y$ and $z$, ($N_{y}$ and $N_{z}$ respectively) for the Galloway-Proctor flow. Since the domain is a square, we set $N_{x}=N_{y}=N$ for the Roberts flow and $N_{y}=N_{z}=N$ for the Galloway-Proctor flow in all examples, so that the total spatial problem size is $N^{2}$. The initial conditions were set to a random perturbation of magnitude $10^{-5}$ for $B_{x}$ and $B_{y}$ in the Roberts flow, and for $B_{y}$ and $B_{z}$ in the Galloway Proctor flow.

In the kinematic approach, the simulation time needs to be long enough that the largest growing mode can be found and a steady growth rate calculated. In this work, this time was found to be around 10-20 turnover time periods. However, in a fully dynamic simulation, much longer simulation times are required to reach a steady state in both the fluid flow and the magnetic field. For example, Smith et al. [45] run simulations until $T=450$ to achieve a statistically steady flow. As the aim of the work is to inform future studies of dynamic simulations, we choose a longer time interval of 50.

Parareal was implemented by splitting the MPI world communicator into a space communicator and a time communicator. The space communicator was utilised by Dedalus to parallelize in space, whilst the time communicator was used to communicate between time slices. Two instances of Dedalus were created on each time slice, one with the coarse resolution ($N_{C}$ spectral modes in each direction) and time step $\Delta t$, and one with the fine resolution ($N_{F}$ spectral modes in each direction) and time step $\delta t$. Interpolation from coarse to fine grids, and restriction from fine to coarse grids, were carried out using the Dedalus set_scales(ratio) method on the field objects, which allows efficient interpolation based on Fourier re-sampling.

4.1 Pseudo-code for the Parareal algorithm in Dedalus

5 Results

To ensure that reported speedups are a like-to-like comparison and meaningful, we need to make sure that the solution provided by Parareal is of the same accuracy as the fine solver used serially. Given the need for highly accurate results in dynamo simulations, we aim for an accuracy of $10^{-5}$. Furthermore, we need to compare Parareal against an efficient fine solver – achieving speedup with Parareal when a much more efficient alternative to the fine solver exists would not provide convincing evidence of Parareal’s usefulness. Therefore, we first find the optimal fine solver in Dedalus to deliver the required accuracy and then parallelise it with Parareal.

5.1 Fixing the fine solver

First, we need to fix the required spatial resolution for $\mathcal{F}$. Convergence in space was tested by running simulations at double the previous spatial resolution until the normalised $L^{2}$ difference between two solutions was $\sim 10^{-15}$. At this point, the error from the spatial discretisation is of the order of machine precision. We denote $U_{N}$ as the solution vector containing $b_{x}$ and $b_{y}$ in real space with resolution $N$. The smaller solution was interpolated onto the same grid size as the fine solution using spectral interpolation so that the difference could be computed. The most highly resolved solution, $U_{Nmax}$, was then used as a reference solution to compute relative error of each $U_{N}$. The results are shown in Figure 2, confirming the expected spectral convergence behaviour. For each $R_{m}$, we set $N_{F}$ to the smallest value that gives a solution with error smaller than $10^{-5}$ (indicated by the dashed red line). Because higher magnetic Reynolds numbers produce smaller scale features, they require better spatial resolution to match our error tolerance.

Next, we fix the time stepping method and time step. Creating a reference solution with a temporal error of the order of machine precision proved to be unfeasible, due to computational constraints, especially in the higher $R_{m}$ cases. Therefore, a result with error lower than $10^{-7}$ was used as a reference solution for setting $\delta t$. This is two orders of magnitude smaller than the desired result of $10^{-5}$ and should provide an accurate estimation of the error due to time-stepping. A comparison of the Runge-Kutta time steppers available in Dedalus is shown in Figure 3. Because RK443 reaches the requested tolerance of $10^{-5}$ with the smallest number of evaluations of the right hand side function, it is the most efficient choice. Similar results were found for other magnetic Reynolds numbers. We therefore use RK443 for the fine method $\mathcal{F}$ throughout this work.

Figure 4 shows the dependence of the solution error on time step size for a range of $R_{m}$ for the Roberts and Galloway-Proctor flows. All simulations were carried out with the optimal $N_{F}$ found from the spatial resolution study. We can see that smaller step sizes are required to meet a given level of accuracy for the Galloway-Proctor flow than for the Roberts flow. This is likely due to the Galloway-Proctor flow depending explicitly on time. In the case of the Roberts flow, a $\delta t$ small enough to satisfy the stability requirements for a given $N_{F}$ is sufficient to also satisfy the accuracy requirement of $10^{-5}$, except for the most simple case of $R_{m}=4$. This affects the performance of Parareal through the ratio of computational run times $R_{f}/R_{c}$ because we have to use essentially the same time step for both the coarse and fine method. In contrast, for the Galloway-Proctor simulations, satisfaction of the stability requirement did not guarantee accuracy within the required tolerance, and a reduced $\delta t$ must be used for the fine solver, leading to a better coarse-to-fine computation time ratio.

5.2 Fixing the coarse solver

The $N_{F}$ and $\delta t$ determined in Section 5.1 are used in the fine solver. We now discuss the different possibilities available for choosing a coarse solver. Using the same spatial resolution with $\Delta t>\delta t$ was not suitable for the Roberts flow, as $\delta t$ was the largest stable time step for a given resolution. It was also unsuitable for the Galloway-Proctor simulations, as the ratio of $\Delta t/\delta t$ was not large enough to give meaningful speedup. Use of a fully implicit coarse solver was rejected, as the spectral spatial discretisation meant that a dense matrix solve would be required at each time step. This large increase in computational complexity would reduce the difference in computations required between the fine and coarse solvers, leading to smaller speed ups. There was little scope to attempt to use reduced physics in this study, as we are already considering the simplest form of dynamo problem. However, this strategy may be useful in further work on a non-linear dynamo. Coarsening in both space and time was found to be the most promising strategy. As time step stability is linked to spatial resolution, reducing spatial resolution allows a larger time step to be taken, even where the fine solver is at the largest stable time step. This means that $N_{C}<N_{F}$ and $\Delta t>\delta t$, opening up the possibility of a large difference in computational complexity between the coarse and fine solvers. However, too aggressive coarsening will lead to a very inaccurate coarse solver and slow convergence. The most efficient amount of spatial coarsening was studied by carrying out Parareal simulations with a wide range of coarse method spatial resolution.

Simulations were carried out for the Roberts flow with $R_{m}=512$. This is moderately high, whilst allowing relatively modest compute resources to be used. The fine solver parameters were fixed, with $N_{F}=160$, and $\delta t=10^{-2}$, while the coarse step $\Delta t$ was set to the highest stable step for the given $N_{C}$. This was found by estimating the error at different time steps for each resolution, as shown in Figure 5(a). A similar study was carried out for the Galloway-Proctor flow (Figure 5(b)). The number of Parareal time slices $N_{P}$ was fixed at 10. Figure 6 shows that the peak speed up is acquired when $N_{C}=0.5N_{F}$. When $N_{C}<0.5N_{F}$, the speed up is reduced by the extra number of Parareal iterations required to converge, and when $N_{C}>0.5N_{F}$, the difference in computational complexity between the coarse and fine solvers is insufficient.

5.3 Scaling Results

Scaling tests were carried out for both the Roberts flow and the Galloway-Proctor flow. Simulations of the Roberts flow were carried out on the ARC 3 HPC facility at the University of Leeds, made up of Intel Xeon E5-2650v4 (Broadwell) CPU’s, with a total of 6,048 cores. Simulations of the Galloway-Proctor flow were carried out on the ARCHER HPC facility, made up of Intel Xeon E5-2697v2 (Ivy Bridge) CPU’s, with a total of 109,056 cores.

A range of $R_{m}$ were simulated, to see the effect on Parareal performance (see Table 1). Scaling performance was compared with pure spatial scaling of the Dedalus solver. Fully parallel in space and in time simulations were also carried out in order to show how Parareal can increase scalability beyond the saturation of spatial scaling. Validation was carried out by comparing computed growth rates with those found in the literature. Growth rates for the Roberts dynamo were found to be consistent with those found by Plunian and Radler [11], which were reported to 2 significant digits, with the peak growth rate for each $R_{m}$ occurring at the same $k_{z}$ wave number reported. Growth rates for the Galloway-Proctor simulations were consistent with those found by Charbonneau et al. [46], with the peak growth rate occurring at a $k_{x}$ of 0.57. The Galloway-Proctor flow had the correct behaviour in terms of growth rate for large $R_{m}$, with the growth rate staying positive, showing the expected fast-dynamo behaviour.

5.4 Roberts Flow

Scaling results for the Roberts flow are shown in Figure 7 for $R_{m}=512$ (upper figures) and $R_{m}=4096$ (lower figures). For both cases, both space parallel scaling and efficiency are superior to Parareal in the beginning. As expected, spatial scaling is better for the $R_{m}=4096$ case with higher spatial resolution, due to higher workload per processor. While Parareal alone is not competitive, in both cases a combined space-time parallelization generates slightly more speedup than a pure spatial parallelization. The theoretical maximum efficiency for Parareal is $1/3$, indicated by a horizontal dashed line, due to the simulation requiring three iterations to converge. However, because of the relatively expensive coarse solver, Parareal’s observed efficiency is mostly substantially lower. As the efficiency of the combined space-time parallelisation is the product of the parallel in space efficiency and the parallel in time efficiency, it is low for high numbers of processors because of the low efficiency of Parareal for the Roberts flow. Despite the larger overall speedup, with efficiencies below 0.1, space-time parallelization using Parareal may not be particularly attractive .

5.5 Galloway Proctor Flow

Results for $Rm$= 3, 300 and 3000 for the Galloway Proctor flow are shown in Figure 8. Performance of Parareal is much better than for the Roberts flow. Parareal speed up is competitive with spatial parallelisation at a relatively low number of processors, and the efficiency of Parareal stays close to $1/3$ over a range of $N_{p}$ and does not fall much as the number of processors increases. The poor performance of Dedalus in parallelising the $Rm$=3 case is attributable to the fact that there are only $16^{2}$ spectral modes. In the $R_{m}=3000$ case, the results show parallel in space and time results, with 32 processors in space. The results show that speed up above that of spatial parallelisation alone is possible, with an efficiency around 0.16. In all cases, Parareal has not yet reached saturation in its scaling performance, and has almost ideal scaling behaviour, except for a constant offset due to the bounds on Parareal scaling. In all of these cases, the number of iterations required to converge was three, so that the efficiency is bounded by $1/3$. This is shown in Fig. 9, where the efficiency of the method is tracked over the different $R_{m}$. Efficiency is close to the bound of $1/3$, and Parareal efficiency does not fall with increasing $R_{m}$. Pure Parareal efficiency was estimated in the $R_{m}=3000$ case by dividing by the efficiency of the spatial parallelisation found for that particular $N_{S}$ (32). In that case, total parallel efficiency is lower than for the other Galloway-Proctor cases due to the combination of the spatial and temporal parallelisation, and is approximately the product of the two, as expected. This reduction in overall efficiency is unavoidable, as parallel in space is more efficient than Parareal for lower numbers of processors, and Parareal only becomes competitive after spatial efficiency falls away. Also shown on this Figure are the efficiencies obtained at different $R_{m}$ for the Roberts flow, highlighting the difference in performance of the method for the two cases.

6 Conclusions

The Parareal algorithm has been found to offer parallel speed up for kinematic dynamo simulations beyond what can be achieved through spatial parallelisation alone. In the case of the simpler steady Roberts dynamo, the speed up is modest and parallel efficiencies are low. Here, owing to the steady nature of the imposed velocity, the difference in computational complexity of the coarse and fine methods was found to be too small for good performance of Parareal. The issue was that the time step size was not a limiting factor on the accuracy of the fine solver; as long as the time step was stable, it was within the accuracy required. Therefore, there was little room to use a coarser resolution for the coarse propagator in Parareal.

Performance for the time-dependent Galloway-Proctor flow was better and the efficiency stayed close to the theoretical limit over a wide range of magnetic Reynolds numbers, while scaling well to large numbers of processors. In this problem, since evolution of the magnetic field depends explicitly on the current time, the accuracy of the solution depends more on the size of the time step. This means that a time step in the coarse solver much larger than that of the fine step is possible, allowing for better Parareal performance. Fully coupled dynamic dynamo simulation is complicated, and has non-linear dependencies, and so the accuracy of the fine solver is expected to behave more like the Galloway Proctor flow. Therefore, the good performance of Parareal for the Galloway-Proctor flow suggests that good performance can be expected also for more complex dynamos. This paper therefore gives an interesting example of how studying a too simple problem can lead to an overly negative assessment of the performance of Parareal.

The parallel efficiency of the Parareal algorithm applied to the Galloway Proctor dynamo is close to the theoretical maximum of $1/k_{con}$. This means that the overheads due to communication are small, in comparison to the serial cost of the coarse method, pointing to an efficient implementation of the algorithm. Performance of the algorithm when applied to this problem does not appear to degrade with increased $R_{m}$, as can be seen in Figure 9. Performance has remained constant, with Parareal efficiency not much lower than $1/3$ for $R_{m}=3$, $R_{m}=300$ and $R_{m}=3000$. This is noteworthy since highly advective problems are thought to cause problems with Parareal convergence, but this has not yet been found in the the highly advective case with $R_{m}$ up to $\sim 10^{3}$.

The results in this paper show that parallel in time methods can speed up dynamo simulations and are therefore worthy of further study. Better parallelization would enable the study of dynamos at larger magnetic Reynolds numbers by reducing simulation times by harnessing the ever growing number of available cores in HPC facilities. Current and future work involves extending our analysis to non-linear dynamic dynamo systems and the interaction of magnetic fields with convection. The development of successful parallel-in-time methods for these problems could allow the integration of geo- and astro- dynamo simulations in parameter regimes closer to reality than have been hitherto possible.

Acknowledgments

A.T.C is supported by the Engineering and Physical Sciences Research Council (EPSRC) Centre for Doctoral Training in Fluid Dynamics (EP/L01615X/1). C.J.D is supported by a Natural Environment Research Council (NERC) Independent Research Fellowship (NE/L011328/1). D. R. thankfully acknowledges support from grants EPSRC EP/P02372X/1 ”A new algorithm to track fast ions in fusion reactors” and NERC NE/R008795/1 ”Parallel Paradigms for Numerical Weather Prediction”. S.M.T is supported by a Levehulme Fellowship and by funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement no. D5S-DLV-786780). This work used the ARCHER UK National Supercomputing Service (http://www.archer.ac.uk) as well as ARC2 and ARC3, part of the High Performance Computing facilities at the University of Leeds. Figures were produced using Matplotlib [47]. All codes used to create the results in this manuscript are publicly available on GitHub at [48].