J. Plasma Fusion Res. SERIES, Vol. 9 (2010)

A new framework for integrated simulation model using
MPMD approach
A.Takayama, K.Shimizu1) , Y.Tomita, and T.Takizuka1)
National Institute for Fusion Science, Toki,Gifu, Japan1)
Japan Atomic Energy Agency, Naka, Ibaraki, Japan(Received: 18 November 2009 / Accepted: 5 April 2010)

We introduce a new eﬀective framework to combine computer codes into integrated simulation
model. The framework employs MPMD (Multiple Program Multiple Data) approach and each computercode is written with MPI (Message Passing Interface) library. Adopting MPMD approach makes eachcomputer code independent, which leads to the capability for maintenance and improvement of theintegrated simulation model. The validity and usability of this approach is shown through a simple
example model, which simulates the integrated divertor code, SONIC .
Keywords: MPMD (Multiple Program Multiple Data), MPI (Message Passing Interface), integrated
simulation

1. Introduction

variables for j-th term of sequence, and identiﬁer of(j)functions, respectively. Sets of variables (xi ) can be(j−1)
determined by sets of the previous outputs (yi),(i ∈ [1, N ]) for well-posed problems. In the single program approach all functions fi are implemented in asingle program and all variables xj are stored into thesingle program. On the other hand, our MPMD approach separates functions fi (xi , pi ), (i ∈ [1, N ]) into
independent program and uses a master program ormaster programs in order to exchange and transform
variables appropriately.

Integrated simulations, such as multiscale simulations, multiple-physics simulations, and so on, havecome up and been performed, recently [1, 2, 3, 4].
In such simulations, at least several computer codes
based on diﬀerent physics or diﬀerent spatio-temporalscales have to be combined with each other. Establishing a simple and useful way to combine computercodes is signiﬁcant for promoting integrated simulations. Though it is possible to combine computercodes simply by restructuring and bundling together
into a single program, such a procedure diminishes thecapability for maintenance and improvement of theintegrated simulation code. The simple and straightforward approach of unifying program is usable forproving validity of integrated simulation. It, however,is not suitable for continuing model development, especially not for group or team development.
Here we introduce a new eﬀective framework tocombine computer codes into integrated simulation
model. The framework employs MPMD (MultipleProgram Multiple Data) approach and each computercode is written with MPI (Message Passing Interface)library [5], [6]. Adopting MPMD approach makes eachcomputer code independent, or at least not strongly
coupled, which leads to the capability for maintenanceand improvement of the integrated simulation model.
Using MPI enables to execute each code on a computer system suited to it, which might enhance theperformance of the integrated simulation execution.

Grouping
In our MPMD approach simulation codes are separated into several programs. MPI library is employed for communication between programs. Adequate grouping of programs is prefereble for using
MPI eﬀectively. Figure 1 is a schematic of grouping
method. In this case there are 14 processes, 4 program groups plus master group. The master group
governs integrated simulation and communicates witheach other program groups. The process with MPIrank 0 in master group is called as ‘grand master’process. First of all, we group programs according
to an identiﬁer which consists of a character string.
Communication between program groups is performed
between MPI rank 0 process in each program groups.
For the purpose of this inter-group communication, weconstruct the other MPI group, which we call ‘data exchange group’.
Construction of General MPI Datatype
Integrated simulations have many data to be exchanged. In order to exchange such data eﬀectively
with MPI library we should make data blocks according to timing of exchanging the data and construct
a General MPI Datatype for each data block. Using

2. MPMD approach
Each computer code in an integrated simulation can be regarded as a function or mapping fi :(j)
(j)
(j)
(j)
(xi , pi ) �→ yi , where xi , pi , yi , and i ∈ [1, N ]are set of input variables, input parameters, output
author’s e-mail: takayama.arimichi@nifs.ac.jp604©2010 by The Japan Society of Plasma
Science and Nuclear Fusion Research

A. Takayama et al., A New Framework for Integrated Simulation Model Using MPMD Approacha integer variables is deﬁned. 3 real variables are forthe current time, the next time interval, and the time
step for calculation, and the integer variable is used
for a ﬂag. Some bits of the integer variable expressedin binary form have the meaning:
TERMINATE indicates termination of the execution.
WAIT means a request for some wait.
FASTMODE indicates that all process are in the ﬁnetime resolution calculation mode described below.
CALCFAST indicates whether the time step control iscarried out.
Fig. 1 Schematic view of grouping processes/programs.
There are 14 processes, which are grouped into 4sub (program) groups and a master group. Intergroup communication is performed within ‘data exchange group’, which consists of process with MPI
rank 0 within each local group.

Control of execution
We introduce two types of execution control. One
is a control of time step. If the control of time stepis on, when an abrupt change of physical quantity or
condition is detected, we rewind time steps and reexecute with ﬁner time step during several time steps.
The other is a sequence control of program execution. Each program group depends on other programgroup(s). If all processes are executed in parallel, thedependency could have been violated. When the sequence control of program execution is switched on,
the grand master process controls the order of program execution. This control, however, could makeprocesses wait for a while and degrades the eﬃciencyof computation.

a routine for easy construction of datatype oﬀered in
our MPMD library helps this step. The usage is thefollowings:
call initMakeTypecall addVariable( time )call addVariable( pwi )call addVariable( pwe )call termMakeType( newtype )In this case a general MPI datatype newtype containing 3 variables of time, pwi, and pwe is constructed.

Structure of Grand Master process
The aim of grand master is governing integrated
simulation and exchanging data. At ﬁrst initialization routines for using MPMD approach is called. In
this step, grouping of programs and construction of
exchange-data structures are carried out. Then thegrand master process receives initial condition datafrom the local master process(es) of the other program group(s). The structure of the main calculation
loop is as follows:

Condition message
In order to organize multiple programs, each program should notify the its own condition or status tothe grand master process at an appropriate timing.
For this purpose, we deﬁne some condition messagesof named integer constant. They are used as the tagof MPI communication.
READY means that the process is ready for calculation
and requests the condition of the next time interval.
DATA INQUIRE means that the process requests input
data for execution.
FAST MODE means that the process has detected anabrupt change of physical quantity or condition
during the calculation.
CALC DONE means that the process has completed thecalculation for the requested time interval andsends results.

mainLoop: docall MPI_Probe(MPI_ANY_SOURCE&TAG..)select case( istat(MPI_TAG) )case( READY )call MPI_Recv(..READY..)! termination condition is satisfied?
call MPI_Send(..READY..)! all processes are being terminated?
! if not, cycle mainLoop
exit mainLoopcase( DATA_INQUIRE )call MPI_Recv(kDataset...)select case( kDataset )case( Prog#1 )call MPI_Send(time,1,prog#1_Itype..)...
end selectOther condition messages can be implemented into
this framework if necessary.
General MPI Datatype for time related information
In order to exchange the condition of the nexttime interval a General MPI Datatype of 3 real and605

A. Takayama et al., A New Framework for Integrated Simulation Model Using MPMD Approachterium neutrals (D and D2 ) with Monte Carlo (MC)method, and the IMPMC code which solves kineticequations for impurity neutrals and ions with MC
method. Absorbed power of NBI heating is determined by OFMC code. The transports of plasma,neutral and impurity are solved iteratively.
The simple model equations simulating theSONIC+OFMC code are introduced as follows. They
consist of ﬁve parts, namely soldor describing the timeevolution of total particle number and stored energy,
neut2d describing the particle source by global particle conﬁnement, impmc for the total radiation power,
and ofmc for the absorbed power of NBI heating.

case( FAST_MODE )call MPI_Recv(..FAST_MODE..)case( CALC_DONE )call MPI_Recv(kDataset..)select case( kDataset )case( Prog#1 )call MPI_Recv(time,1,prog#1_Otype..)! store received data...
end select! all processes have completed?
! if not, cycle mainLoop
call MPI_Barrier(dewld,ierr)! notify whether abrupt change occured
end selectenddo mainLoop

Model of soldor
The number of plasma particles N and the energy
W can be modeled using the following equations:

Structure of Local Master process
Each local master process communicates with thegrand master process and exchanges data. In local
master process initialization routines for using MPMD
approach should be called, too. Then it sends initial condition data to the grand master process. The
structure of the main calculation loop is as follows:

NdN=−+ SN ,dt
τN

(1)

WdW
=−+ POH + PNBI − Prad .
dt
τW(2)

Plasma conﬁnement property described by τN andτW are assumed to bedocall MPI_Send(..READY..)
call MPI_Recv(time,1,timType..)if(is_msg(t_msg,WAIT)) cycle ! wait
call MPI_Send(kProg,..,DATA_INQUIRE..)
call MPI_Recv(time,1,Prog_Itype..)call MPI_Bcast ! broadcast to my groupif(is_msg(t_msg,TERMINATE)) exit
do i=1,iter
call pls_exec( kcnd )! abrupt change is detected? -> exit
enddo
call MPI_Send(kProg,..,CALC_DONE..)
call MPI_Send(time,1,Prog_Otype..)call MPI_Barrier(dewld,ierr)call MPI_Probe(MPI_ANY_SOURCE&TAG)
if( istat(MPI_TAG)==FAST_MODE )thencall MPI_Recv(..FAST_MODE..)! rewind data
endif
enddo

τN = 0.7τW ,0.85[s]τW = 0.30[s]
0.92[s](3)(OH phase)(L-mode phase)(4)(H-mode phase)

L→H and H→L mode transition conditions are taken
to be Te = 0.8[keV] and Te = 0.6[keV] at x = 0.95,
respectively. Here the temperature Te is assumed tohave a given proﬁle and the value is determined by
the energy of plasma W and the number of plasma
particles N .
Model of neut2dParticle source term in equation (1) is given by thesummation of particle ﬂux due to recycling Γrecy , gas
puﬃng Γpuﬀ and neutral beam injection (NBI) ΓNBI .
These three ﬂuxes are given by

3. Example problem

Γrecy = 0.6

In order to verify the validity of the MPMD approach introduced in section 2, we make a simple example problem, which simulates the SONIC code [1]combined with the OFMC code. The SONIC code isa simulation code package designed for analyzing edge
region of 2-dimensional magnetic conﬁnement systems
and consists of the SOLDOR code which solves 2dimensional ﬂuid equations for ion and electron, theNEUT2D code which solves kinetic equations for deu-N,τN
t < 0.25
0.021
Γpuﬀ = 0.75 × 100.25 ≤ t < 3.8
density control 3.8 ≤ t,
ΓNBI =606

PNBI.
80keV

(5)

(6)

(7)

A. Takayama et al., A New Framework for Integrated Simulation Model Using MPMD ApproachFor t ≥ 3.8, Γpuﬀ is adjusted to ﬁx the plasma densityof ne0 = 6.0 × 1019 .
Model of ofmc
inj
Input power due to NBI PNBI
and absorbed power
into plasma PNBI are modeled by
t < 0.80
0.0[MW]
inj(8)
PNBI = 16.0[MW] 0.80 ≤ t < 4.5
10.0[MW] 4.5 ≤ t,
inj,
PNBI = ηPNBI(9)

where absorption rate η, which is a function of ne0determined by N , is given byη={0.4 + (ne0 /1019 − 2.0) · 0.2
1.0

ne0 < 4 × 1019

ne0 ≥ 4 × 1019 .
(10)Model of impmc
Radiation loss term Prad in equation (2) is modeled by the summation of radiation loss due to carbon
and argon impurity;
Prad C = 0.2(POH + PNBI ),(Te0
Prad Ar = PAr (t) 1 + 0.5
10(11)).

(12)

Here a level of the radiation loss power due to argon
is given by{0t < 2.3
PAr (t) =(13)1[MW] t ≥ 2.3.

Fig. 2 Solution of example problem by use of the MPMD
approach with the control of time step and the sequence control of program execution. These ﬁgureshows plasma density, temperature, ﬂuxes, input
and absorbed power, and radiation loss obtained bysoldor, neut2d, ofmc, impmc process, respectively.

4. Solution of example problem with
MPMD approach
We solve the example problem given by equations(1)-(13) by use of the MPMD approach described insection 2. The solution is shown in ﬁgure 2. In thissolution, seven MPI processes are used. Each one process is for master, soldor, neut2d and ofmc, and three
processes are for impmc.
The three impmc processes are for local-master,radiation loss due to carbon and radiation loss due
to argon. The local-master of impmc communicateswith grand-master process, distributes preconditions
to and gathers solutions from two processes used forcarbon and argon.
The soldor process receives N , W , τN and τW ofsoldor previous results, SN determined by the particle
ﬂux from neut2d, PNBI of ofmc, and Prad of impmc,solves eqs.(1-4), and returns N , W , τN and τW . POH
is given as an input parameter.

The neut2d process receives N , τN of soldor andPNBI of ofmc, solves eqs.(5-7), and returns Γrecy , Γpuﬀ ,ΓNBI .
The ofmc process receives ne0 of soldor, solves
eqs.(8-10), and returns PNBI . Input power due to NBI
described by eq.(9) is given as a parameter.
The impmc group receives POH , PNBI , and Te0 ofsoldor, solves eqs.(11-13), and returns Prad .
The time interval of gathering all processes ∆T is0.1[s] for usual case. The time step for each calculation ∆t is 0.02[s] for soldor, and 0.1[s] for the others,which equals to the gathering time interval. Whenabrupt change of physical quantity or condition, ﬁne607

A. Takayama et al., A New Framework for Integrated Simulation Model Using MPMD Approachtime resolution iteration with ∆T = ∆t = 0.01[s] iscarried out ﬁve times. After the ﬁne time resolution
calculation, ∆T and ∆t are turned back.
At t = 0.800 NBI begins as equation (9). ofmc
detects the change and all processes move to the ﬁner
inj
causes jumps oftime step mode. The jump of PNBIPNBI and ΓNBI , which leads to the increase of the temperature Te0 . After the ﬁne time resolution calculation of ﬁve times iteration, all processes return to thenormal time step mode.
At t ∼ 1.40 Tedge exceeds 0.8keV and the system
undergoes a transition to H-mode, which is detectedby soldor and all processes move to the ﬁner timestep mode. Using the ﬁner time step, Tedge is foundto be lower than 0.8keV, which means that the transition to H-mode in the normal time step calculation ismis-predicted one. During ﬁner time step calculation,Tedge exceeds 0.8keV at t = 1.417 and the system undergoes a transition to H-mode. This causes the dropof Γrecy , and then density Ne0 begins to increase.
At t = 2.320 argon injection begins, which is detected by impmc and all processes move to the ﬁnertime step mode. Prad Ar jumps, which degrades tempereture increase. After the ﬁne time resolution calculation of ﬁve times iteration, all processes return tothe normal time step mode.
At t = 3.225, Ne0 > 4.0 × 1019 and this causes
the jump of η. soldor detects this rapid change. Att = 3.830 neut2d detects the beginning of the density
control and jump of Γpuﬀ . At t = 4.535 ofmc detects
inj. At the three timings, all processes
the change of PNBI
move to the ﬁner time step mode once again. After theﬁne time resolution calculation of ﬁve times iteration,
all processes return to the normal time step mode.

is required for the MPMD approach introduced in thisarticle for the accuracy and eﬃciency of calculation.

Acknowledgment
The authors acknowledge Prof. A. Fukuyama at
Kyoto University and Prof. M. Yagi at Kyushu University for intensive discussions and comments.
[1] K. Shimizu, T.Takizuka, K. Ohya, K. Inai, T. Nakano,
A. Takayama, H. Kawashima and K. Hoshino, Nucl.
Fusion 49 (2009) 065028.
[2] A. Fukuyama, S. Murakami, M. Honda, Y. Izumi,M.Yagi, N. Nakajima, Y. Nakamura and T. Ozeki,
Proc. of 20th IAEA Fusion Energy Conf. (Villamoura,Portugal, 2004) IAEA-CSP-25/CD/TH/P2-3.
[3] M. Sato, S. Toda, Y. Nakamura, K.-Y. Watanabe, A.
Fukuyama, S. Murakami, M. Yokoyama, H.Funaba,
H. Yamada and N. Nakajima, Plasma Fusion Res. 3(2008) S1063.
[4] B. Guillerminet, M. Airaj, P. Huynh, G. Huysmans,
F. Iannone, F. Imbeaux, J. B. Lister, G. Manduchi,
P. Strand and the contributors to the European Task
Force on Integrated Modelling Activity, Fusion Eng.
Des. 83 (2008) 442.
[5] Y. Aoyama, ‘Introductory Course in Parallel Programming (MPI version)’ (in Japanese), available fromhttp://accc.riken.jp/hpc/training/text.html[6] P.S. Pacheco, ‘Parallel Programming with MPI’, (Morgan Kaufmann Publishers, 1997); Japanese translation(Baifukan, 2001)

5. Concluding remarks
As shown in the previous section, our MPMD approach successfully carried out the calculation. Independent programs, each of which solves diﬀerent
equations, cooperates through the intermediation of
the grand master process.
The two types of execution control, that is, thecontrol of time step and the sequence control of program execution, also works well. Figures 3 and 4 arethe case without the control of time step or the sequence control of program execution, respectively. In
those cases simulation conditions except for the execution control are same as the case in section 4. Theseﬁgures show that lack of execution control alters ordegrades the solution. Without the sequence controlof program execution, the previous time interval data
is used in some cases, or the synchronization of data
could be lost. This is because the degradation of thesolution occurs.
This indicates that appropriate execution control608

A. Takayama et al., A New Framework for Integrated Simulation Model Using MPMD Approach

Fig. 3 Solution of example problem by use of the MPMD
approach with the sequence control of program execution. The diﬀerence between this ﬁgure and ﬁgure 2 is in that the control of time step is switchedoﬀ in this case.

Fig. 4 Solution of example problem by use of the MPMD
approach with the control of time step. The diﬀerence between this ﬁgure and ﬁgure 2 is in that thesequence control of program execution is switchedoﬀ in this case.

609