Parareal Algorithm Implementation and Simulation in Julia

TL;DR

This paper implements the parareal algorithm in Julia, demonstrating its flexibility with different integrators and providing a simulation tool for performance analysis and educational purposes.

Contribution

The paper provides a complete Julia implementation of the parareal algorithm with customizable integrators and a simulation environment for testing and educational use.

Findings

Implementation supports various integrators like Euler and Runge-Kutta.

Simulation tool aids in understanding and analyzing algorithm performance.

Demonstrates the algorithm's applicability to general first-order initial-value problems.

Abstract

We present a full implementation of the parareal algorithm---an integration technique to solve differential equations in parallel---in the Julia programming language for a fully general, first-order, initial-value problem. We provide a brief overview of Julia---a concurrent programming language for scientific computing. Our implementation of the parareal algorithm accepts both coarse and fine integrators as functional arguments. We use Euler's method and another Runge-Kutta integration technique as the integrators in our experiments. We also present a simulation of the algorithm for purposes of pedagogy and as a tool for investigating the performance of the algorithm.

Figures

Equations

$$y^{\prime }(t)=f(t,y(t))t\in [a,b]$$

$$y(t^{\ast })=y^{\ast }{t}^{\ast }\in [a,b].$$

$$I(\delta ,x_0,y_0,g)=y_0+g(x_0,y_0)\delta$$

$$y_{0,1}$$

$$y_{n+1,1}$$

$$z_{n,1}^0$$

$$z_{n,1}^{m+1}$$

$$y(t_n)=y_{n,1}.$$

$$w_{1,2}$$

$$w_{n,2}$$

$$y_{n,2}$$

$$w_{1,k}$$

$$w_{n,k}$$

$$y_{n,k}$$

$$y^{\prime }(x)=sin(xy),x\in [-20,20]$$

Full text

Concurrent programming, Euler’s method, Julia, Runge-Kutta methods, parareal algorithm, ordinary differential equations.

1. Introduction

The parareal algorithm was first proposed in 2001 by Lions, Maday, and Turinici (Lions et al., 2001) as an integration technique to solve differential equations in parallel. We present a full implementation of the parareal algorithm in the Julia programming language (https://julialang.org) (Moffit and Tate, 2014) for a fully general, first-order, initial-value problem. Furthermore, we present a simulation of the algorithm for purposes of pedagogy and as a tool for investigating the performance of the algorithm. Our implementation accepts both coarse and fine integrators as functional arguments. We use Euler’s method and another Runge-Kutta integration technique as the integrators in our experiments. We start with a brief introduction to the Julia programming language.

2. An Introduction to Julia:

Dynamic, yet Efficient, Scientific/Numerical Programming

Julia is a multi-paradigm language designed for scientific computing; it supports multidimensional arrays, concurrency, and metaprogramming. Due to both Julia’s LLVM-based Just-In-Time compiler and the language design, Julia programs run computationally efficient—approaching and sometimes matching the speed of languages like C. See (web, ) for a graph depicting the relative performance of Julia compared to other common languages for scientific computing on a set of micro-benchmarks.

2.1. Coroutines and Channels in Julia

Coroutines are typically referred to as tasks in Julia, and are not scheduled to run on separate CPU cores. Channels in Julia can be either synchronous or asynchronous, and can be typed. However, if no type is specified in the definition of a channel, then values of any type can be written to that channel, much like unix pipes. Messages are passed between coroutines through channels with the put! and take!() functions. To add tasks to be automatically scheduled, use the schedule() function, or the @schedule and @sync macros. Of course, coroutines have little overhead, but will always run on the same cpu.

The current version of Julia multiplexes all tasks onto a single os thread. Thus, while tasks involving i/o operations benefit from parallel execution, compute bound tasks are effectively executed sequentially on a single os thread. Future versions of Julia may support scheduling of tasks on multiple threads, in which case compute bound tasks will see benefits of parallel execution too(Jul, ).

2.2. Parallel Computing

In addition to tasks, Julia supports parallel computing—functions running on multiple cpus or distributed computers. New processes are spawned with addproc($\mathtt{<}$n$\mathtt{>}$), where $\mathtt{<}$n$\mathtt{>}$ is the number of processes desired. The function addproc returns the pids of the created processes. The function workers returns a list of the processes. Alternatively, the Julia interpreter can be started with the -p $\mathtt{<}$n$\mathtt{>}$ option, where $\mathtt{<}$n$\mathtt{>}$ is the number of processes desired. For instance:

⬇

$ julia

julia> addprocs(3)

3-element Array{Int64,1}:

2

3

4

julia> workers()

3-element Array{Int64,1}:

2

3

4

^D

$

$ julia -p 3

julia> workers()

3-element Array{Int64,1}:

2

3

4

^D

$

Note that the process ids start at 2 because the Julia REPL shell is process 1.

Processes in Julia, which are either locally running or remotely distributed, communicate with each other through message passing.

The function remotecall($\mathtt{<}$Function$\mathtt{>}$, $\mathtt{<}$ProcessID$\mathtt{>}$, $\mathtt{<}$args …$\mathtt{>}$) executes $\mathtt{<}$Function$\mathtt{>}$ on worker $\mathtt{<}$ProcessID$\mathtt{>}$ and returns a value of the Future type, which contains a reference to a location from which the return value can be retrieved, once $\mathtt{<}$Function$\mathtt{>}$ has completed its execution. The Future value can be extracted with the function fetch(), which blocks until the result is available. Thus, the function remotecall is used to send a message while the function fetch is used to receive a message. For instance:

⬇

julia> addprocs(2)

julia> future = remotecall(sqrt, 2, 4)

Future(2,1,3,Nullable{Any}())

julia> fetch(future)

2.0

After the function remotecall is run, the worker process simply waits for the next call to remotecall.

⬇

julia> counter1 = new_counter(3)

(::#1) (generic function with 1 method)

julia> future = remotecall(counter1, 2)

Future(2,1,23,Nullable{Any}())

julia> fetch(future)

4

The Julia macro @spawn simplifies this message-passing protocol for the programmer and obviates the need for explicit use of the low-level remotecall function. Similarly, the macro @parallel can be used to run each iteration of a (for) loop in its own process.

⬇

julia> future = @spawn sqrt(4)

julia> fetch(future)

2.0

julia> addprocs(2)

2-element Array{Int64,1}:

3

4

julia> @everywhere function fib(n)

**if** (n < 2)

    **return** n

**else**

    **return** fib(n-1) + fib(n-2)

**end**

end

julia> @everywhere function fib_parallel(n)

**if** (n < 35)

    **return** fib(n)

**else**

    x = @spawn fib_parallel(n-1)

    y = fib_parallel(n-2)

    **return** fetch(x) + y

**end**

end

julia> @time fib(42)

2.271563 seconds (793 allocations: 40.718 KB)

267914296

julia> @time fib_parallel(42)

3.483601 seconds (344.48 k allocations:

                15.344 MB, 0.25*% gc time)*

There are also remote channels which are writable for more control over synchronizing processes.

2.3. Multidimensional Arrays

Julia supports multidimensional arrays, an important data structure in scientific computing applications, with a simple syntax and their efficient creation and interpretation over many dimensions (Bezanson et al., 2014). The function call ArrayType($\mathtt{<}$dimensions$\mathtt{>}$) creates an array, where the $n^{\mathrm{th}}$ argument in $\mathtt{<}$dimensions$\mathtt{>}$ specifies the size of the $n^{\mathrm{th}}$ dimension of the array. Similarly, the programmer manipulates these arrays using function calls that support infinite-dimensional arrays given only limitations on computational time.

In summary, Julia incorporates concepts and mechanisms—particularly concurrency and multidimensional arrays—which support efficient scientific computing.

3. The Parareal Algorithm

The parareal algorithm is designed to perform parallel-in-time integration for a first-order initial-value problem. The algorithm involves two integration techniques, often known as the ‘coarse’ integrator and the ‘fine’ integrator. For the algorithm to be effective, the coarse integrator must be of substantially lower computational cost than the fine integrator. The reason will become apparent later in this section. Consider the differential equation $\eqref{ode}$ given by

[TABLE]

with its associated initial-value problem $\eqref{ivp}$

[TABLE]

For simplicity, let us assume $t^{*}=a$, so that the solution only extends rightward. To obtain an approximate solution to equation $\eqref{ode}$ satisfying the initial condition $\eqref{ivp}$, we partition our domain into $[t_{0}=a,...,t_{N}=b]$ with uniform step size $\Delta$. We now precisely define an ‘integrator’ as a function from $(0,\infty)\times\mathbb{R}^{2}\times\mathscr{R}$ to $\mathbb{R}$ where $\mathscr{R}$ is the set of all Riemann integrable functions. For example, the integrator $I$ given by

[TABLE]

is the integrator corresponding to Euler’s method with step size $\delta$. Let $\mathscr{C}$ and $\mathscr{F}$ be the coarse and fine integrators, respectively. Define

[TABLE]

Since $y_{n+1,1}$ depends on $y_{n,1}$, this algorithm is inherently sequential. Partition $[t_{n},t_{n+1}]$ into $\{t_{n}^{0}=t_{n},...,t_{n}^{m},...t_{n}^{M}=t_{n+1}\}$ with uniform step size $\delta<\Delta$. Define

[TABLE]

This yields an approximate solution $\{z_{n,1}^{0},...,z_{n,1}^{M}\}$ to $\eqref{ode}$ over $[t_{n},t_{n+1}]$ with initial conditions

[TABLE]

Since $z_{n_{1},1}^{m_{1}}$ does not depend on $z_{n_{2},1}^{m_{2}}$ for $n_{1}\neq n_{2}$, we can compute these approximations in parallel. After the last subproblem is solved, we simply combine the solutions on each subdomain to obtain a solution over the whole interval. However, our values $\{y_{1,1},...,y_{n,1}\}$ are relatively inaccurate. The vertical spikes in the orange graph separating the coarse and fine predictions in Figure 1 illustrate this error. However, $z_{n-1,1}^{M}$ is a better approximation for $\phi(t_{n})$ where $\phi$ is the exact solution to the differential equation. We use this to obtain a better set of points $\{y_{n,2}\}$ for the coarse approximation. We do this by first defining $w_{n,1}=y_{n,1}$ and then defining

[TABLE]

Thus, $w_{n+1,2}$ serves as a new prediction given a more accurate previous prediction from $y_{n,2}$ since $z_{n-1,1}^{M}$ has now been taken into account in calculating $y_{n,2}$. In general, we continue evaluating so that for $k>1$, we have

[TABLE]

Note that since $y_{n,k}$ is dependent on $w_{n,k}$, this step must be done sequentially. As $k$ increases, $w_{n,k}-w_{n,k-1}\to 0$, which means that $y_{n,k}$ converges to the value that the fine integrator would predict if fine integration were simply done sequentially. Thus, each $k$ denotes fine integration over the whole interval. This means that the total computation performed is much greater than if fine integration were performed sequentially. However, the time efficiency of each iteration has the potential to be improved through concurrency. Since fine integration is more computationally intensive, this improvement in the run-time efficiency may compensate for the cumulative computation performed.

Let $K$ be the total number of iterations necessary to achieve a desired accuracy of solution and $P$ be the number of subintervals into which we divide according to the coarse integrator. If $K=1$, then we achieve perfect parallel efficiency. If $K=P$, then we likely slowed the computation down. The parareal algorithm is guaranteed to converge to the solution given by the sequential fine integrator within $P$ iterations. For a more complete treatment of this convergence analysis, we refer the reader to (Gander and Vandewalle, 2007). For fully general pseudocode, we refer the reader to (Aubanel, 2011; Nielsen, 2012).

4. Parareal Algorithm Implementation in Julia

Listing 1 presents an implementation of the parareal algorithm (from the prior section) in Julia. The @async macro within the loop causes the program to evaluate the first expression to its right as a concurrent task (i.e., the for loop assigning values to sub). The @sync macro causes the main program thread to wait until all tasks (spawned in the the first expression to its right with an @async or @parallel macro) complete. Once all concurrent tasks are complete, execution of the program proceeds sequentially. Given the semantics of these macros, the program in Listing 1 correctly perform concurrent integration. The sequential and parallel versions of this implementation have no significant differences in run-time efficiency. However, if a sleep statement is placed in the argument of fineIntegrator, the parallel version runs much faster. This demonstrates that use of those two macros does lead to concurrent program execution.

5. Graphical Algorithm Simulation

The function simulate in Listing 2 creates a graphical simulator of the parareal algorithm. This function can be used to introduce the parareal algorithm to students in a numerical analysis course. The first line gets the sequential solution from the fine integrator (the ‘ideal’ solution) and the second line gets the history of the computations that took place during the parareal execution. The main loop over the variable $k$ then displays the inner workings of the algorithm. The ideal solution is plotted, with a scatter plot of the points obtained from the coarse integrator. To simulate the parallel nature of the algorithm, random progress is made on randomly selected subdomains. Thus, the plot dynamically makes partial progress on different subdomains until all subdomains are finished with the fine integration. After this, the plots are connected into the current iteration’s approximation. During the next iteration, the previous guesses from the coarse integrator are displayed in red and the new guesses from the coarse integrator are displayed in green. As $k$ increases, these guesses converge to the ideal solution.

In addition to the use of this function for pedagogical purposes, it can be used to investigate the types of curves for which the parareal algorithm might be practical. For instance, consider the differential equation

[TABLE]

with $y(-20)=10$, $\Delta=4$ (10 points), and $\delta=0.008$ (500 points). Figure 2 shows the first and ninth iterations. The ninth iteration’s large error on the right end of the interval shows that this is an example where parareal convergence is slow. This is as inefficient as possible, needing as many iterations as subdomains in order for the solution to converge. However, the simulation also shows that if $f(x,y)=sin(x)e^{x}$, then the solution converges after just one iteration. These two examples show that the algorithm’s efficiency can be highly dependent on the integrand. Below the simulation function are Euler’s method and another Runge-Kutta integration technique that can be used as examples to be passed as first-class functions as coarse or fine integration techniques to the ‘parareal’ or ‘simulate’ functions. A Git repository of both the implementation and graphical simulation is available in BitBucket at https://bitbucket.org/sperugin/parareal-implementation-and-simulation-in-julia. All of the graphical plots are generated with the Julia Plots package available at https://juliaplots.github.io/. A video describing this application of Julia is available on YouTube at https://www.youtube.com/watch?v=MtgbeLO6ZM4.