Manual

User Manual:

Open the PDF directly: View PDF .
Page Count: 290 [warning: Documents this large are best viewed by clicking the View PDF Link!]

I Using the Pencil Code
II Programming the Pencil Code
III Appendix
IV Indexes

The Pencil Code:

A High-Order MPI code for MHD Turbulence

User’s and Reference Manual

June 17, 2018

http://www.nordita.org/software/pencil-code/

https://github.com/pencil-code/pencil-code

The PENCIL CODE: multi-purpose and multi-user maintained

http://www.nordita.org/~brandenb/talks/misc/PencilCode04.htm

Figure 1: Check-in patterns as a function of time for different subroutines. The different users are marked

by different symbols and different colors.

Contributors to the code

(in inverse alphabetical order according to their user name)

An up to date list of Pencil Code contributors can be found at GitHub.

wladimir.lyra Wladimir Lyra California State University/JPL

weezy S. Louise Wilkin University of Newcastle

wdobler Wolfgang Dobler Potsdam

vpariev Vladimir Pariev University of Rochester

torkel Ulf Torkelsson Chalmers University

tavo.buk Gustavo Guerrero Stanford University

thomas.gastine Thomas Gastine MPI for Solar System Research

theine Tobias (Tobi) Heinemann IAS Princeton

tarek Tarek A. Yousef University of Trondheim

sven.bingert Sven Bingert MPI for Solar System Research

steveb Steve Berukoff UCLA

snod Andrew Snodin University of Newcastle

pkapyla Petri K¨

apyl¨

a University of Helsinki

nils.e.haugen Nils Erland L. Haugen SINTEF

ngrs Graeme R. Sarson University of Newcastle

NBabkovskaia Natalia Babkovskaia University of Helsinki

mreinhardt Matthias Rheinhardt University of Helsinki

mkorpi Maarit J. K¨

apyl¨

a (n´

ee Korpi, Mantere) University of Helsinki

miikkavaisala Miikka V¨

ais¨

al¨

a University of Helsinki

mee Antony (tOnY) Mee University of Newcastle

mcmillan David McMillan York University, Toronto

mattias Mattias Christensson formerly at Nordita

koenkemel Koen Kemel Nordita, Stockholm

karlsson Torgny Karlsson Nordita

joishi Jeff S. Oishi Kavli Institute for Particle Astrophysics

joern.warnecke J¨

orn Warnecke MPI for Solar System Research, Lindau

Iomsn1 Simon Candelaresi University of Dundee, Dundee

fadiesis Fabio Del Sordo Nordita, Stockholm

dorch Bertil Dorch University of Copenhagen

boris.dintrans Boris Dintrans Observatoire Midi-Pyr´

en´

ees, Toulouse

dhruba.mitra Dhrubaditya Mitra Nordita, Stockholm

ccyang Chao-Chin Yang Lund Observatory

christer Christer Sandin University of Uppsala

Bourdin.KIS Philippe Bourdin MPI for Solar System Research

AxelBrandenburg Axel Brandenburg Nordita

apichat Apichat Neamvonk University of Newcastle

amjed Amjed Mohammed University of Oldenburg

alex.i.hubbard Alex Hubbard Am. Museum Nat. History

michiel.lambrechts Michiel Lambrechts Lund Observatory, Lund University

anders Anders Johansen Lund Observatory, Lund University

mppiyali Piyali Chatterjee University of Oslo

2001–2017 Wolfgang Dobler & Axel Brandenburg

Permission is granted to make and distribute verbatim copies of this manual provided

the copyright notice and this permission notice are preserved on all copies.

iii

Permission is granted to copy and distribute modiﬁed versions of this manual under

the conditions for verbatim copying, provided that the entire resulting derived work is

distributed under the terms of a permission notice identical to this one.

License agreement and giving credit

The content of all ﬁles under :pserver:$USER@svn.nordita.org:/var/cvs/brandenb are

under the GNU General Public License (http://www.gnu.org/licenses/gpl.html).

We, the PENCIL CODE community, ask that in publications and presenta-

tions the use of the code (or parts of it) be acknowledged with reference to

the web site http://www.nordita.org/software/pencil-code/ or (equivalently) to.

https://github.com/pencil-code/pencil-code. As a courtesy to the people involved in

developing particularly important parts of the program (use svn annotate src/*.f90 to

ﬁnd out who did what!) we suggest to give appropriate reference to one or several of the

following (or other appropriate) papers (listed here in temporal order):

Dobler, W., Haugen, N. E. L., Yousef, T. A., & Brandenburg, A.: 2003, “Bottleneck ef-

fect in three-dimensional turbulence simulations,” Phys. Rev. E 68, 026304, 1-8

(astro-ph/0303324)

Haugen, N. E. L., Brandenburg, A., & Dobler, W.: 2003, “Is nonhelical hydromag-

netic turbulence peaked at small scales?” Astrophys. J. Lett. 597, L141-L144

(astro-ph/0303372)

Brandenburg, A., K ¨

apyl¨

a, P., & Mohammed, A.: 2004, “Non-Fickian diffusion and

tau-approximation from numerical turbulence,” Phys. Fluids 16, 1020-1027

(astro-ph/0306521)

Johansen, A., Andersen, A. C., & Brandenburg, A.: 2004, “Simulations of dust-

trapping vortices in protoplanetary discs,” Astron. Astrophys. 417, 361-371

(astro-ph/0310059)

Haugen, N. E. L., Brandenburg, A., & Mee, A. J.: 2004, “Mach number dependence

of the onset of dynamo action,” Monthly Notices Roy. Astron. Soc. 353, 947-952

(astro-ph/0405453)

Brandenburg, A., & Multam ¨

aki, T.: 2004, “How long can left and right handed life forms

coexist?” Int. J. Astrobiol. 3, 209-219 (q-bio/0407008)

McMillan, D. G., & Sarson, G. R.: 2005, “Dynamo simulations in a spherical shell of

ideal gas using a high-order Cartesian magnetohydrodynamics code,” Phys. Earth

Planet. Int.153, 124-135

Heinemann, T., Dobler, W., Nordlund, ˚

A., & Brandenburg, A.: 2006, “Radiative transfer

in decomposed domains,” Astron. Astrophys. 448, 731-737 (astro-ph/0503510)

Dobler, W., Stix, M., & Brandenburg, A.: 2006, “Convection and magnetic ﬁeld genera-

tion in fully convective spheres,” Astrophys. J. 638, 336-347 (astro-ph/0410645)

Snodin, A. P., Brandenburg, A., Mee, A. J., & Shukurov, A.: 2006, “Simulating ﬁeld-

aligned diffusion of a cosmic ray gas,” Monthly Notices Roy. Astron. Soc. 373, 643-

652 (astro-ph/0507176)

Johansen, A., Klahr, H., & Henning, T.: 2006, “Dust sedimentation and self-sustained

Kelvin-Helmholtz turbulence in protoplanetary disc mid-planes,” Astrophys. J.

636, 1121-1134 (astro-ph/0512272)

de Val-Borro, M. and 22 coauthors (incl. Lyra, W.): 2006, “A comparative study

of disc-planet interaction,” Monthly Notices Roy. Astron. Soc. 370, 529-558

(astro-ph/0605237)

Johansen, A., Oishi, J. S., Mac Low, M. M., Klahr, H., Henning, T., & Youdin, A.: 2007,

“Rapid planetesimal formation in turbulent circumstellar disks,” Nature 448,

1022–1025 (arXiv/0708.3890)

Lyra, W., Johansen, A., Klahr, H., & Piskunov, N.: 2008, “Global magnetohydrody-

namical models of turbulence in protoplanetary disks I. A cylindrical potential

on a Cartesian grid and transport of solids,” Astron. Astrophys. 479, 883-901

(arXiv/0705.4090)

Brandenburg, A., R ¨

adler, K.-H., Rheinhardt, M., & K¨

apyl¨

a, P. J.: 2008, “Magnetic diffu-

sivity tensor and dynamo effects in rotating and shearing turbulence,” Astrophys.

J. 676, 740-751 (arXiv/0710.4059)

Lyra, W., Johansen, A., Klahr, H., & Piskunov, N.: 2008, “Embryos grown in the dead

zone. Assembling the ﬁrst protoplanetary cores in low-mass selfgravitating cir-

cumstellar disks of gas and solids,” Astron. Astrophys. 491, L41-L44

Lyra, W., Johansen, A., Klahr, H., & Piskunov, N.: 2009, “Standing on the shoulders of

giants. Trojan Earths and vortex trapping in low-mass selfgravitating protoplan-

etary disks of gas and solids,” Astron. Astrophys. 493, 1125-1139

Lyra, W., Johansen, A., Zsom, A., Klahr, H., & Piskunov, N.: 2009, “Planet formation

bursts at the borders of the dead zone in 2D numerical simulations of circumstel-

lar disks,” Astron. Astrophys. 497, 869-888 (arXiv/0901.1638)

Mitra, D., Tavakol, R., Brandenburg, A., & Moss, D.: 2009, “Turbulent dynamos in spher-

ical shell segments of varying geometrical extent,” Astrophys. J. 697, 923-933

(arXiv/0812.3106)

Haugen, N. E. L., & Kragset, S.: 2010, “Particle impaction on a cylinder in a crossﬂow

as function of Stokes and Reynolds numbers,” J. Fluid Mech. 661, 239-261

Rheinhardt, M., & Brandenburg, A.: 2010, “Test-ﬁeld method for mean-ﬁeld coefﬁcients

with MHD background,” Astron. Astrophys. 520, A28 (arXiv/1004.0689)

Babkovskaia, N., Haugen, N. E. L., Brandenburg, A.: 2011, “A high-order public domain

code for direct numerical simulations of turbulent combustion,” J. Comp. Phys.

230, 1-12 (arXiv/1005.5301)

Johansen, A., Klahr, H., & Henning, Th.: 2011, “High-resolution simulations of planetes-

imal formation in turbulent protoplanetary discs,” Astron. Astrophys. 529, A62

Johansen, A., Youdin, A. N., & Lithwick, Y.: 2012, “Adding particle collisions to the for-

mation of asteroids and Kuiper belt objects via streaming instabilities,” Astron.

Astrophys. 537, A125

Lyra, W. & Kuchner, W. : 2013, “Formation of sharp eccentric rings in debris disks with

gas but without planets,” Nature 499, 184–187

Yang, C.-C., & Johansen, A.: 2016, “Integration of Particle-Gas Systems with Stiff Mu-

tual Drag Interaction,” Astrophys. J. Suppl. Series 224, 39

This list is not always up-to-date. We therefore ask the developers to check in new rele-

vant papers, avoiding however redundancies.

Foreword

This code was originally developed at the Turbulence Summer School of the Helmholtz

Institute in Potsdam (2001). While some SPH and PPM codes for hydrodynamics and

magnetohydrodynamics are publicly available, this does not generally seem to be the

case for higher order ﬁnite-difference or spectral codes. Having been approached by peo-

ple interested in using our code, we decided to make it as ﬂexible as possible and as

user-friendly as seems reasonable, and to put it onto a public

CVS

repository. Since

21 September 2008 it is distributed via https://github.com/pencil-code/pencil-code.

The code can certainly not be treated as a black box (no code can), and in order to solve

a new problem in an optimal way, users will need to ﬁnd their own optimal set of pa-

rameters. In particular, you need to be careful in choosing the right values of viscosity,

magnetic diffusivity, and radiative conductivity.

The PENCIL CODE is primarily designed to deal with weakly compressible turbulent

ﬂows, which is why we use high-order ﬁrst and second derivatives. To achieve good par-

allelization, we use explicit (as opposed to compact) ﬁnite differences. Typical scientiﬁc

targets include driven MHD turbulence in a periodic box, convection in a slab with non-

periodic upper and lower boundaries, a convective star embedded in a fully nonperiodic

box, accretion disc turbulence in the shearing sheet approximation, etc. Furthermore,

nonlocal radiation transport, inertial particles, dust coagulation, self-gravity, chemical

reaction networks, and several other physical components are installed, but this num-

ber increases steadily. In addition to Cartesian coordinates, the code can also deal with

spherical and cylindrical polar coordinates.

Magnetic ﬁelds are implemented in terms of the magnetic vector potential to ensure

that the ﬁeld remains solenoidal (divergence-free). At the same time, having the mag-

netic vector potential readily available is a big advantage if one wants to monitor the

magnetic helicity, for example. The code is therefore particularly well suited for all kinds

of dynamo problems.

The code is normally non-conservative; thus, conserved quantities should only be con-

served up to the discretization error of the scheme (not to machine accuracy). There is

no guarantee that a conservative code is more accurate with respect to quantities that

are not explicitly conserved, such as entropy. Another important quantity that is (to

our knowledge) not strictly conserved by ordinary ﬂux conserving schemes is

magnetic

helicity

There are currently no plans to implement adaptive mesh reﬁnement into the code,

which would cause major technical complications. Given that turbulence is generically

space-ﬁlling, local reﬁnement to smaller scales would often not be very useful anyway.

On the other hand, in some geometries turbulence may well be conﬁned to certain re-

gions in space, so one could indeed gain by solving the outer regions with fewer points.

In order to be cache-efﬁcient, we solve the equations along

pencils

in the xdirection.

One very convenient side-effect is that auxiliary and derived variables use very little

memory, as they are only ever deﬁned on one pencil. The domain can be tiled in the y

and zdirections. On multiprocessor computers, the code can use

MPI

(Message Pass-

ing Interface) calls to communicate between processors. An easy switching mechanism

allows the user to run the code on a machine without MPI libraries (e.g. a notebook

computer). Ghost zones are used to implement boundary conditions on physical and

processor boundaries.

vii

A high level of ﬂexibility is achieved by encapsulating individual physical processes

and variables in individual

modules

, which can be switched on or off in the ﬁle

‘Makefile.local’ in the local ‘src’ directory. This approach avoids the use of difﬁcult-

to-read preprocessor directives, at the price of requiring one dummy module for each

physics module. For nonmagnetic hydrodynamics, for example, one will use the module

‘nomagnetic.f90’ and speciﬁes

MAGNETIC = nomagnetic

in ‘Makefile.local’, while for MHD simulations, ‘magnetic.f90’ will be used:

MAGNETIC = magnetic

Note that the term

module

as used here is only loosely related to Fortran modules: both

‘magnetic.f90’ and ‘nomagnetic.f90’ deﬁne an F90 module named Magnetic — this is the

basis of the switching mechanism we are using.

Input parameters (which are set in the ﬁles ‘start.in’, ‘run.in’) can be changed without

recompilation. Furthermore, one can change the list of variables for monitoring (diag-

nostic) output on the ﬂy, and there are mechanisms for making the code reload new

parameters or exit gracefully at runtime. You may want to check for correctness of these

ﬁles with the command pc_configtest.

The requirements for using the Pencil-MPI code are modest: you can use it on any Linux

or Unix system with a

F95

and

compiler suite, like

GNU gcc

and

gfortran

, together

with the shell

CSH

, and the

Perl

interpreter are mandatory requirements.

Although the PENCIL CODE is mainly designed to run on supercomputers, more than

50% of the users run their code also on Macs, and the other half uses either directly

Linux on their laptops or they use VirtualBox on their Windows machine on which they

install Ubuntu Linux. If you have

IDL

as well, you will be able to visualize the re-

sults (a number of sample procedures are provided), but other tools such as

Python

(OpenDX, data explorer) can also be used and some relevant tools and routines come

with the PENCIL CODE.

If you want to make creative use of the code, this manual will contain far too little in-

formation. Its major aim is to give you an idea of the way the code is organized, so you

can more efﬁciently read the source code, which contains a reasonable amount of com-

ments. You might want to read through the various sample directories that are checked

in. Choose one that is closest to your application and start modifying. For further en-

hancements that you may want to add to the code, you can take as an example the lines

in the code that deal with related variables, functions, diagnostics, equations etc., which

have already been implemented. Just remember: grep is one of your best friends when

you want to understand how certain variables or functions are used in the code.

We will be happy to include user-supplied changes and updates to the code in future

releases and welcome any feedback.

wdobler@gmail.com Potsdam

AxelBrandenburg@gmail.com Stockholm

viii

Acknowledgments

Many people have contributed in different ways to the development of this code. We

thank ﬁrst of all ˚

Ake Nordlund (Copenhagen Observatory) and Bob Stein (University of

Michigan) who introduced us to the idea of using high-order schemes in compressible

ﬂows and who taught us a lot about simulations in general.

The calculation of the power spectra, structure functions, the remeshing procedures,

routines for changing the number of processors, as well as the shearing sheet approxi-

mation and the ﬂux-limited diffusion approximation for radiative transfer were imple-

mented by Nils Erland L. Haugen (University of Trondheim). Tobi Heinemann added

the long characteristics method for radiative transfer as well as hydrogen ionization.

He also added and/or improved shock diffusion for other variables and improved the

resulting timestep control. Anders Johansen, Wladimir (Wlad) Lyra, and Jeff Oishi con-

tributed to the implementation of the dust equations (which now comprises an array of

different components). Antony (Tony) Mee (University of Newcastle) implemented shock

viscosity and added the interstellar module together with Graeme R. Sarson (also Uni-

versity of Newcastle), who also implemented the geodynamo set-up together with David

McMillan (currently also at the University of Newcastle). Tony also included a method

for outputting auxiliary variables and enhanced the overall functionality of the code

and related idl and dx procedures. He also added, together with Andrew Snodin, the

evolution equations for the cosmic ray energy density. Vladimir Pariev (University of

Rochester) contributed to the development and testing of the potential ﬁeld boundary

condition at an early stage. The implementation of spherical and cylindrical coordinates

is due to Dhrubaditya (Dhruba) Mitra and Wladimir Lyra. Wlad also implemented the

global set-up for protoplanetary disks (as opposed to the local shearing sheet formalism).

He also added a N-body code (based on the particle module coded by Anders Johansen

and Tony), and implemented the coupled evolution equations of neutrals and ions for

two-ﬂuid models of ambipolar diffusion. Boris Dintrans is in charge of implementing the

anelastic and Boussinesq modules. Philippe-A. Bourdin implemented HDF5 support and

wrote the optional IO-modules for high-performance computing featuring various com-

munication strategies. He also contributed to the solar-corona module and worked on

the IDL GUI, including the IDL routines for reading and working with large amounts of

data. Again, this list contains other recent items that are not yet fully documented and

acknowledged.

Use of the PPARC supported supercomputers in St Andrews (Mhd) and Leicester (Ukaff)

is acknowledged. We also acknowledge the Danish Center for Scientiﬁc Computing for

granting time on Horseshoe, which is a 512+140 processor Beowulf cluster in Odense

(Horseshoe).

Contents

I Using the PENCIL CODE 1

1 System requirements 1

2 Obtaining the code 2

2.1 Obtaining the code via git or svn . . . . . . . . . . . . . . . . . . . . . . . . 2

2.2 Updating via svn or git . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2

2.3 Getting the last validated version . . . . . . . . . . . . . . . . . . . . . . . . 3

2.4 Getting older versions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4

3 Getting started 5

3.1 Setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

3.1.1 Environment settings . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

3.1.2 Linking scripts and source ﬁles . . . . . . . . . . . . . . . . . . . . . 6

3.1.3 Adapting ‘Makefile.src’ ......................... 6

3.1.4 Running make ............................... 6

3.1.5 Choosing a data directory . . . . . . . . . . . . . . . . . . . . . . . . 7

3.1.6 Running the code . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

3.2 Further tests . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9

4 Code structure 11

4.1 Directory tree . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

4.2 Basic concepts . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12

4.2.1 Data access in pencils . . . . . . . . . . . . . . . . . . . . . . . . . . 12

4.2.2 Modularity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

4.3 Files in the run directories . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

4.3.1 ‘start.in’, ‘run.in’, ‘print.in’ . . . . . . . . . . . . . . . . . . . . . . 14

4.3.2 ‘datadir.in’ ................................ 14

4.3.3 ‘reference.out’ .............................. 14

4.3.4 ‘start.csh’, ‘run.csh’, ‘getconf.csh’ [obsolete; see Sect. 5.1] . . . . . 14

4.3.5 ‘src/ ’.................................... 14

4.3.6 ‘data/ ’ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

5 Using the code 17

5.1 Conﬁguring the code to compile and run on your computer . . . . . . . . . 17

5.1.1 Locating the conﬁguration ﬁle . . . . . . . . . . . . . . . . . . . . . . 17

5.1.2 Structure of conﬁguration ﬁles . . . . . . . . . . . . . . . . . . . . . 18

5.1.3 Compiling the code . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

5.1.4 Running the code . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

5.1.5 Testing the code . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

5.2 Adapting ‘Makefile.src’ [obsolete; see Sect. 5.1] . . . . . . . . . . . . . . . 21

5.3 Changing the resolution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22

5.4 Using a non-equidistant grid . . . . . . . . . . . . . . . . . . . . . . . . . . 23

5.5 Diagnostic output . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25

5.6 Data ﬁles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25

5.6.1 Snapshot ﬁles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25

5.7 Video ﬁles and slices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26

5.8 Averages . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29

5.8.1 One-dimensional output averaged in two dimensions . . . . . . . . 29

5.8.2 Two-dimensional output averaged in one dimension . . . . . . . . . 29

5.8.3 Azimuthal averages . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29

5.8.4 Time averages . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30

5.9 Helper scripts . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31

5.10 RELOAD and STOP ﬁles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33

5.11 RERUN and NEWDIR ﬁles . . . . . . . . . . . . . . . . . . . . . . . . . . . 34

5.12 Start and run parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34

5.13 Physical units . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36

5.14 Minimum amount of viscosity . . . . . . . . . . . . . . . . . . . . . . . . . . 37

5.15 The time step . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38

5.15.1 The usual RK-2N time step . . . . . . . . . . . . . . . . . . . . . . . 38

5.15.2 The Runge-Kutta-Fehlberg time step . . . . . . . . . . . . . . . . . . 38

5.16 Boundary conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39

5.16.1 Where to specify boundary conditions . . . . . . . . . . . . . . . . . 39

5.16.2 How to specify boundary conditions . . . . . . . . . . . . . . . . . . 39

5.17 Restarting a simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40

5.18 One- and two-dimensional runs . . . . . . . . . . . . . . . . . . . . . . . . . 41

5.19 Visualization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41

5.19.1 Gnuplot . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41

5.19.2 Data explorer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41

5.19.3 GDL . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42

5.19.4 IDL . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43

5.19.5 Python . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46

5.20 Running on multi-processor computers . . . . . . . . . . . . . . . . . . . . . 49

5.20.1 How to run a sample problem in parallel . . . . . . . . . . . . . . . 49

5.20.2 Hierarchical networks (e.g. on Beowulf clusters) . . . . . . . . . . . 50

5.20.3 Extra workload caused by the ghost zones . . . . . . . . . . . . . . . 50

5.21 Running in double-precision . . . . . . . . . . . . . . . . . . . . . . . . . . . 51

5.22 Power spectrum . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52

5.23 Structure functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54

5.24 Particles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55

5.24.1 Particles in parallel . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56

5.24.2 Large number of particles . . . . . . . . . . . . . . . . . . . . . . . . 58

5.24.3 Random number generator . . . . . . . . . . . . . . . . . . . . . . . 58

5.25 Non-cartesian coordinate systems . . . . . . . . . . . . . . . . . . . . . . . . 59

6 The equations 60

6.1 Continuity equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60

6.2 Equation of motion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60

6.3 Induction equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61

6.4 Entropy equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61

6.4.1 Viscous heating . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62

6.4.2 Alternative description . . . . . . . . . . . . . . . . . . . . . . . . . . 62

6.5 Transport equation for a passive scalar . . . . . . . . . . . . . . . . . . . . 63

6.6 Bulk viscosity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63

6.6.1 Shock viscosity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63

6.7 Equation of state . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63

6.8 Ionization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64

6.8.1 Ambipolar diffusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65

6.9 Radiative transfer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66

6.10 Self-gravity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67

6.11 Incompressible and anelastic equations . . . . . . . . . . . . . . . . . . . . 67

6.12 Dust equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67

6.13 Cosmic ray pressure in diffusion approximation . . . . . . . . . . . . . . . 68

6.14 Particles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69

6.14.1 Tracer particles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69

6.14.2 Dust particles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69

6.15 N-body solver . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70

6.16 Test-ﬁeld equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71

6.17 Gravitational wave equations . . . . . . . . . . . . . . . . . . . . . . . . . . 71

7 Troubleshooting / Frequently Asked Questions 75

7.1 Download and setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75

7.1.1 Download forbidden . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75

7.1.2 Shell gives error message when sourcing ‘sourceme.X’ . . . . . . . . 75

7.2 Compilation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76

7.2.1 Problems compiling syscalls . . . . . . . . . . . . . . . . . . . . . . . 76

7.2.2 Unable to open include ﬁle: chemistry.h . . . . . . . . . . . . . . . . 76

7.2.3 Compiling with

ifc

under Linux . . . . . . . . . . . . . . . . . . . . . 76

7.2.4 Segmentation fault with

ifort

8.0 under Linux . . . . . . . . . . . . 77

7.2.5 The underscore problem: linking with

MPI

. . . . . . . . . . . . . . 77

7.2.6 Compilation stops with the cryptic error message: . . . . . . . . . . 77

7.2.7 The code doesn’t compile, . . . . . . . . . . . . . . . . . . . . . . . . . 78

7.2.8 Some samples don’t even compile, . . . . . . . . . . . . . . . . . . . 78

7.2.9 Internal compiler error with Compaq/Dec F90 . . . . . . . . . . . . 79

7.2.10 Assertion failure under SunOS . . . . . . . . . . . . . . . . . . . . . 79

7.2.11 After some dirty tricks I got pencil code to compile with MPI, ... . . 80

7.2.12 Error: Symbol ’mpi comm world’ at (1) has no IMPLICIT type . . . 80

7.2.13 Error: Can’t open included ﬁle ’mpif.h’ . . . . . . . . . . . . . . . . . 81

7.3 Pencil check . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81

7.3.1 The pencil check complains for no reason. . . . . . . . . . . . . . . . 81

7.3.2 The pencil check reports MISSING PENCILS and quits . . . . . . 81

7.3.3 The pencil check reports unnecessary pencils . . . . . . . . . . . . . 81

7.3.4 The pencil check reports that most or all pencils are missing . . . . 81

7.3.5 Running the pencil check triggers mathematical errors in the code 82

7.3.6 The pencil check still complains . . . . . . . . . . . . . . . . . . . . . 82

7.3.7 The pencil check is annoying so I turned it off . . . . . . . . . . . . 82

7.4 Running . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82

7.4.1 Periodic boundary conditions in ‘start.x’ . . . . . . . . . . . . . . . 82

7.4.2 csh problem? . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82

7.4.3 ‘run.csh’ doesn’t work: . . . . . . . . . . . . . . . . . . . . . . . . . . 83

7.4.4 Code crashes after restarting . . . . . . . . . . . . . . . . . . . . . . 83

7.4.5 auto-test gone mad...? . . . . . . . . . . . . . . . . . . . . . . . . . . 83

7.4.6 Can I restart with a different number of cpus? . . . . . . . . . . . . 84

7.4.7 Can I restart with a different number of cpus? . . . . . . . . . . . . 84

7.4.8 fft xyz parallel 3D: nygrid needs to be an integer multiple... . . . . 84

7.4.9 Unit-agnostic calculations? . . . . . . . . . . . . . . . . . . . . . . . 85

7.5 Visualization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86

xii

7.5.1 ‘start.pro’ doesn’t work: . . . . . . . . . . . . . . . . . . . . . . . . . 86

7.5.2 ‘start.pro’ doesn’t work: . . . . . . . . . . . . . . . . . . . . . . . . . 86

7.5.3 Something about tag name undeﬁned: . . . . . . . . . . . . . . . . . 86

7.5.4 Something INC in start.pro . . . . . . . . . . . . . . . . . . . . . . . 86

7.5.5 nl2idl problem when reading param2.nml . . . . . . . . . . . . . . . 87

7.5.6 Spurious dots in the time series ﬁle . . . . . . . . . . . . . . . . . . 87

7.5.7 Problems with pc_varcontent.pro . . . . . . . . . . . . . . . . . . . 87

7.6 General questions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88

7.6.1 “Installation” procedure . . . . . . . . . . . . . . . . . . . . . . . . . 88

7.6.2 Small numbers in the code . . . . . . . . . . . . . . . . . . . . . . . . 88

7.6.3 Why do we need a /lphysics/ namelist in the ﬁrst place? . . . . . . 89

7.6.4 Can I run the code on a Mac? . . . . . . . . . . . . . . . . . . . . . . 90

7.6.5 Pencil Code discussion forum . . . . . . . . . . . . . . . . . . . . . . 90

7.6.6 The manual . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90

II Programming the PENCIL CODE 91

8 Understanding the code 95

8.1 Example: how is the continuity equation being solved? . . . . . . . . . . . 95

9 Adapting the code 97

9.1 The PENCIL CODE coding standard . . . . . . . . . . . . . . . . . . . . . . . 97

9.2 Adding new output diagnostics . . . . . . . . . . . . . . . . . . . . . . . . . 98

9.3 The f-array . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100

9.4 The df-array . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100

9.5 The fp-array . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101

9.6 The pencil case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101

9.6.1 Pencil check . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102

9.6.2 Adding new pencils . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103

9.7 Adding new physics: the Special module . . . . . . . . . . . . . . . . . . . . 103

9.8 Adding switchable modules . . . . . . . . . . . . . . . . . . . . . . . . . . . 104

9.9 Adding your initial conditions: the InitialCondition module . . . . . . . . . 104

10 Testing the code 106

10.1 How to set up periodic tests (auto-tests) . . . . . . . . . . . . . . . . . . . . 106

11 Useful internals 108

11.1 Global variables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108

11.2 Subroutines and functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108

III Appendix 111

A Timings 111

A.1 Test case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118

A.2 Running the code . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119

A.3 Triolith . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119

A.4 Lindgren . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119

B Coding standard 122

xiii

B.1 File naming conventions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122

B.2 Fortran Code . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122

B.2.1 Indenting and whitespace . . . . . . . . . . . . . . . . . . . . . . . . 122

B.2.2 Comments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123

B.2.3 Module names . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124

B.2.4 Variable names . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124

B.2.5 Emacs settings . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125

B.3 Other best practices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126

B.4 General changes to the code . . . . . . . . . . . . . . . . . . . . . . . . . . . 126

C Some speciﬁc initial conditions 127

C.1 Random velocity or magnetic ﬁelds . . . . . . . . . . . . . . . . . . . . . . . 127

C.2 Turbulent initial with given spectrum . . . . . . . . . . . . . . . . . . . . . 127

C.3 Beltrami ﬁelds . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128

C.4 Magnetic ﬂux rings: initaa=’fluxrings’ . . . . . . . . . . . . . . . . . . . 128

C.5 Vertical stratiﬁcation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129

C.5.1 Isothermal atmosphere . . . . . . . . . . . . . . . . . . . . . . . . . . 129

C.5.2 Polytropic atmosphere . . . . . . . . . . . . . . . . . . . . . . . . . . 130

C.5.3 Changing the stratiﬁcation . . . . . . . . . . . . . . . . . . . . . . . 131

C.5.4 The Rayleigh number . . . . . . . . . . . . . . . . . . . . . . . . . . . 131

C.5.5 Entropy boundary condition . . . . . . . . . . . . . . . . . . . . . . . 132

C.5.6 Temperature boundary condition at the top . . . . . . . . . . . . . . 132

C.6 Potential-ﬁeld boundary condition . . . . . . . . . . . . . . . . . . . . . . . 132

C.7 Planet solution in the shearing box . . . . . . . . . . . . . . . . . . . . . . . 134

D Some speciﬁc boundary conditions 135

D.1 Perfect-conductor boundary condition . . . . . . . . . . . . . . . . . . . . . 135

D.2 Stress-free boundary condition . . . . . . . . . . . . . . . . . . . . . . . . . 135

D.3 Normal-ﬁeld-radial boundary condition . . . . . . . . . . . . . . . . . . . . 136

E High-frequency ﬁlters 137

E.1 Conservative hyperdissipation . . . . . . . . . . . . . . . . . . . . . . . . . . 137

E.2 Hyperviscosity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 139

E.2.1 Conservative case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 139

E.2.2 Non-conservative cases . . . . . . . . . . . . . . . . . . . . . . . . . . 140

E.2.3 Choosing the coefﬁcient . . . . . . . . . . . . . . . . . . . . . . . . . 141

E.2.4 Turbulence with hyperviscosity . . . . . . . . . . . . . . . . . . . . . 141

E.3 Anisotropic hyperdissipation . . . . . . . . . . . . . . . . . . . . . . . . . . . 142

E.4 Hyperviscosity in Burgers shock . . . . . . . . . . . . . . . . . . . . . . . . 142

F Special techniques 144

F.1 After changing

REAL PRECISION

. . . . . . . . . . . . . . . . . . . . . . . 144

F.2 Remeshing (regridding) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 144

F.3 Restarting from a run with less physics . . . . . . . . . . . . . . . . . . . . 145

F.4 Restarting with particles from a run without them . . . . . . . . . . . . . . 146

G Runs and reference data 148

G.1 Shock tests . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 148

G.1.1 Sod shock tube problem . . . . . . . . . . . . . . . . . . . . . . . . . 148

G.1.2 Temperature jump . . . . . . . . . . . . . . . . . . . . . . . . . . . . 148

G.2 Random forcing function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 148

xiv

G.3 Three-layered convection model . . . . . . . . . . . . . . . . . . . . . . . . . 149

G.4 Magnetic helicity in the shearing sheet . . . . . . . . . . . . . . . . . . . . 150

H Numerical methods 154

H.1 Sixth-order spatial derivatives . . . . . . . . . . . . . . . . . . . . . . . . . 154

H.2 Upwind derivatives to avoid ‘wiggles’ . . . . . . . . . . . . . . . . . . . . . . 155

H.3 The bidiagonal scheme for cross-derivatives . . . . . . . . . . . . . . . . . . 156

H.4 The 2N-scheme for time-stepping . . . . . . . . . . . . . . . . . . . . . . . . 157

H.5 Diffusive error from the time-stepping . . . . . . . . . . . . . . . . . . . . . 158

H.6 Ionization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 159

H.7 Radiative transfer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 160

H.7.1 Solving the radiative transfer equation . . . . . . . . . . . . . . . . 160

H.7.2 Angular integration . . . . . . . . . . . . . . . . . . . . . . . . . . . . 161

I Curvilinear coordinates 163

I.1 Covariant derivatives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 163

I.2 Differential operators . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 163

I.2.1 Gradient . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 163

I.2.2 Divergence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 164

I.2.3 Curl . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 164

I.2.4 Advection operator . . . . . . . . . . . . . . . . . . . . . . . . . . . . 165

I.2.5 Mixed advection operator . . . . . . . . . . . . . . . . . . . . . . . . 165

I.2.6 Shear term . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 166

I.2.7 Another mixed advection operator . . . . . . . . . . . . . . . . . . . 166

I.2.8 Strain Matrix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 166

I.2.9 Lambda effect . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 166

I.2.10 Laplacian of a scalar . . . . . . . . . . . . . . . . . . . . . . . . . . . 167

I.2.11 Hessian of a scalar . . . . . . . . . . . . . . . . . . . . . . . . . . . . 167

I.2.12 Double curl . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 168

I.2.13 Gradient of a divergence . . . . . . . . . . . . . . . . . . . . . . . . . 169

J Switchable modules 171

K Startup and run-time parameters 172

K.1 Startup parameters for ‘start.in’ . . . . . . . . . . . . . . . . . . . . . . . . 172

K.2 Runtime parameters for ‘run.in’ . . . . . . . . . . . . . . . . . . . . . . . . 179

K.3 Parameters for ‘print.in’ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 186

K.4 Parameters for ‘video.in’ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 222

K.5 Parameters for ‘phiaver.in’ . . . . . . . . . . . . . . . . . . . . . . . . . . . 224

K.6 Parameters for ‘xyaver.in’ . . . . . . . . . . . . . . . . . . . . . . . . . . . . 225

K.7 Parameters for ‘xzaver.in’ . . . . . . . . . . . . . . . . . . . . . . . . . . . . 231

K.8 Parameters for ‘yzaver.in’ . . . . . . . . . . . . . . . . . . . . . . . . . . . . 232

K.9 Parameters for ‘yaver.in’ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 235

K.10 Parameters for ‘zaver.in’ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 237

K.11 Boundary conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 240

K.11.1 Boundary condition

bcx

. . . . . . . . . . . . . . . . . . . . . . . . . 240

K.11.2 Boundary condition

bcy

. . . . . . . . . . . . . . . . . . . . . . . . . 243

K.11.3 Boundary condition

bcz

. . . . . . . . . . . . . . . . . . . . . . . . . 245

K.12 Initial condition parameter dependence . . . . . . . . . . . . . . . . . . . . 249

IV Indexes 253

xvi

Part I

Using the PENCIL CODE

1 System requirements

To use the code, you will need the following:

1. Absolutely needed:

•

F95

compiler

•

compiler

2. Used heavily (if you don’t have one of these, you will need to adjust many things

manually):

•a

Unix

Linux

-type system with

make

and

csh

•

Perl

(remember: if it doesn’t run Perl, it’s not a computer)

3. The following are dispensable, but enhance functionality in one way or the other:

•an

MPI

implementation (for parallelization on multiprocessor systems)

•

alias

OpenDX

data explorer

(for 3-D visualization of results)

•

IDL

(for visualization of results; the 7-minute demo license will do for many

applications)

2 THE PENCIL CODE

2 Obtaining the code

The code is now distributed via https://github.com/pencil-code/pencil-code, where

you can either download a tarball, or, preferably, download it via

svn

git

. In Iran and

some other countries, GitHub is not currently available. To alleviate this problem, we

have made a recent copy available on http://www.nordita.org/software/pencil-code/.

If you want us to update this tarball, please contact us.

To ensure at least some level of stability of the

svn/git

versions, a set of test problems

(listed in ‘$PENCIL_HOME/bin/auto-test’) are routinely tested. This includes all problems

in ‘$PENCIL_HOME/samples’. See Sect. 10 for details.

2.1 Obtaining the code via git or svn

1. Many machines have

svn

installed (try svn -v or which svn). On Ubuntu, for ex-

ample,

svn

comes under the package name subversion.

2. The code is now saved under Github, git can be obtained in Linux by typing sudo

apt-get install git

3. Unless you are a privileged users with write access, you can download the code

with the command

git clone https://github.com/pencil-code/pencil-code.git

svn checkout https://github.com/pencil-code/pencil-code/trunk/ ...\\

pencil-code --username MY_GITHUB_USERNAME

In order to push your changes to the repository, you have to ask the maintainer of

pencil code for push access (to become a contributor), or put a pull request to the

maintainer of the code.

Be sure to run auto-test before you check anything back in again. It can be very

annoying for someone else to ﬁgure out what’s wrong, especially if you are just up

to something else. At the very least, you should do

pc_auto-test --level=0 --no-pencil-check -C

This allows you to run just 2 of the most essential tests starting with all the no-

modules and then most-modules.

2.2 Updating via svn or git

Independent of how you installed the code in the ﬁrst place (from tarball or via

svn/git

you can update your version using

svn/git

. If you have done nontrivial alterations to

your version of the code, you ought to be careful about upgrading: although

svn/git

is an

excellent tool for distributed programming, conﬂicts are quite possible, since many of us

are going to touch many parts of the code while we develop it further. Thus, despite the

2.3 Getting the last validated version 3

fact that the code is under

svn/git

, you should probably back up your important changes

before upgrading.

Here is the upgrading procedure for

git

1. Perform a git update of the tree:

unix> git pull

2. Fix any conﬂicts you encounter and make sure the examples in the directory

‘samples/’ are still working.

Here is the upgrading procedure for

svn

1. Perform a svn update of the tree:

unix> pc_svnup

2. Fix any conﬂicts you encounter and make sure the examples in the directory

‘samples/’ are still working.

If you have made useful changes, please contact one of the (currently) 10 “Contributors”

(listed under https://github.com/pencil-code/pencil-code) who can give you push or

check-in permission. Be sure to have sufﬁcient comments in the code and please follow

our standard coding conventions explained in Section 9.1. There is also a script to check

and ﬁx the most common stylebreaks, pc codingstyle.

2.3 Getting the last validated version

The script pc_svnup accepts arguments -val or -validated, which means that the current

changes on a user’s machine will be merged into the last working version. This way

every user can be sure that any problems with the code must be due to the current

changes done by this user since the last check-in.

Examples:

unix> pc_svnup -src -s -validated

brings all ﬁles in ‘$PENCIL_HOME/src’ to the last validated status, and merges all your

changes into this version. This allows you to work with this, but in order to check in

your changes you have to update everything to the most recent status ﬁrst, i.e.

unix> pc_svnup -src

Your own changes will be merged into this latest version as before.

NOTE: The functionality of the head of the trunk should be preserved at all times. How-

ever, accidents do happen. For the beneﬁt of all other developers, any errors should

be corrected within 1-2 hours. This is the reason why the code comes with a ﬁle

‘pencil-code/license/developers.txt’, which should contain contact details of all de-

velopers. The pc_svnup -val option allows all other people to stay away from any trou-

ble.

4 THE PENCIL CODE

2.4 Getting older versions

You may ﬁnd that the latest

svn

version of the code produces errors. If you have reasons

to believe that this is due to changes introduced on 27 November 2008 (to give an ex-

ample), you can check out the version prior to this by specifying a revision number with

svn update -r #####. One reason why one cannot always reproduce exactly the same

situation too far back in time is connected with the fact that processor architecture and

the compiler were different, resulting e.g. in different rounding errors.

3. Getting started 5

3 Getting started

To get yourself started, you should run one or several examples which are provided in

one of the ‘samples/’ subdirectories. Note that you will only be able to fully assess the

numerical solutions if you visualize them with

IDL

or other tools (see Sect. 5.19).

3.1 Setup

3.1.1 Environment settings

The functionality of helper scripts and IDL routines relies on a few environment vari-

ables being set correctly. The simplest way to achieve this is to go to the top directory of

the code and source one of the two scripts ‘sourceme.csh’ or ‘sourceme.sh’ (depending on

the type of shell you are using):

csh> cd pencil-code

csh> source ./sourceme.csh

for

tcsh

csh

users; or

sh> cd pencil-code

sh> . ./sourceme.sh

for users of

bash

Bourne shell

, or similar shells. You should get output similar to

PENCIL_HOME = </home/dobler/f90/pencil-code>

Adding /home/dobler/f90/pencil-code/bin to PATH

Apart from the PATH variable, the environment variable IDL_PATH is set to something

like ./idl:../idl:+$PENCIL_HOME/idl:./data:<IDL_DEFAULT> .

Note 1 The <IDL_DEFAULT> mechanism does not work for IDL versions 5.2 or older. In

this case, you will have to edit the path manually, or adapt the ‘sourceme’ scripts.

Note 2 If you don’t want to rely on the ‘sourceme’ scripts’ (quite heuristic) ability to cor-

rectly identify the code’s main directory, you can set the environment variable PENCIL_-

HOME explicitly before you run the source command.

Note 3 Do not just source the ‘sourceme’ script from your shell startup ﬁle (‘~/.cshrc’

or ‘~/.bashrc’, because it outputs a few lines of diagnostics for each sub-shell, which will

break many applications. To suppress all output, follow the instructions given in the

header documentation of ‘sourceme.csh’ and ‘sourceme.sh’. Likewise, output from other

ﬁles invoked by source should also be suppressed.

Note 4 The second time you source ‘sourceme’, it will not add anything to your PATH

variable. This is on purpose to avoid cluttering of your environment: you can source the

ﬁle as often as you like (in your shell startup script, then manually and in addition in

some script you have written), without thinking twice. If, however, at the ﬁrst sourcing,

the setting of PENCIL_HOME was wrong, this mechanism would keep you from ever adding

6 THE PENCIL CODE

the right directory to your PATH. In this case, you need to ﬁrst undeﬁne the environment

variable PENCIL_HOME:

csh> unsetenv PENCIL_HOME

csh> source ./sourceme.csh

sh> unset PENCIL_HOME

sh> . ./sourceme.sh

3.1.2 Linking scripts and source ﬁles

With your environment set up correctly, you can now go to the directory you want to

work in and set up subdirectories and links. This is accomplished by the script ‘pc_-

setupsrc’, which is located in ‘$PENCIL_HOME/bin’ and is thus now in your executable

path.

For concreteness, let us assume you want to use ‘samples/conv-slab’ as your

run direc-

tory

, i.e. you want to run a three-layer slab model of solar convection. You then do the

following:

unix> cd samples/conv-slab

unix> pc_setupsrc

src already exists

2 files already exist in src

The script has linked a number of scripts from ‘$PENCIL_HOME/bin’, generated a directory

‘src’ for the source code and linked the Fortran source ﬁles (plus a few more ﬁles) from

‘$PENCIL_HOME/src’ to that directory:

unix> ls -F

reference.out src/

start.csh@ run.csh@ getconf.csh@

start.in run.in print.in

3.1.3 Adapting ‘Makefile.src’

This step requires some input from you, but you only have to do this once for each

machine you want to run the code on. See Sect. 5.2 for a description of the steps you

need to take here.

Note: If you are lucky and use compilers similar to the ones we have, you may be able

to skip this step; but blame yourself if things don’t compile, then. If not, you can run

make with explicit ﬂags, see Sect. 5.2 and in particular Table 1.

3.1.4 Running make

Next, you run make in the ‘src’ subdirectory of your run directory. Since you are

using one of the predeﬁned test problems, the settings in ‘src/Makefile.local’ and

‘src/cparam.local’ are all reasonable, and you just do

3.1 Setup 7

unix> make

If you have set up the compiler ﬂags correctly, compilation should complete successfully.

3.1.5 Choosing a data directory

The code will by default write data like snapshot ﬁles to the subdirectory ‘data’ of the

run directory. Since this will involve a large volume of IO-operations (at least for large

grid sizes), one will normally try to avoid writing the data via NFS. The recommended

way to set up a ‘data’ data directory is to generate a corresponding directory on the local

disc of the computer you are running on and (soft-)link it to ‘./data’. Even if the link is

part of an NFS directory, all the IO operations will be local. For example, if you have a

local disc ‘/scratch’, you can do the following:

unix> mkdir -p /scratch/$USER/pencil-data/samples/conv-slab

unix> ln -s /scratch/$USER/pencil-data/samples/conv-slab ./data

This is done automatically by the pc_mkdatadir command which, in turn, is invoked

when making a new run directory with the pc_newrun command, for example.

Even if you don’t have an NFS-mounted directory (say, on your notebook computer), it

is probably still a good idea to have code and data well separated by a scheme like the

one described above.

An alternative to symbolic links, is to provide a ﬁle called ‘datadir.in’ in the root of

the run directory. This ﬁle should contain one line of text specifying the absolute or

relative data directory path to use. This facility is useful if one wishes to switch one run

directory between different data directories. It is suggested that in such cases symbolic

links are again made in the run directory to the various locations, then the ‘datadir.in’

need contain only a short relative path.

3.1.6 Running the code

You are now ready to start the code:

unix> start.csh

Linux cincinnatus 2.4.18-4GB #1 Wed Mar 27 13:57:05 UTC 2002 i686 unknown

Non-MPI version

datadir = data

Fri Aug 8 21:36:43 CEST 2003

src/start.x

CVS: io_dist.f90 v. 1.61 (brandenb ) 2003/08/03 09:26:55

[. . . ]

CVS: start.in v. 1.4 (dobler ) 2002/10/02 20:11:14

nxgrid,nygrid,nzgrid= 32 32 32

thermodynamics: assume cp=1

uu: up-down

piecewise polytropic vertical stratification (lnrho)

init_lnrho: cs2bot,cs2top= 1.450000 0.3333330

e.g. for ionization runs: cs2bot,cs2top not yet set

piecewise polytropic vertical stratification (ss)

8 THE PENCIL CODE

start.x has completed successfully

0.070u 0.020s 0:00.14 64.2% 0+0k 0+0io 180pf+0w

Fri Aug 8 21:36:43 CEST 2003

This runs ‘src/start.x’ to construct an initial condition based on the parameters

set in ‘start.in’. This initial condition is stored in ‘data/proc0/var.dat’ (and in

‘data/proc1/var.dat’, etc. if you run the multiprocessor version). It is fair to say that

this is now a rather primitive routine; see ‘pencil-code/idl/read’ for various reading

routines. You can then visualize the data using standard idl language.

If you visualize the proﬁles using

IDL

(see below), the result should bear some resem-

blance to Fig. 2, but with different values in the ghost zones (the correct values are set

at run-time only) and a simpler velocity proﬁle.

Now we run the code:

unix> run.csh

This executes ‘src/run.x’ and carries out

time steps, where

and other run-time

parameters are speciﬁed in ‘run.in’. On a decent PC (1.7 GHz), 50 time steps take about

10 seconds.

The relevant part of the code’s output looks like

--it----t-------dt-------urms----umax----rhom------ssm-----dtc----dtu---dtnu---dtchi-

0 0.34 6.792E-03 0.0060 0.0452 14.4708 -0.4478 0.978 0.013 0.207 0.346

10 0.41 6.787E-03 0.0062 0.0440 14.4707 -0.4480 0.978 0.013 0.207 0.345

20 0.48 6.781E-03 0.0064 0.0429 14.4705 -0.4481 0.977 0.012 0.207 0.345

30 0.54 6.777E-03 0.0067 0.0408 14.4703 -0.4482 0.977 0.012 0.207 0.345

40 0.61 6.776E-03 0.0069 0.0381 14.4702 -0.4482 0.977 0.011 0.207 0.346

and lists

1. the number

of the current time step;

2. the time,

;

3. the time step,

;

4. the rms velocity,

urms

=phu2i;

5. the maximum velocity,

umax

= max |u|;

6. the mean density,

rhom

=hρi;

7. the mean entropy,

ssm

=hsi/cp;

8. the time step in units of the acoustic Courant step,

dtc

=δt/(cs0δxmin);

9. the time step in units of the advective time step,

dtu

=δt/(cδt δx/ max |u|);

10. the time step in units of viscous time step,

dtnu

=δt/(cδt,vδx2/νmax);

11. the time step in units of the conductive time step,

dtchi

=δt/(cδt,vδx2/χmax).

The entries in this list can be added, removed or reformatted in the ﬁle ‘print.in’, see

Sects 5.5 and K.3. The output is also saved in ‘data/time_series.dat’ and should be

identical to the content of ‘reference.out’.

3.2 Further tests 9

ln ρ

0 1 2 3

−0.5

0.0

0.5

1.0

1.5

−0.02 0.00 0.02 0.04

−0.5

0.0

0.5

1.0

1.5

Entropy s

−0.6−0.4−0.2 0.0 0.2

−0.5

0.0

0.5

1.0

1.5

Temperature T

0.5 1.0 1.5 2.0

−0.5

0.0

0.5

1.0

1.5

Figure 2: Stratiﬁcation of the three-layer convection model in ‘samples/conv-slab’ after 50 timesteps

(t= 0.428). Shown are (from left to right) density ρ, vertical velocity uz, entropy s/cpand temperature

Tas functions of the vertical coordinate zfor about ten different vertical lines in the computational box.

The dashed lines denote domain boundaries: z < −0.68 is the lower ghost zone (points have no physical

signiﬁcance); −0.68 < z < 0is a stably stratiﬁed layer (ds/dz > 0); 0< z < 1is the unstable layer

(ds/dz < 0); 1< z < 1.32 is the isothermal top layer; z > 1.32 is the upper ghost zone (points have no

physical signiﬁcance).

If you have

IDL

, you can visualize the stratiﬁcation with (see Sect. 5.19.4 for details)

unix > idl

IDL > pc_read_var,obj=var,/trimall

IDL > tvscl,var,uu(*,*,0,0)

which shows uxin the xy plane through the ﬁrst meshpoint in the zdirection. There

have been some now outdates speciﬁc routines that produce results like that shown in

Fig. 2.

Note: If you want to run the code with

MPI

, you will probably need to adapt

‘getconf.csh’, which deﬁnes the commands and ﬂags used to run MPI jobs (and which

is sourced by the scripts ‘start.csh’ and ‘run.csh’). Try

csh -v getconf.csh

csh -x getconf.csh

to see how ‘getconf.csh’ makes its decisions. You would add a section for the host name

of your machine with the particular settings. Since ‘getconf.csh’ is linked from the cen-

tral directory ‘pencil-code/bin’, your changes will be useful for all your other runs too.

3.2 Further tests

There is a number of other tests in the ‘samples/’ directory. You can use the script

‘bin/auto-test’ to automatically run these tests and have the output compared to refer-

10 THE PENCIL CODE

ence results.

4. Code structure 11

4 Code structure

4.1 Directory tree

pencil-code

src idl | dx doc bin samples

Makefile conv-slab interlocked-fluxrings ...

srcstart, run

timestep, deriv

equ

hydro | nohydro

density | nodensity

entropy | noentropy

grav_z | grav_r | nograv

magnetic | nomagnetic

Makefile.local

cparam.local

Figure 3: The basic structure of the code

The overall directory structure of the code is shown in Fig. 3. Under ‘pencil-code’, there

are currently the following ﬁles and directories:

bin/ config/ doc/ idl/ license/ perl/ samples/ sourceme.sh utils/

bugs/ dx/ lib/ misc/ README sourceme.csh src/ www/

Almost all of the source code is contained in the directory ‘src/’, but in order to encap-

sulate individual applications, the code is compiled separately for each run in a local

directory ‘src’ below the individual run directory, like e. g. ‘samples/conv-slab/src’.

It may be a good idea to keep your own runs also under

svn

cvs

(which is older than

but similar to

svn

), but this would normally be a different repository. On the machine

where you are running the code, you may want to check them out into a subdirectory of

12 THE PENCIL CODE

‘pencil-code/’. For example, we have our own runs in a repository called ‘pencil-runs’,

so we do

unix> cd $PENCIL_HOME

unix> svn co runs pencil-runs

In this case, ‘runs’ contains individual run directories, grouped in classes (like ‘spher’ for

spherical calculations, or ‘kinematic’ for kinematic dynamo simulations). The current

list of classes in our own ‘pencil-runs’ repository is

1d-tests/ disc/ kinematic/ rings/

2d-tests/ discont/ Misc/ slab_conv/

3d-tests/ discussion/ OLD/ test/

buoy_tube/ forced/ pass_only/

convstar/ interstellar/ radiation/

The directory ‘forced/’ contains some forced turbulence runs (both magnetic and non-

magnetic); ‘gravz/’ contains runs with vertical gravity; ‘rings/’ contains decaying MHD

problems (interlocked ﬂux rings as initial condition, for example); and ‘kinematic/’ con-

tains kinematic dynamo problems where the hydrodynamics is turned off entirely. The

ﬁle ‘samples/README’ should contain an up-to-date list and short description of the indi-

vidual classes.1

The subdirectory ‘src’ of each run directory contains a few local conﬁguration ﬁles (cur-

rently these are ‘Makefile.local’ and ‘cparam.local’) and possibly ‘ctimeavg.local’. To

compile the samples, links the ﬁles ‘.f90’, ‘.c’ and ‘Makefile.src’ need to be linked from

the top ﬁle[src/]src directory to the local directory ‘./src’. These links are set up by the

script pc_setupsrc) when used in the root of a run directory.

General-purpose visualization routines for

IDL

are in the directories ‘idl’ and ‘dx’,

respectively. There are additional and more specialized

IDL

directories in the different

branches under ‘pencil-runs’.

The directory ‘doc’ contains this manual; ‘bin’ contains a number of utility scripts

(mostly written in

csh

and

Perl

), and in particular the ‘start.csh’, ‘run.csh’, and

‘getconf.csh’ scripts. The ‘.svn’ directory is used (you guessed it) by

.svn

, and is not

normally directly accessed by the user; ‘bugs’, ﬁnally is used by us for internal purposes.

The ﬁles ‘sourceme.csh’ and ‘sourceme.sh’ will set up some environment variables — in

particular PATH — and aliases/shell functions for your convenience. If you do not want

to source one of these ﬁles, you need to make sure your

IDL

path is set appropriately

(provided you want to use

IDL

) and you will need to address the scripts from ‘bin’ with

their explicit path name, or adjust your PATH manually.

4.2 Basic concepts

4.2.1 Data access in pencils

Unlike the CRAY computers that dominated supercomputing in the 80s and early 90s,

all modern computers have a cache that constitutes a signiﬁcant bottleneck for many

1Our ‘pencil-runs’ directory also contains runs that were done some time ago. Occasionally, we try to

update these, especially if we have changed names or other input conventions.

4.2 Basic concepts 13

codes. This is the case if large three-dimensional arrays are constantly used within each

time step, which has the obvious advantage of working on long arrays and allows vector-

ization of elementary machine operations. This approach also implies conceptual sim-

plicity of the code and allows extensive use of the intuitive F90 array syntax. However,

a more cache-efﬁcient way of coding is to calculate an entire time step (or substep of a

multi-stage time-stepping scheme) only along a one-dimensional pencil of data within

the numerical grid. This technique is more efﬁcient for modern RISC processors: on

Linux PCs and SGI workstations, for example, we have found a speed-up by about 60%

in some cases. An additional advantage is a drastic reduction in temporary storage for

auxiliary variables within each time step.

4.2.2 Modularity

Each run directory has a ﬁle ‘src/Makefile.local’ in which you choose certain

modules

which tell the code whether or not entropy, magnetic ﬁelds, hydrodynamics, forcing, etc.

should be invoked. For example, the settings for forced turbulent MHD simulations are

HYDRO = hydro

DENSITY = density

ENTROPY = noentropy

MAGNETIC = magnetic

GRAVITY = nogravity

FORCING = forcing

MPICOMM = nompicomm

GLOBAL = noglobal

IO = io_dist

FOURIER = nofourier

This ﬁle will be processed by

make

and the settings are thus assignments of

make

variables. Apart from the physics modules (equation of motion: yes, density [pressure]:

yes, entropy equation: no, magnetic ﬁelds: yes, gravity: no, forcing: yes), a few technical

modules can also be used or deactivated; in the example above, these are

MPI

(switched

off), additional global variables (none), input/output (distributed), and

FFT

(not used).

The table in Sect. J in the Appendix lists all currently available modules.

Note that most of these

make

variables must be set, but they will normally obtain rea-

sonable default values in ‘Makefile’ (so you only need to set the non-standard ones in

‘Makefile.local’). It is by using this switching mechanism through make that we achieve

high ﬂexibility without resorting to excessive amounts of cryptic preprocessor directives

or other switches within the code.

Many possible combinations of modules have already been tested and examples are part

of the distribution, but you may be interested in a combination which was never tried

before and which may not work yet, since the modules are not fully orthogonal. In such

cases, we depend on user feedback for ﬁxing problems and documenting the changes for

others.

2We stress once more that we are not talking about F90 modules here, although there is some connec-

tion, as most of our modules deﬁne F90 modules: For example each of the modules

gravity simple

grav r

and

nogravity

deﬁnes a Fortran module

Gravity

14 THE PENCIL CODE

4.3 Files in the run directories

4.3.1 ‘start.in’, ‘run.in’, ‘print.in’

These ﬁles specify the startup and runtime parameters (see Sects. 5.12 and K.2), and the

list of diagnostic variables to print (see 5.5). They specify the setup of a given simulation

and are kept under

svn

in the individual ‘samples’ directories.

You may want to check for the correctness of these conﬁguration ﬁles by issuing the

command pc_configtest.

4.3.2 ‘datadir.in’

If this ﬁle exists, it must contain the name of an existing directory, which will be used as

data directory

, i. e. the directory where all results are written. If ‘datadir.in’ does not

exist, the data directory is ‘data/’.

4.3.3 ‘reference.out’

If present, ‘reference.out’ contains the output you should obtain in the given run direc-

tory, provided you have not changed any parameters. To see whether the results of your

run are OK, compare ‘time_series.dat’ to ‘reference.out’:

unix> diff data/time_series.dat reference.out

4.3.4 ‘start.csh’, ‘run.csh’, ‘getconf.csh’ [obsolete; see Sect. 5.1]

These are links to ‘$PENCIL_HOME/bin’. You will be constantly using the scripts

‘start.csh’ and ‘run.csh’ to initialize the code. Things that are needed by both (like the

name of the mpirun executable,

MPI

options, or the number of processors) are located in

‘getconf.csh’, which is never directly invoked.

4.3.5 ‘src/ ’

The ‘src’ directory contains two local ﬁles, ‘src/Makefile.local’ and ‘src/cparam.local’,

which allow the user to choose individual modules (see 4.2.2) and to set parameters like

the grid size and the number of processors for each direction. These two ﬁles are part

of the setup of a given simulation and are kept under

svn

in the individual ‘samples’

directories.

The ﬁle ‘src/cparam.inc’ is automatically generated by the script ‘mkcparam’ and contains

the number of fundamental variables for a given setup.

All other ﬁles in ‘src/’ are either links to source ﬁles (and ‘Makefile.src’) in the

‘$PENCIL_HOME/src’ directory, or object and module ﬁles generated by the compiler.

4.3 Files in the run directories 15

4.3.6 ‘data/ ’

This directory (the name of which will actually be overwritten by the ﬁrst line of

‘datadir.in’, if that ﬁle is present; see §4.3.2) contains the output from the code:

‘data/dim.dat’The global array dimensions.

‘data/legend.dat’The header line specifying the names of the diagnostic variables in

‘time_series.dat’.

‘data/time_series.dat’Time series of diagnostic variables (also printed to stdout). You

can use this ﬁle directly for plotting with

Gnuplot

IDL

Xmgrace

or similar tools (see

also §5.19).

‘data/tsnap.dat’, ‘data/tvid.dat’Time when the next snapshot ‘VARN’ or animation

slice should be taken.

‘data/params.log’Keeps a log of all your parameters: ‘start.x’ writes the startup pa-

rameters to this ﬁle, ‘run.x’ appends the runtime parameters and appends them anew,

each time you use the ‘RELOAD’ mechanism (see §5.10).

‘data/param.nml’Complete set of startup parameters, printed as Fortran namelist.

This ﬁle is read in by ‘run.x’ (this is how values of startup parameters are propagated

to ‘run.x’) and by

IDL

(if you use it).

‘data/param2.nml’Complete set of runtime parameters, printed as Fortran namelist.

This ﬁle is read by

IDL

(if you use it).

‘data/index.pro’Can be used as include ﬁle in

IDL

and contains the column in which

certain variables appear in the diagnostics ﬁle (‘time_series.dat’). It also contains the

positions of variables in the ‘VARN’ ﬁles. These positions depend on whether

entropy

noentropy

, etc, are invoked. This is a temporary solution and the ﬁle may disappear in

future releases.

‘data/interstellar.dat’Unformatted ﬁle containing the time at which the next su-

pernova event will occur, under certain supernova schemes. (Only needed by the

inter-

stellar

module.)

‘data/proc0’, ‘data/proc1’, . . . These are the directories containing data from the in-

dividual processors. So after running an

MPI

job on two processors, you will have the

two directories ‘data/proc0’ and ‘data/proc1’. Each of the directories can contain the

following ﬁles:

‘var.dat’binary ﬁle containing the latest snapshot;

16 THE PENCIL CODE

‘VARN’binary ﬁle containing individual snapshot number N;

‘dim.dat’ASCII ﬁle containing the array dimensions as seen by the given processor;

‘time.dat’ASCII ﬁle containing the time corresponding to ‘var.dat’ (not actually used

by the code, unless you use the

io mpiodist.f90

module);

‘grid.dat’binary ﬁle containing the part of the grid seen by the given processor;

‘seed.dat’the random seed for the next time step (saved for reasons of reproducibility).

For multi-processor runs with velocity forcing, the ﬁles ‘procN/seed.dat’ must all

contain the same numbers, because globally coherent waves of given wavenumber

are used;

‘X.xy’, ‘X.xz’, ‘X.yz’two-dimensional sections of variable X, where Xstands for the

corresponding variable. The current list includes

bx.xy bx.xz by.xy by.xz bz.xy bz.xz divu.xy lnrho.xz

ss.xz ux.xy ux.xz uz.xy uz.xz

Each processor writes its own slice, so these need to be reassembled if one wants

to plot a full slice.

5. Using the code 17

5 Using the code

5.1 Conﬁguring the code to compile and run on your computer

Note: We recommend to use the procedure described here, rather than the old method

described in Sects. 5.2 and 4.3.4.

Quick instructions: You may compile with a default compiler-speciﬁc conﬁguration:

1. Single-processor using the GNU compiler collection:

unix> pc_build -f GNU-GCC

2. Multi-processor using GNU with MPI support:

unix> pc_build -f GNU-GCC_MPI

Many compilers are supported already, please refer to the available conﬁg ﬁles in

‘$PENCIL_HOME/config/compilers/*.conf’, e.g. ‘Intel.conf’ and ‘Intel_MPI.conf’.

If you have to set up some compiler options speciﬁc to a certain host system you work on,

or if you like to create a host-speciﬁc conﬁguration ﬁle so that you can simply execute

pc_build without any options, you can clone an existing host-ﬁle, just include an existing

compiler conﬁguration, and simply only add the options you need. A good example of

a host-ﬁle is ‘$PENCIL_HOME/config/hosts/IWF/host-andromeda-GNU_Linux-Linux.conf’.

You may save a clone under ‘$PENCIL_HOME/config/hosts/<ID>.conf’, where ‘<ID>’ is to

be replaced by the output of pc_build -i. This will be the new default for pc_build.

If you don’t know what this was all about, read on.

In essence, conﬁguration, compiling and running the code work like this:

1. Create a conﬁguration ﬁle for your computer’s host ID.

2. Compile the code using pc_build.

3. Run the code using pc_run

In the following, we will discuss the essentials of this scheme. Exhaustive documen-

tation is available with the commands perldoc Pencil::ConfigFinder and perldoc

PENCIL::ConfigParser.

5.1.1 Locating the conﬁguration ﬁle

Commands like pc_build and pc_run use the Perl module ‘Pencil::ConfigFinder’ to lo-

cate an appropriate conﬁguration ﬁle and ‘Pencil::ConfigParser’ to read and interpret

it. When you use ‘ConfigFinder’ on a given computer, it constructs a host ID for the sys-

tem it is running on, and looks for a ﬁle ‘host_ID.conf’ in any subdirectory of ‘$PENCIL_-

HOME/config/hosts’.

E.g., if the host ID is “workhorse.pencil.org”, ‘ConfigFinder’ would consider the ﬁle

‘$PENCIL_HOME/config/hosts/pencil.org/workhorse.pencil.org.conf’.

18 THE PENCIL CODE

Note 1: The location in the tree under ‘hosts/’ is irrelevant, which allows you to group

related hosts by institution, owner, hardware, etc.

Note 2: ‘ConfigFinder’ actually uses the following search path:

1. ‘./config’

2. ‘$PENCIL_HOME/config-local’

3. ‘$HOME/.pencil/config-local’

4. ‘$PENCIL_HOME/config’

This allows you to override part of the ‘config/’ tree globally on the given ﬁle system,

or locally for a particular run directory, or for one given copy of the PENCIL CODE. This

search path is used both, for locating the conﬁguration ﬁle for your host ID, and for

locating included ﬁles (see below).

The host ID is constructed based on information that is easily available for your system.

The algorithm is as follows:

1. Most commands using ‘ConfigFinder’ have a ‘--host-id’ (sometimes abbreviated

as ‘-H’) option to explicitly set the host ID.

2. If the environment variable

PENCIL HOST ID

is set, its value is used.

3. If any of the ﬁles ‘./host-ID’, ‘$PENCIL_HOME/host-ID’, ‘$HOME/.pencil/host-ID’, ex-

ists, its ﬁrst line is used.

4. If ‘ConfigFinder’ can get hold of a

fully qualiﬁed host name

, that is used as host

ID.

5. Else, a combination of host name, operating system name and possibly some other

information characterizing the system is used.

6. If no conﬁg ﬁle for the host ID is found, the current operating system name is tried

as fallback host ID.

To see which host IDs are tried (up to the ﬁrst one for which a conﬁguration ﬁle is found),

run

unix> pc_build --debug-config

This command will tell you the host-ID of the system that you are using:

unix> pc_build -i

5.1.2 Structure of conﬁguration ﬁles

It is strongly recommended to include in a users conﬁguration ﬁle one of the preset

compiler suite conﬁguration ﬁles. Then, only minor options need to be set by a user, e.g.

the optimization ﬂags. One of those user conﬁguration ﬁles looks rather simple:

# Simple config file. Most files don’t need more.

%include compilers/GNU-GCC

or if you prefer a different compiler:

5.1 Conﬁguring the code to compile and run on your computer 19

# Simple Intel compiler suite config file.

%include compilers/Intel

A more complex ﬁle (using MPI with additional options) would look like this:

# More complex config file.

%include compilers/GNU-GCC_MPI

%section Makefile

MAKE_VAR1 = -j4 # joined compilation with four threads

FFLAGS += -O3 -Wall -fbacktrace # don’t redefine, but append with ’+=’

%endsection Makefile

%section runtime

mpiexec = mpirun # some MPI backends need a redefinition of mpiexec

%endsection runtime

%section environment

SCRATCH_DIR=/var/tmp/$USER

%endsection environment

Adding ” MPI” to a compiler suite’s name is usually sufﬁcient to use MPI.

Note 3: We strongly advise not to mix Fortran- and C-compilers from different manu-

facturers or versions by manually including multiple separate compiler conﬁgurations.

Note 4: We strongly advise to use at maximum the optimization levels ’-O2’ for the

Intel compiler and ’-O3’ for all other compilers. Higher optimization levels implicate an

inadequate loss of precision.

The ‘.conf’ ﬁles consist of the following elements:

Comments: A#sign and any text following it on the same line are ignored.

Sections: There are three sections:

Makeﬁle for setting make parameters

runtime for adding compiler ﬂags used by pc_run

environment shell environment variables for compilation and running

Include statements: An %include ... statement is recursively replaced by the con-

tents of the ﬁles it points to.3

The included path gets a .conf sufﬁx appended. Included paths are typically “ab-

solute”, e.g.

%include os/Unix

will include the ﬁle ‘os/Unix.conf’ in the search path listed above (typically from

‘$PENCIL_HOME/config’). However, if the included path starts with a dot, it is a rel-

ative path, so

3However, if the include statement is inside a section, only the ﬁle’s contents inside that section are

inserted.

20 THE PENCIL CODE

%include ./Unix

will only search in the directory where the including ﬁle is located.

Assignments: Statements like FFLAGS += -O3 or mpiexec=mpirun are assignments and

will set parameters that are used by pc_build/make or pc_run.

Lines ending with a backslash ‘\’ are continuation lines.

If possible, one should always use incremental assignments, indicated by using a

+= sign instead of =, instead of redeﬁning certain ﬂags.

Thus,

CFLAGS += -O3

CFLAGS += -I../include -Wall

is the same as

CFLAGS = $(CFLAGS) -O3 -I../include -Wall

5.1.3 Compiling the code

Use the pc_build command to compile the code:

unix> pc_build # use default compiler suite

unix> pc_build -f Intel_MPI # specify a compiler suite

unix> pc_build -f os/GNU_Linux,mpi/open-mpi # explicitly specify config files

unix> pc_build VAR=something # set variables for the makefile

unix> pc_build --cleanall # remove generated files

The third example circumvents the whole host ID mechanism by explicitly instructing

pc_build which conﬁguration ﬁles to use. The fourth example shows how to deﬁne extra

variables (VAR=something) for the execution of the Makefile.

See pc_build --help for a complete list of options.

5.1.4 Running the code

Use the pc_run command to run the code:

unix> pc_run # start if necessary, then run

unix> pc_run start

unix> pc_run run

unix> pc_run start run^3 # start, then run 3 times

unix> pc_run start run run run # start, then run 3 times

unix> pc_run ^3 # start if necessary, then run 3 times

See pc_run --help for a complete list of options.

5.1.5 Testing the code

The pc_auto-test command uses pc_build for compiling and pc_run for running the

tests. Here are a few useful options:

5.2 Adapting ‘Makefile.src’ﬁle]Makeﬁle.src@Makefile.src [obsolete; see Sect. 5.1]21

unix> pc_auto-test

unix> pc_auto-test --no-pencil-check # suppress pencil consistency check

unix> pc_auto-test --max-level=1 # run only tests in level 0 and 1

unix> pc_auto-test --time-limit=2m # kill each test after 2 minutes

See pc_auto-test --help for a complete list of options.

The pencil-test script will use pc_auto-test if given the ‘--use-pc_auto-test’ or ‘-b’

option:

unix> pencil-test --use-pc_auto-test

unix> pencil-test -b # ditto

unix> pencil-test -b -Wa,--max-level=1,--no-pencil-check # quick pencil

See pencil-test --help for a complete list of options.

5.2 Adapting ‘Makefile.src’ [obsolete; see Sect. 5.1]

By default, one should use the above described conﬁguration mechanism for compila-

tion. If for whatever reason one needs to work with a modiﬁed ‘Makefile’, there is a

mechanism for picking the right set of compiler ﬂags based on the hostname.

To give you an idea of how to add your own machines, let us assume you have several

Linux boxes running a compiler f95 that needs the options ‘-O2 -u’, while one of them,

Janus

, runs a compiler f90 which needs the ﬂags ‘-O3’ and requires the additional op-

tions ‘-lmpi -llam’ for

MPI

The ‘Makefile.src’ you need will have the following section:

### Begin machine dependent

## IRIX64:

[. . . ] (leave as it is in the original Makeﬁle)

## OSF1:

[. . . ] (leave as it is in the original Makeﬁle)

## Linux:

[. . . ] (leave everything from original Makeﬁle and add:)

#FC=f95

#FFLAGS=-O2 -u

#FC=f90 #(Janus)

#FFLAGS=-O3 #(Janus)

#LDMPI=-lmpi -llam #(Janus)

## SunOS:

[. . . ] (leave as it is in the original Makeﬁle)

## UNICOS/mk:

[. . . ] (leave as it is in the original Makeﬁle)

## HI-UX/MPP:

[. . . ] (leave as it is in the original Makeﬁle)

## AIX:

[. . . ] (leave as it is in the original Makeﬁle)

22 THE PENCIL CODE

Table 1: Compiler ﬂags for common compilers. Note that some combinations of OS and compiler require

much more elaborate settings; also, if you use MPI, you will have to set LDMPI.

Compiler FC FFLAGS CC CFLAGS

Unix/POSIX:

GNU gfortran -O3 gcc -O3 -DFUNDERSC=1

Intel ifort -O2 icc -O3 -DFUNDERSC=1

PGI pgf95 -O3 pgcc -O3 -DFUNDERSC=1

G95 g95 -O3 -fno-second-underscore gcc -O3 -DFUNDERSC=1

Absoft f95 -O3 -N15 gcc -O3 -DFUNDERSC=1

IBM XL xlf95 -qsuffix=f=f90 -O3 xlc -O3 -DFUNDERSC=1

outdated:

IRIX Mips f90 -64 -O3 -mips4 cc -O3 -64 -DFUNDERSC=1

Compaq f90 -fast -O3 cc -O3 -DFUNDERSC=1

### End machine dependent

Note 1 There is a script for adapting the Makeﬁle: ‘adapt-mkfile’. In the example

above, #(Janus) is not a comment, but marks this line to be activated (uncommented) by

adapt-mkfile if your hostname (‘uname -n‘) matches ‘Janus’ or ‘janus’ (capitalization

is irrelevant). You can combine machine names with a vertical bar: a line containing

#(onsager|Janus) will be activated on both,

Janus

and

Onsager

Note 2 If you want to experiment with compiler ﬂags, or if you want to get things

running without setting up the machine-dependent section of the ‘Makefile’, you can set

make

variables at the command line in the usual manner:

src> make FC=f90 FFLAGS=’-fast -u’

will use the compiler f90 and the ﬂags ‘-fast -u’ for both compilation and linking. Ta-

ble 1 summarizes ﬂags we use for common compilers.

5.3 Changing the resolution

It is advisable to produce a new run directory each time you run a new case. (This does

not include restarts from an old run, of course.) If you have a 323run in some directory

‘runA_32a’, then go to its parent directory, i.e.

runA_32a> cd ..

forced> pc_newrun runA_32a runA_64a

forced> cd runA_64a/src

forced> vi cparam.local

and edit the ‘cparam.local’ for the new resolution.

If you have ever wondered why we don’t do dynamic allocation of the main variable (f)

array, the main reason it that with static allocation the compiler can check whether we

are out of bounds.

5.4 Using a non-equidistant grid 23

5.4 Using a non-equidistant grid

We introduce a non-equidistant grid zi(by now, this is also implemented for the other

directions) as a function z(ζ)of an equidistant grid ζiwith grid spacing ∆ζ= 1.

The way the parameters are handled, the box size and position are not changed when

you switch to a non-equidistant grid, i.e. they are still determined by

xyz0

and

Lxyz

The ﬁrst and second derivatives can be calculated by

dz =df

dζ

dz =1

z′f′(ζ),d2f

dz2=1

z′2f′′(ζ)−z′′

z′3f′(ζ),(1)

which can be written somewhat more compactly using the inverse function ζ(z):

dz =ζ′(z)f′(ζ),d2f

dz2=ζ′2(z)f′′(ζ) + ζ′′(z)ζ′(z)f′(ζ).(2)

Internally, the code uses the quantities

dz 1

≡1/z′=ζ′(z)and

dz tilde

≡ −z′′/z′2=

ζ′′/ζ′, and stores them in ‘data/procN/grid.dat’.

The parameters

lequidist

(a 3-element logical array),

grid func

coeff grid

≥2-element real array) are used to choose a non-equidistant grid and deﬁne

the function z(ζ). So far, one can choose between grid_function=’sinh’,grid_-

function=’linear’ (which produces an equidistant grid for testing purposes), and

grid_function=’step-linear’.

The sinh proﬁle: For grid_function=’sinh’, the function z(ζ)is given by

z(ζ) = z0+Lz

sinh a(ζ−ζ∗) + sinh a(ζ∗−ζ1)

sinh a(ζ2−ζ∗) + sinh a(ζ∗−ζ1),(3)

where z0and z0+Lzare the lowest and uppermost levels, ζ1,ζ2are the ζvalues represent-

ing those levels (normally ζ1= 0, ζ2=Nz−1for a grid of Nzvertical layers [excluding

ghost layers]), and ζ∗is the ζvalue of the inﬂection point of the sinh function. The z

coordinate and ζvalue of the inﬂection point are related via

z∗=z0+Lz

sinh a(ζ∗−ζ1)

sinh a(ζ2−ζ∗) + sinh a(ζ∗−ζ1),(4)

which can be inverted (“after some algebra”) to

ζ∗=ζ1+ζ2

2+1

aartanh 2z∗−z0

Lz−1tanh aζ2−ζ1

2.(5)

General proﬁle: For a general monotonic function ψ() instead of sinh we get,

z(ζ) = z0+Lz

ψ[a(ζ−ζ∗)] + ψ[a(ζ∗−ζ1)]

ψ[a(ζ2−ζ∗)] + ψ[a(ζ∗−ζ1)] ,(6)

and the reference point ζ∗is found by inverting

z∗=z0+Lz

ψ[a(ζ∗−ζ1)]

ψ[a(ζ2−ζ∗)] + ψ[a(ζ∗−ζ1)] ,(7)

numerically.

24 THE PENCIL CODE

Duct ﬂow: The proﬁle function grid_function=’duct’ generates a grid proﬁle for tur-

bulent ﬂow in ducts. For a duct ﬂow, most gradients are steepest close to the walls, and

hence very ﬁne resolution is required near the walls. Here we follow the method of [22]

and use a Chebyshev-type grid with a cosine distribution of the grid points such that in

the ydirection they are located at

yj=hcos θj,(8)

where

θj=(Ny−j)π

Ny−1, j = 1,2,...,Ny(9)

and h=Ly/2is the duct half-width.

Currently this method is only adapted for ducts where xis the stream-wise direction, z

is in the span-wise direction and the walls are at y=y0and y=y0+Ly.

In order to have ﬁne enough resolution, the ﬁrst grid point should be a distance δ=

0.05 lwfrom the wall, where

lw=ν

uτ

, uτ=rτw

ρ,(10)

and τwis the shear wall stress. This is accomplished by using at least

Ny≥N∗

y=π

arccos1−δ

h+ 1 (11)

=πrh

2δ+ 1 −π

24r2δ

h+O"δ

h3/2#(12)

grid points in the y-direction. After rounding up to the next integer value, we ﬁnd that

the truncated condition

Ny≥&πrh

2δ'+ 1 (13)

(where ⌈x⌉denotes the ceiling function, i.e. the smallest integer equal to, or larger than,

x) gives practically identical results.

Example: To apply the sinh proﬁle, you can set the following in ‘start.in’ (this exam-

ple is from ‘samples/sound-spherical-noequi’):

&init_pars

[...]

xyz0 = -2., -2., -2. ! first corner of box

Lxyz = 4., 4., 4. ! box size

lperi = F , F , F ! periodic direction?

lequidist = F, F, T ! non-equidistant grid in z

xyz_star = , , -2. ! position of inflection point

grid_func = , , ’sinh’ ! sinh function: linear for small, but

! exponential for large z

coeff_grid = , , 0.5

5.5 Diagnostic output 25

The parameter array

coeff grid

represents z∗and a; the bottom height z0and the total

box height Lzare taken from

xyz0

and

Lxyz

as in the equidistant case. The inﬂection

point of the sinh proﬁle (the part where it is linear) is not in the middle of the box,

because we have set

xyz star(3)

(i. e. z∗) to −2.

5.5 Diagnostic output

Every

it1

time steps (

it1

is a runtime parameter, see Sect. K.2), the code writes moni-

toring output to

stdout

and, parallel to this, to the ﬁle ‘data/time_series.dat’. The vari-

ables that appear in this listing and the output format are deﬁned in the ﬁle ‘print.in’

and can be changed without touching the code (even while the code is running). A simple

example of ‘print.in’ may look like this:

t(F10.3)

urms(E13.4)

rhom(F10.5)

oum

which means that the output table will contain time

in the ﬁrst column formatted as

(F10.3), followed by the mean squared velocity,

urms

(i.e. hu2i1/2) in the second column

with format (E13.4), the average density

rhom

(hρi, which allows to monitor mass con-

servation) formatted (F10.5) and the kinetic helicity

oum

(that is hω·ui) in the last

column with the default format (E10.2).4The corresponding diagnostic output will look

----t---------urms--------rhom------oum----

0.000 4.9643E-03 14.42457 -8.62E-06

0.032 3.9423E-03 14.42446 -5.25E-06

0.063 6.8399E-03 14.42449 -3.50E-06

0.095 1.1437E-02 14.42455 -2.58E-06

0.126 1.6980E-02 14.42457 -1.93E-06

5.6 Data ﬁles

5.6.1 Snapshot ﬁles

Snapshot ﬁles contain the values of all evolving variables and are sufﬁcient to restart a

run. In the case of an MHD simulation with entropy equation, for example, the snapshot

ﬁles will contain the values of velocity, logarithmic density, entropy and the magnetic

vector potential.

There are two kinds of snapshot ﬁles: the current snapshot and permanent snapshots,

both of which reside in the directory ‘data/procN/’. The parameter

isav

determines the

frequency at which the current snapshot ‘data/procN/var.dat’ is written. If you keep

this frequency too high, the code will spend a lot of time on I/O, in particular for large

jobs; too low a frequency makes it difﬁcult to follow the evolution interactively during

test runs. There is also the

ialive

parameter. Setting this to 1 or 10 gives an updated

4The format speciﬁers are like in Fortran, apart from the fact that the Eformat will use standard

scientiﬁc format, corresponding to the Fortran 1pE syntax. Seasoned Fortran IV programmers may use

formats like (0pE13.4) to enjoy nostalgic feelings, or (1pF10.5) if they depend on getting wrong numbers.

26 THE PENCIL CODE

timestep in the ﬁles ‘data/proc*/alive.info’. You can put

ialive=0

to turn this off to

limit the I/O on your machine.

The permanent snapshots ‘data/proc*/VARN’ are written every

dsnap

time units. These

ﬁles are numbered consecutively from N= 1 upward and for long runs they can occupy

quite some disk space. On the other hand, if after a run you realize that some addi-

tional quantity qwould have been important to print out, these ﬁles are the only way to

reconstruct the time evolution of qwithout re-running the code.

File structure Snapshot ﬁles consist of the following

Fortran records

1. variable vector f[

nvar

]

2. time t[1], coordinate vectors x[

], y[

], z[

], grid spacings δx [1], δy [1], δz

[1], shearing-box shift ∆y[1]

All numbers (apart from the record markers) are single precision (4-byte) ﬂoating point

numbers, unless you use double precision (see §5.21, in which case all numbers are 8-

byte ﬂoating point numbers, while the record markers remain 4-byte integers.

The script pc_tsnap allows you to determine the time tof a snapshot ﬁle:

unix> pc_tsnap data/proc0/var.dat

data/proc0/var.dat: t = 8.32456

unix> pc_tsnap data/proc0/VAR2

data/proc0/VAR2: t = 2.00603

5.7 Video ﬁles and slices

We use the terms

video ﬁles

and

slice ﬁles

interchangeably. These ﬁles contain a time

series of values of one variable in a given plane. The output frequency of these video

snapshots is set by the parameter

dvid

(in code time units).

When output to video ﬁles is activated by some settings in ‘run.in’ (see example below)

and the existence of ‘video.in’, slices are written for four planes:

1. x-zplane (yindex

; ﬁle sufﬁx .xz)

2. y-zplane (yindex

; sufﬁx .yz)

3. x-yplane (yindex

; sufﬁx .xy)

4. another slice parallel to the x-yplane (yindex

iz2

; sufﬁx .xy2)

You can specify the position of the slice planes by setting the parameters

and

iz2

in the namelist

run pars

in ‘run.in’. Alternatively, you can set the input parameter

slice position

to one of ’p’ (periphery of box) or ’m’ (middle of box). Or you can also

specify the z−position in terms of z using the tags

zbot slice

and

ztop slice

. In this case,

the

zbot slice

slice will have the sufﬁx .xy and the

ztop slice

the sufﬁx .xy2

In the ﬁle ‘video.in’ of your run directory, you can choose for which variables you want

to get video ﬁles; valid choices are listed in §K.4.

Fortran record

is marked by the 4-byte integer byte count of the data in the record at the beginning

and the end, i.e. has the form hNbytes,raw data, Nbytesi

5.7 Video ﬁles and slices 27

The

slice ﬁles

are written in each processor directory ‘data/proc*/’ and have a ﬁle name

indicating the individual variable (e. g. ‘slice_uu1.yz’ for a slice of uxin the y-zplane).

Before visualizing slices one normally wants to combine the sub-slices written by each

processor into one global slice (for each plane and variable). This is done by running

‘src/read_videofiles.x’, which will prompt for the variable name, read the individual

sub-slices and write global slices to ‘data/’ Once all global slices have been assembled

you may want to remove the local slices ‘data/proc*/slice*’.

To read all sub-slices demanded in ‘video.in’ at once use ‘src/read_all_videofiles.x’.

This program doesn’t expect any user input and can thus be submitted in computing

queues.

For visualization of slices, you can use the

IDL

routines ‘rvid_box.pro’, ‘rvid_-

plane.pro’, or ‘rvid_line.pro’ which allows the ﬂag ‘/png’ for writing

PNG

images that

can then be combined into an

MPEG

movie using

mpeg encode

. Based on ‘rvid_box’,

you can write your own video routines in

IDL

An example Suppose you have set up a run using

entropy.f90

and

magnetic.f90

(most

probably together with

hydro.f90

and other modules). In order to animate slices of en-

tropy sand the z-component Bzof the magnetic ﬁeld, in planes passing through the

center of the box, do the following:

1. Write the following lines to ‘video.in’ in your run directory:

divu

2. Edit the namelist

run pars

in the ﬁle ‘run.in’. Request slices by setting

write slices

and set

dvid

and

slice position

to reasonable values, say

!lwrite_slices=T !(no longer works; write requested slices into video.in)

dvid=0.05

slice_position=’m’

3. Run the PENCIL CODE:

unix> start.csh

unix> run.csh

4. Run ‘src/read_videofiles.x’ to assemble global slice ﬁles from those scattered

across ‘data/proc*/’:

unix> src/read_videofiles.x

enter name of variable (lnrho, uu1, ..., bb3): ss

unix> src/read_videofiles.x

enter name of variable (lnrho, uu1, ..., bb3): bb3

5. Start

IDL

and run ‘rvid_box’:

unix> idl

IDL> rvid_box,’bb3’

IDL> rvid_box,’ss’,min=-0.3,max=2.

etc.

28 THE PENCIL CODE

Another example Suppose you have set up a run using

magnetic.f90

and some other

modules. This run studies some process in a “surface” layer inside the box. This “surface”

can represent a sharp change in density or turbulence. So you deﬁned your box setting

the z= 0 point at the surface. Therefore, your ‘start.in’ ﬁle will look something similar

to:

&init_pars

lperi=T,T,F

bcz = ’s’,’s’,’a’,’hs’,’s’,’s’,’a’

xyz0 = -3.14159, -3.14159, -3.14159

Lxyz = 6.28319, 6.28319, 9.42478

Now you can analyze quickly the surface of interest and some other xy−slice setting

zbot slice

and

ztop slice

in the ‘run.in’ ﬁle:

&run_pars

slice_position=’c’

zbot_slice=0.

ztop_slice=0.2

In this case, the slices with the sufﬁx .xy will be at the “surface” and the ones with the

sufﬁx .xy2 will be at the position z= 0.2above the surface. And you can visualize this

slices by:

1. Write the following lines to ‘video.in’ in your run directory:

2. Edit the namelist

run pars

in the ﬁle ‘run.in’ to include

zbot slice

and

ztop slice

3. Run the PENCIL CODE:

unix> start.csh

unix> run.csh

4. Run ‘src/read_videofiles.x’ to assemble global slice ﬁles from those scattered

across ‘data/proc*/’:

unix> src/read_videofiles.x

enter name of variable (lnrho, uu1, ..., bb3): bb3

5. Start

IDL

, load the slices with ‘pc_read_video’ and plot them at some time:

unix> idl

IDL> pc_read_video, field=’bb3’,ob=bb3,nt=ntv

IDL> tvscl,bb3.xy(*,*,100)

IDL> tvscl,bb3.xy2(*,*,100)

etc.

File structure

Slice ﬁles

consist of one Fortran record (see footnote on page 26) for

each slice, which contains the data of the variable (without ghost zones), the time tof

the snapshot and the position of the sliced variable (e. g. the xposition for a y-zslice):

1. data1[

], time t1[1], position1[1]

5.8 Averages 29

2. data2[

], time t2[1], position2[1]

3. data3[

], time t3[1], position3[1]

etc.

5.8 Averages

5.8.1 One-dimensional output averaged in two dimensions

In the ﬁle ‘xyaver.in’, z-dependent (horizontal) averages are listed. They are written to

the ﬁle ‘data/xyaverages.dat’. A new line of averages is written every

it1

th time steps.

There is the possibility to output two-dimensional averages. The result then de-

pends on the remaining dimension. The averages are listed in the ﬁles ‘xyaver.in’,

‘xzaver.in’, and ‘yzaver.in’ where the ﬁrst letters indicate the averaging directions.

The output is then stored to the ﬁles ‘data/xyaverages.dat’, ‘data/xzaverages.dat’, and

‘data/yzaverages.dat’. The output is written every

it1d

th time steps.

The rms values of the so deﬁned mean magnetic ﬁelds are referred to as

bmz

bmy

and

bmx

, respectively, and the rms values of the so deﬁned mean velocity ﬁelds are referred

to as

umz

umy

, and

umx

. (The last letter indicates the direction on which the averaged

quantity still depends.)

See Sect. 9.2 on how to add new averages.

In idl such xy-averages can be read using the procedure ‘pc_read_xyaver’.

5.8.2 Two-dimensional output averaged in one dimension

There is the possibility to output one-dimensional averages. The result then depends

on the remaining two dimensions. The averages are listed in the ﬁles ‘yaver.in’,

‘zaver.in’, and ‘phiaver.in’ where the ﬁrst letter indicates the averaging direction.

The output is then stored to the ﬁles ‘data/yaverages.dat’, ‘data/zaverages.dat’, and

‘data/phiaverages.dat’.

See Sect. 9.2 on how to add new averages.

Disadvantage: The output ﬁles, e.g. ‘data/zaverages.dat’, can be rather big because each

average is just appended to the ﬁle.

5.8.3 Azimuthal averages

Azimuthal averages are controlled by the ﬁle ‘phiaver.in’, which currently supports the

quantities listed in Sect. K.5. In addition, one needs to set

lwrite phiaverages

lwrite -

yaverages

, or

lwrite zaverages

.true.

. For example, if ‘phiaver.in’ contains the single

line

b2mphi

then you will get azimuthal averages of the squared magnetic ﬁeld B2.

30 THE PENCIL CODE

Azimuthal averages are written every

d2davg

time units to the ﬁles

‘data/averages/PHIAVGN’. The ﬁle format of azimuthal-average ﬁles consists of the

following

Fortran records

1. number of radial points Nr,φ−avg [1], number of vertical points Nz,φ−avg [1], number

of variables Nvar,φ−avg[1], number of processors in zdirection [1]

2. time t[1], positions of cylindrical radius rcyl [Nr,φ−avg] and z[Nz,φ−avg] for the grid,

radial spacing δrcyl [1], vertical spacing δz [1]

3. averaged data [Nr,φ−avg×Nz,φ−avg]

4. label length [1], labels of averaged variables [Nvar,φ−avg]

All numbers are 4-byte numbers (ﬂoating-point numbers or integers), unless you use

double precision (see §5.21).

To read and visualize azimuthal averages in

IDL

, use ‘$PENCIL_HOME/idl/files/pc_-

read_phiavg.pro’

IDL> avg = pc_read_phiavg(’data/averages/PHIAVG1’)

IDL> contour, avg.b2mphi, avg.rcyl, avg.z, TITLE=’!17B!U2!N!X’

or have a look at ‘$PENCIL_HOME/idl/phiavg.pro’ for a more sophisticated example.

5.8.4 Time averages

Time averages need to be prepared in the ﬁle ‘src/ctimeavg.local’, since they use extra

memory. They are controlled by the averaging time τavg (set by the parameter

tavg

‘run.in’), and by the indices

idx tavg

of variables to average.

Currently, averaging is implemented as exponential (memory-less) average,6

hfit+δt =hfit+δt

τavg

[f(t+δt)− hfit],(16)

which is equivalent to

hfit=

e−(t−t′)/τavg f(t′)dt′.(17)

Here t0is the time of the snapshot the calculation started with, i.e. the snapshot read by

the last run.x command. Note that the implementation (16) will approximate Eq. (17)

only to ﬁrst-order accuracy in δt. In practice, however, δt is small enough to make this

accuracy sufﬁce.

6At some point we may also implement the more straight-forward average

hfit+δt =hfit+δt

t−t0+δt[f(t+δt)− hfit],(14)

which is equivalent to

hfit=1

t−t0

f(t′)dt′,(15)

but we do not expect large differences.

5.9 Helper scripts 31

In ‘src/ctimeavg.local’, you need to set the number of slots used for time averages.

Each of these slots uses mx ×my ×mz ﬂoating-point numbers, i.e. half as much memory

as each fundamental variable.

For example, if you want to get time averages of all variables, set

integer, parameter :: mtavg=mvar

in ‘src/ctimeavg.local’, and don’t set

idx tavg

in ‘run.in’.

If you are only interested in averages of variables 1–3and 6–8(say, the velocity vector

and the magnetic vector potential in a run with ‘hydro.f90’, ‘density.f90’, ‘entropy.f90’

and ‘magnetic.f90’), then set

integer, parameter :: mtavg=6

in ‘src/ctimeavg.local’, and set

idx_tavg = 1,2,3,6,7,8 ! time-average velocity and vector potential

in ‘run.in’.

Permanent snapshots of time averages are written every

tavg

time units to

the ﬁles ‘data/proc*/TAVN’. The current time averages are saved periodically in

‘data/proc*/timeavg.dat’ whenever ‘data/proc*/var.dat’ is written. The ﬁle format for

time averages is equivalent to that of the snapshots; see §5.6.1 above.

5.9 Helper scripts

The ‘bin’ directory contains a collection of utility scripts, some of which are discussed

elsewhere, Here is a list of the more important ones.

adapt-mkﬁle Activate the settings in a ‘Makefile’ that apply to the given computer,

see §5.2.

auto-test Verify that the code compiles and runs in a set of run directories and compare

the results to the reference output. These tests are carried out routinely to ensure

that the

svn

version of the code is in a usable state.

cleanf95 Can be use to clean up the output from the Intel x86 Fortran 95 compiler

(ifc).

copy-proc-to-proc Used for restarting in a different directory. Example

copy-proc-to-proc seed.dat ../hydro256e.

copy-snapshots Copy snapshots from a processor-local directory to the global direc-

tory. To be started in the background before ‘run.x’ is invoked. Used by ‘start.csh’

and ‘run.csh’ on network connected processors.

pc copyvar var1 var2 source dest Copies snapshot ﬁles from one directory (source)

to another (dest). See documentation in ﬁle.

pc copyvar v v dir Copies all ‘var.dat’ ﬁles from current directory to ‘var.dat’ in ’dir’

run directory. Used for restarting in a different directory.

32 THE PENCIL CODE

pc copyvar N v Used to restart a run from a particular snapshot ‘VARN’. Copies a

speciﬁed snapshot ‘VARN’ to ‘var.dat’ where Nand (optionally) the target run di-

rectory are given on the command line.

cvs-add-rundir Add the current run directory to the

svn

repository.

cvsci run Similar to cvs-add-rundir, but it also checks in the ‘*.in’ and ‘src/*.local’

ﬁles. It also checks in the ﬁles ‘data/time_series.dat’, ‘data/dim.dat’ and

‘data/index.pro’ for subsequent processing in

IDL

on another machine. This is

particularly useful if collaborators want to check each others’ runs.

dx *These script perform several data collection or reformatting exercises required to

read particular ﬁles into

. They are called internally by some of the

macros

in the ‘dx/macros/’ directory.

getconf.csh See §4.3.4

gpgrowth Plot simple time evolution with Gnuplot’s ASCII graphics for fast orienta-

tion via a slow modem line.

local Materialize a symbolic link

mkcparam Based on ‘Makefile’ and ‘Makefile.local’, generate ‘src/cparam.inc’,

which speciﬁes the number

mvar

of fundamental variables, and

maux

of auxiliary

variables. Called by the ‘Makefile’.

pc mkdatadir Creates a link to a data directory in a suitable workspace. By default

this is on ‘/var/tmp/’, but different locations are speciﬁed for different machines.

mkdotin Generate minimal ‘start.in’, ‘run.in’ ﬁles based on ‘Makefile’ and

‘Makefile.local’.

mkinpars Wrapper around ‘mkdotin’ — needs proper documentation.

mkproc-tree Generates a multi-processor(‘procN’/), directory structure. Useful when

copying data ﬁles in a processor tree, such as slice ﬁles.

mkwww Generates a template HTML ﬁle for describing a run of the code, showing

input parameters and results.

move-slice Moves all the slice ﬁles from a processor tree structure, ‘procN/’, to a new

target tree creating directories where necessary.

nl2idl Transform a Fortran

namelist

(normally the ﬁles ‘param.nml’, ‘param2.nml’ writ-

ten by the code) into an

IDL

structure. Generates an

IDL

ﬁle that can be sourced

from ‘start.pro’ or ‘run.pro’.

pacx-adapt-makeﬁle Version of adapt-makeﬁle for highly distributed runs using

PACX MPI.

pc newrun Generates a new run directory from an old one. The new one contains a

copy of the old ‘*.local’ ﬁles, runs pc_setupsrc, and makes also a copy of the old

‘*.in’ and ‘k.dat’ ﬁles.

pc newscan Generates a new scan directory from an old one. The new one contains

a copy of the old, e.g. the one given under ‘samples/parameter_scan’. Look in the

‘README’ ﬁle for details.

5.10 RELOAD and STOP ﬁles 33

pc inspectrun Check the execution of the current run: prints legend and the last few

lines of the ‘time_series.dat’ ﬁle. It also appends this result to a ﬁle called ‘SPEED’,

which contains also the current wall clock time, so you can work out the speed of

the code (without being affected by i/o time).

read videoﬁles.csh The script for running read videoﬁles.x.

remote-top Create a ﬁle ‘top.log’ in the relevant ‘procN/’ directory containing the

output of top for the appropriate processor. Used in batch scripts for multi-

processor runs.

run.csh The script for producing restart ﬁles with the initial condition; see §4.3.4

scpdatadir Make a tarball of data directory, ‘data/’ and use scp to secure copy to copy

it to the speciﬁed destination.

pc setupsrc Link ‘start.csh’, ‘run.csh’ and ‘getconf.csh’ from ‘$PENCIL_HOME/bin’.

Generate ‘src/’ if necessary and link the source code ﬁles from ‘$PENCIL_HOME/src’

to that directory.

start.csh The script for initializing the code; see §4.3.4

summarize-history Evaluate ‘params.log’ and print a history of changes.

timestr Generate a unique time string that can be appended to ﬁle names from shell

scripts through the backtick mechanism.

pc tsnap Extract time information from a snapshot ﬁle, ‘VARN’.

There are several additional scripts on ‘pencil-code/utils’. Some are located in sepa-

rate folders according to users. There could be redundancies, but it is often just as easy

to write your own new script than ﬁguring out how something else works.

5.10 RELOAD and STOP ﬁles

The code periodically (every

time steps) checks for the existence of two ﬁles, ‘RELOAD’

and ‘STOP’, which can be used to trigger certain behavior.

Reloading run parameters In the directory where you started the code, create the ﬁle

‘RELOAD’ with

unix> touch RELOAD

34 THE PENCIL CODE

to force the code to re-read the runtime parameters from ‘run.in’. This will happen the

next time the code is writing monitoring output (the frequency of this happening is

controlled by the input parameter

, see Sect. 5.12).

Each time the parameters are reloaded, the new set of parameters is appended (in the

form of

namelists

) to the ﬁle ‘data/params.log’ together with the time t, so you have a

full record of your changes. If ‘RELOAD’ contains any text, its ﬁrst line will be written to

‘data/params.log’ as well, which allows you to annotate changes:

unix> echo "Reduced eta to get fields growing" > RELOAD

Use the command summarize-history to print a history of changes.

Stopping the code In the directory where you started the code, create the ﬁle ‘STOP’

with

unix> touch STOP

to stop the code in a controlled manner (it will write the latest snapshot). Again, the

action will happen the next time the code is writing monitoring output.

5.11 RERUN and NEWDIR ﬁles

After the code ﬁnishes (e.g., when the ﬁnal timestep number is reached or when a ‘STOP’

ﬁle is found), the ‘run.csh’ script checks whether there is a ‘RERUN’ ﬁle. If so, the code will

simply run again, perhaps even after you have recompiled the code. This is useful in the

development phase when you changed something in the code, so you don’t need to wait

for a new slot in the queue!

Even more naughty, as Tony says, is the ‘NEWDIR’ ﬁle, where you can enter a new direc-

tory path (relative path is ok, e.g. ../conv-slab). If nothing is written in this ﬁle (e.g.

via touch NEWDIR) it stays in the same directory. On distributed machines, the ‘NEWDIR’

method will copy all the ‘VAR#’ and ‘var.dat’ ﬁles back to and from the sever. This can be

useful if you want to run with new data ﬁles, but you better do it in a separate directory,

because with ‘NEWDIR’ the latest data from the code are written back to the server before

running again.

Oh, by the way, if you want to be sure that you haven’t messed up the content of the pair

of ‘NEWDIR’ ﬁles, you may want to try out the pc_jobtransfer command. It writes the

decisive ‘STOP’ ﬁle only after the script has checked that the content of the two ‘NEWDIR’

ﬁles points to existing run directory paths, so if the new run crashes, the code returns

safely to the old run directory. I’m not sure what Tony would say now, but this is now

obviously extremely naughty.

5.12 Start and run parameters

All input parameters in ‘start.in’ and ‘run.in’ are grouped in Fortran

namelists

. This

allows arbitrary order of the parameters (within the given namelist; the namelists need

no longer be in the correct order), as well as enhanced readability through inserted For-

tran comments and whitespace. One namelist (

init pars

run pars

) contains general

parameters for initialization/running and is always read in. All other namelists are spe-

ciﬁc to individual modules and will only be read if the corresponding module is used.

5.12 Start and run parameters 35

The syntax of a namelist (in an input ﬁle like ‘start.in’) is

&init_pars

ip=5, Lxyz=2,4,2

— in this example, the name of the namelist is

init pars

, and we read just two variables

(all other variables in the namelist retain their previous value):

, which is set to 5, and

Lxyz

, which is a vector of length three and is set to (2,4,2).

While all parameters from the namelists can be set, in most cases reasonable default

values are preset. Thus, the typical ﬁle ‘start.in’ will only contain a minimum set of

variables or (if you are very minimalistic) none at all. If you want to run a particular

problem, it is best to start by modifying an existing example that is close to your appli-

cation.

Before starting a simulation run, you may want to execute the command pc_configtest

in order to test the correctness of your changes to these conﬁguration ﬁles.

As an example, we give here the start parameters for ‘samples/helical-MHDturb’

&init_pars

cvsid=’$Id:$’, ! identify version of start.in

xyz0 = -3.1416, -3.1416, -3.1416, ! first corner of box

Lxyz = 6.2832, 6.2832, 6.2832, ! box size

lperi = T , T , T , ! periodic in x, y, z

random_gen=’nr_f90’

&hydro_init_pars

&density_init_pars

gamma=1.

&magnetic_init_pars

initaa=’gaussian-noise’, amplaa=1e-4

The three entries specifying the location, size and periodicity of the box are just given

for demonstration purposes here — in fact a periodic box from −πto −πin all three

directions is the default. In this run, for reproducibility, we use a random number gen-

erator from the Numerical Recipes [30], rather than the compiler’s built-in generator.

The adiabatic index γis set explicitly to 1(the default would have been 5/3) to achieve

an isothermal equation of state. The magnetic vector potential is initialized with uncor-

related, normally distributed random noise of amplitude 10−4.

The run parameters for ‘samples/helical-MHDturb’ are

&run_pars

cvsid=’$Id:$’, ! identify version of start.in

nt=10, it1=2, cdt=0.4, cdtv=0.80, isave=10, itorder=3

dsnap=50, dvid=0.5

random_gen=’nr_f90’

&hydro_run_pars

36 THE PENCIL CODE

&density_run_pars

&forcing_run_pars

iforce=’helical’, force=0.07, relhel=1.

&magnetic_run_pars

eta=5e-3

&viscosity_run_pars

nu=5e-3

Here we run for

= 10 timesteps, every second step, we write a line of diagnostic output;

we require the time step to keep the advective

Courant number

≤0.4and the diffusive

Courant number

≤0.8, save ‘var.dat’ every 20 time steps, and use the 3-step time-

stepping scheme described in Appendix H.4 (the Euler scheme

itorder

= 1 is only useful

for tests). We write permanent snapshot ﬁle ‘VARN’ every

dsnap

= 50 time units and 2d

slices for animation every

dvid

= 0.5time units. Again, we use a deterministic random

number generator. Viscosity νand magnetic diffusivity ηare set to 5×10−3(so the mesh

Reynolds number is about urmsδx/ν = 0.3×(2π/32)/5×10−3≈12, which is in fact rather a

bit to high). The parameters in

forcing run pars

specify fully helical forcing of a certain

amplitude.

A full list of input parameters is given in Appendix K.

5.13 Physical units

Many calculations are unit-agnostic, in the sense that all results remain the same in-

dependent of the unit system in which you interpret the numbers. E. g. if you simulate

a simple hydrodynamical ﬂow in a box of length L= 1.and get a maximum velocity of

umax = 0.5after t= 3 time units, then you may interpret this as L= 1 m,umax = 0.5 m/s,

t= 3 s, or as L= 1 pc,umax = 0.5 pc/Myr,t= 3 Myr, depending on the physical system you

have in mind. The units you are using must of course be consistent, thus in the second

example above, the units for diffusivities would be pc2/Myr, etc.

The units of magnetic ﬁeld and temperature are determined by the values µ0= 1 and

cp= 1 used internally by the code7. This means that if your units for density and velocity

are [ρ]and [v], then magnetic ﬁelds will be in

[B] = pµ0[ρ] [v]2,(18)

and temperatures are in

[T] = [v]2

=γ−1

[v]2

R/µ .(19)

For some choices of density and velocity units, Table 2 shows the resulting units of

magnetic ﬁeld and temperature.

7Note that cp= 1 is only assumed if you use the module

noionization.f90

. If you work with

ioniza-

tion.f90

, temperature units are speciﬁed by

unit temperature

as described below.

5.14 Minimum amount of viscosity 37

Table 2: Units of magnetic ﬁeld and temperature for some choices of [ρ]and [v]according to Eqs. (18) and

(19). Values are for a monatomic gas (γ= 5/3) of mean atomic weight ¯µg= ¯µ/1 g in grams.

[ρ] [v] [B] [T]

1 kg/m31 m/s 1.12 mT = 11.2 G ¯µg

0.6×2.89 ×10−5K

1 g/cm31 cm/s 3.54 ×10−4T = 3.54 G ¯µg

0.6×2.89 nK

1 g/cm31 km/s 35.4 T = 354 kG ¯µg

0.6×28.9 K

1 g/cm310 km/s 354 T = 3.54 MG ¯µg

0.6×2 890 K

On the other hand, as soon as material equations are used (e.g. one of the popular pa-

rameterizations for radiative losses, Kramers opacity, Spitzer conductivities or ioniza-

tion, which implies well-deﬁned ionization energies), the corresponding routines in the

code need to know the units you are using. This information is speciﬁed in ‘start.in’ or

‘run.in’ through the parameters

unit system

unit length

unit velocity

unit density

and

unit temperature

8like e. g.

unit_system=’SI’,

unit_length=3.09e16, unit_velocity=978. ! [l]=1pc, [v]=1pc/Myr

Note: The default unit system is unit_system=’cgs’ which is a synonym for unit_-

system=’Babylonian cubits’.

5.14 Minimum amount of viscosity

We emphasize that, by default, the code works with constant diffusion coefﬁcients (vis-

cosity ν, thermal diffusivity χ, magnetic diffusivity η, or passive scalar diffusivity D).

If any of these numbers is too small, you would need to have more meshpoints to get

consistent numerical solutions; otherwise the code develops wiggles (‘ringing’) and will

eventually crash. A useful criterion is given by the mesh Reynolds number based on the

maximum velocity,

Remesh = max(|u|) max(δx, δy, δz)/ν, (20)

which should not exceed a certain value which can be problem-dependent. Often the

largest possible value of Remesh is around 5. Similarly there exist mesh P´

eclet and mesh

magnetic Reynolds numbers that should not be too large.

Note that in some cases, ‘wiggles’ in ln ρwill develop despite sufﬁciently large diffusion

coefﬁcients, essentially because the continuity equation contains no dissipative term.

For convection runs (but not only for these), we have found that this can often be pre-

vented by

upwinding

, see Sect. H.2.

If the Mach number of the code approaches unity, i.e. if the rms velocity becomes com-

parable with the speed of sound, shocks may form. In such a case the mesh Reynolds

number should be smaller. In order to avoid excessive viscosity in the unshocked regions,

8Note: the parameter

unit temperature

is currently only used in connection with

ionization.f90

. If you

are working with

noionization.f90

, the temperature unit is completely determined by Eq. (19) above.

38 THE PENCIL CODE

one can use the so-called shock viscosity (Sect. 6.6.1) to concentrate the effects of a low

mesh Reynolds number to only those areas where it is necessary.

5.15 The time step

5.15.1 The usual RK-2N time step

RK-2N refers to the third order Runge-Kutta scheme by Williamson (1980) with a mem-

ory consumption of two chunks. Therefore the 2N in the name.

The time step is normally speciﬁed as Courant time step through the coefﬁcients

cdt

(cδt),

cdtv

(cδt,v) and

cdts

(cδt,s). The resulting Courant step is given by

δt = min cδt

δxmin

Umax

, cδt,v

δx2

min

Dmax

, cδt,s

Hmax ,(21)

where

δxmin ≡min(δx, δy, δz) ; (22)

Umax ≡max |u|+qc2

s+v2

A,(23)

csand vAdenoting sound speed and Alfv´

en speed, respectively;

Dmax = max(ν, γχ, η, D),(24)

where νdenotes kinematic viscosity, χ=K/(cpρ)thermal diffusivity and ηthe magnetic

diffusivity; and

Hmax = max 2νS2+ζshock(∇·u)2+...

cvT,(25)

where dots indicate the presence of other terms on the rhs of the entropy equation.

To ﬁx the time step δt to a value independent of velocities and diffusivities, explicitly set

the run parameter

, rather than

cdt

cdtv

(see p. 180).

If the time step exceeds the viscous time step the simulation may actually run ok for

quite some time. Inspection of images usually helps to recognize the problem. An exam-

ple is shown in Fig. 4.

Timestepping is accomplished using the Runge-Kutta 2N scheme. Regarding details of

this scheme see Sect. H.4.

5.15.2 The Runge-Kutta-Fehlberg time step

A ﬁfth order Runge-Kutta-Fehlberg time stepping procedure is available. It is used

mostly for chemistry application, often together with the double precision option. In

order to make this work, you need to compile with

TIMESTEP = timestep_rkf

in ‘src/Makefile.local’. In addition, you must put itorder=5 in ‘run.in’. An example

application is ‘samples/1d-tests/H2_flamespeed’. This procedure is still experimental.

5.16 Boundary conditions 39

Figure 4: Example of a velocity slice from a run where the time step is too long. Note the spurious

checkerboard modulation in places, for example near x=−0.5and −2.5< y < −1.5. This is an example

of a hyperviscous turbulence simulations with 5123meshpoints and a third order hyperviscosity of ν3=

5×10−12. Hyperviscosity is explained in the Appendix E.

5.16 Boundary conditions

5.16.1 Where to specify boundary conditions

In most tests that come with the PENCIL CODE, boundary conditions are set in ‘run.in’,

which is a natural choice. However, this may lead to unexpected initial data written by

‘start.x’, since when you start the code (via ‘start.csh’), the boundary conditions are

unknown and ‘start.x’ will then ﬁll the ghost zones assuming periodicity (the default

boundary condition) in all three directions. These ghost data will never be used in a

calculation, as ‘run.x’ will apply the boundary conditions before using any ghost-zone

values.

To avoid these periodic conditions in the initial snapshot, you can set the boundary

conditions in ‘start.in’ already. In this case, they will be inherited by ‘run.x’, unless you

also explicitly set boundary conditions in ‘run.in’.

5.16.2 How to specify boundary conditions

Boundary conditions are implemented through three layers of

ghost points

on either

boundary, which is quite a natural choice for an MPI code that uses ghost zones for

representing values located on the neighboring processors anyway. The desired type of

boundary condition is set through the parameters

{

x,y,z

}in ‘run.in’; the nomenclature

used is as follows. Set

{

x,y,z

}to a sequence of letters like

bcx = ’p’,’p’,’p’, ’p’, ’p’

for periodic boundaries, or

bcz = ’s’,’s’,’a’,’a2’,’c1:c2’

for non-periodic ones. Each element corresponds to one of the variables, which are those

of the variables ux,uy,uz,ln ρ,s/cp,Ax,Ay,Az,ln cthat are actually used in this order.

The following conditions are available:

‘p’periodic boundary condition

‘a’antisymmetric condition w. r. t. the boundary, i. e. vanishing value

40 THE PENCIL CODE

‘s’symmetric condition w. r. t. the boundary, i. e. vanishing ﬁrst derivative

‘a2’antisymmetry w. r. t. the arbitrary value on the boundary, i. e. vanishing second

derivative

‘c1’special boundary condition for ln ρand s: constant heat ﬂux through the boundary

‘c2’special boundary condition for s: constant temperature at the boundary — requires

boundary condition a2 for ln ρ

‘cT’special boundary condition for sor ln T: constant temperature at the boundary (for

arbitrarily set ln ρ)

‘ce’special boundary condition for s: set temperature in ghost points to value on bound-

ary (for arbitrarily set ln ρ)

‘db’low-order one-sided derivatives (“no boundary condition”) for density

‘she’shearing-sheet boundary condition (default when the module

Shear

is used)

‘g’force the value of the corresponding ﬁeld on vertical boundaries (should be used

in combination with the force lower bound and force upper bound ﬂags set in the

namelist init pars)

‘hs’special boundary condition for ln ρand swhich enforces hydrostatic equilibrium on

vertical boundaries

The extended syntax a:b(e. g. ‘c1:c2’) means: use boundary condition aat the left/lower

boundary, but bat the right/upper one.

If you build a new ‘run.in’ ﬁle from another one with a different number of variables

(

noentropy

vs.

entropy

, for example), you need to remember to adjust the length of the

arrays

bcx

bcz

. The advantage of the present approach is that it is very easy to ex-

change all boundary conditions by a new set of conditions in a particular direction (for

example, make everything periodic, or switch off shearing sheet boundary conditions

and have stress-free instead).

5.17 Restarting a simulation

When a run stops at the end of a simulation, you can just resubmit the job again, and

it will start from the latest snapshot saved in ‘data/proc*/var.dat’. The value of the

latest time is saved in a separate ﬁle, ‘data/proc*/time.dat’. On parallel machines it is

possible that some (or just one) of the ‘var.dat’ are corrupt; for example after a system

crash. Check for ﬁle size and date, and restart from a good ‘VAR’Nﬁle instead.

If you want to run on a different machine, you just need to copy the ‘data/proc*/var.dat’

(and, just to be sure) ‘data/proc*/time.dat’) ﬁles into a new directory tree. You may also

need the ‘data/proc*/seed.dat’ ﬁles for the random number generator. The easiest way

to get all these other ﬁles is to run start.csh again on the new machine (or in a new

directory) and then to overwrite the ‘data/proc*/var.dat’ ﬁles with the correct ones.

For restarting from runs that didn’t have magnetic ﬁelds, passive scalar ﬁelds, or test

ﬁelds, see Sect. F.3.

5.18 One- and two-dimensional runs 41

5.18 One- and two-dimensional runs

If you want to run two-dimensional problems, set the number of mesh points in one

direction to unity, e.g.

nygrid

=1 or

nzgrid

=1 in ‘cparam.local’. Remember that the num-

ber of mesh points is still divisible by the number of processors. For 2D-runs, it is also

possible to write only 2D-snapshots (i.e. VAR ﬁles written only in the considered (x, y)

or (x, z)plane, with a size seven times smaller as we do not write the third unused di-

rection). To do that, please add the logical ﬂag ‘lwrite 2d=T’ in the namelist init pars in

‘start.in’.

Similarly, for one-dimensional problems, set, for example,

nygrid

=1 and

nzgrid

=1 in

‘cparam.local’. You can even do a zero-dimensional run, but then you better set

(rather than

cdt

), because there is no Courant condition for the time step.

See 0d, 1d, 2d, and 3d tests with examples.

5.19 Visualization

5.19.1 Gnuplot

Simple visualization can easily be done using

Gnuplot

(http://www.gnuplot.info), an

open-source plotting program suitable for two-dimensional plots.

For example, suppose you have the variables

---it-----t-------dt-------urms-----umax-----rhom-----ssm------dtc---

in ‘time_series.dat’ and want to plot urms(t). Just start gnuplot and type

gnuplot> plot "data/time_series.dat" using 2:4 with lines

If you work over a slow line and want to see both urms(t)and umax(t), use ASCII graphics:

gnuplot> set term dump

gnuplot> set logscale y

gnuplot> plot "data/time_series.dat" using 2:4 title "urms", \

gnuplot> "data/time_series.dat" using 2:5 title "umax"

5.19.2 Data explorer

(

data explorer

;http://www.opendx.org) is an open-source tool for visualization of

three-dimensional data.

The PENCIL CODE provides a few networks for

. It is quite easy to read in a snapshot

ﬁle from

(the only tricky thing is the four extra bytes at the beginning of the ﬁle,

representing a Fortran record marker), and whenever you run ‘start.x’, the code writes

a ﬁle ‘var.general’ that tells

all it needs to know about the data structure.

As a starting point for developing your own

programs or

networks

, you can use a

few generic

scripts provided in the directory ‘dx/basic/’. From the run directory,

start

with

unix> dx -edit $PENCIL_HOME/dx/basic/lnrho

42 THE PENCIL CODE

to load the ﬁle ‘dx/basic/lnrho.net’, and execute it with ✄

✂



✁

Ctl-o or Execute →Execute

Once. You will see a set of iso-surfaces of logarithmic density. If the viewport does not

ﬁt to your data, you can reset it with ✄

✂



✁

Ctl-f. To rotate the object, drag the mouse over

the Image window with the left or right mouse button pressed. Similar networks are

provided for entropy (‘ss.net’), velocity (‘uu.net’) and magnetic ﬁeld (‘bb.net’).

When you expand these simple networks to much more elaborate ones, it is probably

a good idea to separate the different tasks (like Importing and Selecting, visualizing

velocity, visualizing entropy, and Rendering) onto separate pages through Edit →Page.

Note Currently,

can only read in data ﬁles written by one single processor, so from

a multi-processor run, you can only visualize one subregion at a time.

5.19.3 GDL

GDL

, also known as

Gnu Data Language

is a free visualization package that can be

found at http://gnudatalanguage.sourceforge.net/. It aims at replacing the very ex-

pensive

IDL

package (see S. 5.19.4). For the way we use IDL for the Pencil Code, com-

patibility is currently not completely sufﬁcient, but you can use GDL for many of the

visualization tasks. If you get spurious “Error opening ﬁle” messages, you can normally

simply ignore them.

This section tells you how to get started with using GDL for visualization.

Setup As of GDL 0.9 – at least the version packed with Ubuntu Jaunty (9.10) – you

will need to add GDL’s ‘examples/pro/’ directory to your

!PATH

variable. So the ﬁrst call

after starting

GDL

should be

GDL> .run setup_gdl

Starting visualization There are mainly two possibilities for visualization: using a sim-

ple GUI or loading the data with pc_read and work with it interactively. Please note that

the GUI was written and tested only with IDL, see §5.19.4.

Here, the pc_read family of IDL routines to read the data is described. Try

GDL> pc_read

to get an overview.

To plot a time series, use

GDL> pc_read_ts, OBJECT=ts

GDL> help, ts, /STRUCT ;; (to see which slots are available)

GDL> plot, ts.t, ts.umax

GDL> oplot, ts.t, ts.urms

Alternatively, you could simply use the ‘ts.pro’ script:

GDL> .run ts

To work with data from ‘var.dat’ and similar snapshot ﬁles, you can e.g. use the follow-

ing routines:

5.19 Visualization 43

GDL> pc_read_dim, OBJECT=dim

GDL> $$PENCIL_HOME/bin/nl2idl -d ./data/param.nml> ./param.pro

GDL> pc_read_param, OBJECT=par

GDL> pc_read_grid, OBJECT=grid

GDL> pc_read_var, OBJECT=var

Having thus read the data structures, we can have a look at them to see what informa-

tion is available:

GDL> help, dim, /STRUCT

GDL> help, par, /STRUCT

GDL> help, grid, /STRUCT

GDL> help, var, /STRUCT

To visualize data, we can e.g. do9

GDL> plot, grid.x, var.ss[*, dim.ny/2, dim.nz/2]

GDL> contourfill, var.ss[*, *, dim.nz/2], grid.x, grid.y

GDL> ux_slice = var.uu[*, *, dim.nz/2, 0]

GDL> uy_slice = var.uu[*, *, dim.nz/2, 1]

GDL> wdvelovect, ux_slice, uy_slice, grid.x, grid.y

GDL> surface, var.lnrho[*, *, dim.nz/2, 0]

See also Sect. 5.19.4.

5.19.4 IDL

IDL

is a commercial visualization program for two-dimensional and simple three-

dimensional graphics. It allows to access and manipulate numerical data in a fashion

quite similar to how Fortran handles them.

In ‘$PENCIL_HOME/idl’, we provide a number of general-purpose

IDL

scripts that we are

using all the time in connection with the PENCIL CODE. While

IDL

is quite an expensive

software package, it is quite useful for visualizing results from numerical simulations.

In fact, for many applications, the 7-minute demo version of

IDL

is sufﬁcient.

Visualization in IDL The Pencil Code GUI is a data post-processing tool for the usage

on a day-to-day basis. It allows fast inspection of many physical quantities, as well as ad-

vanced features like horizontal averages, streamline tracing, freely orientable 2D-slices,

and extraction of cut images and movies. To use the Pencil Code GUI, it is sufﬁcient to

run:

unix> idl

IDL> .r pc_gui

If you like to load only some subvolume of the data, like any 2D-slices from the given

data snapshots, or 3D-subvolumes, it is possible to choose the corresponding options in

the ﬁle selector dialog. The Pencil Code GUI offers also some options to be set on the

command-line, please refer to their description in the source code.

9If contourfill produces just contour lines instead of a color-coded plot, your version of GDL is too

old. E.g. the version shipped with Ubuntu 9.10 is based on GDL 0.9rc1 and has this problem.

44 THE PENCIL CODE

There are also other small GUIs available, e.g. the ﬁle ‘time-series.dat’ can easily be

analyzed with the command:

unix> idl

IDL> pc_show_ts

The easiest way to derive physical quantities at the command-line is to use one of the

many pc_read_var-variants (pc_read_var_raw is recommended for large datasets be-

cause of its high efﬁciency regarding computation and memory usage) for reading the

data. With that, one can make use of pc_get_quantity to derive any implemented physi-

cal quantity. Packed in a script, this is the recommended way to get reproducible results,

without any chance for accidental errors on the interactive IDL command-line.

Alternatively, by using the command-line to see the time evolution of e.g. velocity and

magnetic ﬁeld (if they are present in you run), start

IDL

10 and run ‘ts.pro’:

unix> idl

IDL> .run ts

The

IDL

script ‘ts.pro’ reads the time series data from ‘data/time_series.dat’ and sorts

the column into the structure

, with the slot names corresponding to the name of the

variables (taken from the header line of ‘data/time_series.dat’). Thus, you can refer to

time as ts.t, to the rms velocity as ts.urms, and in order to plot the mean density as a

function of time, you would simply type

IDL> plot, ts.t, ts.rhom

The basic command sequence for working with a snapshot is:

unix> idl

IDL> .run start

IDL> .run r

IDL>

[speciﬁc commands]

You run ‘start.pro’ once to initialize (or reinitialize, if the mesh size has changed, for

example) the ﬁelds and read in the startup parameters from the code. To read in a new

snapshot, run ‘r.pro’ (or ‘rall.pro’, see below).

If you are running in parallel on several processors, the data are scattered over different

directories. To reassemble everything in

IDL

, use

10 If you run IDL from the command line, you will highly appreciate the following tip: IDL’s command

line editing is broken beyond hope. But you can ﬁx it, by running IDL under rlwrap, a wrapper for the

excellent GNU

readline

library.

Download and install rlwrap from http://utopia.knoware.nl/~hlub/uck/rlwrap/ (on some systems

you just need to run ‘emerge rlwrap’, or ‘apt-get install rlwrap’), and alias your idl command:

csh> alias idl ’rlwrap -a -c idl’

bash> alias idl=’rlwrap -a -c idl’

From now on, you can

•use long command lines that correctly wrap around;

•type the ﬁrst letters of a command and then ✄

✂



✁

PageUp to recall commands starting with these letters;

•capitalize, uppercase or lowercase a word with ✄

✂



✁

Esc-C, ✄

✂



✁

Esc-U, ✄

✂



✁

Esc-L;

•use command line history across IDL sessions (you might need to create ‘~/.idl_history’ for this);

•complete ﬁle names with ✄

✂



✁

Tab (works to some extent);

•. . . use all the other

readline

features that you are using in bash,octave,bc,gnuplot,ftp, etc.

5.19 Visualization 45

IDL> .r rall

instead of .r r (here, .r is a shorthand for .run). The procedure ‘rall.pro’ reads (and

assembles) the data from all processors and correspondingly requires large amounts of

memory for very large runs. If you want to look at just the data from one processor, use

‘r.pro’ instead.

If you need the magnetic ﬁeld or the current density, you can calculate them in IDL by

IDL> bb=curl(aa)

IDL> jj=curl2(aa)

By default one is reading always the current snapshot ‘data/procN/var.dat’; if you want

to read one of the permanent snapshots, use (for example)

IDL> varfile=’VAR2’

IDL> .r r

(or

.r rall

)

See Sect. 5.6.1 for details on permanent snapshots.

With ‘r.pro’, you can switch the part of the domain by changing the variable

datadir

IDL> datadir=’data/proc3’

IDL> .r r

will read the data written by processor 3.

Reading data into IDL arrays or structures As an alternative to the method described

above, there is also the possibility to read the data into structures. This makes some

more operations possible, e.g. reading data from an IDL program where the command

.r is not allowed.

An efﬁcient and still scriptable way would look like the following:

IDL> pc_read_var_raw, obj=var, tags=tags

IDL> bb = pc_get_quantity (’B’, var, tags)

IDL> jj = pc_get_quantity (’j’, var, tags)

This reads the data into an array ‘var’, as well as the array indices of the contained phys-

ical variables into a separate structure ‘tags’. To use a caching mechanism within pc_-

get_quantity, please refer to the documentation and the examples contained in ‘pencil-

code/idl/pc get quantity.pro’, where you can also start adding not yet implemented phys-

ical quantities.

To read a snapshot ’VAR10’ into the IDL structure ff, type the following command

IDL> pc_read_var, obj=ff, varfile=’VAR10’, /trimall

The option /trimall removes ghost zones from the data. A number of other options are

documented in the source code of pc_read_var. You can see what data the structure

contains by using the command tag_names

IDL> print, tag_names(ff)

T X Y Z DX DY DZ UU LNRHO AA

11 Keep in mind that jj=curl(bb) would use iterated ﬁrst derivatives instead of the second derivatives

and thus be numerically less accurate than jj=curl2(aa), particularly at small scales.

46 THE PENCIL CODE

One can access the individual variables by typing ff.varname, e.g.

IDL> help, ff.aa

<Expression> FLOAT = Array[32, 32, 32, 3]

There are a number of ﬁles that read different data into structures. They are placed in

the directory $PENCIL_HOME/idl/files. Here is a list of ﬁles (including suggested options

to call them with)

•pc_read_var_raw, obj=var, tags=tags

Efﬁciently read a snapshot into an array.

•pc_read_var, obj=ff, /trimall

Read a snapshot into a structure.

•pc_read_ts, obj=ts

Read the time series into a structure.

•pc_read_xyaver, obj=xya

pc_read_xzaver, obj=xza

pc_read_yzaver, obj=yza

Read 1-D time series into a structure.

•pc_read_const, obj=cst

Read code constants into a structure.

•pc_read_pvar, obj=fp

Read particle data into a structure.

•pc_read_param, obj=par

Read startup parameters into a structure.

•pc_read_param, obj=par2, /param2

Read runtime parameters into a structure.

Other options are documented in the source code of the ﬁles.

For some examples on how to use these routines, see Sect. 5.19.3.

5.19.5 Python

Pencil supports reading, processing and the visualization of data using

python

. A num-

ber of scripts are placed in the subfolder ‘$PENCIL_HOME/python’. A README ﬁle placed

under that subfolder contains the information needed to read Pencil output data into

python.

Installation For modern operating systems, Python is generally installed together with

the system. If not, it can be installed via your preferred package manager or downloaded

from the website https://www.python.org/. For convenience, it is strongly recommend

to also install IPython, which is a more convenient console for Python. You will also need

the NumPy, matplotlib and Tk library.

Perhaps the easiest way to obtain all the required software mentioned above is to install

either Continuum’s Anaconda or Enthought’s Canopy. These Python distributions also

provide (or indeed are) integrated graphical development environments.

5.19 Visualization 47

In order for Python to ﬁnd the Pencil Code commands you will have to add to your

.bashrc:

export PYTHONPATH=$PENCIL_HOME/python

ipythonrc If you use IPython, for convenience, you should modify your

.ipython/ipythonrc (create it if it doesn’t exist) and add:

import_all pencil

Additional, add to your .ipython/profile_default/startup/init.py the following lines:

import numpy as np

import pylab as plt

import pencil as pc

import matplotlib

from matplotlib import rc

plt.ion()

matplotlib.rcParams[’savefig.directory’] = ’’

.pythonrc In case you are on a cluster and don’t have access to IPython you can edit

your .pythonrc:

#!/usr/bin/python

import numpy as np

import pylab as plt

import pencil as pc

import atexit

#import readline

import rlcompleter

# enables search with CTR+r in the history

try:

import readline

except ImportError:

print "Module readline not available."

else:

import rlcompleter

readline.parse_and_bind("tab: complete")

# enables command history

historyPath = os.path.expanduser("~/.pyhistory")

def save_history(historyPath=historyPath):

import readline

readline.write_history_file(historyPath)

if os.path.exists(historyPath):

readline.read_history_file(historyPath)

atexit.register(save_history)

del os, atexit, readline, rlcompleter, save_history, historyPath

plt.ion()

48 THE PENCIL CODE

create the ﬁle .pythonhistory and add to your .bashrc:

export PYTHONSTARTUP=~/.pythonrc

Pencil Code Commands in General For a list of all Pencil Code commands start

IPython and type pc. <TAB> (as with auto completion). To access the help of any com-

mand just type the command followed by a ’?’ (no spaces), e.g.:

In [1]: pc.dot?

Type: function

String Form:<function dot at 0x7f9d96cb0cf8>

File: ~/pencil-code/python/pencil/math/vector_multiplication.py

Definition: pc.dot(a, b)

Docstring:

take dot product of two pencil-code vectors a & b with shape

a.shape = (3,mz,my,mx)

You can also use help(pc.dot) for a more complete documentation of the command.

There are various reading routines for the Pencil Code data. All of them return an object

with the data. To store the data into a user deﬁned variable type e.g.

In [1]: ts = pc.read_ts()

Most commands take some arguments. For most of them there is a default value, e.g.

In [1]: pc.read_ts(filename=’time_series.dat’, datadir=’data’, plot_data=True)

You can change the values by simply typing e.g.

In [1]: pc.read_ts(plot_data = False)

Reading and Plotting Time Series Reading the time series ﬁle is very easy. Simply

type

In [1]: ts = pc.read_ts()

and Python stores the data in the variable ts. The physical quantities are members of

the object ts and can be accessed accordingly, e.g. ts.t,ts.emag. To check which other

variables are stored simply do the tab auto completion ts. <TAB>.

Plot the data with the matplotlib commands:

In [1]: plt.plot(ts.t, ts.emag)

The standard plots are not perfect and need a little polishing. See

the online wiki for a few examples on how to make pretty plots

(https://github.com/pencil-code/pencil-code/wiki/PythonForPencil). You can

save the plot into a ﬁle using the GUI or with

In [1]: plt.savefig(’plot.eps’)

Reading and Plotting VAR ﬁles and slice ﬁles Read var ﬁles:

In [1]: var = pc.read_var()

5.20 Running on multi-processor computers 49

Read slice ﬁles:

In [1]: slices, t = pc.read_slices(field=’bb1’, extension=’xy’)

This returns the array slices with indices slices[nTimes, my, mx] and the time array

If you want to plot e.g. the x-component of the magnetic ﬁeld at the central plane simply

type:

In [1]: plt.imshow(zip(*var.bb[0, 128, :, :]), origin=’lower’, extent=[-4, 4, -4, 4])

For a complete list of arguments of plt.imshow refer to its documentation.

For a more interactive plot use:

In [1]: pc.animate_interactive(slices, t)

Be aware: arrays from the reading routines are ordered f[nvar, mz, my, mx], i.e. re-

versed to IDL. This affects reading var ﬁles and slice ﬁles.

5.20 Running on multi-processor computers

The code is parallelized using

MPI

(

message passing interface

) for a simple domain

decomposition (data-parallelism), which is a straight-forward and very efﬁcient way of

parallelizing ﬁnite-difference codes. The current version has a few restrictions, which

need to be kept in mind when using the MPI features.

The global number of grid points (but excluding the ghost zones) in a given direction

must be an exact multiple of the number of processors you use in that direction. For

example if you have nprocy=8 processors for the ydirection, you can run a job with

nygrid=64 points in that direction, but if you try to run a problem with nygrid=65 or

nygrid=94, the code will complain about an inconsistency and stop. (So far, this has not

been a serious restriction for us.)

5.20.1 How to run a sample problem in parallel

To run the sample problem in the directory ‘samples/conv-slab’ on 16 CPUs, you need

to do the following (in that directory):

1. Edit ‘src/Makefile.local’ and replace

MPICOMM = nompicomm

MPICOMM = mpicomm

2. Edit ‘src/cparam.local’ and replace

integer, parameter :: ncpus=1, nprocy=1, nprocz=ncpus/nprocy, nprocx=1

integer, parameter :: nxgrid=32, nygrid=nxgrid, nzgrid=nxgrid

integer, parameter :: ncpus=16, nprocy=4, nprocz=ncpus/nprocy, nprocx=1

integer, parameter :: nxgrid=128, nygrid=nxgrid, nzgrid=nxgrid

50 THE PENCIL CODE

The ﬁrst line speciﬁes a 4×4layout of the data in the yand zdirection. The second

line increases the resolution of the run because running a problem as small as 323

on 16 CPUs would be wasteful. Even 1283may still be quite small in that respect.

For performance timings, one should try and keep the size of the problem per CPU

the same, so for example 2563on 16 CPUs should be compared with 1283on 2 CPUs.

3. Recompile the code

unix> (cd src; make)

4. Run it

unix> start.csh

unix> run.csh

Make sure that all CPUs see the same ‘data/’ directory; otherwise things will go wrong.

Remember that in order to visualize the full domain with IDL (rather than just the

domain processed and written by one processor), you need to use ‘rall.pro’ instead of

‘r.pro’.

5.20.2 Hierarchical networks (e.g. on Beowulf clusters)

On big Beowulf clusters, a group of nodes is often connected with a switch of modest

speed, and all these groups are connected via a ntimes faster uplink switch. When

bandwidth-limited, it is important to make sure that consecutive processors are mapped

onto the mesh such that the load on the uplink is .ntimes larger than the load on the

slower switch within each group. On a 512 node cluster, where groups of 24 processors

are linked via fast ethernet switches, which in turn are connected via a Gigabit uplink

(∼10 times faster), we found that

nprocy

=4 is optimal. For 128 processors, for example

we ﬁnd that

nprocy

nprocz

= 4×32 is the optimal layout, while. For comparison, 8×16

is 3 times slower, and 16 ×8is 17 (!) times slower. These results can be understood from

the structure of the network, but the basic message is to watch out for such effects and

to try varying

nprocy

and

nprocz

In cases where nygrid>nzgrid it may be advantageous to swap the ordering of processor

numbers. This can be done by setting

lprocz slowest

=F.

5.20.3 Extra workload caused by the ghost zones

Normally, the workload caused by the ghost zones is negligible. However, if one increases

the number of processors, a signiﬁcant fraction of work is done in the ghost zones. In

other words, the effective mesh size becomes much larger than the actual mesh size.

Consider a mesh of size Nw=Nx×Ny×Nz, and distribute the task over Pw=Px×Py×

Pzprocessors. If no communication were required, the number of points per processor

would be Nw

=Nx×Ny×Nz

Px×Py×Pz

.(26)

However, for ﬁnite difference codes some communication is required, and the amount of

communication depends on spatial order of the scheme, Q. The PENCIL CODE works by

5.21 Running in double-precision 51

default with sixth order ﬁnite differences, so Q= 6, i.e. one needs 6 ghost zones, 3 on

each end of the mesh. With Q6= 0 the number of points per processor is

N(eﬀ)

=Nx

+Q×Ny

+Q×Nz

+Q.(27)

There is efﬁcient scaling only when

min Nx

,Ny

,Nz

Pz≫Q. (28)

In the special case were Nx=Ny=Nz≡N=N1/3

w, with Px= 1 and Py=Pz≡P=P1/2

we have

N(eﬀ)

= (N+Q)×N

P+Q2

.(29)

For N= 128 and Q= 6 the effective mesh size exceeds the actual mesh size by a factor

N(eﬀ)

= (N+Q)×N

P+Q2

×Pw

.(30)

These factors are listed in Table 3.

Table 3: N(eﬀ)

w/Nwversus Nand P.

P\N128 256 512 1024 2048

1 1.15 1.07 1.04 1.02 1.01

2 1.19 1.09 1.05 1.02 1.01

4 1.25 1.12 1.06 1.03 1.01

8 1.34 1.16 1.08 1.04 1.02

16 1.48 1.22 1.11 1.05 1.03

32 1.68 1.31 1.15 1.07 1.04

64 1.98 1.44 1.21 1.10 1.05

128 2.45 1.64 1.30 1.14 1.07

256 3.21 1.93 1.43 1.20 1.10

512 4.45 2.40 1.62 1.29 1.14

Ideally, one wants to keep the work in the ghost zones at a minimum. If one accepts

that 20–25% of work are done in the ghost zones, one should use 4 processors for 1283

mesh points, 16 processors for 2563mesh points, 64 processors for 5123mesh points, 256

processors for 10243mesh points, and 512 processors for 15363mesh points.

5.21 Running in double-precision

With many compilers, you can easily switch to double precision (8-byte ﬂoating point

numbers) as follows.

Add the lines

# Use double precision

REAL_PRECISION = double

52 THE PENCIL CODE

to ‘src/Makefile.local’ and (re-)run pc_setupsrc.

REAL PRECISION

is set to ‘double’, the ﬂag

FFLAGS DOUBLE

is appended to

the Fortran compile ﬂags. The default for

FFLAGS DOUBLE

is -r8, which works for

g95

or ifort; for

gfortran

, you need to make sure that

FFLAGS DOUBLE

is set to

-fdefault-real-8.

You can see the ﬂags in ‘src/Makefile.inc’, which is the ﬁrst place to check if you have

problems compiling for double precision.

Using double precision might be important in turbulence runs where the resolution

is 2563meshpoints and above (although such runs often seem to work ﬁne at single

precision). To continue working in double precision, you just say lread_from_other_-

prec=T in run_pars; see Sect. F.1.

5.22 Power spectrum

Given a real variable u, its Fourier transform ˜uis given by

˜u(kx, ky, kz) = F(u(x, y, z)) = 1

NxNyNz

Nx−1

p=0

Ny−1

q=0

Nz−1

r=0

u(xp, yq, zr)

×exp(−ikxxp) exp(−ikyyq) exp(−ikzzr),(31)

where

|kx|<πNx

,|ky|<πNy

,|kz|<πNz

when Lis the size of the simulation box. The three-dimensional power spectrum P(k)is

deﬁned as

P(k) = 1

2˜u˜u∗,(32)

where

k=qk2

x+k2

y+k2

z.(33)

Note that we can only reasonably calculate P(k)for k < πNx/Lx.

To get power spectra from the code, edit ‘run.in’ and add for example the following lines

dspec=5., ou_spec=T, ab_spec=T !(for energy spectra)

oned=T

under run_pars. The kinetic (vel_spec) and magnetic (mag_spec) power spectra will

now be calculated every 5.0 (dspec) time units. The Fourier spectra is calculated using

fftpack

. In addition to calculating the three-dimensional power spectra also the one-

dimensional power spectra will be calculated (oned).

In addition one must edit ‘src/Makefile.local’ and add the lines

FOURIER = fourier_fftpack

POWER = power_spectrum

Running the code will now create the ﬁles ‘powerhel_mag.dat’ and ‘power_kin.dat’

containing the three-dimensional magnetic and kinetic power spectra respectively. In

addition to these three-dimensional ﬁles we will also ﬁnd the one-dimensional ﬁles

‘powerbx_x.dat’, ‘powerby_x.dat’, ‘powerbz_x.dat’, ‘powerux_x.dat’, ‘poweruy_x.dat’ and

5.22 Power spectrum 53

‘poweruz_x.dat’. In these ﬁles the data are stored such that the ﬁrst line contains the

time of the snapshot, the following

nxgrid

/2numbers represent the power at each

wavenumber, from the smallest to the largest. If several snapshots have been saved,

they are being stored immediately following the preceding snapshot.

You can read the results with the idl procure ‘power’, like this:

power,’_kin’,’_mag’,k=k,spec1=spec1,spec2=spec2,i=n,tt=t,/noplot

power,’hel_kin’,’hel_mag’,k=k,spec1=spec1h,spec2=spec2h,i=n,tt=t,/noplot

If powerhel is invoked, krms is written during the ﬁrst computation. The relevant out-

put ﬁle is ‘power_krms.dat’. This is needed for a correct calculation of kused in the

realizability conditions.

A caveat of the implementation of Fourier transforms in the PENCIL CODE is that, due

to the parallelization, the permitted resolution is limited to the case when one direction

is an integer multiple of the other. So, it can be done for

Nx = n*Ny

Unfortunately, for some applications one wants Nx < Ny. Wlad experimented with ar-

bitrary resolution by interpolating xto the same resolution of yprior to transposing,

then transform the interpolated array and then interpolating it back (check ‘fourier_-

transform_y’ in ‘fourier_fftpack.f90’).

A feature of our current implementation with xparallelization is that fft_xyz_-

parallel_3D requires nygrid to be an integer multiple of nprocy*nprocz. Examples of

good mesh layouts are listed in Table 4.

Table 4: Examples of mesh layouts for which Fourier transforms with xparallelization is possible.

ny nprocyx nprocy nprocz ncpus

256 1 16 16 256

256 2 16 16 512

256 4 16 16 1024

256 8 16 16 2048

288 2 16 18 576

512 2 16 32 1024

512 4 16 16 1024

512 4 16 32 2048

576 4 18 32 2304

576 8 18 32 4608

576 16 18 32 9216

1024 4 32 32 4096

1024 4 16 64 4096

1024 8 16 32 4096

1152 4 36 32 4608

1152 4 32 36 4608

2304 2 32 72 4608

2304 4 36 64 9216

2304 4 32 72 9216

To visualize with

IDL

just type power and you get the last snapshot of the three-

54 THE PENCIL CODE

dimensional power spectrum. See head of ‘$PENCIL_HOME/idl/power.pro’ for options to

power.

5.23 Structure functions

We deﬁne the p-th order longitudinal structure function of uas

long(l) = h|ux(x+l, y, z)−ux(x, y, z)|pi,(34)

while the transverse is

trans(l) = h|uy(x+l, y, z)−uy(x, y, z)|pi+h|uz(x+l, y, z)−uz(x, y, z)|pi.(35)

Edit ‘run.in’ and add for example the following lines

dspec=2.3,

lsfu=T, lsfb=T, lsfz1=T, lsfz2=T

under run_pars. The velocity (lsfu), magnetic (lsfb) and Elsasser (lsfz1 and lsfz2)

structure functions will now be calculated every 2.3 (dspec) time unit.

In addition one must edit ‘src/Makefile.local’ and add the line

STRUCT_FUNC = struct_func

In ‘src/cparam.local’, deﬁne

lb nxgrid

and make sure that

nxgrid = nygrid = nzgrid = 2**lb_nxgrid

E.g.

integer, parameter :: lb_nxgrid=5

integer, parameter :: nxgrid=2**lb_nxgrid,nygrid=nxgrid,nzgrid=nxgrid

Running the code will now create the ﬁles:

‘sfu-1.dat’, ‘sfu-2.dat’, ‘sfu-3.dat’ (velocity),

‘sfb-1.dat’, ‘sfb-2.dat’, ‘sfb-3.dat’ (magnetic ﬁeld),

‘sfz1-1.dat’, ‘sfz1-2.dat’, ‘sfz1-3.dat’ (ﬁrst Elsasser variable),

‘sfz2-1.dat’, ‘sfz2-2.dat’, ‘sfz2-3.dat’ (second Elsasser variable),

which contains the data of interest. The ﬁrst line in each ﬁle contains the time tand

the number

qmax

, such that the largest moment calculated is

qmax

−1. The next

imax

numbers represent the ﬁrst moment structure function for the ﬁrst snapshot, here

imax

= 2ln(

nxgrid

)

ln 2 −2.(36)

The next

imax

numbers contain the second moment structure function, and so on until

qmax

−1. The following

imax

numbers then contain the data of the signed third order

structure function i.e. S3

long(l) = h[ux(x+l, y, z)−ux(x, y, z)]3i.

The following

imax

qmax

×2numbers are zero if

nr directions

= 1 (default), otherwise

they are the same data as above but for the structure functions calculated in the y and

z directions.

5.24 Particles 55

If the code has been run long enough as to calculate several snapshots, these snapshots

will now follow, being stored in the same way as the ﬁrst snapshot.

To visualize with

IDL

just type structure and you get the time-average of the ﬁrst order

longitudinal structure function (be sure that ‘pencil-runs/forced/idl/’ is in IDL_PATH).

See head of ‘pencil-runs/forced/idl/structure.pro’ for options to structure.

5.24 Particles

The PENCIL CODE has modules for tracer particles and for dust particles (see Sect. 6.14).

The particle modules are chosen by setting the value of the variable PARTICLES in

Makefile.local to either particles_dust or particles_tracers. For the former case

each particle has six degrees of freedom, three positions and three velocities. For the

latter it sufﬁces to have only three position variables as the velocity of the particles are

equal to the instantaneous ﬂuid velocity at that point. In addition one can choose to have

several additional internal degrees of freedoms for the particles. For example one can

temporally evolve the particles radius by setting PARTICLES_RADIUS to particles_radius

in Makefile.local.

All particle infrastructure is controlled and organized by the Particles_main module.

This module is automatically selected by Makefile.src if PARTICLES is different from

noparticles. Particle modules are compiled as a separate library. This way the main

part of the Pencil Code only needs to know about the particles_main.a library, but not

of the individual particle modules.

For a simulation with particles one must in addition deﬁne a few parameters in

cparam.local. Here is a sample of cparam.local for a parallel run with 2,000,000 parti-

cles:

integer, parameter :: ncpus=16, nprocy=4, nprocz=4, nprocx=1

integer, parameter :: nxgrid=128, nygrid=256, nzgrid=128

integer, parameter :: npar=2000000, mpar_loc=400000, npar_mig=1000

integer, parameter :: npar_species=2

The parameter npar is the number of particles in the simulation, mpar_loc is the number

of particles that is allowed on each processor and npar_mig is the number of particles

that are allowed to migrate from one processor to another in any time-step. For a non-

parallel run it is enough to specify npar. The number of particle species is set through

npar_species (assumed to be one if not set). The particle input parameters are given in

start.in and run.in. Here is a sample of the particle part of start.in for dust particles:

&particles_init_pars

initxxp=’gaussian-z’, initvvp=’random’

zp0=0.02, delta_vp0=0.01, eps_dtog=0.01, tausp=0.1

lparticlemesh_tsc=T

The initial positions and velocities of the dust particles are set in initxxp and initvvp.

The next four input parameters are further speciﬁcations of the initial condition. Inter-

action between the particles and the mesh, e.g. through drag force or self-gravity, require

a mapping of the particles on the mesh. The PENCIL CODE currently supports NGP

(Nearest Grid Point, default), CIC (Cloud in Cell, set lparticlemesh_cic=T) and TSC

56 THE PENCIL CODE

(Triangular Shaped Cloud, set lparticlemesh_tsc=T). See Youdin & Johansen (2007) for

details.

Here is a sample of the particle part of run.in (also for dust particles):

&particles_run_pars

ldragforce_gas_par=T

cdtp=0.2

The logical ldragforce_gas_par determines whether the dust particles inﬂuence the gas

with a drag force. cdtp tells the code how many friction times should be resolved in a

minimum time-step.

The sample run ‘samples/sedimentation/’ contains the latest setup for dust particles.

5.24.1 Particles in parallel

The particle variables (e.g. xiand vi) are kept in the arrays fp and dfp. For parallel

runs, particles must be able to move from processor to processor as they pass out of the

(x, y, z)-interval of the local processor. Since not all particles are present at the same

processor at the same time (hopefully), there is some memory optimization in making

fp not big enough to contain all the particles at once. This is achieved by setting the code

variable mpar_loc less than npar in cparam.local for parallel runs. When running with

millions of particles, this trick is necessary to keep the memory need of the code down.

The communication of migrating particles between the processors happens as follows

(see the subroutine redist_particles_procs in particles_sub.f90):

1. In the beginning of each time-step all processors check if any of their particles have

crossed the local (x, y, z)-interval. These particles are called migrating particles. A

run can have a maximum of npar_mig migrating particles in each time-step. The

value of npar_mig must be set in cparam.local. The number should (of course) be

slightly larger than the maximum number of migrating particles at any time-step

during the run. The diagnostic variable nmigmax can be used to output the maxi-

mum number of migrating particles at a given time-step. One can set lmigration_-

redo=T in &particles_run_pars to force the code to redo the migration step if more

than npar_mig want to migrate. This does slow the code down somewhat, but has

the beneﬁt that the code does not stop when more than npar_mig particles want to

migrate.

2. The index number of the receiving processor is then calculated. This requires some

assumption about the grid on other processors and will currently not work for

nonequidistant grids. Particles do not always pass to neighboring processors as the

global boundary conditions may send them to the other side of the global domains

(periodic or shear periodic boundary conditions).

3. The migrating particle information is copied to the end of fp, and the empty spot

left behind is ﬁlled up with the particle of the highest index number currently

present at the processor.

4. Once the number of migrating particles is known, this information is shared with

neighboring processors (including neighbors over periodic boundaries) so that they

5.24 Particles 57

all know how many particles they have to receive and from which processors.

5. The communication happens as directed MPI communication. That means that

processors 0 and 1 can share migrating particles at the same time as processors 2

and 3 do it. The communication happens from a chunk at the end of fp (migrating

particles) to a chunk that is present just after the particle of the highest index

number that is currently at the receiving processor. Thus the particles are put

directly at their ﬁnal destination, and the migrating particle information at the

source processor is simply overwritten by other migrating particles at the next

time-step.

6. Each processor keeps track of the number of particles that it is responsible for. This

number is stored in the variable npar_loc. It must never be larger than mpar_loc

(see above). When a particle leaves a processor, npar_loc is reduced by one, and

then increased by one at the processor that receives that particle. The maximum

number of particles at any processor is stored in the diagnostic variable nparmax. If

this value is not close to npar/ncpus, the particles have piled up in such a way that

computations are not evenly shared between the processors. One can then try to

change the parallelization architecture (nprocy and nprocz) to avoid this problem.

In simulations with many particles (comparable to or more than the number of grid

cells), it is crucial that particles are shared relatively evenly among the processors. One

can as a ﬁrst approach attempt to not parallelize directions with strong particle density

variations. However, this is often not enough, especially if particles clump locally.

Alternatively one can use Particle Block Domain Decomposition (PBDD, see Johansen

et al. 2011). The steps in Particle Block Domain Decomposition scheme are as follows:

1. The ﬁxed mesh points are domain-decomposed in the usual way (with

ncpus=nprocx×nprocy×nprocz).

2. Particles on each processor are counted in bricks of size nbx×nby×nbz (typically

nbx=nby=nbz=4).

3. Bricks are distributed among the processors so that each processor has approxi-

mately the same number of particles

4. Adopted bricks are referred to as blocks.

5. The Pencil Code uses a third order Runge-Kutta time-stepping scheme. In the be-

ginning of each sub-time-step particles are counted in blocks and the block counts

communicated to the bricks on the parent processors. The particle density assigned

to ghost cells is folded across the grid, and the ﬁnal particle density (deﬁned on the

bricks) is communicated back to the adopted blocks. This step is necessary because

the drag force time-step depends on the particle density, and each particle assigns

density not just to the nearest grid point, but also to the neighboring grid points.

6. In the beginning of each sub-time-step the gas density and gas velocity ﬁeld is

communicated from the main grid to the adopted particle blocks.

7. Drag forces are added to particles and back to the gas grid points in the adopted

blocks. This partition aims at load balancing the calculation of drag forces.

8. At the end of each sub-time-step the drag force contribution to the gas velocity ﬁeld

is communicated from the adopted blocks back to the main grid.

58 THE PENCIL CODE

Particle Block Domain Decomposition is activated by setting PARTICLES = particles_-

dust_blocks and PARTICLES_MAP = particles_map_blocks in Makefile.local. A sam-

ple of cparam.local for Particle Block Domain Decomposition can be found in

samples/sedimentation/blocks:

integer, parameter :: ncpus=4, nprocx=2, nprocy=2, nprocz=1

integer, parameter :: nxgrid=32, nygrid=32, nzgrid=32

integer, parameter :: npar=10000, mpar_loc=5000, npar_mig=100

integer, parameter :: npar_species=4

integer, parameter :: nbrickx=4, nbricky=4, nbrickz=4, nblockmax=32

The last line deﬁnes the number of bricks in the total domain – here we divide the grid

into 4×4×4bricks each of size 8×8×8grid points. The parameter nblockmax tells the

code the maximum number of blocks any processor may adopt. This should not be so low

that there is not room for all the bricks with particles, nor so high that the code runs out

of memory.

5.24.2 Large number of particles

When dealing with large number of particles, one needs to make sure that the number

of particles npar is less than the maximum integer that the compiler can handle with.

The maximum integer can be checked by the Fortran intrinsic function huge,

program huge_integers

print *, huge(0_4) ! for default Fortran integer (32 Bit)

print *, huge(0_8) ! for 64 Bit integer in Fortran

end program huge_integers

If the number of particles npar is larger than default maximum integer, one can promote

the maximum integer to 64 Bit by setting

integer(kind=8), parameter :: npar=4294967296

in the cparam.local ﬁle. This works because the data type of npar is only set here. It is

worth noting that one should not use the ﬂag

FFLAGS += -integer-size 64

to promote all the integers to 64 Bit. This will break the Fortran-C interface. One

should also make sure that npar_mig<=npar/ncpus. It is also beneﬁcial to set mpar_-

loc=2*npar/ncpus.

5.24.3 Random number generator

There are several methods to generate random number in the code. It is worth noting

that when simulating coagulation with the super-particle approach, one should use the

intrinsic random number generator of FORTRAN instead of the one implemented in the

code. When invoking random_number_wrapper, there will be back-reaction to the gas ﬂow.

This unexpected back-reaction can be tracked by inspecting the power spectra, which

exhibits the oscillation at the tail. To avoid this, one should set luser_random_number_-

wrapper=F under the module particles_coag_run_pars in run.in.

5.25 Non-cartesian coordinate systems 59

5.25 Non-cartesian coordinate systems

Spherical coordinates are invoked by adding the following line in the ﬁle ‘start.in’

&init_pars

coord_system=’spherical_coords’

One can also invoke cylindrical coordinates by saying cylindrical_coords instead. In

practice, the names (x, y, z)are still used, but they refer then to (r, θ, φ)or (r, φ, z)instead.

When working with curvilinear coordinates it is convenient to use differential operators

in the non-coordinate basis, so the derivatives are taken with respect to length, and

not in a mixed fashion with respect to length for ∂/∂r and with respect to angle for

∂/∂θ and ∂/∂φ. The components in the non-coordinate basis are denoted by hats, see,

e.g., [26], p. 213; see also Appendix B of [27]. For spherical polar coordinates the only

nonvanishing Christoffel symbols (or connection coefﬁcients) are

Γˆ

ˆrˆ

θ= Γˆ

ˆrˆ

φ=−Γˆrˆ

θˆ

θ=−Γˆrˆ

φˆ

φ= 1/r, (37)

Γˆ

φˆ

θˆ

φ=−Γˆ

θˆ

φˆ

φ= cot θ/r. (38)

The Christoffel symbols enter as correction terms for the various differential opera-

tors in addition to a term calculated straightforwardly in the non-coordinate basis. The

derivatives of some relevant Christoffel symbols are

Γˆ

ˆrˆ

θ,ˆ

θ= Γˆ

ˆrˆ

φ, ˆ

φ= Γˆ

φˆ

θˆ

φ, ˆ

φ= 0 (39)

Γˆ

ˆrˆ

θ,ˆr= Γˆ

ˆrˆ

φ,ˆr=−r−2(40)

Γˆ

φˆ

θˆ

φ,ˆ

θ=−r−2sin−2θ(41)

Further details are given in Appendix I.

60 THE PENCIL CODE

6 The equations

The equations solved by the PENCIL CODE are basically the standard compressible

MHD equations. However, the modular structure allows some variations of the MHD

equations, as well as switching off some of the equations or individual terms of the

equation (nomagnetic, noentropy, etc.).

In this section the equations are presented in their most complete form. It may be ex-

pected that the code can evolve most subsets or simpliﬁcations of these equations.

6.1 Continuity equation

In the code the continuity equation, ∂ρ/∂t +∇·ρu= 0, is written in terms of ln ρ,

D ln ρ

Dt=−∇·u.(42)

Here ρdenotes density, uthe ﬂuid velocity, tis time and D/Dt≡∂/∂t +u·∇is the

convective derivative.

6.2 Equation of motion

In the equation of motion, using a perfect gas, the pressure term, can be expressed as

−ρ−1∇p=−c2

s(∇s/cp+∇ln ρ), where the squared sound speed is given by

s=γp

ρ=c2

s0 exp γs/cp+ (γ−1) ln ρ

ρ0,(43)

and γ=cp/cvis the ratio of speciﬁc heats, or adiabatic index. Note that c2

sis proportional

to the temperature with c2

s= (γ−1)cpT.

The equation of motion is accordingly

Dt=−c2

s∇s

+ ln ρ−∇Φgrav +j×B

+ν∇2u+1

3∇∇ ·u+ 2S·∇ln ρ+ζ(∇∇ ·u) ; (44)

Here Φgrav is the gravity potential, jthe electric current density, Bthe magnetic ﬂux

density, νis kinematic viscosity, ζdescribes a bulk viscosity, and, in Cartesian coordi-

nates

Sij =1

2∂ui

∂xj

+∂uj

∂xi−2

3δij∇·u(45)

is the rate-of-shear tensor that is traceless, because it can be written as the generic rate-

of-strain tensor minus its trace. In curvilinear coordinates, we have to replace partial

differentiation by covariant differentiation (indicated by semicolons), so we write Sij =

2(ui;j+uj;i)−1

3δij∇·u.

The interpretation of the two viscosity terms varies greatly depending upon the Viscos-

ity module used, and indeed on the parameters given to the module. See §6.6.

For isothermal hydrodynamics, see §6.4 below.

6.3 Induction equation 61

6.3 Induction equation

∂A

∂t =u×B−ηµ0j.(46)

Here Ais the magnetic vector potential, B=∇ × Athe magnetic ﬂux density, η=

1/(µ0σ)is the magnetic diffusivity (σdenoting the electrical conductivity), and µ0the

magnetic vacuum permeability. This form of the induction equation corresponds to the

Weyl gauge

Φ = 0, where Φdenotes the scalar potential.

6.4 Entropy equation

The current thermodynamics module

entropy

formulates the thermal part of the physics

in terms of entropy s, rather than thermal energy e, which you may be more famil-

iar with. Thus the two fundamental thermodynamical variables are ln ρand s. The

reason for this choice of variables is that entropy is the natural physical variable for

(at least) convection processes: the sign of the entropy gradient determines convective

(in)stability, the Rayleigh number is proportional to the entropy gradient of the associ-

ated hydrostatic reference solution, etc. The equation solved is

ρT Ds

Dt=H − C +∇·(K∇T) + ηµ0j2+ 2ρνS⊗S+ζρ (∇·u)2.(47)

Here, Tis temperature, cpthe speciﬁc heat at constant pressure, Hand Care explicit

heating and cooling terms, Kis the radiative (thermal) conductivity, ζdescribes a bulk

viscosity, and Sis the rate-of-shear tensor that is traceless.

In the entropy module we solve for the speciﬁc entropy, s. The heat conduction term on

the right hand side can be written in the form

∇·(K∇T)

ρT (48)

=cpχh∇2ln T+∇ln T·∇(ln T+ ln χ+ ln ρ)i(49)

=cpχγ∇2s/cp+ (γ−1)∇2ln ρ

+cpχ[γ∇s/cp+ (γ−1)∇ln ρ]·[γ(∇s/cp+∇ln ρ) + ∇ln χ],(50)

where χ=K/(ρcp)is the thermal diffusivity. The latter equation shows that the diffu-

sivity for sis γχ, which is what we have used in Eq. (24).

In an alternative formulation for a constant K, the heat conduction term on the right

hand side can also be written in the form

∇·(K∇T)

ρT =K

ρh∇2ln T+ (∇ln T)2i(51)

which is the form used when constant Kis chosen.

Note that by setting γ= 1 and initially s= 0, one obtains an isothermal equation of

state (albeit at some unnecessary expense of memory). Similarly, by switching off the

evolution terms of entropy, one immediately gets polytropic behavior (if swas initially

constant) or generalized polytropic behavior (where sis not uniform, but ∂s/∂t = 0).

A better way to achieve isothermality is to use the

noentropy

module.

62 THE PENCIL CODE

6.4.1 Viscous heating

We can write the viscosity as the divergence of a tensor τij,j ,

ρ∂ui

∂t =... +τij,j ,(52)

where τij = 2νρSij is the stress tensor. The viscous power density Pis

P=uiτij,j (53)

=∂

∂xj

(uiτij)−ui,j τij (54)

The term under the divergence is the viscous energy ﬂux and the other term is the

kinetic energy loss due to heating. The heating term +ui,jτij is positive deﬁnite, because

τij is a symmetric tensor and the term only gives a contribution from the symmetric part

of ui,j, which is 1

2(ui,j +uj,i), so

ui,jτij =1

2νρ(ui,j +uj,i)(2Sij ).(55)

But, because Sij is traceless, we can add anything proportional to δij and, in particular,

we can write

ui,jτij =1

2(ui,j +uj,i)(2νρSij )(56)

2(ui,j +uj,i −1

3δij∇·u)(2νρSij)(57)

= 2νρS2,(58)

which is positive deﬁnite.

6.4.2 Alternative description

By setting pretend_lnTT=T in init_pars or run_pars (i.e. the general part of the name

list) the logarithmic temperature is used instead of the entropy. This has computational

advantages when heat conduction (proportional to K∇T) is important. Another alterna-

tive is to use another module, i.e. set ENTROPY=temperature_idealgas in ‘Makefile.local’.

When pretend_lnTT=T is set, the entropy equation

∂s

∂t =−u·∇s+1

ρT RHS (59)

is replaced by

∂ln T

∂t =−u·∇ln T+1

ρcvTRHS −(γ−1) ∇·u,(60)

where RHS is the right hand side of equation (47).

6.5 Transport equation for a passive scalar 63

6.5 Transport equation for a passive scalar

In conservative form, the equation for a passive scalar is

∂

∂t(ρc) + ∇·[ρcu−ρD∇c] = 0.(61)

Here cdenotes the concentration (per unit mass) of the passive scalar and Dits diffusion

constant (assumed constant). In the code this equation is solved in terms of ln c,

D ln c

Dt=D∇2ln c+ (∇ln ρ+∇ln c)·∇ln c.(62)

Using ln cinstead of chas the advantage that it enforces c > 0for all times. However,

the disadvantage is that one cannot have c= 0. For this reason we ended up using the

non-logarithmic version by invoking PSCALAR=pscalar_nolog.

6.6 Bulk viscosity

For a monatomic gas it can be shown that the bulk viscosity vanishes. We therefore don’t

use it in most of our runs. However, for supersonic ﬂows, or even otherwise, one might

want to add a shock viscosity which, in its simplest formulation, take the form of a bulk

viscosity.

6.6.1 Shock viscosity

Shock viscosity, as it is used here and also in the Stagger Code of ˚

Ake Nordlund, is

proportional to positive ﬂow convergence, maximum over ﬁve zones, and smoothed to

second order,

ζshock =cshock Dmax

5[(−∇·u)+]E(min(δx, δy, δz))2,(63)

where cshock is a constant deﬁning the strength of the shock viscosity. In the code this

dimensionless coefﬁcient is called nu_shock, and it is usually chosen to be around unity.

Assume that the shock viscosity only enters as a bulk viscosity, so the whole stress

tensor is then

τij = 2ρνSij +ρζshockδij ∇·u.(64)

Assume ν=const, but ζ6=const, so

ρ−1Fvisc =ν∇2u+1

3∇∇ ·u+ 2S·∇ln ρ+ζshock [∇∇ ·u+ (∇ln ρ+∇ln ζshock)∇·u].

(65)

and

ρ−1Γvisc = 2νS2+ζshock(∇·u)2.(66)

In the special case with periodic boundary conditions, we have 2hS2i=hω2i+4

3h(∇·u)2i.

6.7 Equation of state

In its present conﬁguration only hydrogen ionization is explicitly included. Other con-

stituents (currently He and H2) can have ﬁxed values. The pressure is proportional to

the total number of particles, i.e.

p= (nHI +nHII +nH2+ne+nHe +...)kBT. (67)

64 THE PENCIL CODE

It is convenient to normalize to the total number of H both in atomic and in molecular

hydrogen, nHtot ≡nH+ 2nH2, where nHI +nHII =nH, and deﬁne xe≡ne/nHtot,xHe ≡

nHe/nHtot, and xH2≡nH2/nHtot. Substituting nH=nHtot −2nH2, we have

p= (1 −xH2+xe+xHe +...)nHtotkBT. (68)

This can be written in the more familiar form

p=R

µρT, (69)

where R=kB/muand muis the atomic mass unit (which is for all practical purposes the

same as mHtot) and

µ=nH+ 2nH2+ne+ 4nHe

nH+nH2+ne+nHe

=1 + 4xHe

1−xH2+xe+xHe

(70)

is the mean molecular weight (which is here dimensionless; see Kippenhahn & Weigert

1990, p. 102). The factor 4 is really to be substituted for 3.97153. Some of the familiar

relations take still the usual form, in particular e=cvTand h=cpTwith cv=3

2R/µ and

cp=5

2R/µ.

The number ratio, xHe, is more commonly expressed as the mass ratio, Y=

mHenHe/(mHnHtot +mHenHe), or Y= 4xHe/(1 + 4xHe), or 4xHe = (1/Y −1)−1. For exam-

ple, Y= 0.27 corresponds to xHe = 0.092 and Y= 0.25 to xHe = 0.083. Note also that for

100% H2abundance, xH2= 1/2.

In the following, the ionization fraction is given as y=ne/nH, which can be different

from xeif there is H2. Substituting for nHin terms of nHtot yields y=xe/(1 −2xH2).

6.8 Ionization

This part of the code can be invoked by setting EOS=eos_ionization (or EOS=eos_-

temperature_ionization) in the ‘Makefile.local’ ﬁle. The equation of state described

below works for variable ionization, and the entropy offset is different from that used in

Eq. (43), which is now no longer valid. As a replacement, one can use EOS=eos_fixed_-

ionization, where the degree of ionization can be given by hand. Here the normalization

of the entropy is the same as for EOS=eos_ionization. This case is described in more de-

tail below.12

We treat the gas as being composed of partially ionized hydrogen and neutral helium.

These are four different particle species, each of which regarded as a perfect gas.

The ionization fraction y, which gives the ratio of ionized hydrogen to the total amount

of hydrogen nH, is obtained from the Saha equation which, in this case, may be written

1−y=1

nHmekBT

2π~23/2

exp −χH

kBT.(71)

The temperature Tcannot be obtained directly from the PENCIL CODE’s independent

variables ln ρand s, but is itself dependent on y. Hence, the calculation of yessentially

becomes a root ﬁnding problem.

12We omit here the contribution of H2.

6.8 Ionization 65

The entropy of a perfect gas consisting of particles of type iis known from the Sackur-

Tetrode equation

Si=kBNi ln "1

ntot mikBT

2π~23/2#+5

2!.(72)

Here Niis the number of particles of a single species and ntot is the total number density

of all particle species.

In addition to the individual entropies we also have to take the entropy of mixing,

Smix =−NtotkBPipiln pi, into account. Summing up everything, we can get a closed ex-

pression for the speciﬁc entropy sin terms of y,ln ρand T, which may be solved for

Figure 5: Dependence of temperature on entropy for different values of the density.

For given ln ρand swe are then able to calculate the ionization fraction yby ﬁnding the

root of

f(y) = ln "1−y

nHmekBT(y)

2π~23/2#−χH

kBT(y).(73)

In the ionized case, several thermodynamic quantities of the gas become dependent on

the ionization fraction ysuch as its pressure, P= (1+ y+xHe)nHkBT, and its internal

energy, E=3

2(1 + y+xHe)nHkBT+yχH, where xHe gives the ratio of neutral helium to the

total amount of hydrogen. The dependence of temperature on entropy is shown in Fig. 5

for different values of the density.

For further details regarding the procedure of solving for the entropy see Sect. H.6 in

the appendix.

6.8.1 Ambipolar diffusion

Another way of dealing with ionization in the PENCIL CODE is through use of the neu-

trals module. That module solves the coupled equations of neutral and ionized gas, in a

two-ﬂuid model

66 THE PENCIL CODE

∂ρi

∂t =−∇·(ρiui) + G(74)

∂ρn

∂t =−∇·(ρnun)− G (75)

∂(ρiui)

∂t =−∇·(ρiui:ui)−∇pi+pe+B2

2µ0+F(76)

∂(ρnun)

∂t =−∇·(ρnun:un)−∇pn− F (77)

∂A

∂t =ui×B(78)

where the subscripts nand iare for neutral and ionized, respectively. The terms Gand

F, through which the two ﬂuids exchange mass and momentum, are given by

G=ζρn−αρ2

i(79)

F=ζρnun−αρ2

iui+γρiρn(un−ui).(80)

In the above equations, ζis the ionization coefﬁcient, αis the recombination coefﬁcient,

and γthe collisional drag strength. By the time of writing (spring 2009), these three

quantities are supposed constant. The electron pressure peis also assumed equal to the

ion pressure. Only isothermal neutrals are supported so far.

In the code, Eq. (74) and Eq. (76) are solved in ‘density.f90’ and ‘hydro.f90’ respectively.

Equation 75 is solved in ‘neutraldensity.f90’ and Eq. (77) in ‘neutralvelocity.f90’. The

sample ‘1d-test/ambipolar-diffusion’ has the current setup for a two-ﬂuid simulation

with ions and neutrals.

6.9 Radiative transfer

Here we only state the basic equations. A full description about the implementation is

given in Sect. H.7 and in the original paper by Heinemann et al. (2006).

The basic equation for radiative transfer is

dτ =−I+S , (81)

where

τ≡

χ(s′)ds′(82)

is the optical depth (sis the geometrical coordinate along the ray).

Note that radiative transfer is called also in ‘start.csh’, and again each time a snapshot

is being written, provided the output of auxiliary variables is being requested lwrite_-

aux=T. (Also, of course, the pencil check runs radiative transfer 7 times, unless you put

pencil_check_small=F.)

6.10 Self-gravity 67

6.10 Self-gravity

The PENCIL CODE can consider the self-gravity of the ﬂuid in the simulation box by

adding the term

∂u

∂t =...−∇φself (83)

to the equation of motion. The self-potential φself (or just φfor simplicity) satisﬁes Pois-

son’s equation

∇2φ= 4πGρ . (84)

The solution for a single Fourier component at scale kis

φk=−4πGρk

k2.(85)

Here we have assumed periodic boundary conditions. The potential is obtained by

Fourier-transforming the density, then ﬁnding the corresponding potential at that scale,

and ﬁnally Fourier-transforming back to real space.

The x-direction in the shearing sheet is not strictly periodic, but is rather shear periodic

with two connected points at the inner and outer boundary separated by the distance

∆y(t) = mod[(3/2)Ω0Lxt, Ly]in the y-direction. We follow here the method from [17]

to allow for shear-periodic boundaries in the Fourier method for self-gravity. First we

take the Fourier transform along the periodic y-direction. We then shift the entire y-

direction by the amount δy(x) = ∆y(t)x/Lxto make the x-direction periodic. Then we

proceed with Fourier transforms along xand then z. After solving the Poisson equation

in Fourier space, we transform back to real space in the opposite order. We differ here

from the method by [17] in that we shift in Fourier space rather than in real space13.

The Fourier interpolation formula has the advantage over polynomial interpolation in

that it is continuous and smooth in all its derivatives.

6.11 Incompressible and anelastic equations

This part has not yet been documented and is still under development.

6.12 Dust equations

The code treats gas and dust as two separate ﬂuids14. The dust and the gas interact

through a drag force. This force can most generally be written as an additional term to

the equation of motion as

Dud

Dt=...−1

τs

(ud−u).(86)

Here τsis the so-called stopping time of the considered dust species. This measures the

coupling strength between dust and gas. In the Epstein drag regime

τs=adρs

csρ,(87)

13We were kindly made aware of the possibility of interpolating in Fourier space by C. McNally on his

website.

14See master’s thesis of A. Johansen (can be downloaded from

http://www.mpia.de/homes/johansen/research_en.php)

68 THE PENCIL CODE

where adis the radius of the dust grain and ρsis the solid density of the dust grain.

Two other important effects work on the dust. The ﬁrst is coagulation controlled

by the discrete coagulation equation

dnk

dt=1

i+j=k

Aijninj−nk

∞

i=1

Aikni.(88)

In the code Ndiscrete dust species are considered. Also the bins are logarithmically

spaced in order to give better mass resolution. It is also possible to keep track of both

number density and mass density of each bin, corresponding to having a variable grain

mass in each bin.

Dust condensation is controlled by the equation

dt=1

τcond

Nd−1

d.(89)

Here Nis the number of monomers in the dust grain (such as water molecules) and dis

the physical dimension of the dust grain. The condensation time τcond is calculated from

τcond

=A1vthαnmon 1−1

Smon ,(90)

where A1is the surface area of a monomer, αis the condensation efﬁciency, nmon is the

number density of monomers in the gas and Smon is the saturation level of the monomer

given by

Smon =Pmon

Psat

.(91)

Here Psat is the saturated vapor pressure of the monomer. Currently only water ice has

been implemented in the code.

All dust species fulﬁll the continuity equation

∂ρd

∂t +∇·(ρdud) = 0.(92)

6.13 Cosmic ray pressure in diffusion approximation

Cosmic rays are treated in the diffusion approximation. The equation of state is pc=

(γc)ecwhere the value of γcis usually somewhere between 14/9and 4/3. In the momen-

tum equation (44) the cosmic ray pressure force, −ρ−1∇pcis added on the right hand

side, and ecsatisﬁes the evolution equation

∂ec

∂t +∇·(ecu) + pc∇·u=∂i(Kij ∂jec) + Qc,(93)

where Qcis a source term and

Kij =K⊥δij + (Kk−K⊥)ˆ

Biˆ

Bj(94)

is an anisotropic diffusivity tensor.

6.14 Particles 69

In the non-conservative formulation of this code it is advantageous to expand the diffu-

sion term using the product rule, i.e.

∂i(Kij∂jec) = −Uc·∇ec+Kij ∂i∂jec.(95)

where Uci=−∂Kij /∂xjacts like an extra velocity trying to straighten magnetic ﬁeld

lines. We can write this term also as Uc=−(Kk−K⊥)∇·(ˆ

Bˆ

B), where the last term

is a divergence of the dyadic product of unit vectors.15 However, near magnetic nulls,

this term can becomes inﬁnite. In order to avoid this problem we are forced to limit

∇·(ˆ

Bˆ

B), and hence |Uc|, to the maximum possible value that can be resolved at a given

resolution.

A physically appealing way of limiting the maximum propagation speed is to restore

an explicit time dependence in the equation for the cosmic ray ﬂux, and to replace the

diffusion term in Eq. (93) by a divergence of a ﬂux that in turn obeys the equation

∂Fci

∂t =−˜

Kij∇jec−Fci

τ(non-Fickian diffusion),(96)

where Kij =τ˜

Kij would be the original diffusion tensor of Eq. (94), if the time derivative

were negligible. Further details are described in Snodin et al. (2006).

6.14 Particles

Particles are entities that each have a space coordinate and a velocity vector, where

a ﬂuid only has a velocity vector ﬁeld (the continuity equation of a ﬂuid in some way

corresponds to the space coordinate of particles). In the code particles are present either

as tracer particles or as dust particles

6.14.1 Tracer particles

Tracer particles always have the local velocity of the gas. The dynamical equations are

thus ∂xi

∂t =u,(97)

where the index iruns over all particles. Here uis the gas velocity at the position of

the particle. One can choose between a ﬁrst order (default) and a second order spline

interpolation scheme (set lquadratic_interpolation=T in &particles_init_pars) to cal-

culate the gas velocity at the position of a tracer particle.

The sample run ‘samples/dust-vortex’ contains the latest setup for tracer particles.

6.14.2 Dust particles

Dust particles are allowed to have a velocity that is not similar to the gas,

dxi

dt=vi.(98)

15In practice, we calculate ∂j(ˆ

Biˆ

Bj) = (δij −2ˆ

Biˆ

Bk)ˆ

BjBk,j /|B|, where derivatives of Bare calculated as

Bi,j =ǫiklAl,jk.

70 THE PENCIL CODE

The particle velocity follows an equation of motion similar to a ﬂuid, only there is no

advection term. Dust particles also experience a drag force from the gas (proportional to

the velocity difference between a particle and the gas).

dvi

dt=...−1

τs

(vi−u).(99)

Here τsis the stopping time of the dust particle. The interpolation of the gas velocity to

the position of a particle is done using one of three possible particle-mesh schemes,

•NGP (Nearest Grid Point, default)

The gas velocity at the nearest grid point is used.

•CIC (Cloud in Cell, set lparticlemesh_cic=T)

A ﬁrst order interpolation is used to obtain the gas velocity ﬁeld at the position of

a particle. Affects 8 grid points.

•TSC (Triangular Shaped Cloud, set lparticlemesh_tsc=T)

A second order spline interpolation is used to obtain the gas velocity ﬁeld at the

position of a particle. Affects 27 grid points.

The particle description is the proper description of dust grains, since they do not feel

any pressure forces (too low number density). Thus there is no guarantee that the grains

present within a given volume will be equilibrated with each other, although drag force

may work for small grains to achieve that. Larger grains (meter-sized in protoplanetary

discs) must be treated as individual particles.

To conserve momentum the dust particles must affect the gas with a friction force as

well. The strength of this force depends on the dust-to-gas ratio ǫd, and it can be safely

ignored when there is much more gas than there is dust, e.g. when ǫd= 0.01. The friction

force on the gas appears in the equation of motion as

∂u

∂t =...−ρ(i)

ρ∂v(i)

∂t drag

(100)

Here ρ(i)

pis the dust density that particle irepresents. This can be set through the pa-

rameter eps_todt in &particle_init_pars. The drag force is assigned from the particles

onto the mesh using either NGP, CIC or TSC assignment. The same scheme is used both

for interpolation and for assignment to avoid any risk of a particle accelerating itself (see

Hockney & Eastwood 1981).

6.15 N-body solver

The N-body code takes advantage of the existing Particles module, which was coded with

the initial intent of treating solid particles whose radius a•is comparable to the mean

free path λof the gas, for which a ﬂuid description is not valid. A N-body implementation

based on that module only needed to include mass as extra state for the particles, solve

for the N2gravitational pair interactions and distinguish between the N-body and the

small bodies that are mapped into the grid as a ρpdensity ﬁeld.

The particles of the N-body ensemble evolve due to their mutual gravity and by inter-

acting with the gas and the swarm of small bodies. The equation of motion for particle i

6.16 Test-ﬁeld equations 71

dvpi

dt =Fgi−

j6=i

GMj

Rij (101)

where Rij =|rpi−rpj|is the distance between particles iand j, and ˆ

Rij is the unit vector

pointing from particle jto particle i. The ﬁrst term of the R.H.S. is the combined gravity

of the gas and of the dust particles onto the particle i, solved via

Fgi=−GZV

[ρg(r) + ρp(r)]Ri

(R2

i+b2

i)3/2dV, (102)

where the integration is carried out over the whole disk. The smoothing distance biis

taken to be as small as possible (a few grid cells). For few particles (<10), calculating the

integral for every particle is practical. For larger ensembles one would prefer to solve

the Poisson equation to calculate their combined gravitational potential.

The evolution of the particles is done with the same third-order Runge-Kutta time-

stepping routine used for the gas. The particles deﬁne the timestep also by the Courant

condition that they should not move more than one cell at a time. For pure particle runs,

where the grid is absent, one can adopt a ﬁxed time-step tp≪2πΩ−1

fp where Ωfp is the

angular frequency of the fastest particle.

By now (spring 2009), no inertial accelerations are included in the N-body module, so

only the inertial frame - with origin at the barycenter of the N-body ensemble - is avail-

able. For a simulation of the circular restricted three-body problem with mass ratio

q=10−3, the Jacobi constant of a test particle initially placed at position (x, y)=(2,0) was

found to be conserved up to one part in 105within the time span of 100 orbits.

We stress that the level of conservation is poor when compared to integrators designed

to speciﬁcally deal with long-term N-body problems. These integrators are usually sym-

plectic, unlike the Runge-Kutta scheme of the PENCIL CODE. As such, PENCIL should

not be used to deal with evolution over millions of years. But for the time-span typical

of astrophysical hydrodynamical simulations, this degree of conservation of the Jacobi

constant can be deemed acceptable.

As an extension of the particle’s module, the N-body is fully compatible with the par-

allel optimization of PENCIL, which further speeds up the calculations. Parallelization,

however, is not yet possible for pure particle runs, since it relies on splitting the grid

between the processors. At the time of writing (spring 2009), the N-body code does not

allow the particles to have a time-evolving mass.

6.16 Test-ﬁeld equations

The test-ﬁeld method is used to calculate turbulent transport coefﬁcients for magneto-

hydrodynamics. This is a rapidly evolving ﬁeld and we refer the interested reader to

recent papers in this ﬁeld, e.g. by Sur et al. (2008) or Brandenburg et al. (2008). For

technical details; see also Sect. F.3.

6.17 Gravitational wave equations

The expansion of the universe with time is described by the scale factor a(τ), where

τis the physical time. Using conformal time, t(τ) = Rτ

0dτ′/a(τ′), and dependent vari-

72 THE PENCIL CODE

Table 5: Scale factor and conformal Hubble parameter for different values of n.

n a HH

0 1 0 0

1/2η/2 1/η 1/η

2/3η2/3 2/η 6/η2

ables that are appropriately scaled with powers of a, the hydromagnetic equations can

be expressed completely without scale factor [9, 16]. This is not true, however, for the

gravitational wave (GW) equations, where a dependence on aremains [16]. The time

dependence of acan be modeled as a power law, a∝τn, where n= 1/2applies to the

radiation-dominated era; see Table 5 the general relationship. To compare with cases

where the expansion is ignored, we put n= 0.

In the transverse traceless (TT) gauge, the six components of the spatial part of the

symmetric tensor characterizing the linearized evolution of the metric perturbations

hij, reduce to two components which, in the linear polarization basis, are the +and

×polarizations. However, the projection onto that basis is computationally intensive,

because it requires nonlocal operations involving Fourier transformations. It is therefore

advantageous to evolve instead the perturbation of the metric tensor, hij , in an arbitrary

gauge, compute then hTT

ij in the TT gauge, and perform then the decomposition into the

linear polarization basis whenever we compute diagnostic quantities such as averages

or spectra. Thus, we solve the linearized GW equation in the form [16]

∂2hij

∂t2=−2H∂hij

∂t +c2∇2hij +16πG

a2c2Tij (103)

for the six components 1≤i≤j≤3, where tis comoving time, ais the scale factor,

H= ˙a/a is the comoving Hubble parameter, Tij is the source term, cis the speed of light,

and Gis Newton’s constant. For n= 0, when the cosmic expansion is ignored, we have

a= 1 and H= 0. We use the PENCIL CODE; for the numerical treatment of Eq. (103)

and equations (105)–(107). For most of the simulations, we use 11523meshpoints on 1152

cores of a Cray XC40 system with 2.3 GHz processors.

The source term is chosen to be the traceless part of the stress tensor,

Tij(x, t) = ρuiuj−BiBj−1

3δij(ρu2−B2).(104)

The removal of the trace is in principle not necessary, but it helps preventing a contin-

uous build-up of a large trace, which would be numerically disadvantageous. We have

ignored here the viscous stress, which is usually small.

We compute Tij by solving the energy, momentum, and induction equations for an ultra-

relativistic gas in the form [9, 11]

∂ln ρ

∂t =−4

3(∇·u+u·∇ln ρ) + 1

ρu·(J×B) + ηJ2,(105)

Dt =u

3(∇·u+u·∇ln ρ)−u

ρu·(J×B) + ηJ2

−1

4∇ln ρ+3

4ρJ×B+2

ρ∇·(ρνS) + f,(106)