Object-oriented implementations of the MPDATA advection equation solver in C++, Python and ...

19 Pages·2013·0.38 MB·English

Object-oriented implementations of the MPDATA advection equation solver
in C++, Python and Fortran
SylwesterArabasa,DorotaJareckaa,AnnaJarugaa,MaciejFijałkowskib
aInstituteofGeophysics,FacultyofPhysics,UniversityofWarsaw
bPyPyTeam
3
1
0
Abstract
2
r Three object-oriented implementations of a prototype solver of the advection equation are introduced. The presented programs
aare based on Blitz++ (C++), NumPy (Python), and Fortran’s built-in array containers. The solvers include an implementation
M
of the MultidimensionalPositive-Definite Advective TransportAlgorithm (MPDATA). The introduced codes exemplify how the
applicationofobject-orientedprogramming(OOP)techniquesallowstoreproducethemathematicalnotationusedintheliterature
9
within the program code. A discussion on the tradeoffs of the programming language choice is presented. The main angles
1
of comparison are code brevity and syntax clarity (and hence maintainability and auditability) as well as performance. In the

]case of Python, a significant performancegain is observed when switching from the standard interpreter (CPython) to the PyPy
h
implementationofPython. Entiresourcecodeofallthreeimplementationsisembeddedinthetextandislicensedundertheterms
p
-oftheGNUGPLlicense.
p
mKeywords: object-orientedprogramming,advectionequation,MPDATA,C++,Fortran,Python
o
c
.Contents 1. Introduction
s
c
i1 Introduction 1 Object oriented programming (OOP) ”has become recog-
s
y nised as the almost unique successful paradigm for creating
h2 Implementation 2 complexsoftware”[1,Sec.1.3]. Itisintriguingthat,whilethe
p 2.1 Arraycontainers. . . . . . . . . . . . . . . . . 2 quotedstatementcomesfromtheverybooksubtitledTheArtof
[
2.2 Containersforsequencesofarrays . . . . . . . 3 ScientificComputing,hardlyany(ifnotnone)ofthecurrently

2 2.3 Staggeredgrid . . . . . . . . . . . . . . . . . . 4 operational weather and climate prediction systems - flagship
v 2.4 Haloregions. . . . . . . . . . . . . . . . . . . 5 examplesof complex scientific software - make extensive use
4
2.5 Arrayindexpermutations . . . . . . . . . . . . 5 ofOOPtechniques. Fortranhasbeenthelanguageofchoicein
3
3 2.6 Prototypesolver . . . . . . . . . . . . . . . . . 6 oceanic[2], weather-prediction[3]andEarthsystem[4]mod-
1 2.7 Periodicboundaries(C++) . . . . . . . . . . . 8 elling,andnoneofits20-centuryeditionswereobject-oriented
1. 2.8 Donor-cellformulæ(C++) . . . . . . . . . . . 8 languages[seee.g.5,fordiscussion].
ApplicationofOOPtechniquesindevelopmentofnumerical
0 2.9 Donor-cellsolver(C++) . . . . . . . . . . . . 9
3 2.10 MPDATAformulæ(C++) . . . . . . . . . . . 9 modellingsoftwaremayhelpto:
1
2.11 MPDATAsolver(C++) . . . . . . . . . . . . . 10
: (i) maintain modularity and separation of programlogic lay-
v 2.12 Usageexample(C++) . . . . . . . . . . . . . 10
ers(e.g.separationofnumericalalgorithms,parallelisation
i
X mechanisms,datainput/output,errorhandlingandthede-
3 Performanceevaluation 12
r scriptionofphysicalprocesses);and
a
4 Discussiononthetradeoffsoflanguagechoice 13 (ii) shorten and simplify the source code and improve its
4.1 OOPforblackboardabstractions . . . . . . . . 14 readabilitybyreproducingwithintheprogramlogicthe
4.2 Performance . . . . . . . . . . . . . . . . . . . 14 mathematicalnotationusedintheliterature.
4.3 Easeofuseandabuse . . . . . . . . . . . . . . 14
4.4 Addedvalues . . . . . . . . . . . . . . . . . . 14 The first application is attainable, yet arguably cumbersome,
withproceduralprogramming. Thelatter,virtuallyimpossible
5 Summaryandoutlook 15 toobtainwithproceduralprogramming,isthefocusofthispa-
per. It also enables the compiler or library authors to relieve
Appendix P Pythoncodeforsections2.7–2.11 16 theuser(i.e. scientificprogrammer)fromhand-codingoptimi-
sations,apracticelongrecognisedashavingastrongnegative
Appendix F Fortrancodeforsections2.7–2.11 17 impactwhendebuggingandmaintenanceareconsidered[6].
PreprintsubmittedtoComputerPhysicsCommunications December11,2013
MPDATA [7] stands for MultidimensionalPositive Definite Programminglanguageconstructswheninlinedin thetextare
AdvectiveTransportAlgorithmandisanexampleofanumer- typesetinbold,e.g. GOTO2.
ical procedure used in weather, climate and ocean simulation
systems [e.g. 8, 9, 10, respectively]. MPDATA is a solver for
2. Implementation
systemsofadvectionequationsofthefollowingform:
Doubleprecisionfloating-pointformatisusedinallthreeim-
∂ψ=−∇·(~vψ) (1)
t
plementations.Thecodesbeginwiththefollowingdefinitions:
that describe evolution of a scalar field ψ transported by the listing C.1 (C++)
3 typedef double real_t;
fluid flow with velocity ~v. Quoting Numerical Recipes once
more,developmentofmethodstonumericallysolvesuchprob- listing P.1 (Python)
3 real_t = ’float64’
lems”isanartasmuchasascience”[1,Sec.20.1],andMP-
DATA is an example of the state-of-the art in this field. MP- listing F.1 (Fortran)
3 module real_m
DATA is designed to accurately solve equation (1) in an ar- 4 implicit none
bitrarynumberofdimensionsassuringpositive-definitenessof 5 integer, parameter :: real_t = kind(0.d0)
6 end module
scalarfieldψandincurringsmallnumericaldiffusion.Allrele-
vantMPDATAformulæaregiveninthetextbutarepresented
whichprovideaconvenientwayofswitchingtodifferentpreci-
without derivationor detailed discussion. For a recent review
sion.
ofMPDATA-basedtechniquesseeSmolarkiewicz[11,andref-
Allcodesarestructuredinawayallowingcompilationofthe
erencestherein].
code in exactly the same order as presented in the text within
Inthispaperweintroduceanddiscussobject-orientedimple-
one sourcefile, hence everyFortranlisting containsdefinition
mentationsofanMPDATA-basedtwo-dimensional(2D)advec-
ofaseparatemodule.
tionequationsolverwritteninC++11(ISO/IEC14882:2011),
Python [13] and Fortran 2008 (ISO/IEC 1539-1:2010). In
2.1. Arraycontainers
the following section we introduce the three implementations
brieflydescribingthealgorithmitselfanddiscussingwhereand Solutionof equation(1) using MPDATA impliesdiscretisa-
how the OOP techniques may be applied in its implementa- tion onto a grid of the ψ and the Courant numberC~ = ~v· ∆t
∆x
tion. ThesyntaxandnomenclatureofOOPtechniquesareused fields, where∆t isthesolvertimestepand∆x isthegridspac-
without introduction, for an overview of OOP in context of ing.
C++, Python and Fortran, consult for example [15, Part II], Presented C++ implementation of MPDATA is built upon
[16, Chapter 5] and [17, Chapter 11], respectively. The third the Blitz++ library1. Blitz offers object-oriented representa-
sectionofthispapercoversperformanceevaluationofthethree tion of n-dimensional arrays, and array-valued mathematical
implementations. The fourth section covers discussion of the expressions. Inparticular,itoffersloop-freenotationforarray
tradeoffs of the programminglanguage choice. The fifth sec- arithmeticsthatdoesnotincurcreationofintermediatetempo-
tionclosesthearticlewithabriefsummary. rary objects. Blitz++ is a header-only library2 – to use it, it
Throughoutthe paperwepresentthethreeimplementations isenoughtoincludetheappropriateheaderfile,andoptionally
by discussing source code listings which coverthe entire pro- exposetherequiredclassestothepresentnamespace:
gram code. Subsections 2.1-2.6 describe all three implemen- listing C.2 (C++)
4 #include <blitz/array.h>
tations,whilesubsequentsections2.7-2.12coverdiscussionof 5 using arr_t = blitz::Array<real_t, 2>;
C++codeonly.TherelevantpartsofPythonandFortrancodes 6 using rng_t = blitz::Range;
7 using idx_t = blitz::RectDomain<2>;
do not differ significantly, and for readability reasons are pre-
sentedinAppendix PandAppendix F,respectively.
Here arr_t, rng_t and idx_t serve as alias identifiers and are
TheentirecodeislicensedunderthetermsoftheGNUGen-
introducedinordertoshortenthecode.
eralPublicLicenselicenseversion3[18].
ThepowerofBlitz++ comesfromthe abilitytoexpressar-
All listings include line numbers printed to the left of the
ray expressions as objects. In particular, it is possible to de-
source code, with separate numbering for C++ (listings pre-
fine a function that returns an array expression; i.e. not the
fixedwithC,blackframe),
resultant array, but an object representing a „recipe” defining
listing C.0 (C++)
1 // code licensed under the terms of GNU GPL v3 theoperationstobeperformedonthearguments. Asa conse-
2 // copyright holder: University of Warsaw quence,the returntypesof suchfunctionsbecomeunintelligi-
ble. Luckily,theautoreturntypedeclarationfromtheC++11
Python(listingsprefixedwithP,blueframe)and
standardallowstosimplifythecodesignificantly,evenmoreif
listing P.0 (Python)
1 # code licensed under the terms of GNU GPL v3 usedthroughthefollowingpreprocessormacro:
2 # copyright holder: University of Warsaw
Fortran(listingsprefixedwithF,redframe). 1Blitz++ is a C++ class library for scientific computing which uses
listing F.0 (Fortran) the expression templates technique to achieve high performance, see
1 ! code licensed under the terms of GNU GPL v3 http://sf.net/projects/blitz/
2 ! copyright holder: University of Warsaw 2Blitz++requireslinkingwithlibblitzifdebuggingmodeisused
2
listing C.3 (C++) C[x] and C[y] arrays constituting the sequence have different
8 #define return_macro(expr) \ sizes (see discussion of the Arakawa-C grid in section 2.3).
9 -> decltype(safeToReturn(expr)) \
10 { return safeToReturn(expr); } Second,theorderofdimensionswouldneedtobedifferentfor
differentlanguagesto assure that the contiguousdimension is
The call to blitz::safeToReturn() function is included in or- usedforoneofthespacedimensionsandnotfortimelevels.
der to ensure that all arrays involved in the expression being In the C++ implementation the Boost5 ptr_vector class
returnedcontinuetoexistinthecallerscope.Forexample,def- is used to represent sequences of Blitz++ arrays and at the
inition of a function returning its array-valued argumentdou- same time to handle automatic freeing of dynamically allo-
bled, reads: auto f(arr_t x) return_macro(2*x). This is the cated memory. The ptr_vectorclass is furthercustomisedby
onlypreprocessormacrodefinedherein. defininga derivedstructure which element-access[ ] operator
For the Python implementation of MPDATA the NumPy3 isoverloadedwithamodulovariant:
package is used. In order to make the code compatible with listing C.4 (C++)
11 #include <boost/ptr_container/ptr_vector.hpp>
boththe standardCPythonaswellasthe alternativePyPyim- 12 struct arrvec_t : boost::ptr_vector<arr_t>
plementationofPython[19],thePythoncodeincludesthefol- 13 {
14 const arr_t &operator[](const int i) const
lowingsequenceofimportstatements: 15 {
listing P.2 (Python) 16 return this->at((i + this->size()) % this->size());
4 try: 17 }
5 import numpypy 18 };
6 from _numpypy.pypy import set_invalidation
7 set_invalidation(False) Consequently the last element of any such sequence may be
8 except ImportError:
9 pass accessedatindex-1,thelastbutoneat-2,andsoon.
10 import numpy InthePythonimplementationthebuilt-intupletypeisused
11 try:
12 numpy.seterr(all=’ignore’) tostoresequencesofNumPyarrays. Employmentofnegative
13 except AttributeError: indices for handling from-the-endaddressing of elements is a
14 pass
built-infeatureofallsequencecontainersinPython.
First, the PyPy’s built-in NumPy implementation named Fortrandoesnotfeatureanybuilt-insequencecontainerca-
numpypy is imported if applicable (i.e. if running PyPy), pableofstoringarrays,henceacustomarrvec_ttypeisintro-
and the lazy evaluation mode is turned on through the duced:
set_invalidation(False) call. PyPy’s lazy evaluation obtained listing F.2 (Fortran)
7 module arrvec_m
with the help of a just-in-time compiler enables to achieve an 8 use real_m
analogous to Blitz++ temporary-array-freehandling of array- 9 implicit none
10
valued expressions (see discussion in section 3). Second, to 11 type :: arr_t
match the settingsof C++ and Fortrancompilersused herein, 12 real(real_t), allocatable :: a(:,:)
13 end type
the NumPy package is instructed to ignore any floating-point 14
errors, if such an option is available in the interpreter4. The 15 type :: arrptr_t
16 class(arr_t), pointer :: p
above lines conclude all code modificationsthat needed to be 17 end type
addedinordertorunthecodewithPyPy. 18
19 type :: arrvec_t
Among the three considered languages only Fortran is 20 class(arr_t), allocatable :: arrs(:)
equippedwithbuilt-inarrayhandlingfacilitiesofpracticaluse 21 class(arrptr_t), allocatable :: at(:)
22 integer :: length
in high-performance computing. Therefore, there is no need 23 contains
forusinganexternalpackageaswithC++andPython.Fortran 24 procedure :: ctor => arrvec_ctor
25 procedure :: init => arrvec_init
array-handlingfeaturesarenotobject-oriented,though. 26 end type
27
28 contains
2.2. Containersforsequencesofarrays 29
30 subroutine arrvec_ctor(this, n)
Asdiscussedabove,discretisationinspaceofthescalarfield 31 class(arrvec_t) :: this
ψ(x,y) into its ψ grid representation requires floating-point 32 integer, intent(in) :: n
[i,j] 33
arraycontainers. Inturn,discretisationintimerequiresacon- 34 this%length = n
tainer class for storing sequences of such arrays, i.e. {ψ[n], 35 allocate(this%at( -n : n-1 ))
36 allocate(this%arrs( 0 : n-1 ))
ψ[n+1]}. Similarly the componentsof the vector fieldC~ are in 37 end subroutine
facta{C[x],C[y]}arraysequence. 38
39 subroutine arrvec_init(this, n, i, j)
Using an additional array dimension to represent the se- 40 class(arrvec_t), target :: this
quence elements is not considered for two reasons. First, the 41 integer, intent(in) :: n
42 integer, intent(in) :: i(2), j(2)
43
3NumPy is a Python package for scientific computing offering support
for multi-dimensional arrays and a library of numerical algorithms, see 5Boostisafreeandopen-sourcecollectionofpeer-reviewedC++libraries
http://numpy.org/ availableathttp://www.boost.org/.SeveralpartsofBoosthavebeeninte-
4numpy.seterr()isnotsupportedinPyPyasofversion1.9 gratedintoorinspirednewadditionstotheC++standard.
3
4445 atlhlios%caatt(e(nt)h%ips=%>artrhsi(sn%)a%rar(si((n1)) : i(2), j(1) : j(2) )) Cve[[icx+]t1o/2,rj]c,oCm[[iyp,]jo−1/n2]enantsdsCur[[iry,o]j+u1/n2]dtiongdψepict.thHeowgreivdevr,alfuraecstioofntahleinC~-
46 this%at(n - this%length)%p => this%arrs(n) [i,j]
47 end subroutine dexing does not have a built-in counterpart in any of the em-
48 end module ployedprogramminglanguages.Adesiredsyntaxwouldtrans-
latei−1/2toi−1andi+1/2toi. OOPoffersaconvenientwayto
The arr_t type is defined solely for the purpose of overcom- implementsuchnotationbyoverloadingthe+and-operators
ing the limitation of lack of an array-of-arraysconstruct, and
forobjectsrepresentingarrayindices.
its only memberfield is a two-dimensionalarray. An arrayof
In the C++ implementation first a global instance h of an
arr_tisusedhereinafterasacontainerforsequencesofarrays.
empty structure hlf_t is defined, and then the plus and minus
Thearrptr_ttypeisdefinedsolelyforthepurposeofover-
operatorsforhlf_tandrng_tareoverloaded:
coming Fortran’s limitation of not supporting allocatables of
listing C.5 (C++)
pointers. arrptr_t’ssingle memberfield is a pointerto an in- 19 struct hlf_t {} h;
20
stance of arr_t. Creating an allocatable of arrptr_t, instead 21 inline rng_t operator+(const rng_t &i, const hlf_t &)
ofamulti-elementpointerofarr_t,ensuresautomaticmemory 22 {
23 return i;
deallocation. 24 }
Type arrptr_t is used to implement the from-the-end ad- 25
26 inline rng_t operator-(const rng_t &i, const hlf_t &)
dressing of elements in arrvec_t. The array data is stored in 27 {
thearrsmemberfield(oftypearr_t).Theatmemberfield(of 28 return i-1;
29 }
type arrptr_t) stores pointers to the elements of arrs. at has
double the length of arrs and is initialised in a cyclic manner
Thisway, the arraysrepresentingvectorfield componentscan
so that the -1 elementof at pointsto the last elementof arrs, be indexed using (i+h,j), (i-h,j) etc. where h represents the
and so on. Assuming psi is an instance of arrptr_t, the (i,j)
half.
elementofthen-tharrayinpsimaybeaccessedwith
In NumPy in order to preventcopying of array data during
psi%at(n)%p%a(i,j). slicingoneneedstooperateontheso-calledarrayviews. Array
Thector(n)methodinitialisesthecontainerforagivennum-
viewsareobtainedwhenindexingthearrayswithobjectsofthe
berofelementsn.Theinit(n,i,j)methodinitialisesthen-thele- Python’sbuilt-itslicetype(ortuplesofsuchobjectsincaseof
mentofthecontainerwithanewlyallocated2Darrayspanning
multi-dimensionalarrays).Pythonforbidsoverloadingofoper-
indices i(1):i(2), and j(1):j(2)in the first, and last dimensions
atorsofbuilt-intypessuchasslices, anddoesnotdefineaddi-
respectively6. tion/subtractionoperatorsforsliceandintpairs.Consequently,
a customlogichasto bedefinednotonlyforfractionalindex-
2.3. Staggeredgrid
ing,butalsoforshiftingtheslicesbyintegerintervals(i±1).It
isimplementedherebydeclaringashiftclasswiththeadequate
operatoroverloads:
listing P.3 (Python)
15 class shift():
ψ[i,j+1] 16 def __init__(self, plus, mnus):
17 self.plus = plus
18 self.mnus = mnus
19 def __radd__(self, arg):
C[[ix−]1/2,j] C[[ix+]1/2,j] 222012 reaatrruggr..nsstttyaoprpet(++arssgee)ll(ff..pplluuss,
ψ[i−1,j] ψ[i,j] 23 )
24 def __rsub__(self, arg):
25 return type(arg)(
C[y] 26 arg.start - self.mnus,
[i,j−1/2] 27 arg.stop - self.mnus
28 )
andtwoinstancesofittorepresentunityandhalfinexpressions
likei+one,i+hlf,whereiisaninstanceofslice7:
Figure1:AschematicoftheArakawa-Cgrid. listing P.4 (Python)
29 one = shift(1,1)
30 hlf = shift(0,1)
The so-called Arakawa-C staggered grid [20] depicted in
Figure 1 is a natural choice for MPDATA. As a consequence, InFortranfractionalarrayindexingisobtainedthroughdef-
the discretised representations of the ψ scalar field, and each initionandinstantiationofanobjectrepresentingthehalf,and
componentof theC~ =~v· ∆t vectorfield in eq. (1) are defined havingappropriateoperatoroverloads:
∆x
over different grid point locations. In mathematical notation
thiscanbeindicatedbyusageoffractionalindices,e.g.C[x] ,
[i−1/2,j] 7One could argue that not using an own implementation of a slice-
representingclassinNumPyisadesignflaw–beingabletomodifybehaviour
ofahypotheticalnumpy.sliceclassthroughinheritancewouldallowtoimple-
6InFortran,whenanarrayispassedasafunctionargumentitsbaseislocally mentthesamebehaviourasobtainedinlistingP.3withouttheneedtorepresent
settounity,regardlessofthesettingatthecallerscope. theunityasaseparateobject
4
listing F.3 (Fortran) 95 integer :: return(2)
49 module arakawa_c_m 96
50 implicit none 97 return = (/ r(1) - n, r(2) + n /)
51 98 end function
52 type :: half_t 99
53 end type 100 function ext_h(r, h) result (return)
54 101 integer, intent(in) :: r(2)
55 type(half_t) :: h 102 type(half_t), intent(in) :: h
56 103 integer :: return(2)
57 interface operator (+) 104
58 module procedure ph 105 return = (/ r(1) - h, r(2) + h /)
59 end interface 106 end function
60 107 end module
61 interface operator (-)
62 module procedure mh
63 end interface Consequently, a range depicted by i± 1/2 may be expressed
64
65 contains in the codeasext(i, h). In allthreeimplementationsthe ext()
66 functionacceptthesecondargumenttobeanintegerora”half”
67 elemental function ph(i, h) result (return)
68 integer, intent(in) :: i (cf. section2.3).
69 type(half_t), intent(in) :: h
70 integer :: return
71 return = i 2.5. Arrayindexpermutations
72 end function Hereinafter,theπd symbolisusedtodenoteacyclicpermu-
73 a,b
74 elemental function mh(i, h) result (return) tation of an orderd of a set {a,b}. It is used to generalise the
75 integer, intent(in) :: i MPDATAformulæintomultipledimensionsusingthefollow-
76 type(half_t), intent(in) :: h
77 integer :: return ingnotation:
78 return = i - 1
79 end function
80 end module 1
ψ ≡ψ +ψ
[i,j]+πd [i+1,j] [i,j+1]
Xd=0 1,0
2.4. Haloregions
TheMPDATAformulædefiningψ[n+1] asafunctionofψ[n] Blitz++ ships with the RectDomain class (aliased here as
[i,j] [i,j] idx_t)forspecifyingarrayrangesinmultipledimensions. The
(discussed in the following sections) feature terms such as
πpermutationisimplementedinC++asafunctionpi()return-
ψ . Onewayofassuringvalidityoftheseformulæonthe
[i−1,j−1]
inganinstanceofidx_t. Inordertoensurecompile-timeeval-
edgesofthedomain(e.g. fori=0)istointroducetheso-called
uation,thepermutationorderispassedviathetemplateparam-
haloregionsurroundingthedomain.Themethodofpopulating
eterd(notethedifferentorderofiandjargumentsinthetwo
the halo region with data depends on the boundary condition
templatespecialisations):
type.Employmentofthehalo-regionlogicimpliesrepeatedus-
ageofarrayrangeextensionsinthecodesuchasi{i±halo. listing C.7 (C++)
37 template<int d>
An ext() function is defined in all three implementation, in 38 inline idx_t pi(const rng_t &i, const rng_t &j);
ordertosimplifycodingofarrayrangeextensions: 39
40 template<>
listing C.6 (C++) 41 inline idx_t pi<0>(const rng_t &i, const rng_t &j)
30 template<class n_t> 42 {
31 inline rng_t ext(const rng_t &r, const n_t &n) { 43 return idx_t({i,j});
32 return rng_t( 44 };
33 (r - n).first(), 45
34 (r + n).last() 46 template<>
35 ); 47 inline idx_t pi<1>(const rng_t &j, const rng_t &i)
36 } 48 {
49 return idx_t({i,j});
listing P.5 (Python) 50 };
31 def ext(r, n):
32 if (type(n) == int) & (n == 1):
33 n = one NumPy uses tuples of slices for addressing multi-
34 return slice( dimensionalarraywithasingleobject.Therefore,thefollowing
35 (r - n).start,
36 (r + n).stop definitionoffunctionpi()sufficestorepresentπ:
37 ) listing P.6 (Python)
listing F.4 (Fortran) 38 def pi(d, *idx):
81 module halo_m 39 return (idx[d], idx[d-1])
82 use arakawa_c_m
83 implicit none In the Fortran implementation pi() returns a pointer to the
84
85 interface ext array elements specified by i and j interpreted as (i,j) or (j,i)
86 module procedure ext_n dependingonthevalueoftheargumentd. Inadditiontopi(),a
87 module procedure ext_h
88 end interface helperspan()functionreturningthelengthofoneofthevectors
89 passedasargumentisdefined:
90 contains
91 listing F.5 (Fortran)
92 function ext_n(r, n) result (return) 108 module pi_m
93 integer, intent(in) :: r(2) 109 use real_m
94 integer, intent(in) :: n 110 implicit none
5
111 contains 67 bcx(i, j, hlo),
112 function pi(d, arr, i, j) result(return) 68 bcy(j, i, hlo)
113 integer, intent(in) :: d 69 {
114 real(real_t), allocatable, target :: arr(:,:) 70 for (int l = 0; l < 2; ++l)
115 real(real_t), pointer :: return(:,:) 71 psi.push_back(new arr_t(ext(i, hlo), ext(j, hlo)));
116 integer, intent(in) :: i(2), j(2) 72 C.push_back(new arr_t(ext(i, h), ext(j, hlo)));
117 select case (d) 73 C.push_back(new arr_t(ext(i, hlo), ext(j, h)));
118 case (0) 74 }
119 return => arr( i(1) : i(2), j(1) : j(2) ) 75
120 case (1) 76 // accessor methods
121 return => arr( j(1) : j(2), i(1) : i(2) ) 77 arr_t state() {
122 end select 78 return psi[n](i,j).reindex({0,0});
123 end function 79 }
124 80
125 pure function span(d, i, j) result(return) 81 arr_t courant(int d)
126 integer, intent(in) :: i(2), j(2) 82 {
127 integer, intent(in) :: d 83 return C[d];
128 integer :: return 84 }
129 select case (d) 85
130 case (0) 86 // helper methods invoked by solve()
131 return = i(2) - i(1) + 1 87 virtual void advop() = 0;
132 case (1) 88
133 return = j(2) - j(1) + 1 89 void cycle()
134 end select 90 {
135 end function 91 n = (n + 1) % 2 - 2;
136 end module 92 }
93
94 // integration logic
The span() function is used to shorten the declarations of ar- 95 void solve(const int nt)
rays to be returnedfromfunctionsin the Fortranimplementa- 96 {
97 for (int t = 0; t < nt; ++t)
tion(seelistingsF.11andF.17–F.20). 98 {
It isworthnotingherethatthe C++ implementationofpi() 99 bcx.fill_halos(psi[n], ext(j, hlo));
100 bcy.fill_halos(psi[n], ext(i, hlo));
isbranchlessthankstoemploymentoftemplatespecialisation. 101 advop();
With Fortran one needs to rely on compiler optimisations to 102 cycle();
103 }
eliminate the conditional expression within the pi() that de- 104 }
pendsonvalueofdwhichisalwaysknownatcompiletime. 105 };
Thesolverstructureisanabstractdefinition(containingapure
2.6. Prototypesolver
virtualmethod)requiringitsdescendantstoimplementatleast
The tasks to be handled by a prototype advection equation theadvop()methodwhichisexpectedtofillpsi[n+1]withan
solverproposedhereinare:
updated(advected)valuesofpsi[n]. Thetwotemplateparame-
(i) storingarraysrepresentingtheψandC~ fieldsandanyre- tersbcx_tandbcy_tallowthesolvertooperatewithanykind
ofboundaryconditionstructuresthatfulfiltherequirementsim-
quiredhousekeepingdata,
pliedbythecallstothemethodsofbcxandbcy,respectively.
(ii) allocatinganddeallocatingtherequiredmemory, Thedonor-cellandMPDATAschemesbothrequireonlythe
previousstate of an advectedfield in order to advancethe so-
(iii) providingaccesstothesolverstate,
lution. Consequently, memory for two time levels (ψ[n] and
ψ[n+1]) is allocated in the constructor. The sizes of the arrays
(iv) performing the integration by invoking the advection-
representingthetwotimelevelsofψaredefinedbythedomain
operatorandboundary-conditionhandlingroutines.
size(nx×ny)plusthehaloregion. Thesizeofthehaloregion
In the following C++ definition of the solver structure, task is an argumentof the constructor. The cycle()methodis used
(i) is represented with the definition of the structure member toswapthetimelevelswithoutcopyinganydata.
fields;task(ii)issplitbetweenthesolver’sconstructorandthe The arrays representingthe C[x] andC[y] componentsofC~,
destructors of arrvec_t; task (iii) is handled by the accessor require(nx+1)×nyandnx×(ny+1)elements,respectively(be-
methods;task(iv)ishandledwithinthesolvemethod: inglaidoutontheArakawa-Cstaggeredgrid).
listing C.8 (C++) PythondefinitionofthesolverclassfollowscloselytheC++
51 template<class bcx_t, class bcy_t> structuredefinition:
52 struct solver
53 { listing P.7 (Python)
54 // member fields 40 class solver(object):
55 arrvec_t psi, C; 41 # ctor-like method
56 int n, hlo; 42 def __init__(self, bcx, bcy, nx, ny, hlo):
57 rng_t i, j; 43 self.n = 0
58 bcx_t bcx; 44 self.hlo = hlo
59 bcy_t bcy; 45 self.i = slice(hlo, nx + hlo)
60 46 self.j = slice(hlo, ny + hlo)
61 // ctor 47
62 solver(int nx, int ny, int hlo) : 48 self.bcx = bcx(0, self.i, hlo)
63 hlo(hlo), 49 self.bcy = bcy(1, self.j, hlo)
64 n(0), 50
65 i(0, nx-1), 51 self.psi = (
66 j(0, ny-1), 52 numpy.empty((
6
53 ext(self.i, self.hlo).stop, 139 implicit none
54 ext(self.j, self.hlo).stop 140
55 ), real_t), 141 type, abstract :: bcd_t
56 numpy.empty(( 142 contains
57 ext(self.i, self.hlo).stop, 143 procedure(bcd_fill_halos), deferred :: fill_halos
58 ext(self.j, self.hlo).stop 144 procedure(bcd_init), deferred :: init
59 ), real_t) 145 end type
60 ) 146
61 147 abstract interface
62 self.C = ( 148 subroutine bcd_fill_halos(this, a, j)
63 numpy.empty(( 149 import :: bcd_t, real_t
64 ext(self.i, hlf).stop, 150 class(bcd_t ) :: this
65 ext(self.j, self.hlo).stop 151 real(real_t), allocatable :: a(:,:)
66 ), real_t), 152 integer :: j(2)
67 numpy.empty(( 153 end subroutine
68 ext(self.i, self.hlo).stop, 154
69 ext(self.j, hlf).stop 155 subroutine bcd_init(this, d, n, hlo)
70 ), real_t) 156 import :: bcd_t
71 ) 157 class(bcd_t) :: this
72 158 integer :: d, n, hlo
73 # accessor methods 159 end subroutine
74 def state(self): 160 end interface
75 return self.psi[self.n][self.i, self.j] 161 end module
76
77 # helper methods invoked by solve()
78 def courant(self,d): Having defined the abstract type for boundary-condition ob-
79 return self.C[d][:] jects, a definition of a solver class following closely the C++
80
81 def cycle(self): andPythoncounterpartsmaybeprovided:
82 self.n = (self.n + 1) % 2 - 2 listing F.7 (Fortran)
83 162 module solver_m
84 # integration logic 163 use arrvec_m
85 def solve(self, nt): 164 use bcd_m
86 for t in range(nt): 165 use arakawa_c_m
87 self.bcx.fill_halos( 166 use halo_m
88 self.psi[self.n], ext(self.j, self.hlo) 167 implicit none
89 ) 168
90 self.bcy.fill_halos( 169 type, abstract :: solver_t
91 self.psi[self.n], ext(self.i, self.hlo) 170 class(arrvec_t), allocatable :: psi, C
92 ) 171 integer :: n, hlo
93 self.advop() 172 integer :: i(2), j(2)
94 self.cycle() 173 class(bcd_t), pointer :: bcx, bcy
95 174 contains
175 procedure :: solve => solver_solve
The key difference stems from the fact that, unlike Blitz++, 176 procedure :: state => solver_state
177 procedure :: courant => solver_courant
NumPydoesnotallowanarrayto havearbitraryindexbase– 178 procedure :: cycle => solver_cycle
inNumPythefirstelementisalwaysaddressedwith0. Conse- 179 procedure(solver_advop), deferred :: advop
180 end type
quently,whileinC++(andFortran)thecomputationaldomain 181
ischosentostartat(i=0,j=0)andhenceapartofthehalore- 182 abstract interface
183 subroutine solver_advop(this)
gion to have negativeindices, in Pythonthe halo regionstarts 184 import solver_t
at(0,0)8. However,sincethewholehalologicishiddenwithin 185 class(solver_t), target :: this
186 end subroutine
thesolver,suchdetailsarenotexposedtotheuser. Thebcxand 187 end interface
bcyboundary-conditionspecificationsarepassedtothesolver 188
189 contains
throughconstructor-like__init__()methodasopposedtotem- 190
plateparametersinC++. 191 subroutine solver_ctor(this, bcx, bcy, nx, ny, hlo)
192 use arakawa_c_m
The above C++ and Python prototype solvers in principle 193 use halo_m
allow to operatewith anyboundaryconditionobjectsthatim- 194 class(solver_t) :: this
195 class(bcd_t), intent(in), target :: bcx, bcy
plement methods called from within the solver. This require- 196 integer, intent(in) :: nx, ny, hlo
ment is checked at compile-time in the case of C++, and at 197
198 this%n = 0
run-time in the case of Python. In order to obtain an analo- 199 this%hlo = hlo
gous behaviour with Fortran, it is required to define, prior to 200 this%bcx => bcx
201 this%bcy => bcy
definition of a solvertype, an abstracttype with deferredpro- 202
cedureshavingabstractinterfaces[sic!,seeTable2.1in21,for 203 this%i = (/ 0, nx - 1 /)
204 this%j = (/ 0, ny - 1 /)
asummaryofapproximatecorrespondenceofOOPnomencla- 205
turebetweenFortranandC++]: 206 call bcx%init(0, nx, hlo)
207 call bcy%init(1, ny, hlo)
listing F.6 (Fortran)
113378 moudsueleabrcrdv_emc_m 222001890 aclallolctahtie(st%phsiis%%pcstio)r(2)
211 block
8Thereasontoallowthedomaintobeginatanarbitraryindexismainlyto 212 integer :: n
easedebuggingincasethecodewouldbeusedinparallelcomputationsusing 213 do n=0, 1
214 call this%psi%init( &
domaindecompositionwhereeachsubdomaincouldhaveitsownindexbase 215 n, ext(this%i, hlo), ext(this%j, hlo) &
correspondingtothelocationwithinthecomputationaldomain 216 )
7
217 end do 124 void fill_halos(const arr_t &a, const rng_t &j)
218 end block 125 {
219 126 a(pi<d>(left_halo, j)) = a(pi<d>(rght_edge, j));
220 allocate(this%C) 127 a(pi<d>(rght_halo, j)) = a(pi<d>(left_edge, j));
221 call this%C%ctor(2) 128 }
222 call this%C%init( & 129 };
223 0, ext(this%i, h), ext(this%j, hlo) &
224 )
225 call this%C%init( & Ashintedbythememberfieldnames,thefill_halos()meth-
226 1, ext(this%i, hlo), ext(this%j, h) & ods fill the left/righthalo regionswith data from the right/left
227 )
228 end subroutine edges of the domain. Thanks to employment of the function
229 pi() describedin section 2.5 the same code may be applied in
230 function solver_state(this) result (return)
231 class(solver_t) :: this anydimension(herebeingatemplateparameter).
232 real(real_t), pointer :: return(:,:) ListingsP.8andF.8containthePythonandFortrancounter-
233 return => this%psi%at(this%n)%p%a( &
234 this%i(1) : this%i(2), & partstolistingC.9.
235 this%j(1) : this%j(2) &
223367 en)d function 2.8. Donor-cellformulæ(C++)
238 MPDATA is an iterative algorithm in which each iteration
239 function solver_courant(this, d) result (return)
240 class(solver_t) :: this takestheformoftheso-calleddonor-cellformula(whichitself
241 integer :: d isafirst-orderadvectionscheme).
242 real(real_t), pointer :: return(:,:)
243 return => this%C%at(d)%p%a MPDATA and donor-cell are explicit forward-in-time algo-
244 end function rithms–theyallowtopredictψ[n+1] asafunctionofψ[n] where
245
246 subroutine solver_cycle(this) n and n + 1 denote two adjacent time levels. The donor-cell
247 class(solver_t) :: this schememaybewrittenas[eq. 2in7]:
248 this%n = mod(this%n + 1 + 2, 2) - 2
249 end subroutine
250 N−1
251 subroutine solver_solve(this, nt) ψ[n+1]=ψ[n] − F ψ[n] ,ψ[n] ,C[d]
225523 cilnatsesg(esr,olivnetre_ntt)(:i:n)th:i:snt [i,j] [i,j] Xd=0(cid:18) (cid:20) [i,j] [i,j]+πd1,0 [i,j]+πd1/2,0(cid:21) (2)
254 integer :: t
255 −F ψ[n] ,ψ[n] ,C[d]
225567 doctal=l t0,hinst%-b1cx%fill_halos( & " [i,j]+πd−1,0 [i,j] [i,j]+πd−1/2,0#!
258 this%psi%at(this%n)%p%a, ext(this%j, this%hlo) & where N is the number of dimensions, and F is the so-called
259 )
260 call this%bcy%fill_halos( & fluxfunction[7,eq.3]:
261 this%psi%at(this%n)%p%a, ext(this%i, this%hlo) &
262 ) F(ψ ,ψ ,C)=max(C,0)·ψ +min(C,0)·ψ
263 call this%advop() L R L R
264 call this%cycle() C+|C| C−|C| (3)
265 end do = ·ψ + ·ψ
266 end subroutine 2 L 2 R
267 end module
ThefluxfunctiontakesthefollowingforminC++:
listing C.10 (C++)
2.7. Periodicboundaries(C++) 130 template<class T1, class T2, class T3>
131 inline auto F(
From this point, only C++ implementation is explained in 132 const T1 &psi_l, const T2 &psi_r, const T3 &C
themaintext. ThePythonandFortranimplementationsarein- 133 ) return_macro(
134 (
cludedinappendicesPandF. 135 (C + abs(C)) * psi_l +
Thesolverdefinitiondescribedinsection2.6requiresagiven 136 (C - abs(C)) * psi_r
137 ) / 2
boundaryconditionobjecttoimplementafill_halos()method. 138 )
AnimplementationofperiodicboundaryconditionsinC++is
providedinthefollowinglisting: Equation2 is split into the termsunder the summation(effec-
listing C.9 (C++) tivelythe1-dimensionaldonor-cellformula):
106 template<int d> listing C.11 (C++)
107 struct cyclic 139 template<int d>
108 { 140 inline auto donorcell(
109 // member fields 141 const arr_t &psi, const arr_t &C,
110 rng_t left_halo, rght_halo; 142 const rng_t &i, const rng_t &j
111 rng_t left_edge, rght_edge;; 143 ) return_macro(
112 144 F(
113 // ctor 145 psi(pi<d>(i, j)),
114 cyclic( 146 psi(pi<d>(i+1, j)),
115 const rng_t &i, const rng_t &j, int hlo 147 C(pi<d>(i+h, j))
116 ) : 148 ) -
117 left_halo(i.first()-hlo, i.first()-1), 149 F(
118 rght_edge(i.last()-hlo+1, i.last() ), 150 psi(pi<d>(i-1, j)),
119 rght_halo(i.last()+1, i.last()+hlo ), 151 psi(pi<d>(i, j)),
120 left_edge(i.first(), i.first()+hlo-1) 152 C(pi<d>(i-h, j))
121 {} 153 )
122 154 )
123 // method invoked by the solver
8
andtheactualtwo-dimensionaldonor-cellformula: For positive-definite ψ, the A and B terms take the following
form9:
listing C.12 (C++) ψ −ψ
115556 vociodndsotnaorrrcveelcl__top&(psi, const int n, A[[di,]j]= ψ[i,j]+πd1,0+ψ[i,j] (6)
157 const arrvec_t &C, [i,j]+πd [i,j]
158 const rng_t &i, const rng_t &j 1,0
115690 ) p{si[n+1](i,j) = psi[n](i,j) B[d] = 1ψ[i,j]+πd1,1+ψ[i,j]+πd0,1−ψ[i,j]+πd1,−1−ψ[i,j]+πd0,−1 (7)
116612 -- ddoonnoorrcceellll<<01>>((ppssii[[nn]],, CC[[01]],, ij,, ji)); [i,j] 2ψ[i,j]+πd1,1+ψ[i,j]+πd0,1+ψ[i,j]+πd1,−1+ψ[i,j]+πd0,−1
163 } If the denominator in equations 6 or 7 equals zero for a
given i and j, the corresponding A and B are set to zero
[i,j] [i,j]
ListingsP.9-P11andF.9-F.13containthePythonandFortran whatmaybeconvenientlyrepresentedwiththewhereconstruct
counterpartstolistingsC.12-C.15. (availableinallthreeconsideredlanguages):
listing C.14 (C++)
179 template<class nom_t, class den_t>
2.9. Donor-cellsolver(C++) 180 inline auto mpdata_frac(
181 const nom_t &nom, const den_t &den
182 ) return_macro(
Asmentionedintheprevioussection,thedonor-cellformula 183 where(den > 0, nom / den, 0)
constitutesanadvectionscheme,hencewemayuseittocreate 184 )
asolver_donorcellimplementationoftheabstractsolverclass:
TheAtermdefinedinequation6takesthefollowingform:
listing C.13 (C++)
164 template<class bcx_t, class bcy_t> listing C.15 (C++)
165 struct solver_donorcell : solver<bcx_t, bcy_t> 185 template<int d>
166 { 186 inline auto mpdata_A(const arr_t &psi,
167 solver_donorcell(int nx, int ny) : 187 const rng_t &i, const rng_t &j
168 solver<bcx_t, bcy_t>(nx, ny, 1) 188 ) return_macro(
169 {} 189 mpdata_frac(
170 190 psi(pi<d>(i+1, j)) - psi(pi<d>(i,j)),
171 void advop() 191 psi(pi<d>(i+1, j)) + psi(pi<d>(i,j))
172 { 192 )
173 donorcell_op( 193 )
174 this->psi, this->n, this->C,
175 this->i, this->j TheBtermdefinedinequation7takesthefollowingform:
176 );
177 } listing C.16 (C++)
178 }; 194 template<int d>
195 inline auto mpdata_B(const arr_t &psi,
196 const rng_t &i, const rng_t &j
The above definition is given as an example only. In the fol- 197 ) return_macro(
198 mpdata_frac(
lowing sections an MPDATA solver of the same structure is 199 psi(pi<d>(i+1, j+1)) + psi(pi<d>(i, j+1)) -
defined. 200 psi(pi<d>(i+1, j-1)) - psi(pi<d>(i, j-1)),
201 psi(pi<d>(i+1, j+1)) + psi(pi<d>(i, j+1)) +
ListingsP.12andF.14containthePythonandFortrancoun- 202 psi(pi<d>(i+1, j-1)) + psi(pi<d>(i, j-1))
terpartstolistingC.16. 203 ) / 2
204 )
2.10. MPDATAformulæ(C++) Equation5takesthefollowingform:
listing C.17 (C++)
MPDATA introduces corrective steps to the algorithm de- 205 template<int d>
206 inline auto mpdata_C_bar(
fined by equation 2 and 3. Each corrective step is a donor- 207 const arr_t &C,
cellstep(eq.2)withtheCourantnumberfieldscorresponding 208 const rng_t &i,
209 const rng_t &j
to the MPDATA antidiffusive velocities of the followingform 210 ) return_macro(
[eqs13,14in7]: 211 (
212 C(pi<d>(i+1, j+h)) + C(pi<d>(i, j+h)) +
213 C(pi<d>(i+1, j-h)) + C(pi<d>(i, j-h))
C[′i[,dj]]+πd1/2,0=(cid:12)(cid:12)C[[id,]j]+πd1/2,0(cid:12)(cid:12)·(cid:20)1−(cid:12)(cid:12)C[[id,]j]+πd1/2,0(cid:12)(cid:12)(cid:21)·A[[di,]j](ψ) 221145 ) ) / 4
(cid:12)(cid:12) N (cid:12)(cid:12) (cid:12)(cid:12) (cid:12)(cid:12) (4)
−(cid:12) C[(cid:12)d] (cid:12)·C[q] (cid:12)·B[d] (ψ) Equation4takethefollowingform:
q=X0,q,d [i,j]+πd1/2,0 [i,j]+πd1/2,0 [i,j] 216 template<int d> listing C.18 (C++)
217 inline auto mpdata_C_adf(
whereψandCrepresentvaluesfromthepreviousiterationand 218 const arr_t &psi,
219 const rng_t &i, const rng_t &j,
where: 220 const arrvec_t &C
221 ) return_macro(
C[q] = 1 · C[q] +C[q] + 222 abs(C[d](pi<d>(i+h, j)))
[i,j]+πd1/2,0 4 (cid:18) [i,j]+πd1,1/2 [i,j]+πd0,1/2 (5)
C[q] +C[q] 9Sinceψ ≥0,|A|≤ 1and|B|≤ 1. SeeSmolarkiewicz[11,Sec. 4.2]for
[i,j]+πd1,−1/2 [i,j]+πd0,−1/2! descriptionofadaptationoftheformulæforadvectionoffieldsofvariablesign
9
223 * (1 - abs(C[d](pi<d>(i+h, j)))) 296 donorcell_op(
224 * mpdata_A<d>(psi, i, j) 297 this->psi, this->n, C_corr, this->i, this->j
225 - C[d](pi<d>(i+h, j)) 298 );
226 * mpdata_C_bar<d>(C[d-1], i, j) 299 }
227 * mpdata_B<d>(psi, i, j) 300 }
228 ) 301 }
302 };
ListingsP.13-P.17andF.15-F.21containthePythonandFor-
The array of sequences of temporary arrays tmp allocated in
trancounterpartstolistingC.16-C.22.
theconstructorisusedtostoretheantidiffusivevelocitiesfrom
2.11. MPDATAsolver(C++) thepresentandoptionallyprevioustimestep(ifusingmorethan
twoiterations).
An MPDATA solver may be now constructed by inheriting
Theadvop()methodcontrollstheMPDATAiterationswithin
fromsolverclasswiththefollowingdefinitioninC++:
one timestep. The first (step = 0) iteration of MPDATA is an
listing C.19 (C++)
229 template<int n_iters, class bcx_t, class bcy_t> unmodifieddonor-cellstep(comparelistingC.15).Subsequent
230 struct solver_mpdata : solver<bcx_t, bcy_t> iterations begin with calculation of the antidiffusive Courant
231 {
232 // member fields fields using formula 4. In order to calculate values spanning
223334 arrnrgv_etci_mt,tjmmp;[2]; an (i−1⁄2 ... i+1⁄2) range using a formula for C[i+1/2,...] only, the
formula is evaluated using extended index ranges im and jm.
235
236 // ctor In the second (step=1) iteration the uncorrectedCourant field
237 solver_mpdata(int nx, int ny) :
238 solver<bcx_t, bcy_t>(nx, ny, 1), (C_unco) points to the original C field, and the antidiffusive
239 im(this->i.first() - 1, this->i.last()), Courant field is written into C_corr which points to tmp[1].
240 jm(this->j.first() - 1, this->j.last())
241 { In the third (step=2) iteration C_unco points to tmp[1] while
242 int n_tmp = n_iters > 2 ? 2 : 1; C_corrpointstotmp[0]. Insubsequentiterationstmp[0]and
243 for (int n = 0; n < n_tmp; ++n)
244 { tmp[1]arealternatelyswapped.
245 tmp[n].push_back(new arr_t( ListingsP.18andF.22containthePythonandFortrancoun-
246 this->C[0].domain()[0], this->C[0].domain()[1])
247 ); terpartstolistingC.23.
248 tmp[n].push_back(new arr_t(
224590 );this->C[1].domain()[0], this->C[1].domain()[1]) 2.12. Usageexample(C++)
251 } The following listing providesan example of how the MP-
252 }
253 DATAsolverdefinedinsection2.11maybeusedtogetherwith
254 // method invoked by the solver thecyclicboundaryconditionsdefinedinsection2.7. Intheex-
255 void advop()
256 { ample a Gaussian signal is advected in a 2D domain defined
257 for (int step = 0; step < n_iters; ++step) over a grid of 24×24 cells. The program first plots the ini-
258 {
259 if (step == 0) tial condition, then performs the integration for 75 timesteps
260 donorcell_op( with three different settings of the number of iterations used
261 this->psi, this->n, this->C, this->i, this->j
262 ); in MPDATA. The velocity field is constant in time and space
263 else (althoughitisnotassumedinthepresentedimplementations).
264 {
265 this->cycle(); Thesignalshapeattheendofeachsimulationisplottedaswell.
266 this->bcx.fill_halos( Plottingisdonewiththehelpofthegnuplot-iostreamlibrary10.
267 this->psi[this->n], ext(this->j, this->hlo)
268 ); The resultant plot is presented herein as Figure 2. The top
269 this->bcy.fill_halos( paneldepictstheinitialcondition.Thethreeotherpanelsshow
270 this->psi[this->n], ext(this->i, this->hlo)
271 ); a snapshot of the field after 75 timesteps. The donor-cell so-
272 lution is characterisedby strongestnumericaldiffusion result-
273 // choosing input/output for antidiff C
274 const arrvec_t inginsignificantdropinthesignalamplitude. Thesignalsad-
275 &C_unco = (step == 1) vectedusingMPDATA showsmallernumericaldiffusionwith
276 ? this->C
277 : (step % 2) the solution obtained with more iterations preserving the sig-
278 ? tmp[1] // odd steps nal altitude more accurately. In all of the simulations the sig-
279 : tmp[0], // even steps
280 &C_corr = (step % 2) nalmaintainsitspositivedefiniteness. Thedomainperiodicity
281 ? tmp[0] // odd steps is apparentin the plots as the maximumof the signalafter 75
282 : tmp[1]; // even steps
283 timestepsislocatednearthedomainwalls.
284 // calculating the antidiffusive C Listings P.19 and F.23-F.24 contain the Python and Fortran
285 C_corr[0](im+h, this->j) = mpdata_C_adf<0>(
286 this->psi[this->n], im, this->j, C_unco counterpartstolisting C.24(with theset-up andplottinglogic
287 ); omitted).
288 this->bcy.fill_halos(C_corr[0], ext(this->i,h));
289
290 C_corr[1](this->i, jm+h) = mpdata_C_adf<1>(
291 this->psi[this->n], jm, this->i, C_unco 10gnuplot-iostreamisaheader-onlyC++libraryallowinggnuplottobecon-
292 ); trolled from C++, see http://stahlke.org/dan/gnuplot-iostream/.
293 this->bcx.fill_halos(C_corr[1], ext(this->j,h));
Gnuplot is a portable command-line driven graphing utility, see
294
295 // donor-cell step http://gnuplot.info/
10