arXiv:1511.04203v2 [physics.plasm-ph] 25 May 2016

GPEC, A REAL-TIME CAPABLE TOKAMAK EQUILIBRIUM
CODE
M. Rampp1 , R. Preuss2 , R. Fischer2 and the ASDEX Upgrade Team2
1

Max Planck Computing and Data Facility
Gießenbachstraße 2, 85748 Garching, Germany
2

Max-Planck-Institut für Plasmaphysik,
Boltzmannstraße 2, 85748 Garching, Germany
May 26, 2016

Abstract
A new parallel equilibrium reconstruction code for tokamak plasmas is presented.
GPEC allows to compute equilibrium flux distributions sufficiently accurate to derive parameters for plasma control within 1 ms of runtime which enables real-time
applications at the ASDEX Upgrade experiment (AUG) and other machines with a
control cycle of at least this size. The underlying algorithms are based on the wellestablished offline-analysis code CLISTE, following the classical concept of iteratively
solving the Grad-Shafranov equation and feeding in diagnostic signals from the experiment. The new code adopts a hybrid parallelization scheme for computing the
equilibrium flux distribution and extends the fast, shared-memory-parallel Poisson
solver which we have described previously by a distributed computation of the individual Poisson problems corresponding to different basis functions. The code is based
entirely on open-source software components and runs on standard server hardware
and software environments. The real-time capability of GPEC is demonstrated by
performing an offline-computation of a sequence of 1000 flux distributions which are
taken from one second of operation of a typical AUG discharge and deriving the relevant control parameters with a time resolution of a millisecond. On current server
hardware the new code allows employing a grid size of 32 × 64 zones for the spatial
discretization and up to 15 basis functions. It takes into account about 90 diagnostic signals while using up to 4 equilibrium iterations and computing more than 20
plasma-control parameters, including the computationally expensive safety-factor q
on at least 4 different levels of the normalized flux.

1

I

INTRODUCTION

or other simplifying assumptions like function
parametrization 11 are adopted e.g. 12 . Lately,
three codes based on fast Grad-Shafranov
solvers have been able to demonstrate real-time
capability: P-EFIT 13 , a GPU-based variant of
the EFIT equilibrium-reconstruction code 9,7 , is
able to compute one iteration within 0.22 ms on
a 65 × 65 grid with 3 basis functions. With a
similar resolution the JANET software 14,4 has
reached runtimes of about 0.5 ms per iteration
on a few CPU cores within the LabVIEW-based
system at AUG. Using a single CPU core
overclocked at 5 GHz, the LIUQE code 15 which
is deployed at the Swiss TCV experiment,
achieves a cycle time of about 0.2 ms with
a spatial resolution of 28 × 65 points. In all
these cases, however, only a single iteration is
performed and only a very small number of basis
functions can be afforded for the parametrization of the source term of the Grad-Shafranov
equation. While such restrictions may well be
justifiable under specific circumstances 16 the
motivation for our work is to substantially push
the limits of numerical accuracy in terms of
spatial resolution, number of basis functions,
and number of iterations and thus the quality of
the equilibrium solution. This is expected to aid
the generality and robustness of the application,
in particular with respect to variations in the
plasma conditions.
In addition to the reliability and accuracy of
real-time equilibrium reconstructions, the sustainability and maintainability of the codes is
considered an important aspect. Typical implementations of analysis tools in fusion science
have a rather long life time which often exceeds
life times of commercial hardware and software
product cycles. The long usability required for
analysis tools poses special requirements on the
software chosen and the maintenance capabilities

Reconstruction of the plasma equilibrium shape
is a key requirement for the operation of current
and forthcoming tokamak experiments e.g. 1,2,3,4 .
The most commonly employed method numerically reconstructs the plasma equilibrium by
iteratively solving the two-dimensional GradShafranov equation 5,6 for the poloidal flux function (ψ), and using diagnostic signals from the
experiment as constraints 7,8 . In short, the
source term of the Grad-Shafranov equation is
expanded into a linear combination of basis functions which allows formulating a linear regression problem for the expansion coefficients, using
the actual diagnostic signals and their forwardmodelling based on the numerical solution for ψ.
The problem can be tackled numerically by solving a number of independent Grad-Shafranovtype equations — one for each basis function
of the expansion — and employing a Picarditeration procedure that allows evaluating the
source term using the solution from the previous iteration step.
A variety of equilibrium-reconstruction
codes based on this strategy 8 or similar approaches 9,7,10 exist and have been routinely used
for diagnostics and data analysis at the various
tokamak experiments. The codes EFIT 9,7 and
CLISTE 8 , for example, have been employed at
DIII-D, EAST and ASDEX Upgrade (AUG),
respectively, a class of medium-sized tokamaks
which are characterized by a control cycle
on the order of a millisecond. Until recently,
however, computing times in the millisecond
range were inaccessible to such ”first-principles”
codes, and hence their applicability for real-time
diagnostic analysis and plasma-control is very
limited, unless the numerical resolution or the
number of iteration steps is drastically reduced
2

The new code enables computing the equilibrium flux distribution and the derived diagnostics and control parameters within 1 ms of runtime, given a grid size of 32 × 64 zones with up
to 15 spline basis functions for the discretization,
and about 90 diagnostic signals. With a runtime
of less than 0.2 ms for a single equilibrium iteration and about 0.2 ms required for computing a
variety of more than 20 relevant control parameters, one can afford 4 iterations for converging
the equilibrium solution. The latter was checked
to be sufficient for reaching accuracies of better
than a percent for the relevant quantities. As a
proof of principle we shall demonstrate the real
time capability of the new code by performing
an offline analysis using data from an AUG discharge. To our knowledge the new GPEC code
to date is one of the fastest (at a given numerical
accuracy) and most accurate (at a given runtime
constraint) of its kind.

necessary to adapt to changing environments.
Open-source software together with off-the-shelf
computer hardware is thought of being perfectly suited for these demands 17 . This in particular applies to all non-standardized components like, e.g., the implementation of complex,
tailored, and evolving numerical algorithms for
equilibrium reconstruction. Commercial, closedsource software (like, e.g., compilers or numerical libraries) and hardware components may well
comply with such demands, provided that they
can be modularly interchanged using standardized interfaces.
To this end we have implemented a new,
parallel equilibrium-reconstruction code, GPEC
(Garching Parallel Equilibrium Code), which
builds on the fast, shared-memory-parallel GradShafranov solver we have developed previously 18
and a two-level hybrid parallelization scheme
which was pointed out in the same paper. The
basic numerical model and functionality of the
new code originate from the well-established and
validated algorithms implemented by the equilibrium codes CLISTE 8 and specifically its descendant IDE 19 , both of which are being used for
comprehensive offline diagnostics and data analysis at the ASDEX Upgrade experiment. GPEC
is based entirely on open-source software components, runs on standard server hardware and
uses the same code base and computer architecture as employed for performing offline analyses
with the IDE code. Such a strategy is considered a major advantage, in particular concerning verification and validation of the code 7 and
also its evolution: On the one hand, algorithmic
and functional innovations in the offline physics
modeling can gradually be taken over by the realtime version. On the other hand, high-resolution
offline analysis can directly take advantage of optimizations achieved for the real-time variant.

The paper is organized as follows: in Section II
we recall the basic equations and describe our
new, hybrid-parallel implementation of the equilibrium solver and its verification and validation.
Its computational performance and in particular real-time capability is demonstrated in Section III on an example application using real
AUG data. Section IV provides a summary and
conclusions.

II
II.A

ALGORITHMS AND IMPLEMENTATION
Equilibrium Reconstruction

Basic equations The Grad-Shafranov equation 5,6 describing ideal magneto-hydrodynamic
equilibrium in two-dimensional tokamak geome3

try reads
∆∗ ψ(r, z) = −2πµ0 r jφ

,

the right-hand side. The updated flux distribution is then given by

(1)

with cylindrical coordinates (r, z), the elliptic
differential operator
∂ 1 ∂
∂2
+ 2
∂r r ∂r ∂z

ψ=

N
X

ck ψk

,

(7)

k=1

(2) which follows from Eqs. (1,4) and the linearity of the operator ∆∗ . The unknown coeffiand the poloidal flux function (in units of Vs), cients ck are determined by experimental data:
Using the N distributions ψk obtained from
ψ(r, z). The toroidal current density profile,
Eqs.
(5,6) the so-called response matrix B :=


∂p(r, ψ) F (ψ) dF (ψ)
+
,
(3) {bl (ψk )}l=1,...,M,k=1,...N , consisting of predictions
jφ := 2π r
∂ψ
µ0 r dψ
for a set of M diagnostic signals {ml }l=1...,M is
consists of two terms, where p(r, ψ) ≡ p(ψ) is calculated. The prediction bl (ψ) is some (linear)
the plasma pressure (isotropic case) and F (ψ) = function employing the flux distribution ψ to
produce the forward-modelled signal. All mearBφ = (µ0 /2π)Ipol is proportional to the total
surements are located within the poloidal flux
poloidal current, Ipol .
grid such that the poloidal flux outside the grid
is not needed. Linear regression of the response
Numerical solution The classical stratmatrix with the measured diagnostic signals fiegy 8,20,16 for numerically solving Eq. (1) is based
nally yields the coefficients ck for Eq. (7). This
on an expansion of the current density profile jφ
procedure is iterated until certain convergence
on the right-hand-side into a linear combination
criteria are fulfilled.
of a number N = Np + NF of basis functions,
The private flux region close to the divertor
is treated differently, although the same poloidal
jφ (r, z) =


flux as inside the plasma occurs. Since there
Np
NF
X
X
1
are no measurements of the current distribution
cNp +k ϕk (ψ) .
ck πk (ψ) +
2π r
µ0 r
in the private flux region and the chosen funck=1
k=1
tional form for the current decay was proven to
(4)
be of minor importance for the equilibrium reIn each cycle of a Picard-iteration scheme, a construction, the current is chosen to decay approximately exponentially with a decay length
number of N Poisson-type equations,
of about 5 mm, starting at the last closed flux
∆∗ ψk = −4π 2 µ0 r 2 πk (ψ) (k = 1 . . . Np ), (5) surface.
∆∗ := r

,

and

The algorithm can be summarized as follows:
∗

2

∆ ψNp +k = −4π ϕk (ψ)

(k = 1 . . . NF ),

(6)

1. evaluate right-hand sides gk (r, z) :=
4π 2 µ0 r 2 πk (ψ i−1 ) and gNp +k (r, z)
:=
2
i−1
4π ϕk (ψ ) of Eqs. (5,6), using the solu-

are solved individually, where the solution ψ
from the last iteration step is used for evaluating
4

tion, ψ i−1 , from the previous iteration step measurements are used in the IDE code in the
(or initialization).
off-line mode, but are not subject of the present
work due to the lack of reliable on-line availabil2. solve the Poisson-type equations ∆∗ ψk = ity. The relatively large number of basis func−gk (r, z) individually for each k = 1, . . . , N . tions (and hence fit parameters), which is moAn additional Poisson problem arises from tivated by the need to allow for equilibria sufa convergence-stabilization procedure which ficiently flexible to address all occurring plasma
introduces the vertical plasma position as an scenarios, requires application of regularization
auxiliary free parameter 8,15 .
constraints. While in the IDE code additional
3. evaluate {bl (ψk )}l=1,...,M,k=1,...,N to con- smoothness (curvature and amplitude) requireare applied to the source profiles p′ and
struct the response matrix B.
Ad- ments
′
ditional columns {bl,N +1 }l=1,...,M and F F (which effectively adds additional columns
{bl,N +1...N +1+Next }l=1,...,M arise from the to the response matrix), GPEC currently adopts
21
abovementioned convergence-stabilization a simpler ridge-regression procedure .
procedure and from deviations in the
currents measured at a number of Next
Implementation and parallelization The
external field coils to account for wall
implementation of GPEC starts out from the seshielding and plasma-induced wall currents,
rial offline-analysis code IDE 19 , thus closely folrespectively.
lowing the concepts of the well-known CLISTE
4. perform linear regression on B · c = m to code 8 . Specifically, for each of the two sets of
basis functions, {πk }k=1...NP and {ϕk }k=1...NF ,
obtain the set of coefficients {ck }k=1,...,N .
GPEC employs a basis of cubic spline polyPN
5. evaluate new solution ψ i =
k=1 ck ψk + nomials which are defined by a number of Np
ψ ext , adding contributions from the mea- and NF knots at positions xi , respectively, with
sured and fitted deviating currents in ex- πk (xi ) = δki , ϕk (xi ) = δki and natural boundary
ternal coils, ψ ext .
conditions. Bicubic Hermite interpolation is employed for evaluating the flux distribution ψk and
6. go back to step 1.
until convergence
its gradients (which are required for computing
is reached, e.g. in terms of the maxithe magnetic field) at the positions of the magmum norm evaluated over all grid points,
netic probes and the flux loops in order to calcukψ i − ψ i−1 k∞ .
late the corresponding forward-modelled signals
The response matrix B is already prepared for for the response matrix. Thus, besides achieving
additional measurements of the motional Stark higher-order interpolation accuracy, smoothness
effect (MSE), the Faraday rotation, the pressure of the magnetic field is guaranteed by construc22
profile (electron, ion and fast particle pressure), tion .
GPEC utilizes the fast, thread-parallel Poisthe q-profile (e.g. from MHD modes), the divertor tile currents constraining F (ψ) on open son solver we have developed previously 18 . It
flux surfaces, the measurements of loop voltage is based on the two-step scheme for solving
and of iso-flux constraints 19 . These additional the Poisson equation with Dirichlet boundary
5

conditions 23 and replaces the cyclic-reduction
scheme which is traditionally employed in the
codes GEC and CLISTE 24,8 by a parallelizable
Fourier-analysis method for decoupling the linear system into tridiagonal blocks 22 . For details
on the implementation and computational performance we refer to refs. 25,18 .
As briefly sketched in ref. 18 the main idea of
the new code is to exploit an additional level of
parallelism by distributing the individual Poisson problems for determining each of the ψk
to different MPIa -processes. The algorithm, together with notes on the parallelization (relevant
MPI communication routines are noted in brackets) is summarized as follows.

5. performs linear regression on m = B · c
thread-parallelization of linear algebra
routines
6. gathers all ck′ ψk′ computed by the other
processes Pk′6=k and computes ψ =
PN
k′=1 ck′ ψk′
collective communication (MPI Allreduce)
The distribution of the Poisson problems to
different processes comes at the price of collective communication (steps 4 and 6), in which
all processes combine their data with all others, using MPI routines MPI Allgather and
MPI Allreduce, respectively.

Each process Pk=1,...,N is assigned to a single
II.B Computation of Plasma-Control
basis functionb , and process PN +1 is dedicated to
Parameters
the additional Poisson problem that arises from
the introduction of the vertical shift parameter 8 . Given the solution for the poloidal flux function
Accordingly, each process Pk :
ψ(r, z), GPEC computes more than 20 derived
1. computes gk (r, z) using the solution ψ(r, z) quantities which are relevant for plasma confrom the previous iteration step (or initial- trol. Two of them, the (r, z)-coordinates of the
magnetic axis, Rmag and Zmag can be obtained
ization)
thread-parallelization over (r, z)-grid immediately by searching the maximum of the
equilibrium flux. For the rest of the quantities
2. employs fast Poisson solver 18 to solve four major computational tasks (labelled A–D
∆∗ ψk = −gk (r, z)
in the following compilation) have to be carried
thread-parallelization 18 out. Their results allow evaluating the individ3. computes column b(ψk ) of the response ma- ual control parameters (summarized under the
bullet points) in a numerically straightforward
trix
thread-parallelization over (r, z)-grid and inexpensive way. For details on the physical definition and derivation of these quantities,
4. gathers columns b(ψk′ ) computed by the see, e.g. Ref. 27 . Implementation details are given
other processes Pk′6=k and assembles re- further below.
sponse matrix B
collective communication (MPI Allgather) (A) determination of grid points which are ena
closed by the separatrix
The ”Message Passing Interface” standard 26 .
b

In general the implementation supports an even distribution of the N + 1 Poisson problems to a number of
N ′ processes, provided N ′ divides N + 1.

• coordinates (Rgeo , Zgeo ) of the geometric axis, defined as the area6

weighted
integral within the separatrix (D) computation of the pressure distribution
R
(r, z)dS/S.
p(r, z) and the poloidal magnetic field

• the horizontal and vertical minor plasma qradius is given by
R
2
ahor
= 2
r 2 dS/S − Rgeo
and
qR
2 , respecbver = 2
z 2 dS/S − Zgeo
tively, which determines the elongation
k = bver /ahor .

• the total energy content of
R the plasma
is given by WMHD = 3/2 pdV .

• the same integration can be used
for computing theR poloidal beta parameter, βpol =
pdV /Zβ , employing the normalization constant Zβ =
2 /(2 · l2 ) (as used in
µ0 · Vplasma · Itor
pol
CLISTE) with lpol and Itor as specified
under (B) and (C), respectively.
• the plasma
R 2self-inductivity li =
−1
(2µ0 Zβ )
Bpol dV is given by an integral over the squared poloidal mag2 = (2πr)−2 · ((∂ ψ)2 +
netic field Bpol
r
(∂z ψ)2 ).
• dRXP , the difference of the poloidal
flux at the two x-points divided by
(2πRaus ) and by the poloidal magnetic field at (Raus , Zmag ): (ψXP1 −
ψXP2 )/|∇ψ|(Raus ,Zmag ) .

(B) identification of contour lines (r(t), z(t))
with ψ(r(t), z(t)) = ψc for a number of particular values of the normalized flux ψc

• from the separatrix curve (with ψc
taken as the value at the innermost
x-point which defines also the lastclosed flux surface LCFS) a number
of geometric properties can be derived
straightforwardly such as the (r, z)coordinates of the uppermost and lowermost point (in z-direction) on the
plasma surface, or Rin (Raus ) as the
R-coordinate of the innermost (outerImplementation and parallelization The
most) point on the plasma surface.
• the length of the poloidal perimeter lpol algorithms for computing the plasma-control paof the plasma is computed by integra- rameters and their basic implementation are
taken directly from the IDE code. We emphasize
tion over the LCFS.
that unlike other real-time codes 15 GPEC makes
• the safety factor q is computed by inte- no algorithmic simplifications compared with the
gration over contours defined by levels offline variant of the code. In particular, GPEC
ψc ∈ [0.25, 0.50, 0.75, 0.95].
uses the same spatial grid resolution as adopted
(C) computation
of the toroidal plasma current for equilibrium reconstruction and employs the
R
high-order interpolation schemes (cubic splines,
Itor = jφ dS (Eq. 3)
bicubic Hermite interpolation) of IDE. In the fol• the indicators for the current-center lowing we briefly summarize the main concepts
are given by the following current- adopted for the computationally most expensive
weighted integrals of
qR the current routines:
R
r 2 jφ dS/Itor ,
density Rsquad =
(A) for performing the integrations . . . dS/S,
R
Zsquad = zjφ dS/Itor .
grid points located outside the separatrix
7

are masked out according to a comparison
with the value of ψ at the x-point.

lation polynomials for intra-grid interpolation, cf. Sect. II.A).

(B) a custom contour-finding algorithm is used
which employs quadratic interpolation to
obtain the values at the intersections of
the spatial grid with the contour curve.
Quantities on the LCFS, like lpol , are determined by quadratic extrapolation from
three equidistant contour levels close to the
plasma surface.

As shall be shown in Sect. III.B the computing time is dominated by finding and integrating along contours as well as by evaluating the
pressure distribution. Accordingly, the computations are grouped into the following mutually
independent tasks which can be assigned to different MPI tasks, with only a few scalar quantities to be communicated for evaluating the final
results:
(C) for the evaluation
R of the toroidal plasma
current Itor = jφ dS the current density
• computation of the pressure distribution
jφ (r, z) is computed according to Eq. (4).
p(r, z) and the poloidal magnetic field
The summands ck πk (ψ) and cNp +k ϕk (ψ) in
principle would be available from the last
• handling of 3 contour levels for extrapolaequilibrium iteration but are distributed
tion to the last closed flux surface
among processors (cf. Sect. II.A). It turns
• handling of 4 contour levels corresponding
out that locally recomputing πk (ψ) and
to the set of safety factors q25 , q50 , q75 , q95
ϕk (ψ) by using the converged equilibrium
solution ψ (a copy of which is available on
• determination of the x-point and handling
each processor) is faster than collecting the
of 1 contour level corresponding to the sepindividual summands by MPI communicaaratrix curve
tion.
Within each MPI task, OpenMP threads are
(D) the computation of the pressure distriused, e.g. for parallelizing the computation of
bution p(r, z) is based on the identity
the pressure over the individual points of the
P Np
∂ψ p(ψ(r, z)) =
k=1 ck πk (ψ) (cf. Eq. 4). (r, z)-grid or for parallelizing over different conFor fast evaluation of p(r, z) on the poloidal tour levels.
mesh, an equidistant grid of ψl , l = 1, . . . , 12
points is constructed. For each ψl the corII.C Verification and Validation
responding pressure contribution πk (ψl ) is
computed by numerical integration starting The codes CLISTE and IDE have been extenat the plasma boundary. A cubic spline in- sively validated and verified by application to
terpolation is used to evaluate πk (ψ(r, z)) on ASDEX Upgrade data and by means of code
this grid. The gradients for computing the comparisons. Thus, only those parts of the new
poloidal magnetic field are obtained using GPEC code need to be verified which are treated
centered finite differences of the neighboring differently from the IDE code. The procedure is
grid points (consistent with the derivatives greatly facilitated by the fact that both codes
used to define the bicubic Hermite interpo- implement the same algorithms and share the
8

same code base. The only differences in the III APPLICATION PERFORGPEC implementation with respect to IDE are
MANCE
1) the prescription of a fixed number of equilibrium iterations and 2) the use of only 12 instead III.A Demonstration of Real-Time
of 30 integration points for computing the presCapability
sure distribution p(r, z) from its gradient ∂ψ p(ψ)
In order to demonstrate the real-time capabil(see item D above) on every grid point (r, z).
ity of the new code a post-processing run is perFigure 1 shows the sensitivity of a number of formed using data from a typical AUG discharge.
selected parameters to the number of equilibrium We chose AUG shot number #23221, and focus
iterations for an arbitrarily chosen time point on a time window between t1 = 3.014 s and
(t = 3.015 ms) of AUG discharge #23221. By t1000 = 4.013 s with 1000 time points, correusing the solution ψ from the last time point sponding to a time resolution of exactly 1 ms.
as an initial value, the plasma geometry (pa- As shall be shown below this simulation of 1000
rameters Rin , Raus , and the total volume of the time points can be performed on standard server
plasma, Vplasma ) is already very well determined. hardware within less than a second of computing
Within only a few iteration steps accuracies of time.
The GPEC code computes a fixed number of
better than a percent compared with the converged solution (obtained with 200 iterations) four equilibrium iteration steps for every time
are obtained also for the quantities which de- point and uses the equilibrium solution ψ(tn ) at
pend on the pressure distribution (βpol , WMHD ), time point tn as the initial value for the next
or on the toroidal current (Itor ). For the real- time point tn+1 . A standard spatial resolution
time demonstration presented in the subsequent of (Nr , Nz ) = (32, 64) grid points is employed,
section we shall fix the number of iterations to and N = 7 basis functions are used. The refour and show that the accuracy level of 1 % is sponse matrix B contains signals from 61 magmaintained over a time sequence of 1 s with 1000 netic probes and 18 flux loops and in addition
time points.
takes into account 10 external coils. Figure 3
provides an overview of the time evolution of a
Figure 2 shows that a number of 12 integration number of characteristic parameters, such as sepoints for the numerical integration of ∂ψ p(ψ) lected plasma-shape parameters, energy content
is more than sufficient in order to reach sub- and safety factor q95 (black, solid lines). For refpercent accuracies for the control parameters erence, the figure shows an extended time winβpolR and WMHD , both of which are proportional dow between t = 0.4 s and t = 4.4 s, which
to pdV (see item D above). The time for com- includes the dynamic start-up phase of the disputing p(r, z) is proportional to the number of charge with significant variations in the plasma
integration points and scales almost ideally with parameters, and also shows results from a run
the number of threads which motivates our spe- which converges the solution ψ until the maxcific choice of 12 points as a multiple of 6 threads imum changes on the grid, kψ i − ψ i−1 k∞ , are
that shall be employed for the majority of com- below 10−4 (blue, dashed lines). Compared with
putations (cf. Sect. III.B).
the converged reference run one notices good
9

Fig. 1. Sensitivity of selected parameters to the number of equilibrium iterations for computing ψ(r, z).
For each quantity the relative deviation from its converged value is shown as a function of the number
of iterations. Note that the deviations for the geometric quantities, Rin , Raus , and the total volume of
the plasma, Vplasma are multiplied by a factor of 10. Data was taken from AUG shot #23221 at time
t = 3.015 ms.

overall agreement of the ”real-time” simulation
with a fixed number of 4 equilibrium iterations
throughout the extended time interval, even during the initial start-up phase which lasts until
approx. 1.5 s. It is hence justified to confine the
analysis on the — arbitrarily chosen — time interval, [3.014 s, 4.013 s] in the ”plateau” phase

of the discharge. Focusing on this time window,
Figure 4 shows that the relative errors for most
quantities are at smaller than one percent. The
comparably large variations of Zmag are due to
its small absolute value. The absolute scatter of
Zmag is in the order of mm. Note that the cycle
time of ∆t = 1 ms is much shorter than the

10

Fig. 2. Sensitivity of βpol and of the total energy content WMHD to the number of integration points, Npsi ,
spent for computing the pressure from its gradient on every grid point (r, z). The relative deviation from
the converged value is shown as a function of the number of integration points, Npsi . Data was taken from
AUG shot #23221 at time t = 3.015 ms.

timescale of the changes of plasma conditions
(cf. the increase of WMHD , βpol at t ≈ 3.5 s).
Thus the equilibrium solution, even if not fully
converged, can easily follow such secular trends,
which is reflected by the fact that there is no
visible change in the magnitudes of the relative
error at t ≈ 3.5 s.
The employed spatial resolution of (Nr , Nz ) =
(32, 64) grid points is a common choice
for real-time analysis in medium-sized tokamaks 9,7,14,4,15 . For the chosen AUG discharge
a comparison with a run using a twofold finer
spatial grid spacing, (Nr , Nz ) = (64, 128), shows
good agreement for the majority of quantities,
with q95 and Router exhibiting the largest sensitivity to the resolution.
We conclude that by performing four iterations per time point and using (Nr , Nz ) =

(32, 64) spatial grid points, sufficiently accurate
equilibrium solutions for deriving real-time control parameters are obtained, at least in the sense
of a proof-of-principle using the chosen example
data from AUG shot #23221. A systematic analysis of the accuracy and a comprehensive validation of the code under various tokamak operational scenarios is beyond the scope of this
paper.

III.B

Computational Performance

Overview The computational performance of
GPEC is assessed on a standard compute cluster with x86 64 CPUs and InfiniBand (FDR 14)
interconnect. The individual compute nodes
are equipped with two Intel Xeon E5-2680v3
”Haswell” CPUs (2.5 GHz, 2x12 cores). We use a

11

resolution
(Nr , Nz , N+1)
(32, 64, 6)
(32, 64, 8)
(32, 64, 8)
(32, 64, 12)
(32, 64, 16)
(32, 64, 16)
(32, 64, 16)
(64, 128, 6)
(64, 128, 8)
(64, 128, 12)
(64, 128, 16)

Ttot
[ms]
0.85
1.08
0.88
0.97
1.55
1.26
0.99
2.07
1.88
2.02
2.04

Teq
[ms]
0.62
0.87
0.67
0.74
1.32
1.02
0.73
1.55
1.28
1.39
1.40

Tctl
[ms]
0.18
0.17
0.17
0.19
0.18
0.18
0.21
0.39
0.48
0.49
0.51

cores
(nodes)
24 (1)
24 (1)
48 (2)
72 (3)
24 (1)
48 (2)
96 (4)
48 (2)
48 (2)
72 (3)
96 (4)

MPI
tasks
6
4
8
12
4
8
16
6
8
12
16

threads
per task
4
6
6
6
6
6
6
8
6
6
6

TABLE I. Overview of the total runtime Ttot per time point (2nd column) for different combinations of
the numerical resolution (1st column) and compute resources. The latter are given in terms of the total
number of cores (nodes) (5th column), of MPI tasks (6th column) and of OpenMP threads per MPI task
(7th column). Ttot is the sum of the runtime Teq (3rd column) for computing the equilibrium distribution
ψ with 4 iteration steps (cf. Sect. II.A), and Tctl (4th column) for deriving plasma-control parameters (cf.
Sect. II.B). The benchmarks were performed on Xeon E5-2680v3 ”Haswell” CPUs (2.5 GHz, 24 cores per
node) and used data from AUG shot #23221.

standard Intel software stack (FORTRAN compiler v14.0, MPI v5.0) on top of the Linux operating system (SLES11 SP3). For the required
linear algebra functionality the OpenBLAS 28 library (version 0.2.13) is employed using the defacto-standard interfaces from BLAS 29,30 and
LAPACK 31,30 , and the FFTW 32 library (version 3.3.4) is used for the discrete sine transforms in the Poisson solver. Both libraries are
open-source software released under the BSD or
the GPL license, respectively. If desired, opensource alternatives to the Intel compilers (e.g.
GCC 33 ) and MPI libraries (e.g. OpenMPI 34 ,
MVAPICH 35 ) could be readily utilized, albeit
possibly with a certain performance impact.
Table I shows an overview of the total runtime per time point and of the individual contributions from computing the equilibrium solution
and from deriving the set of plasma-control parameters. We chose combinations of different nu-

merical resolutions (in terms of the spatial grid
resolution Nr × Nz and number of basis functions N ) and computing resources (in terms of
the number of CPU cores).
With two compute nodes (each with 24 CPU
cores) and using a resolution of (Nr , Nz , N+1) =
(32, 64, 8) a runtime well below 1 ms and hence
real-time capability is reached. When doubling
the number of basis functions (N +1 = 16) realtime calculation is still possible on four compute
nodes, indicating overall good weak scalability of
the code.
Roughly three quarters of the runtime is required for computing the equilibrium solution
and the rest goes into evaluating the set of
plasma control parameters. Hence, with a reduction of the number of equilibrium iterations
to only one or two 13,14,4,15 runtimes of 0.5 ms are
reached on just a single compute server.
With a higher spatial resolution of (Nr , Nz ) =

12

Fig. 3. Time evolution of selected parameters for AUG shot #23221 as computed with the new real-time
code during an extended time window between t = 0.4 s and t = 4.4 s. The time interval [3.014 s, 4.014 s]
which is chosen for performance analysis is indicated by vertical lines. The individual plots show the
coordinates (Rmag , Zmag ) and (RX , ZX ) of the magnetic axis, and of the x-point, respectively, the horizontal
and vertical minor plasma radius, ahor and bver , the plasma self-inductivity, li , the poloidal beta parameter,
13 the safety factor q95 (see Sect. II.B for the definition
βpol , the total energy content of the plasma, WMHD , and
of the quantities). The solid, black lines correspond to a ”real-time” application which is restricted to 4
equilibrium iterations. For comparison, values derived from converged equilibrium solutions (convergence
criterion kψ i − ψ i−1 k∞ < 10−4 ) are shown as dashed, blue lines. A standard spatial resolution of (32, 64)
grid points and 7 basis functions were employed.

Fig. 4. Evolution of the relative deviation (denoted by square brackets) of the quantities which were
computed after 4 equilibrium iterations, from the values computed with converged equilibrium solutions
(convergence criterion kψ i − ψ i−1 k∞ < 10−4 ) during the time interval chosen for performance analysis (as
indicated by vertical lines in Fig. 3). A standard spatial resolution of (32, 64) grid points and 7 basis functions
were employed.

14

(64, 128) a cycle time of 2 ms is possible with 4
iterations, although from the linear algorithmic
complexity of the algorithm 18 , and the amount
of data which is communicated between the processes one would expect the runtime to increase
by a factor of four compared to (Nr , Nz ) =
(32, 64). The reason is a higher parallel efficiency of the Poisson solver which is limited by
OpenMP overhead at lower resolutions 18 as well
as lower communication overhead due to larger
message sizes in the expensive MPI Allreduce
operation which collects the distributed
P fields
ψk (r, z) and computes the sum ψ = N
k=1 ck ψk
(see the rows labelled ”Poisson solver” and ”ψ
sum” in Table II).

tions would be somewhat faster in this case the
resultant doubling of the effective runtime for the
Poisson solver cannot be compensated for.
Plasma-control parameters For the example of (Nr , Nz , N + 1) = (32, 64, 8) computed
with 48 CPU cores, Figure 5 shows how the runtime of Tctl = 0.17 ms (cf. Table I) is composed and illustrates the benefits of exploiting
task parallelism. Due to the final MPI Allreduce
operation in the equilibrium construction algorithm, which comes at little extra effective runtime cost as compared to a single MPI Reduce,
all MPI tasks hold a copy of the equilibrium
flux distribution ψ(r, z). Hence, independent
parts of the time-consuming computations for
the plasma-control parameters like, e.g. evaluating the pressure distribution or determining the
contours for the safety factor q, can be handled
in parallel by different MPI tasks and only scalar
quantities need to be communicated (within a
few microseconds) between the MPI tasks at the
end. If necessary, the remaining load-imbalance
indicated by Figure 5, namely task 0 and task
3 becoming idle after about 0.1 ms, and tasks
4, . . . , 7 not computing anything at all at this
stage, could be exploited for computing additional quantities and/or for further reducing the
effective runtime. Without task-parallelization,
on the contrary, the runtimes would add up to
more than 0.4 ms despite the thorough OpenMP
parallelization (e.g. across different contour levels) within the task.

Equilibrium reconstruction Table II further shows that the runtime for the equilibrium
iterations is dominated
P by computing the distributed sum, ψ = N
k=1 ck ψk , which is also an
implicit synchronization point for all processes,
and by the Poisson solver. Due to the use of
high-order interpolation schemes (cf. Sect.II.A)
a significant fraction of about 25% is spent on the
evaluation of the right-hand sides of Eqs. (5,6)
(”RHS update”) and of the response matrix B
(”RM evaluation”).
In general, the chosen products of number
of MPI tasks times OpenMP threads per task
turned out as the most efficient ones for this
hardware platform with 24 cores per node: using
more than 6 threads (or 8 threads in the case of
(Nr , Nz ) = (64, 128)) results only in a modest
speedup of the Poisson solver 18 but doubles the
required number of cores and network resources. IV CONCLUSIONS
Keeping those constant and reducing the number of MPI tasks instead by a factor of two would With the motivation of achieving subrequire each MPI task to handle two basis func- millisecond runtimes a new, parallel equilibriumtions. Although the collective MPI communica- reconstruction code, GPEC, was presented
15

resolution
(Nr , Nz , N+1)
cores (MPI tasks
· threads/task)
RHS update
Poisson solver
RM evaluation
RM gather∗
lin. regression
ψ sum∗∗
other
total

(32, 64, 8)

(32, 64, 16)

(64, 128, 8)

(64, 128, 16)

48(8 · 6)

96(16 · 6)

48(8 · 6)

96(16 · 6)

T [ms]
0.022
0.036
0.024
0.015
0.013
0.048
0.008
0.167

%
13.2
21.6
14.5
8.9
7.9
28.9
5.0
100.0

T [ms]
0.021
0.036
0.024
0.021
0.017
0.056
0.008
0.182

%
11.7
19.5
13.3
11.4
9.5
30.5
4.2
100.0

T [ms]
0.050
0.085
0.025
0.019
0.014
0.110
0.011
0.319

%
15.7
26.6
7.8
5.9
4.4
36.2
3.4
100.0

T [ms]
0.053
0.086
0.025
0.029
0.018
0.128
0.014
0.353

%
15.0
24.5
7.1
8.1
5.2
36.2
3.9
100.0

TABLE II. Runtime breakdown of a single Picard equilibrium-iteration step (average over 4 iteration
steps times 1000 time points) for different combinations (subset of Table I) of numerical resolution and
computational resources. The individual algorithmic steps from top to bottom correspond to the sequence
in Sect. II.A. Routines labelled with asterisks are dominated by MPI communication (∗ MPI Allgather,
∗∗
MPI Allreduce).

which is suitable for real-time applications in
medium-sized tokamaks like ASDEX Upgrade
(AUG). GPEC implements the classical concept of iteratively solving the Grad-Shafranov
equation and feeding in diagnostic signals
from the experiment 7,8 . Specifically, GPEC is
implemented as a variant of the IDE code 19 ,
which is a descendant of CLISTE 8 . Compared
with these well-established and validated offlineanalysis codes no algorithmic simplifications
are necessary for achieving the desired cycle
times of less than a millisecond, besides limiting
the number of equilibrium iterations to four.
In addition to real-time applications the new
code enables fast and highly accurate offline
analyses soon after the tokamak discharge:
Here, tolerable runtimes are in the range of
0.1–1 s per cycle which allows computing fully
converged equilibria employing highly resolved
spatial and basis functions grids.

a fast shared-memory-parallel Grad-Shafranov
solver 18 together with the MPI-distributed solution of the individual Poisson-type problems and
a thorough parallelization of the post-processing
algorithms for computing the relevant plasmacontrol parameters from the equilibrium flux distribution.
Using data from a typical AUG discharge the
real-time capability of the new code was demonstrated by the offline computation of a sequence
of 1000 time points within less than a second of
runtime. The relative accuracy was ascertained
by comparing the relevant plasma parameters
with a converged run. By allowing four iterations for computing the equilibrium solution the
majority of control parameters can be computed
with an accuracy of a percent or better.

The adopted two-level, hybrid parallelization
scheme allow efficient utilization available compute resources (in terms of the numbers of
The parallelization of the new code builds on nodes, of CPU sockets per node and of cores
16

Fig. 5. Runtime of the major groups (A-D) of computations for the plasma-control parameters (cf.
Sect. II.B) and their assignment to four MPI tasks which execute in parallel. Within each task the individual parts of the algorithm (indicated by different colours) execute in the order from bottom to top.
Accordingly, the dashed horizontal line marks the effective runtime in the application. The data corresponds
to AUG shot #23221 computed with a resolution of (Nr , Nz , N +1) = (32, 64, 8) using in total 8 MPI tasks,
each with 6 OpenMP threads (cf. Table I).

per socket) for a given numerical resolution (in
terms of the spatial grid and number of basis
functions). Moreover, foreseeable advances in
computer technology, in particular the increasing number of CPU cores per compute node, are
expected to push the limits of real-time applications with GPEC towards even higher numerical
accuracies (in terms of affordable resolution and
number of equilibrium iterations). We note that
the benchmarks figures reported in this work can
be considered rather conservative as the computations were performed on a standard com-

pute cluster which comprises hundreds of nodes.
For energy-budget reasons such clusters are typically configured not with CPUs of the highest
clock-frequency. In our case, for example, CPUs
with 2.5 GHz and only 24 cores per node were
available. For deployment in the control system of a tokamak experiment such as AUG, by
contrast, we envision dedicated server hardware
with higher clock frequencies and at least four
CPU sockets, both of which is expected to further boost the computational performance of the
real-time application. Thus it should be rather

17

straightforward to save enough time for the communication of the code with the control system,
which, depending of the specifics of the system,
requires another few hundred microseconds per
cycle.
Finally, it is worth mentioning that GPEC
is based entirely on open-source software components, relies on established industry (FORTRAN, MPI, OpenMP) or de-facto software
standards (BLAS, LAPACK, FFTW), runs on
standard server hardware and software environments and hence can be released to the community and utilized without legal or commercial
restrictions.

ACKNOWLEDGEMENT
We acknowledge stimulating discussions with K.
Lackner and W. Treutterer. Thanks to P. Martin who has provided a subversion of the CLISTE
code with which this study started and to S. Gori
who has developed an initial FORTRAN 90version.
This work has been carried out within the framework of the EUROfusion Consortium and has
received funding from the Euratom research
and training programme 2014-2018 under grant
agreement No 633053. The views and opinions
expressed herein do not necessarily reflect those
of the European Commission.

18

REFERENCES

[1] C. GORMEZANO, A. SIPS, T. LUCE,
S. IDE, A. BECOULET, X. LITAUDON,
A. ISAYAMA, J. HOBIRK, M. WADE,
T. OIKAWA, R. PRATER, A. ZVONKOV,
B. LLOYD, T. SUZUKI, E. BARBATO, P. BONOLI, C. PHILLIPS,
V. VDOVIN, E. JOFFRIN, T. CASPER,
J. FERRON, D. MAZON, D. MOREAU,
R. BUNDY, C. KESSEL, A. FUKUYAMA,
N. HAYASHI, F. IMBEAUX, M. MURAKAMI, A. POLEVOI, & H. S. JOHN.
Chapter 6: Steady state operation. Nuclear
Fusion, 47(6), S285 (2007).
[2] Y. GRIBOV, D. HUMPHREYS, K. KAJIWARA, E. LAZARUS, J. LISTER,
T. OZEKI, A. PORTONE, M. SHIMADA,
A. SIPS, & J. WESLEY.
Chapter 8:
Plasma operation and control. Nuclear Fusion, 47(6), S385 (2007).
[3] M. MARASCHEK, G. GANTENBEIN,
T. P. GOODMAN, S. GÜNTER, D. F.
HOWELL, F. LEUTERER, A. MÜCK,
O. SAUTER, & H. ZOHM. Active control
of MHD instabilities by ECCD in ASDEX
upgrade. Nuclear Fusion, 45, 1369 (2005).
[4] L.
GIANNONE,
M.
REICH,
M.
MARASCHEK,
C.
RAPSON,
W. TREUTTERER, L. BARRERA,
E. POLI, A. BOCK, G. CONWAY,
F. FISCHER, J. FUCHS, K. LACKNER, P. MCCARTHY, A. MLYNEK,
19

R. PREUSS, K.-H. SCHUHBECK, J. STOBER, M. RAMPP, R. MCDERMOTT,
Q. RUAN, L. WENZEL, & ASDEX
UPGRADE TEAM. A data acquisition
system for real-time magnetic equilibrium
reconstruction on ASDEX Upgrade and its
application to NTM stabilization experiments. Fusion Engineering and Design, 88,
3299 (2013).
[5] V. D. SHAFRANOV. About MHD equilibrium configurations. Zh. Eksp. Teor. Fiz.,
33, 710 (1957).
[6] R. LÜST & A. SCHLÜTER. Axisymmetrische magnetohydrodynamische Gleichgewichtskonfigurationen.
Z. Naturforsch., 129, 850 (1957).
[7] L. L. LAO, J. R. FERRON, R. J. GROEBNER, W. HOWL, H. ST. JOHN, E. J.
STRAIT, & T. S. TAYLOR. Equilibrium
analysis of current profiles in tokamaks. Nuclear Fusion, 30, 1035 (1990).
[8] P. J. MCCARTHY, P. MARTIN, &
W. SCHNEIDER. The CLISTE interpretive equilibrium code.
Tech. Rep.
IPP 5/85, Max-Planck-Institut für Plasmaphysik (1999).
[9] L. LAO, H. S. JOHN, R. STAMBAUGH,
A. KELLMAN, & W. PFEIFFER. Reconstruction of current profile parameters and

plasma shapes in tokamaks. Nuclear Fusion, [14] A. BARP, M. CERNA, S. CONCEZZI,
25(11), 1611 (1985).
L. GIANNONE, G. MORROW, Q. RUAN,
A. VEERAMANI, & L. WENZEL. A
[10] J. BLUM, C. BOULBE, & B. FAUGERAS.
real-time Grad-Shafranov PDE solver and
Reconstruction of the equilibrium of the
MIMO controller. Fusion Engineering and
plasma in a tokamak and identification of
Design, 87(12), 2112 (2012). Proceedings
the current density profile in real time.
of the 8th IAEA Technical Meeting on ConJournal of Computational Physics, 231(3),
trol, Data Acquisition, and Remote Partic960 (2012).
ipation for Fusion Research.
[11] B. BRAAMS, W. JILGE, & K. LACKNER. [15] J.-M. MORET, B. DUVAL, H. LE,
Fast determination of plasma parameters
S. CODA, F. FELICI, & H. REIMERDES.
through function parametrization. Nuclear
Tokamak equilibrium reconstruction code
Fusion, 26(6), 699 (1986).
LIUQE and its real time implementation.
Fusion Engineering and Design, 91, 1
[12] L. GIANNONE, M. CERNA, R. COLE,
(2015).
M. FITZEK, A. KALLENBACH, K. LDDECKE, P. MCCARTHY, A. SCARA- [16] J. FERRON, M. WALKER, L. LAO, H. S.
BOSIO, W. SCHNEIDER, A. SIPS,
JOHN, D. HUMPHREYS, & J. LEUER.
W. TREUTTERER, A. VRANCIC,
Real time equilibrium reconstruction for
L. WENZEL, H. YI, K. BEHLER,
tokamak discharge control. Nuclear Fusion,
T. EICH, H. EIXENBERGER, J. FUCHS,
38, 1055 (1998).
G. HAAS, G. LEXA, M. MARQUARDT,
A. MLYNEK, G. NEU, G. RAUPP, M. RE- [17] I. YONEKAWA, A. V. FERNANDEZ, J.ICH, J. SACHTLEBEN, K. SCHUHBECK,
M. FOURNERON, J.-Y. JOURNEAUX,
T. ZEHETBAUER, S. CONCEZZI, T. DEW.-D. KLOTZ, A. WALLANDER, & COBELLE, B. MARKER, M. MUNROE,
DAC TEAM. ITER instrumentation and
N. PETERSEN, & D. SCHMIDT. Data
control system towards long pulse operaacquisition and real-time signal processing
tion. Plasma Fusion Research, 7, 2505047
of plasma diagnostics on ASDEX upgrade
(2012).
using labview RT. Fusion Engineering and
Design, 85(34), 303 (2010). Proceedings [18] M. RAMPP, R. PREUSS, R. FISCHER,
K. HALLATSCHEK, & L. GIANNONE. A
of the 7th IAEA Technical Meeting on
parallel Grad-Shafranov solver for real-time
Control, Data Acquisition, and Remote
control of tokamak plasmas. Fusion Science
Participation for Fusion Research.
and Technology, 62(3), 409 (2012).
[13] X. N. YUE, B. J. XIAO, Z. P. LUO, &
Y. GUO. Fast equilibrium reconstruction [19] R. FISCHER, J. HOBIRK, L. B. ORTE,
A. BOCK, A. BURCKHART, I. CLASSEN,
for tokamak discharge control based on
M. DUNNE, J. FUCHS, L. GIANGPU. Plasma Physics and Controlled FuNONE, K. LACKNER, P. MCCARTHY,
sion, 55(8), 085016 (2013).
20

R. PREUSS, M. RAMPP, S. RATHGEBER, M. REICH, B. SIEGLIN, W. SUTTROP, E. WOLFRUM, & ASDEX
UPGRADE TEAM. Magnetic equilibrium
reconstruction using geometric information
from temperature measurements at asdex
upgrade. In V. NAULIN, C. ANGIONI,
M. BORGHESI, S. RATYNSKAIA,
S. POEDTS, & T. DONN (eds.), 40th EPS
Conference on Plasma Physics, vol. 37D
of Europhysics conference abstracts, page
P2.139. European Physical Society, Geneva
(2013).

[25] R. PREUSS, R. FISCHER, M. RAMPP,
K. HALLATSCHEK, U. VON TOUSSAINT, L. GIANNONE, & P. MCCARTHY. Parallel equilibrium algorithm
for real-time control of tokamak plasmas.
Tech. Rep. IPP R/47, Max-Planck-Institut
für Plasmaphysik (2012).
[26] http://www.mpi-forum.org/.
[27] W. M. STACEY. Fusion Plasma Physics.
Wiley, 2 edn. (2012).
[28] http://www.openblas.net.

[20] P. J. MCCARTHY & ASDEX UPGRADE [29] BLAST FORUM. Basic linear algebra subprograms technical forum standard. InTEAM. Identification of edge-localized moternational Journal of High Performance
ments of the current density profile in a
Applications and Supercomputing, 16(1)
tokamak equilibrium from external mag(2002).
netic measurements. Plasma Physics and
Controlled Fusion, 54(1), 015010 (2012).
[30] http://www.netlib.org.
[21] A. E. HOERL & R. W. KENNARD. [31] E. ANDERSON, Z. BAI, C. BISCHOF,
Ridge regression: Biased estimation for
S. BLACKFORD, J. DEMMEL, J. DONnonorthogonal problems. Technometrics,
GARRA, J. DU CROZ, A. GREENBAUM,
12(1), 55 (1970).
S. HAMMARLING, A. MCKENNEY, &
D. SORENSEN. LAPACK Users’ Guide.
[22] W. H. PRESS, S. A. TEUKOLSKY, W. T.
Society for Industrial and Applied MatheVETTERLING, & B. P. FLANNERY. Numatics, Philadelphia, PA, third edn. (1999).
merical Recipes: The Art of Scientific ComISBN 0-89871-447-8 (paperback).
puting. Cambridge University Press, 3rd
edn. (2007).
[32] M. FRIGO & S. G. JOHNSON. The design
[23] K. V. HAGENOW & K. LACKNER. Computation of axisymmetric MHD equilibria.
In Proc. 7th Conf. on the Numerical Simulation of Plasmas, page 140. New York, US
(1975).
[33]

and implementation of FFTW3. Proceedings of the IEEE, 93(2), 216 (2005). Special
issue on “Program Generation, Optimization, and Platform Adaptation”.
https://gcc.gnu.org/.

[24] K. LACKNER. Computation of ideal MHD [34] http://www.open-mpi.de/.
equilibria. Comp. Phys. Comm., 12, 33
[35] http://mvapich.cse.ohio-state.edu/.
(1976).
21

