A shared memory implementation of pipelined Parareal
Name\corref{cor1
%% Text of abstract The paper introduces an OpenMP implementation of pipelined Parareal and compares it to a standard MPI-based implementation. Both versions yield essentially identical runtimes, but, depending on the compiler, the OpenMP variant consumes about $7$\% less energy. However, its key advantage is a significantly smaller memory footprint. The higher implementation complexity, including manual control of locks, might make it difficult to use in legacy codes, though.
implementation.tex
\section{PararealF90: Pipelined Parareal in MPI and OpenMP}\label{sec:implementation}
The code used here for benchmarking is written in Fortran 90 and available under an open-source BSD license~\cite{PararealF90}.
It is special-purpose and tailored to solve a single benchmark problem, 3D Burgers' equation
\begin{equation}
	u_{t} + u \cdot \nabla u = \nu \Delta u
\end{equation}
on $[0,1]^{3} \subset \mathbb{R}^{3}$ with periodic boundary conditions.
While this a rather specific setup, finite difference stencils are a widely used motif in computational science: even though specifics and runtimes will be different if more complex problems are solved, the general concepts are used in complex application codes, e.g. by means of domain specific embedded languages (DSEL)~\cite{GysiEtAl2015} and the possibility to use Parareal with such a DSEL has been shown~\cite{ArteagaEtAl2015}.
Thus, conclusions in terms of performance of different implementations of Parareal can be generalised to some extent, in particular because the Parareal routines only operate on linear arrays and do not see any specifics of the underlying time steppers or spatial discretisation.
By exchanging the modules implementing the spatial discretisation, \textsc{PararealF90} could be used to solve and benchmark different equations using stencil-based discretisations.

The time integration module provides two different methods, a strong stability preserving Runge-Kutta method (RK3-SSP)~\cite{ShuOsher1989} and a first order forward Euler.
The module for the spatial discretisation provides a fifth order WENO finite difference discretisation~\cite{ShuOsher1989} and a simple first order upwind stencil for the advection term as well as a second and fourth order centred stencil for the diffusive term.
Two modules provide the different implementations of Parareal, either using MPI as described in Subsection~\ref{subsubsec:parareal_alone_mpi} or OpenMP as introduced in Subsection~\ref{subsubsec:parareal_alone_pipe_openmp}. 

For $\mathcal{G}_{\Delta t}$, the forward Euler with upwind and second order centred stencils is used, while $\mathcal{F}_{\delta t}$ uses the RK3-SSP integrator, a 5th order WENO for advection and a fourth order centred stencil for diffusion.
Note that both implementations of Parareal use the same modules to provide the coarse and fine integrator and spatial discretisation.
Since all three versions rely on the same implementation for the actual integrators and spatial discretisation, differences in performance should therefore solely emerge from the different Parareal routines wrapped around the compute routines.
To run coarse and fine serial reference simulations, driver functions are used that directly call the same time stepper routines that are used within Parareal with the same configurations.

The code also comes with a test harness: in particular, the tests guarantee that all three implementations of Parareal produce results that are identical up to a tolerance of $\varepsilon = 10^{-14}$ and thus essentially to round-off error.
To detect possible race conditions, the comparison test can be performed multiple times: up to 100 instances of the test were run and passed on both used architectures described in Subsection~\ref{subsec:hardware}.
The tests also use various randomised parameters (randomised in a small range to avoid too large problems with very long runtimes): diffusion parameter $\nu$, number of time slices and processors $P$, number of finite difference nodes $N$ (each direction has a different value in the tests) and number of fine and coarse steps per time slice.
Furthermore, both implementations of Parareal use three auxiliary buffers per time-slice: $q$ to store the fine value and for communication, $\delta q$ to store the difference $\mathcal{F}_{\delta t}(q) - \mathcal{G}_{\Delta t}(q)$ needed in the correction step and $q_c$ to store the coarse value from the previous iteration.
It seems that without introducing additional communication, three copies per time slice is the minimum storage required for Parareal.

Both used systems have nodes with two multi-core CPUs, see the description in Subsection~\ref{subsec:hardware}.
For the OpenMP version to be efficient, it has take care not only of thread safety but also take into account \emph{non-uniform memory access} (NUMA): fetching data from memory associated with the other CPU requires communication through the intra-node interconnect and incurs overhead.
Thus, in~\textsc{PararealF90}, solution buffers are extended by one dimension, e.g. $\texttt{Q(i,j,k,p)}$, where $i,j,k$ refer to the three spatial coordinates and $p$ to the thread number.
In the discussion below, the first three dimensions are omitted and the buffer corresponding to thread $p$ is indicated simply by $\texttt{Q}(p)$.
By using first touch initialisation, each thread ensures its respective part of the solution data is allocated locally.
This strategy provides thread safety, as long as each thread is only accessing its own part of the buffers, as well as data affinity.

\subsection{Parareal in MPI}\label{subsubsec:parareal_alone_mpi}
\begin{algorithm2e}[t!]
        \caption{Parareal using MPI}\label{alg:mpi}
        \SetKwInOut{Input}{input}
        \Input{Initial value $q_0$; number of iterations $K$}
        \SetCommentSty{textit}
        \DontPrintSemicolon       
        	$q \leftarrow q_0$\\
	$p = \texttt{MPI\_COMM\_RANK}()$\\
	$q \leftarrow \mathcal{G}_{\Delta t}(q, t_{p}, 0)$\label{alg:parareal_mpi_predict_start}\\
	$q_c \leftarrow\mathcal{G}_{\Delta t}(q, t_{p+1}, t_{p})$\label{alg:parareal_mpi_predict_end}\\
	\For{$k=1,K$}{
		$q \leftarrow \mathcal{F}_{\delta t}(q, t_{p+1}, t_{p})$\label{alg:parareal_mpi_fine_start}\\
		$\delta q \leftarrow q - q_c$\label{alg:parareal_mpi_fine_end}\\
		\If{Process not first}{\label{alg:parareal_mpi_correction_start}
			$\texttt{MPI\_RECV}(q, \texttt{source}=p-1)$\\
		}
		\Else{
			$q \leftarrow q_0$
		}
		$q_c \leftarrow \mathcal{G}_{\Delta t}(q, t_{p+1}, t_{p})$\\
		$q \leftarrow q_c + \delta q$\\
		\If{Process not last}{
			$\texttt{MPI\_SEND}(q, \texttt{target}=p+1)$\\
		}\label{alg:parareal_mpi_correction_end}
	}
\end{algorithm2e}
The implementation of Parareal with MPI is straightforward and has been illustrated before~\cite{ArteagaEtAl2015}.
However, it is repeated here for the sake of completeness and sketched in Algorithm~\ref{alg:mpi}.
\begin{itemize}
	\item \textbf{Prediction phase:} lines~\ref{alg:parareal_mpi_predict_start} and~\ref{alg:parareal_mpi_predict_end}.
	Every process generates its own coarse starting value $\vv{q}^0_p = \mathcal{G}_{\Delta t}(\vv{q}_0, t_{p}, 0)$ and computes $\mathcal{G}_{\Delta t}(\vv{q}^0_p, t_{p+1}, t_{p})$ for use in the correction.
Later processes have to perform more time steps and thus take longer to complete the prediction phase, leading to the pipelined execution model sketched in Figure~\ref{fig:pipelining}, since no global synchronisation points exist.
	\item \textbf{Fine integrator:} lines~\ref{alg:parareal_mpi_fine_start} and \ref{alg:parareal_mpi_fine_end}.
Each process  computes the fine value $\mathcal{F}_{\delta t}(\vv{q}^{k}_{p}, t_{p+1}, t_{p})$ in line~\ref{alg:parareal_mpi_fine_start} and, in order to free the buffer $q$ for reuse in the receive, computes the difference $\delta q := \mathcal{F}_{\delta t}( \vv{q}^{k}_{p}, t_{p+1}, t_{p}) - \mathcal{G}_{\Delta t}(\vv{q}^{k}_{p}, t_{p+1}, t_{p})$ between coarse and fine value in line~\ref{alg:parareal_mpi_fine_end} and stores it in $\delta q$.
	\item \textbf{Update phase:} lines \ref{alg:parareal_mpi_correction_start}--\ref{alg:parareal_mpi_correction_end}.
To receive the updated value from the previous time slice, every process but the first posts an \texttt{MPI\_RECV}.
The first process simply copies the initial value $\vv{q}_0$ into the buffer $q$ instead.
After the receive is completed, the coarse value $\mathcal{G}_{\Delta t}(\vv{q}^{k+1}_{p}, t_{p+1}, t_{p})$ of the new starting value is computed and the updated end value $\vv{q}^{k+1}_{p+1} = \mathcal{G}_{\Delta t}(\vv{q}^{k+1}_{p}, t_{p+1}, t_{p}) + \delta q$ is computed and send to the process handling the following time slice using \texttt{MPI\_SEND}.
\end{itemize}

Note how this implementation naturally features pipelining as sketched in Figure~\ref{fig:pipelining} (right).
In order to obtain the non-pipelined execution flow sketched in Figure~\ref{fig:pipelining} (left), \texttt{MPI\_BARRIER} directives would have to be added artificially after the prediction step, that is after line~\ref{alg:parareal_mpi_predict_end}, and before the correction, that is before line~\ref{alg:parareal_mpi_correction_start}.
Parareal without pipelining in MPI would therefore be more complex and deliver reduced performance, so that this strategy seems rather senseless.
Previous benchmarks did not show a significant difference between blocking and non-blocking communication~\cite{ArteagaEtAl2015} and for the sake of simplicity, blocking communication is used here.
%
%
%
\subsection{Parareal in OpenMp}\label{subsubsec:parareal_alone_pipe_openmp}
The implementation of Parareal with pipelining in OpenMP introduced here is sketched in Algorithm~\ref{alg:parareal_openmp_pipe}.
To implement pipelined Parareal with OpenMP, essentially the whole algorithm is enclosed in one parallel region: threads are spawned by the \texttt{OMP PARALLEL} directive in line~\ref{alg:openmp-pipe-parastart} and terminated by \texttt{OMP END PARALLEL} in line~\ref{alg:openmp-pipe-paraend}.
Because this version requires manual synchronisation, a number of OpenMP locks is created using \texttt{OMP\_INIT\_LOCK} (not shown), one for each thread.
During the fine integrator and update step, these locks are set and unset using \texttt{OMP\_SET\_LOCK} and \texttt{OMP\_UNSET\_LOCK} to protect buffers during writes and avoid race conditions.
\begin{algorithm2e}[t!]
%\tiny
        \caption{Parareal with pipelining using OpenMP}\label{alg:parareal_openmp_pipe}
        \SetKwInOut{Input}{input}
        \Input{Initial value $q_0$; number of iterations $K$}
        \SetCommentSty{textit}
        \DontPrintSemicolon       
	$q(0) \leftarrow q_0$\\
	$\texttt{OMP PARALLEL}$\label{alg:openmp-pipe-parastart}\\
	$N = \texttt{OMP\_GET\_MAX\_THREADS}()$\\
	$\texttt{OMP DO}$\label{alg:openmp-pipe-predict-start}\\
	\For{$p=0,N-1$}{
		$q(p) \leftarrow q(0)$\\
		\If{Thread not first}{
			$q(p) \leftarrow \mathcal{G}(q(p), t_{p}, 0)$
		}
		$g_c(p) \leftarrow \mathcal{G}(q(p), t_{p+1}, t_{p})$\\
	}
	$\texttt{OMP END DO NOWAIT}$\label{alg:openmp-pipe-predict-end}\\ 
	\For{$k=1,K$}{\label{alg:openmp-pipe-parareal-start}
		$\texttt{OMP DO ORDERED}$\label{alg:openmp-pipe-ordered-do}\\
		\For{$p=0,N-1$}{
			$\texttt{OMP\_SET\_LOCK(p)}$\label{alg:openmp-pipe-fine-start}\\
			$q(p) \leftarrow \mathcal{F}_{\delta t}(q(p), t_{p+1}, t_p)$\\
			$\delta q(p) \leftarrow q(p) - q_c(p)$\\
			$\texttt{OMP\_UNSET\_LOCK(p)}$\label{alg:openmp-pipe-fine-end}\\
			$\texttt{OMP\_ORDERERD}$\label{alg:openmp-pipe-correction-start}\\
			\If{Thread is first}{
				$\texttt{OMP\_SET\_LOCK(0)}$\\
				$q(0) \leftarrow q_0$\\
				$\texttt{OMP\_UNSET\_LOCK(0)}$\\
			}
			$q_c(p) \leftarrow \mathcal{G}_{\Delta t}(q(p), t_{p+1}, t_{p})$\\
			\If{Thread not last}{
				$\texttt{OMP\_SET\_LOCK(p+1)}$\\
				$q(p+1) \leftarrow q_c(p) + \delta q(p)$\label{alg:openmp-pipe-correction}\\
				$\texttt{OMP\_UNSET\_LOCK(p+1)}$\\
			}
			$\texttt{OMP END ORDERED}\label{alg:openmp-pipe-correction-end}$
		}
		$\texttt{OMP END DO NOWAIT}$\label{alg:openmp-pipe-last-no-wait}
	}\label{alg:openmp-pipe-parareal-end}
	$\texttt{OMP END PARALLEL}$\label{alg:openmp-pipe-paraend}
\end{algorithm2e}

\begin{itemize}
	\item \textbf{Prediction step:} lines~\ref{alg:openmp-pipe-predict-start}--\ref{alg:openmp-pipe-predict-end}.
Just as in the MPI example, each thread is computing its own coarse prediction of its starting value $\vv{q}^{0}_{p}$ in a parallelised loop.
The coarse value $\mathcal{G}_{\Delta t}(\vv{q}^{0}_p, t_{p+1}, t_p)$ is also computed and stored for use in the first iteration.
The later the time slice (indicated by a higher thread number $p$), the more steps the thread must compute and thus the larger its workload.
Therefore, at the end of the coarse prediction loop, the \texttt{NOWAIT} clause is required to avoid implicit synchronisation and enable pipelining.
\footnote{Note that the prediction step could be implemented without a loop, by just using the thread number. 
However, within the Parareal iteration the \texttt{ORDERED} directive is required to guarantee a serialised execution of the update step and this directive is only available as part of a parallelised loop.
Therefore, to keep the code consistent, the prediction step is also written as a parallel loop.}	
	\item \textbf{Parareal iteration:} lines~\ref{alg:openmp-pipe-parareal-start}--\ref{alg:openmp-pipe-parareal-end}.
Here, both the fine integrator and update step are performed inside a single loop over all time slices, parallelised by \texttt{OMP DO} directives.
Because parts of the loop (the update step) have to be executed in serialised order, the \texttt{ORDERERD} directive has to be used in line~\ref{alg:openmp-pipe-ordered-do}.
Again, to avoid implicit synchronisation at the end of the loop, the \texttt{NOWAIT} clause is required in line~\ref{alg:openmp-pipe-last-no-wait}.
		\begin{itemize}
			\item \textbf{Fine integrator:} lines~\ref{alg:openmp-pipe-fine-start}--\ref{alg:openmp-pipe-fine-end}. Before the fine integrator is executed, an \texttt{OMP\_LOCK} is set to indicate that the thread will start writing into buffer $q(p)$. 
Because thread $p-1$ accesses this buffer in its update step, locks are necessary to prevent race conditions and incorrect solutions.
After the lock is set, the thread proceeds with the computation of $\mathcal{F}_{\delta t}(\vv{q}^k_p, t_{p+1}, t_p)$ and computation of the difference between coarse and fine value $\delta q$.
Then, since $q(p)$ is now up to date and $\delta q$ ready, the lock can be released.
		\item \textbf{Update step:} lines~\ref{alg:openmp-pipe-correction-start}--\ref{alg:openmp-pipe-correction-end}. 
The update step has to be performed in proper order, from first to last time slice. 
Therefore, it is enclosed in \texttt{ORDERED} directives, indicating that this part of the loop is to be executed in serial order.
Then, as in the two other versions, the update step is initialised with $\vv{q}^{k+1}_0 = \vv{q}_0$.
For every time slice, the coarse value of the updated initial guess is computed and the update performed.
The updated end value is written into buffer $q(p+1)$ to serve as the new starting value for the following time slice.
However, to prevent thread $p$ from writing into $q(p+1)$ while thread $p+1$ is still running the fine integrator, thread $p$ sets \texttt{OMP\_LOCK} number $p+1$ while performing the update.
		\end{itemize}
\end{itemize}
This implements a pipelined Parareal in OpenMP and achieves runtimes that are essentially identical to the one provided by the MPI version.
The necessity to manually control synchronisation between threads, however, makes it more complex.
As shown in Section~\ref{sec:results}, it can outperform the MPI implementation in terms of memory requirements and energy consumption.
%


intro.tex
\section{Introduction}

Computational science faces a variety of challenges stemming from the massive increase in parallelism in state-of-the-art high-performance computing systems: projections suggest that exascale machines will require 100-million way concurrency~\cite{IESR}.
This development mandates rethinking algorithms with respect to concurrency, fault-tolerance and energy efficiency but also the development of new and inherently parallel numerical methods.
As a novel direction for the parallelisation of the solution of initial value problems, parallel-in-time (or time-parallel) integration schemes are attracting a quickly growing amount of attention~\cite{Gander2015_review}.
One widely studied and used parallel-in-time method is \emph{Parareal}~\cite{LionsEtAl2001}, but other methods are frequently being introduced: RIDC~\cite{ChristliebEtAl2010}, PFASST~\cite{EmmettMinion2012}, MGRIT~\cite{FalgoutEtAl2014_MGRIT} or a DG-based time multi-grid method~\cite{Neumueller2014}.
Parareal is based on ideas from multiple shooting for boundary value problems, which have been adopted for the time dimension~\cite{ChartierPhilippe1993,Kiehl1994}.
It has been applied problems from finance~\cite{BalMaday2002} over plasma physics~\cite{SamaddarEtAl2010} to fluid flow~\cite{Celledoni2009} and neutron kinetics~\cite{BaudronEtAl2014_DDM}.
The principle efficacy of parallelisation-in-time for extreme-scale parallel computations has been demonstrated~\cite{SpeckEtAl2012} and 
recently the first software packages have also started to emerge.

Parareal's iterates between an expensive fine integrator run in parallel and a cheap coarse method which runs in serial and propagates corrections forward in time.
While the unavoidable serial part limits parallel efficiency according to Amdahl's law, some improvements are possible by using a so-called \emph{pipelined} implementation~\cite{Aubanel2011,Minion2010}.
Pipelining reduces the effective cost of the serial correction step in Parareal and therefore improves speedup.
Even more optimisation is possible by using an event-based approach~\cite{BerryEtAl2012}, but this requires a suitable execution framework that is not available on all machines.
In contrast, pipelining comes naturally when implementing Parareal in MPI.
It is, however, not straightforward in OpenMP and so far no shared memory version of Parareal with pipelining has been described.
Studies therefore use almost exclusively MPI to implement Parareal and the very few using OpenMP implementations of Parareal use only the non-pipelined version~\cite{RuprechtKrause2012,RuprechtKrause2014_DDM}.
However, using shared memory can have advantages, since it avoids e.g. the need to allocate buffers for message passing. 
The disadvantage is that naturally OpenMP is limited to a shared memory unit, typically a single compute node.
Since Parareal's convergence tends to deteriorate of too many parallel time slices are computed~\cite{GanderVandewalle2007_SISC} and given the trend to compute nodes with large numbers of cores, shared memory implementations might nevertheless be an attractive choice - as long as the benefit from pipelining does not have to be relinquished.

This paper introduces an OpenMP-based version of pipelined Parareal and compares it to a standard MPI-based implementation.
For the comparison, a special-purpose Fortran code is used~\cite{PararealF90}, which is freely available under a BSD license.
The code is deliberately designed to not use external libraries other than MPI and OpenMP to avoid interference and facilitate rerunning the benchmarks on other architectures:
since it only requires an MPI \textsc{fortran} compiler built with thread support, it can be compiled and run on almost any platform.
outro.tex
\section{Summary}
The paper introduces and analyses an OpenMP implementation of the parallel-in-time method Parareal with pipelining.
Pipelining allows to hide some of the cost of the serial coarse correction step in Parareal and is important to optimise its efficiency (even though it cannot relax the inherent limit on parallel efficiency given by the inverse of the number of required iterations).
Pipelining comes naturally in a distributed memory MPI implementation, but is not straightforward when using OpenMP.
The new OpenMP implementation is compared to a standard MPI variant in terms of runtime, memory footprint and energy consumption for both a Cray compiler and the GCC.
Both versions produce essentially identical runtimes.
For both compilers, using OpenMP leads to significant reductions in memory footprint, but the effect is more pronounced for the Cray compiler.
In terms of energy-to-solution, the results strongly depend on the compiler: while for GCC the OpenMP version is more energy efficient than the MPI version, there is no difference for the Cray compiler.

The results show that contemplating a shared memory strategy to implement ``parallel-across-the-steps'' methods like Parareal can be worthwhile.
Even though it is more complicated, it can reduce memory requirements.
A potential caveat is whether the benefits carry over to the full space-time parallel case, where a parallel-in-time method is combined with spatial decomposition.
For Parareal without pipelining the potential of such a hybrid space-time parallel approach has been illustrated~\cite{RuprechtKrause2014_DDM} but whether applies to the pipelined version introduced here remains to be seen.
paper.tex
%%
%% Copyright 2007, 2008, 2009 Elsevier Ltd
%%
%% This file is part of the 'Elsarticle Bundle'.
%% ---------------------------------------------
%%
%% It may be distributed under the conditions of the LaTeX Project Public
%% License, either version 1.2 of this license or (at your option) any
%% later version.  The latest version of this license is in
%%    http://www.latex-project.org/lppl.txt
%% and version 1.2 or later is part of all distributions of LaTeX
%% version 1999/12/01 or later.
%%
%% The list of all files belonging to the 'Elsarticle Bundle' is
%% given in the file `manifest.txt'.
%%

%% Template article for Elsevier's document class `elsarticle'
%% with numbered style bibliographic references
%% SP 2008/03/01
%%
%%
%%
%% $Id: elsarticle-template-num.tex 4 2009-10-24 08:22:58Z rishi $
%%
%%
%\documentclass[preprint,12pt]{elsarticle}

%% Use the option review to obtain double line spacing
%\documentclass[preprint,review,12pt]{elsarticle}

%% Use the options 1p,twocolumn; 3p; 3p,twocolumn; 5p; or 5p,twocolumn
%% for a journal layout:
%\documentclass[final,1p,times]{elsarticle}
%% \documentclass[final,1p,times,twocolumn]{elsarticle}
\documentclass[final,3p,times]{elsarticle}
%% \documentclass[final,3p,times,twocolumn]{elsarticle}
%% \documentclass[final,5p,times]{elsarticle}
%% \documentclass[final,5p,times,twocolumn]{elsarticle}

%% if you use PostScript figures in your article
%% use the graphics package for simple commands
%% \usepackage{graphics}
%% or use the graphicx package for more complicated commands
%% \usepackage{graphicx}
%% or use the epsfig package if you prefer to use the old commands
%% \usepackage{epsfig}

%% The amssymb package provides various useful mathematical symbols
\usepackage{amssymb}
%% The amsthm package provides extended theorem environments
%% \usepackage{amsthm}

\usepackage{amsmath}
\usepackage[ruled,algo2e,linesnumbered,algonl]{algorithm2e}
\usepackage[free-standing-units,group-separator={,},binary-units=true]{siunitx}

%% The lineno packages adds line numbers. Start line numbering with
%% \begin{linenumbers}, end it with \end{linenumbers}. Or switch it on
%% for the whole article with \linenumbers after \end{frontmatter}.
%% \usepackage{lineno}

%% natbib.sty is loaded by default. However, natbib options can be
%% provided with \biboptions{...} command. Following options are
%% valid:

%%   round  -  round parentheses are used (default)
%%   square -  square brackets are used   [option]
%%   curly  -  curly braces are used      {option}
%%   angle  -  angle brackets are used    <option>
%%   semicolon  -  multiple citations separated by semi-colon
%%   colon  - same as semicolon, an earlier confusion
%%   comma  -  separated by comma
%%   numbers-  selects numerical citations
%%   super  -  numerical citations as superscripts
%%   sort   -  sorts multiple citations according to order in ref. list
%%   sort&compress   -  like sort, but also compresses numerical citations
%%   compress - compresses without sorting
%%
%% \biboptions{comma,round}

% \biboptions{}


%
% OWN COMMANDS
\usepackage{color}
\newcommand{\todo}[1]{\textcolor{red}{[TODO -- #1]}}
\newcommand{\vv}[1]{\mathbf{#1}}
\usepackage[colorlinks = true,
            linkcolor = blue,
            urlcolor  = blue,
            citecolor = blue,
            anchorcolor = blue]{hyperref}
\usepackage{caption}
\usepackage{subcaption}            
            
\journal{}

\begin{document}

\begin{frontmatter}

%% Title, authors and addresses

%% use the tnoteref command within \title for footnotes;
%% use the tnotetext command for the associated footnote;
%% use the fnref command within \author or \address for footnotes;
%% use the fntext command for the associated footnote;
%% use the corref command within \author for corresponding author footnotes;
%% use the cortext command for the associated footnote;
%% use the ead command for the email address,
%% and the form \ead[url] for the home page:
%%
%% \title{Title\tnoteref{label1}}
%% \tnotetext[label1]{}
%% \author{Name\corref{cor1}\fnref{label2}}
%% \ead{email address}
%% \ead[url]{home page}
%% \fntext[label2]{}
%% \cortext[cor1]{}
%% \address{Address\fnref{label3}}
%% \fntext[label3]{}

\title{A shared memory implementation of pipelined Parareal}

%% use optional labels to link authors explicitly to addresses:
%% \author[label1,label2]{<author name>}
%% \address[label1]{<address>}
%% \address[label2]{<address>}

\author[lab1,lab2]{Daniel Ruprecht}
\ead{d.ruprecht@leeds.ac.uk}
\address[lab1]{School of Mechanical Engineering, University of Leeds, Woodhouse Lane, Leeds LS2 9JT, UK}
\address[lab2]{Centre for Computational Medicine in Cardiology, Institute of Computational Science, Universit{\`a} della Svizzera italiana, Via Giuseppe Buffi 13, CH-6900 Lugano, Switzerland}

\begin{abstract}
%% Text of abstract
The paper introduces an OpenMP implementation of pipelined Parareal and compares it to a standard MPI-based implementation.
Both versions yield essentially identical runtimes, but, depending on the compiler, the OpenMP variant consumes about $7$\% less energy.
However, its key advantage is a significantly smaller memory footprint.
The higher implementation complexity, including manual control of locks, might make it difficult to use in legacy codes, though.
\end{abstract}

\begin{keyword}
%% keywords here, in the form: keyword \sep keyword
Parareal \sep parallel-in-time integration \sep pipelining \sep shared memory \sep memory footprint \sep energy-to-solution
%% MSC codes here, in the form: \MSC code \sep code
\MSC[2010] 68W10 \sep 65Y05 \sep 68N19
%% or \MSC[2008] code \sep code (2000 is the default)
\end{keyword}

\end{frontmatter}

%%
%% Start line numbering here if you want
%%
% \linenumbers
\input{intro}

\input{parareal}

\input{implementation}

\input{results}

\input{outro}

\section*{Acknowledgments}
I would like to thank the Centre for Interdisciplinary Research (ZIF) at the University of Bielefeld, Germany, for inviting me for a research visit in August 2015. 
The writing of this article greatly benefited from the tranquil and productive atmosphere at ZIF.
I gratefully acknowledge Andrea Arteaga's support with the energy measurements.
%\bibliographystyle{model1-num-names}
\bibliographystyle{elsarticle}
\bibliography{pint_mod,hpc_mod,MyCodes,other}

%% Authors are advised to submit their bibtex database files. They are
%% requested to list a bibtex style file in the manuscript if they do
%% not want to use elsarticle-num.bst.

%% References without bibTeX database:

% \begin{thebibliography}{00}

%% \bibitem must have the following form:
%%   \bibitem{key}...
%%

% \bibitem{}

% \end{thebibliography}


\end{document}

%%
%% End of file `elsarticle-template-num.tex'.

parareal.tex
\section{The method: Parareal}
The starting point for Parareal is an initial value problem of the form
\begin{equation}
	\label{eq:ivp}
	\dot{\vv{q}} = \mathbf{f} (\mathbf{q}(t),t), \ \vv{q}(0) = \vv{q}_{0}, \ t \in (0, T],
\end{equation}
which, in the numerical examples below, arises from the spatial discretisation of a PDE ("method-of-lines") with $\vv{q} \in \mathbb{R}^{N_{\rm dof}}$ being a vector containing all degrees-of-freedom.
Let 
$\mathcal{F}_{\delta t}$ denote a numerical procedure for the approximate solution of \eqref{eq:ivp}, for example a Runge-Kutta method. 
Denote further by
\begin{equation}
	\vv{q} = \mathcal{F}_{\delta t}\left( \vv{\tilde{q}}, t_{2}, t_{1} \right)
\end{equation}
the result of approximately integrating \eqref{eq:ivp} forward in time from some starting value\footnote{Here, the prescribed value $\vv{q}_0$ at $t=0$ is referred to as initial value, while the value $\vv{q}_p$ the beginning of a time slice is called starting value.} $\mathbf{\tilde{q}}$ at a time $t_{1}$ to a time $t_{2} > t_{1}$ using $\mathcal{F}_{\delta t}$.
That is, 
\begin{equation}
	\label{eq:fine_serial}
	\vv{q} = \mathcal{F}_{\delta t}( \vv{q}_0, T, 0)
\end{equation}
would indicate sequentially solving the full initial value problem in serial using the time stepper denoted by $\mathcal{F}_{\delta t}$.
In order to parallelise the numerical integration of~\eqref{eq:ivp}, Parareal and other so-called \emph{parallel-across-the-steps} methods~\cite{Burrage1997} introduce a decomposition of the time interval $[0,T]$ into time-slices $[t_{p}, t_{p+1}]$, $p=0, \ldots, P-1$ where $P$ is the number of cores, equal to the number of processes or threads, to be used in time.
In the implementations sketched in Section~\ref{sec:implementation}, time slice $[t_{p}, t_{p+1}]$ is assigned to the core running either thread number $p$ (in the OpenMP variant) or the process with rank $p$ (in the MPI version).
For simplicity, assume here that all time slices have the same length and that the whole interval $[0,T]$ is covered with a single set of time slices - otherwise, some form of restarting or \emph{moving window}~\cite{SchreiberEtAl2015} has to be used.

Introduce now a second time integrator denoted $\mathcal{G}_{\Delta t}$, which has to be much cheaper to compute but can also be much less accurate (commonly referred to as the ''coarse propagator'').
Typically, a lower order method with a larger time step is used here, but a lower order spatial discretisation can be used as well.
Using fewer degrees of freedom on the coarse level is also possible, but necessitates the introduction of suitable transfer operators (interpolation and restriction) and can degrade convergence.
Parareal starts off with a \emph{prediction step}, computing a rough guess of the starting value $\vv{q}^0_{p}$ at the beginning of each time slice by
\begin{equation}
	\vv{q}^0_{p} = \mathcal{G}_{\Delta t}(\vv{q}_0, t_{p}, 0), \quad p=0,\ldots,P-1.
\end{equation}
Here, subscript $p$ indicates an approximation of the solution at time $t_p$.
These approximate starting values are used to start the following iteration, run concurrently on each time slice,
\begin{equation}
	\label{eq:parareal}
	\vv{q}^{k}_{p+1} = \mathcal{G}_{\Delta t}(\vv{q}^{k}_{p}, t_{p+1}, t_{p}) + \mathcal{F}_{\delta t}(\vv{q}^{k-1}_{p}, t_{p+1}, t_{p}) - \mathcal{G}_{\Delta t}(\vv{q}^{k-1}_{p}, t_{p+1}, t_{p}), \quad p=0,\ldots,P-1, \quad k=1,\ldots,K.
\end{equation}
As $k$ increases, this iteration converges at the endpoints of the slices to the same solution provided by~\eqref{eq:fine_serial}, that is $\vv{q}^{k}_p \to \mathcal{F}_{\delta t}(\vv{q}_0, t_p, 0)$ for $p=0,\ldots,P$.
Because the computationally expensive evaluation of the fine propagator (referred to as the \emph{fine integrator step} in Section~\ref{sec:implementation}) can be parallelised across time slices, iteration~\eqref{eq:parareal} can run in less wall clock time than the direct time-serial integration~\eqref{eq:fine_serial} -- provided the coarse method is cheap enough and the number of required iterations $K$ is small.
Computation of the fine values is followed by the \emph{correction step}: the updated value $\vv{q}^{k}_{p}$ from the previous time slice is obtained, the coarse integrator is applied to it and~\eqref{eq:parareal} is evaluated to obtain $\vv{q}^{k}_{p+1}$ which is then used to correct the value on the next time slice and so on.
Note that, in contrast to the fine integrator, the correction step has to be performed in correct order, going step by step from the first time slice to the last.
Also, when using a distributed memory parallelisation, the value $\vv{q}_{p}^{k}$, required in~\eqref{eq:parareal} to compute $\mathcal{G}_{\Delta t}(\vv{q}_{p}^{k}, t_{p+1}, t_p)$, has to be communicated from the process handling time slice $[t_{p-1}, t_{p}]$ to the one handling $[t_{p},t_{p+1}]$.

\subsection{Performance model}
The expected performance of Parareal, referred to here as \emph{projected speedup}, can be described by a simple theoretical model~\cite{Minion2010}.
Here, the model is briefly repeated to more clearly illustrate the effect and importance of pipelining Parareal, see Subsection~\ref{subsec:pipelining}.
It is also used as a baseline to compare against the measured speedups reported in Section~\ref{sec:results}.

Denote, as above, by $P$ the number of cores in time, equal to the number of time slices.
Further, denote by $c_{\rm c}$ the cost of integrating over one time slice using $\mathcal{G}_{\Delta t}$ and by $c_{\rm f}$ the cost when using $\mathcal{F}_{\delta t}$.
Because all time slices are assumed to consist of the same number of steps and an explicit method is used here, it can be assumed that $c_{\rm f}$ and $c_{\rm c}$ are identical for all time slices.
This does not necessarily hold for an implicit method where differences in the number of iterations required by the nonlinear solver to converge can introduce load imbalances.
A simple theoretical model for the speedup of Parareal using $K$ iterations against running the fine method in serial now reads
\begin{equation}
	\label{eq:speedup_nopipe}
	s_{\rm np}(P) = \frac{ P c_{\rm f} }{ \left(1 + K \right) P c_{\rm c} + K c_{\rm f}} = \frac{1}{ \left(1 + K \right) \frac{c_{\rm c}}{c_{\rm f}} + \frac{K}{P}}.
\end{equation}
Equation~\eqref{eq:speedup_nopipe} illustrates the necessity for a cheap coarse method to minimise the ratio $c_{\rm c}/c_{\rm f}$ while still guaranteeing rapid convergence to minimise the terms $K/P$ and $1+K$.

\subsection{Pipelining Parareal}\label{subsec:pipelining}
\begin{figure}[t]
	\centering
	\begin{minipage}[t]{0.4\textwidth}
    	        \vspace{0pt} 
		\includegraphics[width=0.99\textwidth]{images/NonpipelinedParareal.pdf}
	\end{minipage}\hspace*{0.1\textwidth}
	\begin{minipage}[t]{0.4\textwidth}
     	        \vspace{0pt} 
		\includegraphics[width=0.99\textwidth]{images/PipelinedParareal.pdf}
	\end{minipage}
	\caption{Execution diagram for Parareal without (left) and with (right) pipelining. Red blocks correspond to $c_{\text{f}}$, the time needed to run the fine integrator $\mathcal{F}_{\delta t}$ over one time slice $[t_{p}, t_{p+1}]$ while blue blocks correspond to $c_{\text{c}}$, the time required by the coarse method $\mathcal{G}_{\Delta t}$. Pipelining allows to hide some of the cost of the coarse method. While it comes naturally when using MPI to implement Parareal (there, \texttt{Thread} would refer to a \texttt{Process}), using simple loop-based parallelism with OpenMP results in the non-pipelined version shown on the left. A pipelined OpenMP Parareal is introduced in Section~\ref{sec:implementation}.}	
	\label{fig:pipelining}	
\end{figure}
It has been pointed out that a proper implementation of Parareal allows to hide some of the cost of the coarse propagator and the name \emph{pipelining} has been coined for this approach~\cite{Minion2010}.
Figure~\ref{fig:pipelining} sketches the execution diagrams of both a non-pipelined (left) and pipelined (right) implementation for four time slices.
As can be seen, pipelining reduces the effective cost of the coarse correction step in each iteration from $P \times c_{\textrm{c}}$ to $c_{\textrm{c}}$ -- but note that the initial prediction step still has cost $P \times c_{\rm c}$ as before.
For pipelined Parareal, estimate~\eqref{eq:speedup_nopipe} changes to
\begin{equation}
	\label{eq:speedup_pipe}
	s_{\rm p}(P) = \frac{P c_{\rm f}}{P c_{\rm c} + K c_{\rm c} + K c_{\rm f}} = \frac{1}{\left(1 + \frac{K}{P} \right) \frac{c_{\rm c}}{c_{\rm f}} + \frac{K}{P}}.
\end{equation}
Because $K/P \ll K$, the pipelined version allows for better speedup, that is $s_{\rm np}(P) \leq s_{\rm p}(P)$.
However, because pipelining only hides cost from the coarse integrator, the effect is smaller when the coarse method is very cheap and $c_{\rm c}/c_{\rm s} \ll 1$. In that case, the term $K/P$ dominates the estimate which is not affected by pipelining.

results.tex
\section{Numerical results}\label{sec:results}
In this section, the OpenMP and MPI implementation are compared with respect to runtime in Subsection~\ref{subsec:runtime}, memory footprint in Subsection~\ref{subsec:memory} and energy consumption in Subsection~\ref{subsec:energy}.
The parameters for the simulation are a viscosity parameter of $\nu = 0.02$ and a spatial discretisation on both levels with $N_x = N_y = N_z = 40$ grid points in every direction.
The simulation is run until $T=1.0$ with a coarse time step of $\Delta t = 1/192$ and a fine step of $\delta t = 1/240$.
Because of the quite high computational cost of the WENO-5 method in comparison to a cheap first order upwind scheme and the fact that RK3SSP needs three evaluations of the right hand side per step while the Euler method needs only one, the coarse propagator is about a factor of forty faster, despite the fact that the coarse step is only a factor of $1.25$ larger than the fine.

To fix the number of iterations to a meaningful value which guarantees comparable accuracy from Parareal and serial fine integrator, we estimate the discretisation error of $\mathcal{F}_{\delta t}$ by comparing against a reference solution with time step $\delta t / 10$.
This gives estimates for the fine relative error at $T=1$ of about $e_{\rm fine} \approx 5.9 \times 10^{-5}$ and for the coarse error of about $e_{\rm coarse} \approx 7.3 \times 10^{-2}$.
For $P=24$ time slices, after three iterations, the defect between Parareal and the fine solution is approximately $1.4 \times 10^{-4}$, after four iterations $1.5 \times 10^{-5}$.
We therefore fix the number of iterations to $K=4$ so that for all values of $P$ Parareal produces a solution with the same accuracy as the fine integrator.

\subsection{Hardware}\label{subsec:hardware}
Benchmarks are run on a single node of two different systems, respectively.
The first one is \textsc{cub}, a commodity Linux cluster at the Institute of Computational Science in Lugano, consisting of 3x14 IBM Blade nodes.
Each node has two quad-core Opteron (Barcelona) CPUs, for a total of 8 cores per node, and 16 GigaByte main memory.
Nodes are connected through an Infiniband interconnect.
The used MPI implementation is OpenMPI-1.4.2\footnote{\url{http://www.open-mpi.org/}}, compiled with GCC-4.5.0\footnote{\url{https://gcc.gnu.org/}} with flags \texttt{-O3} for maximum optimisation and \texttt{-fopenmp} to enable OpenMP.
Note that, as the code is stand alone, no external libraries have to be linked.

The second system is \textsc{Piz Dora} at the Swiss National Supercomputing Centre.\footnote{Detailed specification can be found here \url{http://www.cscs.ch/computers/piz_daint/index.html}} \textsc{Dora} is a Cray XC40 with a total of $1,256$ compute nodes.
Each node contains two $12$-core Intel Haswell CPUs and has 64 or 128 GigaByte of RAM and nodes are connected through a Cray Aries interconnect, using a dragonfly network topology.
Two compilers are tested, the GCC-4.9.2 and the Cray Fortan compiler version 8.3.12.
Both use the MPICH MPI library version 7.2.2.
As on \textsc{cub}, compiler flags \texttt{-O3} and \texttt{-fopenmp} (GCC) or \texttt{omp} (Cray compiler) are used.
Performance data for each completed job is generated using the Cray \emph{Resource Utilisation Reporting} tool RUR~\cite{Rur2013}.
RUR collects compute node statistics before and after each job and provides data on user and system time, maximum memory used, amount of I/O operations, consumed energy and other metrics.
However, it only collects data for a full node and not for individual CPUs or cores.

\subsection{Wall clock time and speedup}\label{subsec:runtime}
%
\begin{figure}[t]
	\centering
	\begin{subfigure}[t]{0.49\textwidth}
		\centering
		\caption{\textsc{dora}}\label{fig:runtime_dora}
		\includegraphics[scale=1]{images/runtime_dora.pdf}
	\end{subfigure}
	\begin{subfigure}[t]{0.49\textwidth}
		\centering
		\caption{\textsc{cub}}\label{fig:runtime_cub}
		\includegraphics[scale=1]{images/runtime_cub.pdf}
	\end{subfigure}
    \caption{Five run averages of runtime with relative standard deviation below $0.01$ on \textsc{dora} (\ref{fig:runtime_dora}) and $0.025$ on \textsc{cub} (\ref{fig:runtime_cub}).}
        \label{fig:runtimes}      
\end{figure}
\begin{figure}[h!t]
	\centering
	\begin{subfigure}[t]{0.49\textwidth}
		\centering
		\caption{\textsc{dora}}\label{fig:speedup_dora}
		\includegraphics[scale=1]{images/speedup_dora.pdf}
	\end{subfigure}
	\begin{subfigure}[t]{0.49\textwidth}
		\centering
		\caption{\textsc{cub}}\label{fig:speedup_cub}
		\includegraphics[scale=1]{images/speedup_cub.pdf}
	\end{subfigure}
    \caption{Speedup computed from average runtimes shown in Figure~\ref{fig:runtimes} on \textsc{dora} (\ref{fig:speedup_dora}) and \textsc{cub} (\ref{fig:speedup_cub}).}
        \label{fig:speedup}        
\end{figure}
%
At first, runtime and speedup compared to the serial execution of the fine integrator are assessed.
On both \textsc{Cub} and \textsc{Dora}, for each variant of Parareal and each value of $P$, five runs are performed and the average runtime is reported here.
Measured runtimes are quite stable across different runs: the largest relative standard deviation of all performed five-run ensembles is smaller than $0.025$ on \textsc{cub} and smaller than $0.01$ on \textsc{dora}.
Therefore, plots show only the average values without error bars, because those are hardly recognisable and clutter the figure.

Figure~\ref{fig:runtimes} shows runtimes in seconds depending on the number of cores on \textsc{Dora} (left) and \textsc{cub} (right).
For \textsc{Dora}, only results from the Cray compiler are shown, which tends to generate slightly faster code than the GCC.
The runtime of the serial fine integrator is indicated by a horizontal black line.
Note that the y axis is scaled logarithmically, so the distance from $5$ to $10$ seconds is the same as from $10$ to $20$.
Because the more modern CPUs on \textsc{dora} are faster, runtimes are generally smaller on \textsc{dora} than on \textsc{cub}.
For $P=8$, for example, Parareal runtimes on \textsc{cub} are around $10~\text{seconds}$ but only around $7~\text{seconds}$ on \textsc{dora}.
Both versions give almost identical performance: on both \textsc{dora} and \textsc{cub}, OpenMP is marginally faster than the MPI version.

In addition, Figure~\ref{fig:speedup} shows the speedup relative to the fine integrator run serially.
Both versions fall short of the theoretically possible speedup indicated by the black line, but differences between MPI and OpenMP are small.
As far as runtimes and speedup are concerned, there is no indication that using the more complex OpenMP version provides benefits.


\subsection{Memory footprint}\label{subsec:memory}
\begin{figure}[t!]
	\centering
	\begin{subfigure}[t]{0.49\textwidth}
		\centering
		\caption{GCC compiler}
		\label{fig:memory_gcc}
		\includegraphics[scale=1.0]{images/memory_dora_gcc.pdf}	
	\end{subfigure}
	\begin{subfigure}[t]{0.49\textwidth}
		\centering
		\caption{Cray compiler}
		\label{fig:memory_cray}
		\includegraphics[scale=1.0]{images/memory_dora_cray.pdf}
	\end{subfigure}
	\caption{Maximum memory allocated in MegaByte for GCC (\ref{fig:memory_gcc}) and Cray compiler (\ref{fig:memory_cray}) for the three different versions of Parareal depending on the number of used cores $P$. The black line indicates expected memory consumption computed as number of cores time memory footprint of serial fine integrator.}
	\label{fig:memory}
\end{figure}
The memory footprint of the code is measured only on \textsc{dora} where RUR is available.
In contrast to runtime and energy, the memory footprint, as expected, does not vary between runs. 
Therefore~Figure~\ref{fig:memory} shows a visualisation of the data from a single run with no averaging.
The bars indicate the maximum required memory in MegaByte (\SI{}{\mega\byte}) while the black line indicates the expected memory consumption using $P$ cores computed as
\begin{equation}
	\label{eq:memory}
	m(P) = P \times m_{\rm serial}
\end{equation}
where $m_{\rm serial}$ is the value measured for a reference run of the fine integrator.
Because copies of the solution have to be stored for every time slice, the total memory required for Parareal can be expected to increase linearly with the number of cores in time.
Note, however, that memory required \emph{per core} stays constant if it follows~\eqref{eq:memory}.

For the OpenMP variant compiled with GCC, the memory footprint shown in Figure~\ref{fig:memory_gcc} exactly matches the expected values.
The Cray compiler, shown in Figure~\ref{fig:memory_cray}, leads to a smaller than expected memory footprint, but memory requirements still increase linearly with the number of time slices.
For both compilers, the MPI version causes a noticeable overhead in terms of memory footprint, most likely because of internal allocation of additional buffers for sending and receiving~\cite{Rabenseifner2009}.
For both OpenMP and MPI, the total memory footprint is \emph{larger} for the Cray than for the GCC compiler, but the effect is much more pronounced for the MPI implementation (\SI{329}{\mega\byte} versus \si{261}{\mega\byte}) than for the OpenMP variant (\si{200}{\mega\byte} versus \si{193}{\mega\byte}).

It is important to note that both implementations allocate three auxiliary buffers per core.
The overhead in terms of memory in MPI does thus not simply stem from allocating an additional buffer for communication, but comes from within the MPI library.
The OpenMP implementation avoids this overhead.
Given that memory will be a much more precious resource on future supercomputers and that the memory requirements of ''across-the-steps'' parallel-in-time methods have raised concerns~\cite{Keyes2011}, these savings might be important.
Since the memory overhead grows as more and more cores are used in time, savings from OpenMP will be especially pronounced when Parareal is used with many time slices on nodes with a large numbers of cores.

\subsection{Energy-to-solution}\label{subsec:energy}
%
\begin{figure}[t!]
	\centering
	\begin{subfigure}[t]{0.49\textwidth}
		\centering
		\caption{GCC compiler}\label{fig:energy_gcc}
		\includegraphics[scale=1.0]{images/energy_dora_gcc.pdf}
	\end{subfigure}
	\begin{subfigure}[t]{0.49\textwidth}
		\centering
		\caption{Cray compiler}\label{fig:energy_cray}
		\includegraphics[scale=1.0]{images/energy_dora_cray.pdf}
	\end{subfigure}
	\caption{Energy-to-solution in Kilojoule for GCC in (\ref{fig:energy_gcc}) and Cray compiler (\ref{fig:energy_cray}) for the three versions of Parareal, averaged across $50$ samples, each using all $24$ cores of the node. The largest relative standard deviation is $0.07$, error bars indicate $95\%$ confidence intervals.}
	\label{fig:energy}
\end{figure}
The RUR tool reports the energy-to-solution for every completed job.
Because RUR can only measure energy usage for a full node, results are reported here for runs using the full number of $P=24$ cores available on a \textsc{dora} node.
In contrast to runtimes and memory footprint, energy measurements show significant variations between runs due to random fluctuations: thus, the presented values are averages over ensembles of $50$ runs for each version of Parareal.
This number of runs has been sufficient to reduce the relative standard deviation to below $0.09$ in both configurations and therefore gives a robust indication of actual energy requirements.

When using the GCC compiler, the MPI version consumes more energy than OpenMP.
The resulting $95\%$ confidence intervals (assuming energy-to-solution is normally distributed), are $844.04 \pm 15.89$\si{}{\joule} for MPI and $783.72 \pm 11.53$\si{}{\joule} for OpenMP.
Figure~\ref{fig:energy} gives a graphical representation of these values including the confidence intervals.
The average for OpenMP is well outside the confidence interval for the MPI version, so this is very unlikely just a chance result.
Moreover, because runtimes are almost identical, the differences in energy-to-solution cannot simply be attributed to differences in time-to-solution.
It is also noteworthy that for the GCC compiler OpenMP without pipelining (results not shown here), even though it is \emph{slower} than the MPI version, still requires less energy ($844.04 \pm 15.89$\si{}{\joule} versus $801.14 \pm 13.18$\si{}{\joule}).
Tracking down the precise reason for the differences in energy-to-solution and power requirement will require detailed tracing of power uptake, which is only possible on specially prepared machines~\cite{Canturk2003} and thus left for future work.

Note that the energy consumption of Parareal has previously been studied~\cite{ArteagaEtAl2015}.
By comparing against a simple theoretical model, it has been shown that the energy overhead of Parareal (defined as energy-to-solution of Parareal divided by energy-to-solution of the fine serial integrator), is mostly due to Parareal's intrinsic suboptimal parallel efficiency.
While improving parallel efficiency of parallel-in-time integration remains the main avenue for improving energy efficiency, the results here suggest that in some cases a shared memory approach can provide non-trivial additional savings.

For code generated with the Cray compiler, both OpenMP and MPI lead to almost identical energy requirements:
here, confidence intervals are $784.24 \pm 11.93$\si{}{\joule} for MPI and $784.24 \pm 12.18$\si{}{\joule} for OpenMP.
It seems likely that the compiler optimises the message passing to take advantage of the shared memory on the single node.
Supposedly, the MPI version handles communication in a way that is similar to what is explicitly coded in the OpenMP version.
However, as shown in Subsection~\ref{subsec:memory}, this automatic optimisation comes at the expense of a significantly larger memory footprint.

%
%\begin{table}[h]
%\centering
%\begin{tabular}{|c|ccc|} \hline
%		        \multicolumn{4}{|c|}{\textbf{GCC compiler}} \\ \hline
%		       &  $\varnothing$energy (Joule) &  $\varnothing$runtime (sec.) & power (Watt) \\ \hline
%MPI                 & 844.04              & 2.455               & 343.80 \\
%OpenMP         & 801.14              & 2.930               & 273.38 \\
%OpenMP(pipe) & 783.72             & 2.448               & 320.15 \\ \hline
%\end{tabular}
%\begin{tabular}{|c|ccc|} \hline
%		        \multicolumn{4}{|c|}{\textbf{Cray compiler}} \\ \hline
%		       &  $\varnothing$energy (Joule) &  $\varnothing$runtime (sec.) & power (Watt) \\ \hline
%MPI                 & 784.24              & 2.146               & 365.44  \\
%OpenMP         &  833.12             &  2.400              &  347.13 \\
%OpenMP(pipe) & 784.72             &  2.088              &  375.82 \\ \hline
%\end{tabular}
%\caption{A rough estimate of power consumption obtained by dividing averaged energy-to-solution by averaged runtimes for the three different variants run on $P=24$ cores. The differences in estimated power uptake strongly suggest that the differences in consumed energy are not due to different runtimes alone.}
%\label{tab:power}
%\end{table}
paper.bbl
\begin{thebibliography}{10}
\expandafter\ifx\csname url\endcsname\relax
  \def\url#1{\texttt{#1}}\fi
\expandafter\ifx\csname urlprefix\endcsname\relax\def\urlprefix{URL }\fi
\expandafter\ifx\csname href\endcsname\relax
  \def\href#1#2{#2} \def\path#1{#1}\fi

\bibitem{IESR}
J.~Dongarra, P.~Beckman, al., The international exascale software roadmap, Int.
  Journal of High Performance Computer Applications 25~(1).
\newblock \href {http://dx.doi.org/10.1177/1094342010391989}
  {\path{doi:10.1177/1094342010391989}}.

\bibitem{Gander2015_review}
M.~J. Gander,
  \href{{http://www.unige.ch/%7Egander/Preprints/50YearsTimeParallel.pdf}}{{50
  years of Time Parallel Time Integration}}, in: Multiple Shooting and Time
  Domain Decomposition, Springer, 2015.
\newline\urlprefix\url{{http://www.unige.ch/%7Egander/Preprints/50YearsTimeParallel.pdf}}

\bibitem{LionsEtAl2001}
J.-L. Lions, Y.~Maday, G.~Turinici, {A "parareal" in time discretization of
  {PDE}'s}, Comptes Rendus de l'Académie des Sciences - Series I - Mathematics
  332 (2001) 661--668.
\newblock \href {http://dx.doi.org/10.1016/S0764-4442(00)01793-6}
  {\path{doi:10.1016/S0764-4442(00)01793-6}}.

\bibitem{ChristliebEtAl2010}
A.~J. Christlieb, C.~B. Macdonald, B.~W. Ong, {Parallel high-order
  integrators}, SIAM Journal on Scientific Computing 32~(2) (2010) 818--835.
\newblock \href {http://dx.doi.org/10.1137/09075740X}
  {\path{doi:10.1137/09075740X}}.

\bibitem{EmmettMinion2012}
M.~Emmett, M.~L. Minion, {Toward an Efficient Parallel in Time Method for
  Partial Differential Equations}, Communications in Applied Mathematics and
  Computational Science 7 (2012) 105--132.
\newblock \href {http://dx.doi.org/10.2140/camcos.2012.7.105}
  {\path{doi:10.2140/camcos.2012.7.105}}.

\bibitem{FalgoutEtAl2014_MGRIT}
R.~D. Falgout, S.~Friedhoff, T.~V. Kolev, S.~P. MacLachlan, J.~B. Schroder,
  {Parallel time integration with multigrid}, SIAM Journal on Scientific
  Computing 36 (2014) C635--C661.
\newblock \href {http://dx.doi.org/10.1137/130944230}
  {\path{doi:10.1137/130944230}}.

\bibitem{Neumueller2014}
M.~J. Gander, M.~Neumueller, \href{http://arxiv.org/abs/1409.5254}{{Analysis of
  a Time Multigrid Algorithm for {DG}-Discretizations in Time}} (2014).
\newline\urlprefix\url{http://arxiv.org/abs/1409.5254}

\bibitem{ChartierPhilippe1993}
P.~Chartier, B.~Philippe, {A parallel shooting technique for solving
  dissipative {ODE}'s}, Computing 51~(3-4) (1993) 209--236.
\newblock \href {http://dx.doi.org/10.1007/BF02238534}
  {\path{doi:10.1007/BF02238534}}.

\bibitem{Kiehl1994}
M.~Kiehl, {Parallel multiple shooting for the solution of initial value
  problems}, Parallel Computing 20~(3) (1994) 275--295.
\newblock \href {http://dx.doi.org/10.1016/S0167-8191(06)80013-X}
  {\path{doi:10.1016/S0167-8191(06)80013-X}}.

\bibitem{BalMaday2002}
G.~Bal, Y.~Maday, {A "Parareal" time discretization for non-linear {PDE}'s with
  application to the pricing of an American Put}, in: L.~Pavarino, A.~Toselli
  (Eds.), {Recent Developments in Domain Decomposition Methods}, Vol.~23 of
  {Lecture Notes in Computational Science and Engineering}, Springer Berlin,
  2002, pp. 189--202.
\newblock \href {http://dx.doi.org/10.1007/978-3-642-56118-4_12}
  {\path{doi:10.1007/978-3-642-56118-4_12}}.

\bibitem{SamaddarEtAl2010}
D.~Samaddar, D.~E. Newman, R.~S. Sánchez, {Parallelization in time of
  numerical simulations of fully-developed plasma turbulence using the parareal
  algorithm}, Journal of Computational Physics 229 (2010) 6558--6573.
\newblock \href {http://dx.doi.org/10.1016/j.jcp.2010.05.012}
  {\path{doi:10.1016/j.jcp.2010.05.012}}.

\bibitem{Celledoni2009}
E.~Celledoni, T.~Kvamsdal, {Parallelization in time for thermo-viscoplastic
  problems in extrusion of aluminium}, International Journal for Numerical
  Methods in Engineering 79~(5) (2009) 576--598.
\newblock \href {http://dx.doi.org/10.1002/nme.2585}
  {\path{doi:10.1002/nme.2585}}.

\bibitem{BaudronEtAl2014_DDM}
A.-M. Baudron, J.-J. Lautard, Y.~Maday, O.~Mula, {The parareal in time
  algorithm applied to the kinetic neutron diffusion equation}, in: {Domain
  Decomposition Methods in Science and Engineering XXI}, {Lecture Notes in
  Computational Science and Engineering}, Springer International Publishing,
  2014, pp. 437--445.
\newblock \href {http://dx.doi.org/10.1007/978-3-319-05789-7_41}
  {\path{doi:10.1007/978-3-319-05789-7_41}}.

\bibitem{SpeckEtAl2012}
R.~Speck, D.~Ruprecht, R.~Krause, M.~Emmett, M.~L. Minion, M.~Winkel,
  P.~Gibbon, {A massively space-time parallel {N}-body solver}, in:
  {Proceedings of the International Conference on High Performance Computing,
  Networking, Storage and Analysis}, {SC '12}, IEEE Computer Society Press, Los
  Alamitos, CA, USA, 2012, pp. 92:1--92:11.
\newblock \href {http://dx.doi.org/10.1109/SC.2012.6}
  {\path{doi:10.1109/SC.2012.6}}.

\bibitem{Aubanel2011}
E.~Aubanel, {Scheduling of Tasks in the Parareal Algorithm}, Parallel Computing
  37 (2011) 172--182.
\newblock \href {http://dx.doi.org/10.1016/j.parco.2010.10.004}
  {\path{doi:10.1016/j.parco.2010.10.004}}.

\bibitem{Minion2010}
M.~L. Minion, {A Hybrid Parareal Spectral Deferred Corrections Method},
  Communications in Applied Mathematics and Computational Science 5~(2) (2010)
  265--301.
\newblock \href {http://dx.doi.org/10.2140/camcos.2010.5.265}
  {\path{doi:10.2140/camcos.2010.5.265}}.

\bibitem{BerryEtAl2012}
L.~A. Berry, W.~R. Elwasif, J.~M. Reynolds-Barredo, D.~Samaddar, R.~S.
  Sánchez, D.~E. Newman, {Event-based parareal: A data-flow based
  implementation of parareal}, Journal of Computational Physics 231~(17) (2012)
  5945--5954.
\newblock \href {http://dx.doi.org/10.1016/j.jcp.2012.05.016}
  {\path{doi:10.1016/j.jcp.2012.05.016}}.

\bibitem{RuprechtKrause2012}
D.~Ruprecht, R.~Krause, {Explicit parallel-in-time integration of a linear
  acoustic-advection system}, Computers \& Fluids 59~(0) (2012) 72--83.
\newblock \href {http://dx.doi.org/10.1016/j.compfluid.2012.02.015}
  {\path{doi:10.1016/j.compfluid.2012.02.015}}.

\bibitem{RuprechtKrause2014_DDM}
R.~Krause, D.~Ruprecht, {Hybrid Space-Time Parallel Solution of {B}urgers'
  Equation}, in: {Domain Decomposition Methods in Science and Engineering
  {XXI}}, Vol.~98 of {Lecture Notes in Computational Science and Engineering},
  Springer International Publishing, 2014, pp. 647--655.
\newblock \href {http://dx.doi.org/10.1007/978-3-319-05789-7_62}
  {\path{doi:10.1007/978-3-319-05789-7_62}}.

\bibitem{GanderVandewalle2007_SISC}
M.~J. Gander, S.~Vandewalle, {Analysis of the Parareal Time-Parallel
  Time-Integration Method}, SIAM Journal on Scientific Computing 29~(2) (2007)
  556--578.
\newblock \href {http://dx.doi.org/10.1137/05064607X}
  {\path{doi:10.1137/05064607X}}.

\bibitem{PararealF90}
D.~Ruprecht, Parareal{F90}: Implementing {P}arareal -- {OpenMP} or {MPI}?,
  release v1.0 (2015).
\newblock \href {http://dx.doi.org/10.5281/zenodo.31288}
  {\path{doi:10.5281/zenodo.31288}}.

\bibitem{Burrage1997}
K.~Burrage, {Parallel methods for {ODE}s}, Advances in Computational
  Mathematics 7 (1997) 1--3.
\newblock \href {http://dx.doi.org/10.1023/A:1018997130884}
  {\path{doi:10.1023/A:1018997130884}}.

\bibitem{SchreiberEtAl2015}
M.~Schreiber, A.~Peddle, T.~Haut, B.~Wingate,
  \href{http://arxiv.org/abs/1506.05157}{A decentralized
  parallelization-in-time approach with parareal} (2015).
\newline\urlprefix\url{http://arxiv.org/abs/1506.05157}

\bibitem{GysiEtAl2015}
T.~Gysi, C.~Osuna, O.~Fuhrer, M.~Bianco, T.~C. Schulthess, Stella: A
  domain-specific tool for structured grid methods in weather and climate
  models, in: Proceedings of the International Conference on High Performance
  Computing, Networking, Storage and Analysis, 2015.
\newblock \href {http://dx.doi.org/10.1145/2807591.2807627}
  {\path{doi:10.1145/2807591.2807627}}.

\bibitem{ArteagaEtAl2015}
A.~Arteaga, D.~Ruprecht, R.~Krause, {A stencil-based implementation of Parareal
  in the {C++} domain specific embedded language {STELLA}}, Applied Mathematics
  and Computation 267 (2015) 727--741.
\newblock \href {http://dx.doi.org/10.1016/j.amc.2014.12.055}
  {\path{doi:10.1016/j.amc.2014.12.055}}.

\bibitem{ShuOsher1989}
C.-W. Shu, S.~Osher, Efficient implementation of essentially non-oscillatory
  shock-capturing schemes {II}, Journal of Computational Physics 83 (1989)
  32--78.
\newblock \href {http://dx.doi.org/10.1016/0021-9991(89)90222-2}
  {\path{doi:10.1016/0021-9991(89)90222-2}}.

\bibitem{Rur2013}
A.~Barry,
  \href{https://cug.org/proceedings/cug2013_proceedings/includes/files/pap103-file2.pdf}{Resource
  utilization reporting: Gathering and evaluating {HPC} system usage}, in:
  CUG2013 Proceedings, 2013.
\newline\urlprefix\url{https://cug.org/proceedings/cug2013_proceedings/includes/files/pap103-file2.pdf}

\bibitem{Rabenseifner2009}
R.~Rabenseifner, G.~Hager, G.~Jost, Hybrid {M}{P}{I}/{O}pen{M}{P} parallel
  programming on clusters of multi-core {S}{M}{P} nodes, in: 17th Euromicro
  International Conference on Parallel, Distributed and Network-based
  processing, 2009, pp. 427--436.
\newblock \href {http://dx.doi.org/10.1109/PDP.2009.43}
  {\path{doi:10.1109/PDP.2009.43}}.

\bibitem{Keyes2011}
D.~E. Keyes, Exaflop/s: The why and the how, Comptes Rendus M\'{e}canique
  339~(2–3) (2011) 70 -- 77.
\newblock \href {http://dx.doi.org/10.1016/j.crme.2010.11.002}
  {\path{doi:10.1016/j.crme.2010.11.002}}.

\bibitem{Canturk2003}
C.~Isci, M.~Martonosi,
  \href{http://dl.acm.org/citation.cfm?id=956417.956567}{Runtime power
  monitoring in high-end processors: Methodology and empirical data}, in:
  Proceedings of the 36th Annual IEEE/ACM International Symposium on
  Microarchitecture, MICRO 36, IEEE Computer Society, Washington, DC, USA,
  2003, pp. 93--.
\newline\urlprefix\url{http://dl.acm.org/citation.cfm?id=956417.956567}

\end{thebibliography}
