Manual

User Manual: Pdf

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

The Pencil Code:
A High-Order MPI code for MHD Turbulence
User’s and Reference Manual
February 20, 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.
ii
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
Copyright c
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 modified 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.
iv
License agreement and giving credit
The content of all files 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
find 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 field 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 field-
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
v
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 first 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 crossflow
as function of Stokes and Reynolds numbers,J. Fluid Mech. 661, 239-261
Rheinhardt, M., & Brandenburg, A.: 2010, “Test-field method for mean-field coefficients
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.
vi
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 finite-difference or spectral codes. Having been approached by peo-
ple interested in using our code, we decided to make it as flexible 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 find 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
flows, which is why we use high-order first and second derivatives. To achieve good par-
allelization, we use explicit (as opposed to compact) finite differences. Typical scientific
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 fields are implemented in terms of the magnetic vector potential to ensure
that the field 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 flux conserving schemes is
magnetic
helicity
.
There are currently no plans to implement adaptive mesh refinement into the code,
which would cause major technical complications. Given that turbulence is generically
space-filling, local refinement to smaller scales would often not be very useful anyway.
On the other hand, in some geometries turbulence may well be confined to certain re-
gions in space, so one could indeed gain by solving the outer regions with fewer points.
In order to be cache-efficient, 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 defined 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 flexibility is achieved by encapsulating individual physical processes
and variables in individual
modules
, which can be switched on or off in the file
Makefile.local’ in the local ‘src’ directory. This approach avoids the use of difficult-
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 specifies
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’ define an F90 module named Magnetic — this is the
basis of the switching mechanism we are using.
Input parameters (which are set in the files ‘start.in’, ‘run.in’) can be changed without
recompilation. Furthermore, one can change the list of variables for monitoring (diag-
nostic) output on the fly, 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
files 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
C
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
,
DX
(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 efficiently 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 first 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
flows 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 flux-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 field 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-fluid 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 Scientific Computing for
granting time on Horseshoe, which is a 512+140 processor Beowulf cluster in Odense
(Horseshoe).
ix
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 files . . . . . . . . . . . . . . . . . . . . . 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 Configuring the code to compile and run on your computer . . . . . . . . . 17
5.1.1 Locating the configuration file . . . . . . . . . . . . . . . . . . . . . . 17
5.1.2 Structure of configuration files . . . . . . . . . . . . . . . . . . . . . 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 files . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
5.6.1 Snapshot files . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
5.7 Video files and slices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
x
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 files . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
5.11 RERUN and NEWDIR files . . . . . . . . . . . . . . . . . . . . . . . . . . . 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
xi
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-field equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
7 Troubleshooting / Frequently Asked Questions 72
7.1 Download and setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
7.1.1 Download forbidden . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
7.1.2 Shell gives error message when sourcing ‘sourceme.X . . . . . . . . 72
7.2 Compilation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
7.2.1 Problems compiling syscalls . . . . . . . . . . . . . . . . . . . . . . . 73
7.2.2 Unable to open include file: chemistry.h . . . . . . . . . . . . . . . . 73
7.2.3 Compiling with
ifc
under Linux . . . . . . . . . . . . . . . . . . . . . 73
7.2.4 Segmentation fault with
ifort
8.0 under Linux . . . . . . . . . . . . 74
7.2.5 The underscore problem: linking with
MPI
. . . . . . . . . . . . . . 74
7.2.6 Compilation stops with the cryptic error message: . . . . . . . . . . 74
7.2.7 The code doesn’t compile, . . . . . . . . . . . . . . . . . . . . . . . . . 75
7.2.8 Some samples don’t even compile, . . . . . . . . . . . . . . . . . . . 75
7.2.9 Internal compiler error with Compaq/Dec F90 . . . . . . . . . . . . 76
7.2.10 Assertion failure under SunOS . . . . . . . . . . . . . . . . . . . . . 76
7.2.11 After some dirty tricks I got pencil code to compile with MPI, ... . . 77
7.2.12 Error: Symbol ’mpi comm world’ at (1) has no IMPLICIT type . . . 77
7.2.13 Error: Can’t open included file ’mpif.h’ . . . . . . . . . . . . . . . . . 78
7.3 Pencil check . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
7.3.1 The pencil check complains for no reason. . . . . . . . . . . . . . . . 78
7.3.2 The pencil check reports MISSING PENCILS and quits . . . . . . 78
7.3.3 The pencil check reports unnecessary pencils . . . . . . . . . . . . . 78
7.3.4 The pencil check reports that most or all pencils are missing . . . . 78
7.3.5 Running the pencil check triggers mathematical errors in the code 79
7.3.6 The pencil check still complains . . . . . . . . . . . . . . . . . . . . . 79
7.3.7 The pencil check is annoying so I turned it off . . . . . . . . . . . . 79
7.4 Running . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
7.4.1 Periodic boundary conditions in ‘start.x . . . . . . . . . . . . . . . 79
7.4.2 csh problem? . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
7.4.3 ‘run.csh doesn’t work: . . . . . . . . . . . . . . . . . . . . . . . . . . 80
7.4.4 Code crashes after restarting . . . . . . . . . . . . . . . . . . . . . . 80
7.4.5 auto-test gone mad...? . . . . . . . . . . . . . . . . . . . . . . . . . . 80
7.4.6 Can I restart with a different number of cpus? . . . . . . . . . . . . 81
7.4.7 Can I restart with a different number of cpus? . . . . . . . . . . . . 81
7.4.8 fft xyz parallel 3D: nygrid needs to be an integer multiple... . . . . 81
7.4.9 Unit-agnostic calculations? . . . . . . . . . . . . . . . . . . . . . . . 82
7.5 Visualization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
7.5.1 ‘start.pro doesn’t work: . . . . . . . . . . . . . . . . . . . . . . . . . 83
xii
7.5.2 ‘start.pro doesn’t work: . . . . . . . . . . . . . . . . . . . . . . . . . 83
7.5.3 Something about tag name undefined: . . . . . . . . . . . . . . . . . 83
7.5.4 Something INC in start.pro . . . . . . . . . . . . . . . . . . . . . . . 83
7.5.5 nl2idl problem when reading param2.nml . . . . . . . . . . . . . . . 84
7.5.6 Spurious dots in the time series file . . . . . . . . . . . . . . . . . . 84
7.5.7 Problems with pc_varcontent.pro . . . . . . . . . . . . . . . . . . . 84
7.6 General questions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85
7.6.1 “Installation” procedure . . . . . . . . . . . . . . . . . . . . . . . . . 85
7.6.2 Small numbers in the code . . . . . . . . . . . . . . . . . . . . . . . . 85
7.6.3 Why do we need a /lphysics/ namelist in the first place? . . . . . . 86
7.6.4 Can I run the code on a Mac? . . . . . . . . . . . . . . . . . . . . . . 87
7.6.5 Pencil Code discussion forum . . . . . . . . . . . . . . . . . . . . . . 87
7.6.6 The manual . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
II Programming the PENCIL CODE 89
8 Understanding the code 93
8.1 Example: how is the continuity equation being solved? . . . . . . . . . . . 93
9 Adapting the code 95
9.1 The PENCIL CODE coding standard . . . . . . . . . . . . . . . . . . . . . . . 95
9.2 Adding new output diagnostics . . . . . . . . . . . . . . . . . . . . . . . . . 96
9.3 The f-array . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98
9.4 The df-array . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98
9.5 The fp-array . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99
9.6 The pencil case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99
9.6.1 Pencil check . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100
9.6.2 Adding new pencils . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101
9.7 Adding new physics: the Special module . . . . . . . . . . . . . . . . . . . . 101
9.8 Adding switchable modules . . . . . . . . . . . . . . . . . . . . . . . . . . . 102
9.9 Adding your initial conditions: the InitialCondition module . . . . . . . . . 102
10 Testing the code 104
10.1 How to set up periodic tests . . . . . . . . . . . . . . . . . . . . . . . . . . . 104
11 Useful internals 106
11.1 Global variables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106
11.2 Subroutines and functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106
III Appendix 109
A Timings 109
A.1 Test case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116
A.2 Running the code . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117
A.3 Triolith . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117
A.4 Lindgren . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117
B Coding standard 120
B.1 File naming conventions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120
xiii
B.2 Fortran Code . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120
B.2.1 Indenting and whitespace . . . . . . . . . . . . . . . . . . . . . . . . 120
B.2.2 Comments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121
B.2.3 Module names . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122
B.2.4 Variable names . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122
B.2.5 Emacs settings . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123
B.3 Other best practices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124
B.4 General changes to the code . . . . . . . . . . . . . . . . . . . . . . . . . . . 124
C Some specific initial conditions 125
C.1 Random velocity or magnetic fields . . . . . . . . . . . . . . . . . . . . . . . 125
C.2 Turbulent initial with given spectrum . . . . . . . . . . . . . . . . . . . . . 125
C.3 Beltrami fields . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126
C.4 Magnetic flux rings: initaa=’fluxrings’ . . . . . . . . . . . . . . . . . . . 126
C.5 Vertical stratification . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127
C.5.1 Isothermal atmosphere . . . . . . . . . . . . . . . . . . . . . . . . . . 127
C.5.2 Polytropic atmosphere . . . . . . . . . . . . . . . . . . . . . . . . . . 128
C.5.3 Changing the stratification . . . . . . . . . . . . . . . . . . . . . . . 129
C.5.4 The Rayleigh number . . . . . . . . . . . . . . . . . . . . . . . . . . . 129
C.5.5 Entropy boundary condition . . . . . . . . . . . . . . . . . . . . . . . 130
C.5.6 Temperature boundary condition at the top . . . . . . . . . . . . . . 130
C.6 Potential-field boundary condition . . . . . . . . . . . . . . . . . . . . . . . 130
C.7 Planet solution in the shearing box . . . . . . . . . . . . . . . . . . . . . . . 132
D Some specific boundary conditions 133
D.1 Perfect-conductor boundary condition . . . . . . . . . . . . . . . . . . . . . 133
D.2 Stress-free boundary condition . . . . . . . . . . . . . . . . . . . . . . . . . 133
D.3 Normal-field-radial boundary condition . . . . . . . . . . . . . . . . . . . . 134
E High-frequency filters 135
E.1 Conservative hyperdissipation . . . . . . . . . . . . . . . . . . . . . . . . . . 135
E.2 Hyperviscosity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137
E.2.1 Conservative case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137
E.2.2 Non-conservative cases . . . . . . . . . . . . . . . . . . . . . . . . . . 138
E.2.3 Choosing the coefficient . . . . . . . . . . . . . . . . . . . . . . . . . 139
E.2.4 Turbulence with hyperviscosity . . . . . . . . . . . . . . . . . . . . . 139
E.3 Anisotropic hyperdissipation . . . . . . . . . . . . . . . . . . . . . . . . . . . 140
E.4 Hyperviscosity in Burgers shock . . . . . . . . . . . . . . . . . . . . . . . . 140
F Special techniques 142
F.1 After changing
REAL PRECISION
. . . . . . . . . . . . . . . . . . . . . . . 142
F.2 Remeshing (regridding) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142
F.3 Restarting from a run with less physics . . . . . . . . . . . . . . . . . . . . 143
G Runs and reference data 145
G.1 Shock tests . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 145
G.1.1 Sod shock tube problem . . . . . . . . . . . . . . . . . . . . . . . . . 145
G.1.2 Temperature jump . . . . . . . . . . . . . . . . . . . . . . . . . . . . 145
G.2 Random forcing function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 145
G.3 Three-layered convection model . . . . . . . . . . . . . . . . . . . . . . . . . 146
G.4 Magnetic helicity in the shearing sheet . . . . . . . . . . . . . . . . . . . . 147
xiv
H Numerical methods 151
H.1 Sixth-order spatial derivatives . . . . . . . . . . . . . . . . . . . . . . . . . 151
H.2 Upwind derivatives to avoid ‘wiggles’ . . . . . . . . . . . . . . . . . . . . . . 152
H.3 The bidiagonal scheme for cross-derivatives . . . . . . . . . . . . . . . . . . 153
H.4 The 2N-scheme for time-stepping . . . . . . . . . . . . . . . . . . . . . . . . 154
H.5 Diffusive error from the time-stepping . . . . . . . . . . . . . . . . . . . . . 155
H.6 Ionization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 156
H.7 Radiative transfer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 157
H.7.1 Solving the radiative transfer equation . . . . . . . . . . . . . . . . 157
H.7.2 Angular integration . . . . . . . . . . . . . . . . . . . . . . . . . . . . 158
I Switchable modules 161
J Startup and run-time parameters 161
J.1 Startup parameters for start.in . . . . . . . . . . . . . . . . . . . . . . . . 161
J.2 Runtime parameters for run.in . . . . . . . . . . . . . . . . . . . . . . . . 169
J.3 Parameters for ‘print.in . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 176
J.4 Parameters for ‘video.in . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 212
J.5 Parameters for ‘phiaver.in . . . . . . . . . . . . . . . . . . . . . . . . . . . 213
J.6 Parameters for ‘xyaver.in . . . . . . . . . . . . . . . . . . . . . . . . . . . . 214
J.7 Parameters for ‘xzaver.in . . . . . . . . . . . . . . . . . . . . . . . . . . . . 221
J.8 Parameters for ‘yzaver.in . . . . . . . . . . . . . . . . . . . . . . . . . . . . 222
J.9 Parameters for ‘yaver.in . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 224
J.10 Parameters for ‘zaver.in . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 226
J.11 Boundary conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 229
J.11.1 Boundary condition
bcx
. . . . . . . . . . . . . . . . . . . . . . . . . 229
J.11.2 Boundary condition
bcy
. . . . . . . . . . . . . . . . . . . . . . . . . 232
J.11.3 Boundary condition
bcz
. . . . . . . . . . . . . . . . . . . . . . . . . 234
J.12 Initial condition parameter dependence . . . . . . . . . . . . . . . . . . . . 238
IV Indexes 243
xv
xvi
1
Part I
Using the PENCIL CODE
1 System requirements
To use the code, you will need the following:
1. Absolutely needed:
F95
compiler
C
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)
DX
alias
OpenDX
or
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
or
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
or
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 figure 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 first 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, conflicts 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 conflicts 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 conflicts 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 sufficient comments in the code and please follow
our standard coding conventions explained in Section 9.1. There is also a script to check
and fix 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 files 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 first, 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 benefit of all other developers, any errors should
be corrected within 1-2 hours. This is the reason why the code comes with a file
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 find 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
,
DX
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
or
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 file (‘~/.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
files 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
file 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 first 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 first undefine the environment
variable PENCIL_HOME:
csh> unsetenv PENCIL_HOME
csh> source ./sourceme.csh
or
sh> unset PENCIL_HOME
sh> . ./sourceme.sh
3.1.2 Linking scripts and source files
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 files (plus a few more files) 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 flags, 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 predefined 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 flags correctly, compilation should complete successfully.
3.1.5 Choosing a data directory
The code will by default write data like snapshot files 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 file called ‘datadir.in’ in the root of
the run directory. This file 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 profiles 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 profile.
Now we run the code:
unix> run.csh
This executes ‘src/run.x’ and carries out
nt
time steps, where
nt
and other run-time
parameters are specified 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
it
of the current time step;
2. the time,
t
;
3. the time step,
dt
;
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δx2max);
11. the time step in units of the conductive time step,
dtchi
=δt/(cδt,vδx2max).
The entries in this list can be added, removed or reformatted in the file ‘print.in’, see
Sects 5.5 and J.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
z
uz
−0.02 0.00 0.02 0.04
−0.5
0.0
0.5
1.0
1.5
z
Entropy s
−0.6−0.4−0.2 0.0 0.2
−0.5
0.0
0.5
1.0
1.5
z
Temperature T
0.5 1.0 1.5 2.0
−0.5
0.0
0.5
1.0
1.5
z
Figure 2: Stratification 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
significance); 0.68 < z < 0is a stably stratified 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 significance).
If you have
IDL
, you can visualize the stratification 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 first meshpoint in the zdirection. There
have been some now outdates specific 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 defines the commands and flags used to run MPI jobs (and which
is sourced by the scripts ‘start.csh’ and ‘run.csh’). Try
csh -v getconf.csh
or
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 files 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
or
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 flux rings as initial condition, for example); and ‘kinematic/’ con-
tains kinematic dynamo problems where the hydrodynamics is turned off entirely. The
file ‘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 configuration files (cur-
rently these are ‘Makefile.local’ and ‘cparam.local’) and possibly ‘ctimeavg.local’. To
compile the samples, links the files ‘.f90’, ‘.c’ and ‘Makefile.src’ need to be linked from
the top file[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
or
DX
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’, finally is used by us for internal purposes.
The files ‘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 files, 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 significant 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-efficient 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 efficient 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 file ‘src/Makefile.local’ in which you choose certain
modules
2,
which tell the code whether or not entropy, magnetic fields, 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 file 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 fields: 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. I 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 flexibility 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 fixing 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 define F90 modules: For example each of the modules
gravity simple
,
grav r
and
nogravity
defines 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 files specify the startup and runtime parameters (see Sects. 5.12 and J.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 configuration files by issuing the
command pc_configtest.
4.3.2 ‘datadir.in
If this file 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 files, ‘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 files are part
of the setup of a given simulation and are kept under
svn
in the individual ‘samples
directories.
The file ‘src/cparam.inc’ is automatically generated by the script ‘mkcparam’ and contains
the number of fundamental variables for a given setup.
All other files in ‘src/’ are either links to source files (and ‘Makefile.src’) in the
$PENCIL_HOME/src’ directory, or object and module files 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 first line of
datadir.in’, if that file is present; see §4.3.2) contains the output from the code:
data/dim.datThe global array dimensions.
data/legend.datThe header line specifying the names of the diagnostic variables in
time_series.dat’.
data/time_series.datTime series of diagnostic variables (also printed to stdout). You
can use this file directly for plotting with
Gnuplot
,
IDL
,
Xmgrace
or similar tools (see
also §5.19).
data/tsnap.dat’, ‘data/tvid.datTime when the next snapshot ‘VARN’ or animation
slice should be taken.
data/params.logKeeps a log of all your parameters: ‘start.x’ writes the startup pa-
rameters to this file, ‘run.x’ appends the runtime parameters and appends them anew,
each time you use the ‘RELOAD’ mechanism (see §5.10).
data/param.nmlComplete set of startup parameters, printed as Fortran namelist.
This file 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.nmlComplete set of runtime parameters, printed as Fortran namelist.
This file is read by
IDL
(if you use it).
data/index.proCan be used as include file in
IDL
and contains the column in which
certain variables appear in the diagnostics file (‘time_series.dat’). It also contains the
positions of variables in the ‘VARN’ files. These positions depend on whether
entropy
or
noentropy
, etc, are invoked. This is a temporary solution and the file may disappear in
future releases.
data/interstellar.datUnformatted file 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 files:
var.datbinary file containing the latest snapshot;
16 THE PENCIL CODE
VARNbinary file containing individual snapshot number N;
dim.datASCII file containing the array dimensions as seen by the given processor;
time.datASCII file containing the time corresponding to ‘var.dat’ (not actually used
by the code, unless you use the
io mpiodist.f90
module);
grid.datbinary file containing the part of the grid seen by the given processor;
seed.datthe random seed for the next time step (saved for reasons of reproducibility).
For multi-processor runs with velocity forcing, the files ‘procN/seed.dat’ must all
contain the same numbers, because globally coherent waves of given wavenumber
are used;
X.xy’, ‘X.xz’, ‘X.yztwo-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 Configuring 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-specific configuration:
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 config files in
$PENCIL_HOME/config/compilers/*.conf’, e.g. ‘Intel.conf’ and ‘Intel_MPI.conf’.
If you have to set up some compiler options specific to a certain host system you work on,
or if you like to create a host-specific configuration file so that you can simply execute
pc_build without any options, you can clone an existing host-file, just include an existing
compiler configuration, and simply only add the options you need. A good example of
a host-file 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, configuration, compiling and running the code work like this:
1. Create a configuration file 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 configuration file
Commands like pc_build and pc_run use the Perl module ‘Pencil::ConfigFinder’ to lo-
cate an appropriate configuration file 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 file ‘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 file
$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 file 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 configuration file for your host ID, and for
locating included files (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 files ‘./host-ID’, ‘$PENCIL_HOME/host-ID’, ‘$HOME/.pencil/host-ID’, ex-
ists, its first line is used.
4. If ‘ConfigFinder’ can get hold of a
fully qualified 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 config file 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 first one for which a configuration file 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 configuration files
It is strongly recommended to include in a users configuration file one of the preset
compiler suite configuration files. Then, only minor options need to be set by a user, e.g.
the optimization flags. One of those user configuration files 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 Configuring the code to compile and run on your computer 19
# Simple Intel compiler suite config file.
%include compilers/Intel
A more complex file (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 sufficient 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 configurations.
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’ files consist of the following elements:
Comments: A#sign and any text following it on the same line are ignored.
Sections: There are three sections:
Makefile for setting make parameters
runtime for adding compiler flags 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 files it points to.3
The included path gets a .conf suffix appended. Included paths are typically “ab-
solute”, e.g.
%include os/Unix
will include the file ‘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 file’s contents inside that section are
inserted.
20 THE PENCIL CODE
%include ./Unix
will only search in the directory where the including file 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 redefining certain flags.
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 configuration files to use. The fourth example shows how to define 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’file]Makefile.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 configuration mechanism for compila-
tion. If for whatever reason one needs to work with a modified ‘Makefile’, there is a
mechanism for picking the right set of compiler flags 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 flags ‘-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 Makefile)
## OSF1:
[. . . ] (leave as it is in the original Makefile)
## Linux:
[. . . ] (leave everything from original Makefile 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 Makefile)
## UNICOS/mk:
[. . . ] (leave as it is in the original Makefile)
## HI-UX/MPP:
[. . . ] (leave as it is in the original Makefile)
## AIX:
[. . . ] (leave as it is in the original Makefile)
22 THE PENCIL CODE
Table 1: Compiler flags 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 Makefile: ‘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 flags, 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 flags ‘-fast -u’ for both compilation and linking. Ta-
ble 1 summarizes flags 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 first and second derivatives can be calculated by
df
dz =df
dz =1
zf(ζ),d2f
dz2=1
z2f′′(ζ)z′′
z3f(ζ),(1)
which can be written somewhat more compactly using the inverse function ζ(z):
df
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′′/z2=
ζ′′, and stores them in ‘data/procN/grid.dat’.
The parameters
lequidist
(a 3-element logical array),
grid func
,
coeff grid
(a
2-element real array) are used to choose a non-equidistant grid and define
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 profile: 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=Nz1for a grid of Nzvertical layers [excluding
ghost layers]), and ζis the ζvalue of the inflection point of the sinh function. The z
coordinate and ζvalue of the inflection 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 2zz0
Lz1tanh aζ2ζ1
2.(5)
General profile: 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 flow: The profile function grid_function=’duct’ generates a grid profile for tur-
bulent flow in ducts. For a duct flow, most gradients are steepest close to the walls, and
hence very fine resolution is required near the walls. Here we follow the method of [18]
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=(Nyj)π
Ny1, 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 fine enough resolution, the first 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
NyN
y=π
arccos1δ
h+ 1 (11)
=πrh
2δ+ 1 π
24r2δ
h+O"δ
h3/2#(12)
grid points in the y-direction. After rounding up to the next integer value, we find that
the truncated condition
Ny&πrh
2δ'+ 1 (13)
(where xdenotes the ceiling function, i.e. the smallest integer equal to, or larger than,
x) gives practically identical results.
Example: To apply the sinh profile, 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 zand a; the bottom height z0and the total
box height Lzare taken from
xyz0
and
Lxyz
as in the equidistant case. The inflection
point of the sinh profile (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. J.2), the code writes monitor-
ing output to
stdout
and, parallel to this, to the file ‘data/time_series.dat’. The vari-
ables that appear in this listing and the output format are defined in the file ‘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
t
in the first 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
like
----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 files
5.6.1 Snapshot files
Snapshot files contain the values of all evolving variables and are sufficient to restart a
run. In the case of an MHD simulation with entropy equation, for example, the snapshot
files will contain the values of velocity, logarithmic density, entropy and the magnetic
vector potential.
There are two kinds of snapshot files: 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 difficult 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 specifiers are like in Fortran, apart from the fact that the Eformat will use standard
scientific 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 files ‘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
files 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 files are the only way to
reconstruct the time evolution of qwithout re-running the code.
File structure Snapshot files consist of the following
Fortran records
5:
1. variable vector f[
mx
×
my
×
mz
×
nvar
]
2. time t[1], coordinate vectors x[
mx
], y[
my
], z[
mz
], 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) floating point
numbers, unless you use double precision (see §5.21, in which case all numbers are 8-
byte floating point numbers, while the record markers remain 4-byte integers.
The script pc_tsnap allows you to determine the time tof a snapshot file:
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 files and slices
We use the terms
video files
and
slice files
interchangeably. These files 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 files 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
iy
; file suffix .xz)
2. y-zplane (yindex
ix
; suffix .yz)
3. x-yplane (yindex
iz
; suffix .xy)
4. another slice parallel to the x-yplane (yindex
iz2
; suffix .xy2)
You can specify the position of the slice planes by setting the parameters
ix
,
iy
,
iz
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 zposition in terms of z using the tags
zbot slice
and
ztop slice
. In this case,
the
zbot slice
slice will have the suffix .xy and the
ztop slice
the suffix .xy2
In the file ‘video.in’ of your run directory, you can choose for which variables you want
to get video files; valid choices are listed in §J.4.
5A
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 files and slices 27
The
slice files
are written in each processor directory ‘data/proc*/’ and have a file 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 flag ‘/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 field, in planes passing through the
center of the box, do the following:
1. Write the following lines to ‘video.in’ in your run directory:
ss
bb
divu
2. Edit the namelist
run pars
in the file ‘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 files 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 defined your box setting
the z= 0 point at the surface. Therefore, your ‘start.in’ file 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 xyslice setting
zbot slice
and
ztop slice
in the ‘run.in’ file:
&run_pars
slice_position=’c’
zbot_slice=0.
ztop_slice=0.2
In this case, the slices with the suffix .xy will be at the “surface” and the ones with the
suffix .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:
bb
2. Edit the namelist
run pars
in the file ‘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 files 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 files
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[
nx
×
ny
×
nz
], time t1[1], position1[1]
5.8 Averages 29
2. data2[
nx
×
ny
×
nz
], time t2[1], position2[1]
3. data3[
nx
×
ny
×
nz
], time t3[1], position3[1]
etc.
5.8 Averages
5.8.1 One-dimensional output averaged in two dimensions
In the file ‘xyaver.in’, z-dependent (horizontal) averages are listed. They are written to
the file ‘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 files ‘xyaver.in’,
xzaver.in’, and ‘yzaver.in’ where the first letters indicate the averaging directions.
The output is then stored to the files ‘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 defined mean magnetic fields are referred to as
bmz
,
bmy
and
bmx
, respectively, and the rms values of the so defined mean velocity fields 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 files ‘yaver.in’,
zaver.in’, and ‘phiaver.in’ where the first letter indicates the averaging direction.
The output is then stored to the files ‘data/yaverages.dat’, ‘data/zaverages.dat’, and
data/phiaverages.dat’.
See Sect. 9.2 on how to add new averages.
Disadvantage: The output files, e.g. ‘data/zaverages.dat’, can be rather big because each
average is just appended to the file.
5.8.3 Azimuthal averages
Azimuthal averages are controlled by the file ‘phiaver.in’, which currently supports the
quantities listed in Sect. J.5. In addition, one needs to set
lwrite phiaverages
,
lwrite -
yaverages
, or
lwrite zaverages
to
.true.
. For example, if ‘phiaver.in’ contains the single
line
b2mphi
then you will get azimuthal averages of the squared magnetic field B2.
30 THE PENCIL CODE
Azimuthal averages are written every
d2davg
time units to the files
data/averages/PHIAVGN’. The file format of azimuthal-average files consists of the
following
Fortran records
:
1. number of radial points Nr,φavg [1], number of vertical points Nz,φavg [1], number
of variables Nvaravg[1], number of processors in zdirection [1]
2. time t[1], positions of cylindrical radius rcyl [Nr,φavg] and z[Nzavg] 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 [Nvaravg]
All numbers are 4-byte numbers (floating-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 file ‘src/ctimeavg.local’, since they use extra
memory. They are controlled by the averaging time τavg (set by the parameter
tavg
in
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=
t
Z
t0
e(tt)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 first-order accuracy in δt. In practice, however, δt is small enough to make this
accuracy suffice.
6At some point we may also implement the more straight-forward average
hfit+δt =hfit+δt
tt0+δt[f(t+δt)− hfit],(14)
which is equivalent to
hfit=1
tt0
t
Z
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 floating-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 13and 68(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 files ‘data/proc*/TAVN’. The current time averages are saved periodically in
data/proc*/timeavg.dat’ whenever ‘data/proc*/var.dat’ is written. The file 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-mkfile 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 files from one directory (source)
to another (dest). See documentation in file.
pc copyvar v v dir Copies all ‘var.dat’ files 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
specified 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
files. It also checks in the files ‘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 files into
DX
. They are called internally by some of the
DX
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 specifies 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 specified for different machines.
mkdotin Generate minimal ‘start.in’, ‘run.in’ files 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 files in a processor tree, such as slice files.
mkwww Generates a template HTML file for describing a run of the code, showing
input parameters and results.
move-slice Moves all the slice files from a processor tree structure, ‘procN/’, to a new
target tree creating directories where necessary.
nl2idl Transform a Fortran
namelist
(normally the files ‘param.nml’, ‘param2.nml’ writ-
ten by the code) into an
IDL
structure. Generates an
IDL
file that can be sourced
from ‘start.pro’ or ‘run.pro’.
pacx-adapt-makefile Version of adapt-makefile 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’ files, runs pc_setupsrc, and makes also a copy of the old
*.in’ and ‘k.dat’ files.
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’ file for details.
5.10 RELOAD and STOP files 33
pc inspectrun Check the execution of the current run: prints legend and the last few
lines of the ‘time_series.dat’ file. It also appends this result to a file 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 videofiles.csh The script for running read videofiles.x.
remote-top Create a file ‘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 files 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 specified destination.
pc setupsrc Link ‘start.csh’, ‘run.csh’ and ‘getconf.csh’ from ‘$PENCIL_HOME/bin’.
Generate ‘src/’ if necessary and link the source code files 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 file names from shell
scripts through the backtick mechanism.
pc tsnap Extract time information from a snapshot file, ‘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 figuring out how something else works.
5.10 RELOAD and STOP files
The code periodically (every
it
time steps) checks for the existence of two files, ‘RELOAD
and ‘STOP’, which can be used to trigger certain behavior.
Reloading run parameters In the directory where you started the code, create the file
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
it
, see Sect. 5.12).
Each time the parameters are reloaded, the new set of parameters is appended (in the
form of
namelists
) to the file ‘data/params.log’ together with the time t, so you have a
full record of your changes. If ‘RELOAD’ contains any text, its first 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 file ‘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 files
After the code finishes (e.g., when the final timestep number is reached or when a ‘STOP
file is found), the ‘run.csh’ script checks whether there is a ‘RERUN’ file. 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’ file, where you can enter a new direc-
tory path (relative path is ok, e.g. ../conv-slab). If nothing is written in this file (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’ files back to and from the sever. This can be
useful if you want to run with new data files, 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’ files, you may want to try out the pc_jobtransfer command. It writes the
decisive ‘STOP’ file only after the script has checked that the content of the two ‘NEWDIR
files 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-
cific 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 file 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):
ip
, 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 file ‘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 configuration files.
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 [24], 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 104.
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
nt
= 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 file ‘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×103(so the mesh
Reynolds number is about urmsδx/ν = 0.3×(2π/32)/5×10312, 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 J.
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 flow 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 field 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 fields will be in
[B] = pµ0[ρ] [v]2,(18)
and temperatures are in
[T] = [v]2
cp
=γ1
γ
[v]2
R.(19)
For some choices of density and velocity units, Table 2 shows the resulting units of
magnetic field 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 specified by
unit temperature
as described below.
5.14 Minimum amount of viscosity 37
Table 2: Units of magnetic field 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 ×105K
1 g/cm31 cm/s 3.54 ×104T = 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-defined ionization energies), the corresponding routines in the
code need to know the units you are using. This information is specified 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 coefficients (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 sufficiently large diffusion
coefficients, 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 specified as Courant time step through the coefficients
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
1
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 fix the time step δt to a value independent of velocities and diffusivities, explicitly set
the run parameter
dt
, rather than
cdt
or
cdtv
(see p. 170).
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 fifth 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×1012. 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 fill 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
bc
{
x,y,z
}in ‘run.in’; the nomenclature
used is as follows. Set
bc
{
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:
pperiodic boundary condition
aantisymmetric condition w. r. t. the boundary, i. e. vanishing value
40 THE PENCIL CODE
ssymmetric condition w. r. t. the boundary, i. e. vanishing first derivative
a2antisymmetry w. r. t. the arbitrary value on the boundary, i. e. vanishing second
derivative
c1special boundary condition for ln ρand s: constant heat flux through the boundary
c2special boundary condition for s: constant temperature at the boundary — requires
boundary condition a2 for ln ρ
cTspecial boundary condition for sor ln T: constant temperature at the boundary (for
arbitrarily set ln ρ)
cespecial boundary condition for s: set temperature in ghost points to value on bound-
ary (for arbitrarily set ln ρ)
dblow-order one-sided derivatives (“no boundary condition”) for density
sheshearing-sheet boundary condition (default when the module
Shear
is used)
gforce the value of the corresponding field on vertical boundaries (should be used
in combination with the force lower bound and force upper bound flags set in the
namelist init pars)
hsspecial 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’ file 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
to
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 file, ‘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 file size and date, and restart from a good ‘VARNfile 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’) files into a new directory tree. You may also
need the ‘data/proc*/seed.dat’ files for the random number generator. The easiest way
to get all these other files is to run start.csh again on the new machine (or in a new
directory) and then to overwrite the ‘data/proc*/var.dat’ files with the correct ones.
For restarting from runs that didn’t have magnetic fields, passive scalar fields, or test
fields, 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 files 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 flag ‘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
dt
(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
DX
(
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
DX
. It is quite easy to read in a snapshot
file from
DX
(the only tricky thing is the four extra bytes at the beginning of the file,
representing a Fortran record marker), and whenever you run ‘start.x’, the code writes
a file ‘var.general’ that tells
DX
all it needs to know about the data structure.
As a starting point for developing your own
DX
programs or
networks
, you can use a
few generic
DX
scripts provided in the directory ‘dx/basic/’. From the run directory,
start
DX
with
unix> dx -edit $PENCIL_HOME/dx/basic/lnrho
42 THE PENCIL CODE
to load the file ‘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
fit 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 field (‘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,
DX
can only read in data files 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 sufficient, but you can use GDL for many of the
visualization tasks. If you get spurious “Error opening file” 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 GDLs ‘examples/pro/’ directory to your
!PATH
variable. So the first 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 files, 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 sufficient.
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 sufficient 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 file 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 file ‘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 efficiency 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 field (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
ts
, 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>
[specific commands]
You run ‘start.pro’ once to initialize (or reinitialize, if the mesh size has changed, for
example) the fields 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 fix 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 first 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 file 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 field or the current density, you can calculate them in IDL by
11
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 efficient 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 first 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 files that read different data into structures. They are placed in
the directory $PENCIL_HOME/idl/files. Here is a list of files (including suggested options
to call them with)
pc_read_var_raw, obj=var, tags=tags
Efficiently 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 files.
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 file 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 find 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 file .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 defined 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 file 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 file using the GUI or with
In [1]: plt.savefig(’plot.eps’)
Reading and Plotting VAR files and slice files Read var files:
In [1]: var = pc.read_var()
5.20 Running on multi-processor computers 49
Read slice files:
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
t.
If you want to plot e.g. the x-component of the magnetic field 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 files and slice files.
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 efficient way of
parallelizing finite-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
by
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
by
integer, parameter :: ncpus=16, nprocy=4, nprocz=ncpus/nprocy, nprocx=1
integer, parameter :: nxgrid=128, nygrid=nxgrid, nzgrid=nxgrid
50 THE PENCIL CODE
The first line specifies 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 find 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 significant 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
Pw
=Nx×Ny×Nz
Px×Py×Pz
.(26)
However, for finite 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 finite 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(eff)
w
Pw
=Nx
Px
+Q×Ny
Py
+Q×Nz
Pz
+Q.(27)
There is efficient scaling only when
min Nx
Px
,Ny
Py
,Nz
PzQ. (28)
In the special case were Nx=Ny=NzN=N1/3
w, with Px= 1 and Py=PzP=P1/2
w,
we have
N(eff)
w
Pw
= (N+Q)×N
P+Q2
.(29)
For N= 128 and Q= 6 the effective mesh size exceeds the actual mesh size by a factor
N(eff)
w
Nw
= (N+Q)×N
P+Q2
×Pw
Nw
.(30)
These factors are listed in Table 3.
Table 3: N(eff)
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 floating 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.
If
REAL PRECISION
is set to ‘double’, the flag
FFLAGS DOUBLE
is appended to
the Fortran compile flags. 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 flags in ‘src/Makefile.inc’, which is the first 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 fine at single
precision). To continue working in double precision, you just say lread_from_other_-
prec=T in run_pars; see Section 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
Nx1
X
p=0
Ny1
X
q=0
Nz1
X
r=0
u(xp, yq, zr)
×exp(ikxxp) exp(ikyyq) exp(ikzzr),(31)
where
|kx|<πNx
Lx
,|ky|<πNy
Ly
,|kz|<πNz
Lz
,
when Lis the size of the simulation box. The three-dimensional power spectrum P(k)is
defined 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 files ‘powerhel_mag.dat’ and ‘power_kin.dat
containing the three-dimensional magnetic and kinetic power spectra respectively. In
addition to these three-dimensional files we will also find the one-dimensional files
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 files the data are stored such that the first 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 first computation. The relevant out-
put file 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 define the p-th order longitudinal structure function of uas
Sp
long(l) = h|ux(x+l, y, z)ux(x, y, z)|pi,(34)
while the transverse is
Sp
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’, define
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 files:
sfu-1.dat’, ‘sfu-2.dat’, ‘sfu-3.dat’ (velocity),
sfb-1.dat’, ‘sfb-2.dat’, ‘sfb-3.dat’ (magnetic field),
sfz1-1.dat’, ‘sfz1-2.dat’, ‘sfz1-3.dat’ (first Elsasser variable),
sfz2-1.dat’, ‘sfz2-2.dat’, ‘sfz2-3.dat’ (second Elsasser variable),
which contains the data of interest. The first line in each file contains the time tand
the number
qmax
, such that the largest moment calculated is
qmax
1. The next
imax
numbers represent the first moment structure function for the first 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 first snapshot.
To visualize with
IDL
just type structure and you get the time-average of the first 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 suffices to have only three position variables as the velocity of the particles are
equal to the instantaneous fluid 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 define 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 specifications 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 influence 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 benefit 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 filled 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 final 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 first 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 fixed 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 final particle density (defined 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 field 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 field
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 defines 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 file. This works because the data type of npar is only set here. It is
worth noting that one should not use the flag
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 beneficial 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 flow.
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
Since the spring of 2007 spherical and cylindrical polar coordinates have been imple-
mented, although this development is not yet completed. Spherical coordinates are in-
voked by adding the following line in the file ‘start.in
&init_pars
coord_system=’spherical_coords’
Another possibility is to put cylindrical_coords instead. In practice, the names (x, y, z)
are still used, but they refer then to (r, θ, φ)or (r, φ, z)instead.
Bug reports, corrections, and improvements on these are appreciated.
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 simplifications 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.(37)
Here ρdenotes density, uthe fluid 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
ρ1p=c2
s(s/cp+ln ρ), where the squared sound speed is given by
c2
s=γp
ρ=c2
s0 exp γs/cp+ (γ1) ln ρ
ρ0,(38)
and γ=cp/cvis the ratio of specific heats, or adiabatic index. Note that c2
sis proportional
to the temperature with c2
s= (γ1)cpT.
The equation of motion is accordingly
Du
Dt=c2
ss
cp
+ ln ρΦgrav +j×B
ρ
+ν2u+1
3∇∇ ·u+ 2S·ln ρ+ζ(∇∇ ·u) ; (39)
Here Φgrav is the gravity potential, jthe electric current density, Bthe magnetic flux
density, νis kinematic viscosity, ζdescribes a bulk viscosity, and, in Cartesian coordi-
nates
Sij =1
2ui
xj
+uj
xi2
3δij·u(40)
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 =
1
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.(41)
Here Ais the magnetic vector potential, B=∇ × Athe magnetic flux 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 +·(KT) + ηµ0j2+ 2ρνSS+ζρ (·u)2.(42)
Here, Tis temperature, cpthe specific 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 specific entropy, s. The heat conduction term on
the right hand side can be written in the form
·(KT)
ρT (43)
=cpχh2ln T+ln T·(ln T+ ln χ+ ln ρ)i(44)
=cpχγ2s/cp+ (γ1)2ln ρ
+cpχ[γs/cp+ (γ1)ln ρ]·[γ(s/cp+ln ρ) + ln χ],(45)
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
·(KT)
ρT =K
ρh2ln T+ (ln T)2i(46)
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 ,(47)
where τij = 2νρSij is the stress tensor. The viscous power density Pis
P=uiτij,j (48)
=
xj
(uiτij)ui,j τij (49)
The term under the divergence is the viscous energy flux and the other term is the
kinetic energy loss due to heating. The heating term +ui,jτij is positive definite, 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).(50)
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)(51)
=1
2(ui,j +uj,i 1
3δij·u)(2νρSij )(52)
= 2νρS2,(53)
which is positive definite.
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 KT) 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 (54)
is replaced by
ln T
t =u·ln T+1
ρcvTRHS (γ1) ·u,(55)
where RHS is the right hand side of equation (42).
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ρDc] = 0.(56)
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=D2ln c+ (ln ρ+ln c)·ln c.(57)
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 flows, 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 flow convergence, maximum over five zones, and smoothed to
second order,
ζshock =cshock Dmax
5[(·u)+]E(min(δx, δy, δz))2,(58)
where cshock is a constant defining the strength of the shock viscosity. In the code this
dimensionless coefficient 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.(59)
Assume ν=const, but ζ6=const, so
ρ1Fvisc =ν2u+1
3∇∇ ·u+ 2S·ln ρ+ζshock [∇∇ ·u+ (ln ρ+ln ζshock)·u].
(60)
and
ρ1Γvisc = 2νS2+ζshock(·u)2.(61)
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 configuration only hydrogen ionization is explicitly included. Other con-
stituents (currently He and H2) can have fixed values. The pressure is proportional to
the total number of particles, i.e.
p= (nHI +nHII +nH2+ne+nHe +...)kBT. (62)
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 define xene/nHtot,xHe
nHe/nHtot, and xH2nH2/nHtot. Substituting nH=nHtot 2nH2, we have
p= (1 xH2+xe+xHe +...)nHtotkBT. (63)
This can be written in the more familiar form
p=R
µρT, (64)
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
1xH2+xe+xHe
(65)
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
2Rand
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’ file. The equation of state described
below works for variable ionization, and the entropy offset is different from that used in
Eq. (38), 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
as
y2
1y=1
nHmekBT
2π~23/2
exp χH
kBT.(66)
The temperature Tcannot be obtained directly from the PENCIL CODEs independent
variables ln ρand s, but is itself dependent on y. Hence, the calculation of yessentially
becomes a root finding 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π~23/2#+5
2!.(67)
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 specific entropy sin terms of y,ln ρand T, which may be solved for
T.
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 finding the
root of
f(y) = ln "1y
y2
1
nHmekBT(y)
2π~23/2#χH
kBT(y).(68)
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-fluid model
66 THE PENCIL CODE
ρi
t =·(ρiui) + G(69)
ρn
t =·(ρnun)− G (70)
(ρiui)
t =·(ρiui:ui)pi+pe+B2
2µ0+F(71)
(ρnun)
t =·(ρnun:un)pn− F (72)
A
t =ui×B(73)
where the subscripts nand iare for neutral and ionized, respectively. The terms Gand
F, through which the two fluids exchange mass and momentum, are given by
G=ζρnαρ2
i(74)
F=ζρnunαρ2
iui+γρiρn(unui).(75)
In the above equations, ζis the ionization coefficient, αis the recombination coefficient,
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. (69) and Eq. (71) are solved in ‘density.f90’ and ‘hydro.f90’ respectively.
Equation 70 is solved in ‘neutraldensity.f90’ and Eq. (72) in neutralvelocity.f90’. The
sample ‘1d-test/ambipolar-diffusion’ has the current setup for a two-fluid 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
dI
=I+S , (76)
where
τ
s
Z
0
χ(s)ds(77)
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 fluid in the simulation box by
adding the term
u
t =...φself (78)
to the equation of motion. The self-potential φself (or just φfor simplicity) satisfies Pois-
son’s equation
2φ= 4πGρ . (79)
The solution for a single Fourier component at scale kis
φk=4πk
k2.(80)
Here we have assumed periodic boundary conditions. The potential is obtained by
Fourier-transforming the density, then finding the corresponding potential at that scale,
and finally 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 [13]
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 [13] 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 fluids14. 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
(udu).(81)
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ρ,(82)
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 first is coagulation controlled
by the discrete coagulation equation
dnk
dt=1
2X
i+j=k
Aijninjnk
X
i=1
Aikni.(83)
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
dN
dt=1
τcond
Nd1
d.(84)
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
1
τcond
=A1vthαnmon 11
Smon ,(85)
where A1is the surface area of a monomer, αis the condensation efficiency, nmon is the
number density of monomers in the gas and Smon is the saturation level of the monomer
given by
Smon =Pmon
Psat
.(86)
Here Psat is the saturated vapor pressure of the monomer. Currently only water ice has
been implemented in the code.
All dust species fulfill the continuity equation
ρd
t +·(ρdud) = 0.(87)
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 (39) the cosmic ray pressure force, ρ1pcis added on the right hand
side, and ecsatisfies the evolution equation
ec
t +·(ecu) + pc·u=i(Kij jec) + Qc,(88)
where Qcis a source term and
Kij =Kδij + (KkK)ˆ
Biˆ
Bj(89)
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(Kijjec) = Uc·ec+Kij ijec.(90)
where Uci=Kij /∂xjacts like an extra velocity trying to straighten magnetic field
lines. We can write this term also as Uc=(KkK)·(ˆ
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 infinite. 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 flux, and to replace the
diffusion term in Eq. (88) by a divergence of a flux that in turn obeys the equation
Fci
t =˜
KijjecFci
τ(non-Fickian diffusion),(91)
where Kij =τ˜
Kij would be the original diffusion tensor of Eq. (89), 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 fluid only has a velocity vector field (the continuity equation of a fluid 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,(92)
where the index iruns over all particles. Here uis the gas velocity at the position of
the particle. One can choose between a first 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.(93)
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 fluid, 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
(viu).(94)
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 first order interpolation is used to obtain the gas velocity field 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 field 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)
p
ρv(i)
t drag
(95)
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 ais comparable to the mean
free path λof the gas, for which a fluid 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 field.
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
is
6.16 Test-field equations 71
dvpi
dt =Fgi
N
X
j6=i
GMj
R2
ij
ˆ
Rij (96)
where Rij =|rpirpj|is the distance between particles iand j, and ˆ
Rij is the unit vector
pointing from particle jto particle i. The first 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, (97)
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 define 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 fixed time-step tp2π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=103, 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 specifically 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-field equations
The test-field method is used to calculate turbulent transport coefficients for magneto-
hydrodynamics. This is a rapidly evolving field and we refer the interested reader to
recent papers in this field, e.g. by Sur et al. (2008) or Brandenburg et al. (2008). For
technical details, see also Sect. F.3.
72 THE PENCIL CODE
7 Troubleshooting / Frequently Asked Questions
7.1 Download and setup
7.1.1 Download forbidden
A: Both GitHub and SourceForge are banned from countries on the United
States Office of Foreign Assets Control sanction list, including Cuba, Iran, Libya,
North Korea, Sudan and Syria; see http://de.wikipedia.org/wiki/GitHub and
http://en.wikipedia.org/wiki/SourceForge. As a remedy, you might download a tar-
ball from http://pencil-code.nordita.org/; see also Section 2.
7.1.2 When sourcing the ‘sourceme.sh’/‘sourceme.csh’ file or running pc_setupsrc, I get
error messages from the shell, like ‘if: Expression Syntax.’ or ‘set: Variable name
must begin with a letter.
A: This sounds like a buggy shell setup, either by yourself or your system administrator
— or a shell that is even more idiosyncratic than the ones we have been working with.
To better diagnose the problem, collect the following information before filing a bug
report to us:
1. uname -a
2. /bin/csh -v
3. echo $version
4. echo $SHELL
5. ps -p $$
6. If you have problems while sourcing the ‘sourceme’ script,
(a) unset the PENCIL_HOME variable:
for
csh
and similar: unsetenv PENCIL_HOME
for
bash
and similar: unexport PENCIL_HOME; unset PENCIL_HOME
(b) switch your shell in verbose mode,
for
csh
and similar: set verbose; set echo
for
bash
and similar: set -v; set -x
then source again.
7. If you have problems with pc_setupsrc, run it with csh in verbose mode:
/bin/csh -v -x $PENCIL_HOME/bin/pc_setupsrc
7.2 Compilation 73
7.2 Compilation
7.2.1 Linker can’t find the syscalls functions:
ld: 0711-317 ERROR: Undefined symbol: .is_nan_c
ld: 0711-317 ERROR: Undefined symbol: .sizeof_real_c
ld: 0711-317 ERROR: Undefined symbol: .system_c
ld: 0711-317 ERROR: Undefined symbol: .get_env_var_c
ld: 0711-317 ERROR: Undefined symbol: .get_pid_c
ld: 0711-317 ERROR: Undefined symbol: .file_size_c
A: The Pencil Code needs a working combination of a Fortran- and a C-compiler. If this
is not correctly set up, usually the linker won’t find the functions inside the
syscalls
module. If that happens, either the combination of C- and Fortran-compiler is inappro-
priate (e.g.
ifort
needs
icc
), or the compiler needs additional flags, like
g95
might need
the option ‘-fno-second-underscore’ and
xlf
might need the option ‘-qextname’. Please
refer to Sect. 5.2, Table 1.
7.2.2 Make gives the following error now:
PGF90-S-0017-Unable to open include file: chemistry.h (nochemistry.f90: 43)
0 inform, 0 warnings, 1 severes, 0 fatal for chemistry
Line 43 of the nochemistry routine, only has ’contains’.
A: This is because somebody added a new module (together with a corresponding
nomod-
ule.f90
and a
module.h
file (chemistry in this case). These files didn’t exist before, so you
need to say:
pc_setupsrc
If this does not help, say first make clean and then pc_setupsrc.
7.2.3 How do I compile the PENCIL CODE with the Intel (
ifc
) compiler under
Linux
?
A: The PENCIL CODE should compile successfully with
ifc
6.x,
ifc
7.0, sufficiently recent
versions of
ifc
7.1 (you should get the latest version; if yours is too old, you will typically
get an ‘internal compiler error’ during compilation of ‘src/hydro.f90’), as well as with
recent versions of
ifort
8.1 (8.0 may also work).
You can find the
ifort
compiler at ftp://download.intel.com/software/products/compilers/downlo
On many current (as of November 2003) Linux systems, there is a mismatch between
the
glibc
versions used by the compiler and the linker. To work around this, use the
following flag for compiling
FC=ifc -i_dynamic
and set the environment variable
LD_ASSUME_KERNEL=2.4.1; export LD_ASSUME_KERNEL
or
74 THE PENCIL CODE
setenv LD_ASSUME_KERNEL 2.4.1
This has solved the problems e.g. on a system with glibc-2.3.2 and kernel 2.4.22.
Thanks to Leonardo J. Milano (http://udel.edu/~lmilano/) for part of this info.
7.2.4 I keep getting segmentation faults with ‘start.x’ when compiling with
ifort
8.0
A: There was/is a number of issues with
ifort
8.0. Make sure you have the latest patches
applied to the compiler. A number of things to consider or try are:
1. Compile with the the ‘-static -nothreads’ flags.
2. Set your stacksize to a large value (but a far too large value may be problematic,
too), e. g.
limit stacksize 256m
ulimit -s 256000
3. Set the environment variable KMP STACKSIZE to a large value (like 100M)
See also http://softwareforums.intel.com/ids/board/message?board.id=11&message.id=1375
7.2.5 When compiling with MPI on a Linux system, the linker complains:
mpicomm.o: In function ‘mpicomm_mpicomm_init_’:
mpicomm.o(.text+0x36): undefined reference to ‘mpi_init_’
mpicomm.o(.text+0x55): undefined reference to ‘mpi_comm_size_’
mpicomm.o(.text+0x6f): undefined reference to ‘mpi_comm_rank_’
[...]
A: This is the infamous
underscore problem
. Your
MPI
libraries have been compiled
with
G77
without the option ‘-fno-second-underscore’, which makes the
MPI
symbol
names incompatible with other Fortran compilers.
As a workaround, use
MPICOMM = mpicomm_
in ‘Makefile.local’. Or, even better, you can set this globally (for the given computer)
by inserting that line into the file ‘~/.adapt-mkfile.inc’ (see perldoc adapt-mkfile for
more details).
7.2.6 Compilation stops with the cryptic error message:
f95 -O3 -u -c .f90.f90
Error : Could not open sourcefile .f90.f90
compilation aborted for .f90.f90 (code 1)
make[1]: *** [.f90.o] Error 1
What is the problem?
A: There are two possibilities:
7.2 Compilation 75
1. One of the variables for
make
has not been set, so
make
expands it to the empty
string. Most probably you forgot to specify a module in ‘src/Makefile.local’. One
possibility is that you have upgraded from an older version of the code that did not
have some of the modules the new version has.
Compare your ‘src/Makefile.local’ to one of the examples that work.
2. One of the variables for
make
has a space appended to it, e.g. if you use the line
MPICOMM = mpicomm_
(see §7.2.5) with a trailing blank, you will encounter this error message. Remove
the blank. This problem can also occur if you added a new module (and have an
empty space after the module name in ‘src/Makefile.src’, i.e. CHIRAL=nochiral ),
in which case the compiler will talk about “circular dependence” for the file
nochiral’.
7.2.7 The code doesn’t compile,
. . . there is a problem with
mvar
:
make start.x run.x
f95 -O3 -u -c cdata.f90
Error: cdata.f90, line 71: Implicit type for MVAR
detected at MVAR@)
[f95 terminated - errors found by pass 1]
make[1]: *** [cdata.o] Error 2
A: Check and make sure that ‘mkcparam’ (directory ‘$PENCIL_HOME/bin’) is in your path. If
this doesn’t help, there may be an empty cparam.inc’ file in your ‘src’ directory. Remove
cparam.inc’ and try again (Note that ‘cparam.inc’ is automatically generated from the
Makefile’).
7.2.8 Some samples don’t even compile,
as you can see on the web, http://www.nordita.org/software/pencil-code/tests.html.
samples/helical-MHDturb:
Compiling.. not ok:
make start.x run.x read_videofiles.x
make[1]: Entering directory ‘/home/dobler/f90/pencil-code/samples/helical-MHDturb/src’
/usr/lib/lam/bin/mpif95 -O3 -c initcond.f90
/usr/lib/lam/bin/mpif95 -O3 -c density.f90
use Gravity, only: gravz, nu_epicycle
^
Error 208 at (467:density.f90) : No such entity in the module
Error 355 : In procedure INIT_LNRHO variable NU_EPICYCLE has not been given a type
Error 355 : In procedure POLYTROPIC_LNRHO_DISC variable NU_EPICYCLE has not been given
3 Errors
compilation aborted for density.f90 (code 1)
make[1]: *** [density.o] Error 1
76 THE PENCIL CODE
make[1]: Leaving directory ‘/home/dobler/f90/pencil-code/samples/helical-MHDturb/src’
make: *** [code] Error 2
A: Somebody may have checked in something without having run auto-test beforehand.
The problem here is that something has been added in one module, but not in the corre-
sponding no-module. You can of course check with
svn
who it was. . .
7.2.9 Internal compiler error with Compaq/Dec F90
The Dec Fortran optimizer has occasional problems with nompicomm.f90’:
make start.x run.x read_videofiles.x
f90 -fast -O3 -tune ev6 -arch ev6 -c cparam.f90
[...]
f90 -fast -O3 -tune ev6 -arch ev6 -c nompicomm.f90
otal vm 2755568 otal vm 2765296 otal vm 2775024
otal vm 2784752 otal...
Assertion failure: Compiler internal error - please submit problem r...
GEM ASSERTION, Compiler internal error - please submit problem report
Fatal error in: /usr/lib/cmplrs/fort90_540/decfort90 Terminated
*** Exit 3
Stop.
*** Exit 1
Stop.
A: The occurrence of this problem depends upon the grid size; and the problem never
seems to occur with ‘mpicomm.f90’, except when ncpus=1. The problem can be avoided by
switching off the loop transformation optimization (part of the ‘-O3’ optimization), via:
#OPTFLAGS=-fast -O3 -notransform_loops
This is currently the default compiler setting in ‘Makefile’, although it has a measurable
performance impact (some 8% slowdown).
7.2.10 Assertion failure under SunOS
Under SunOS, I get an error message like
user@sun> f90 -c param_io.f90
Assertion failed: at_handle_table[at_idx].tag == VAR_TAG,
file ../srcfw/FWcvrt.c, line 4018
f90: Fatal error in f90comp: Abort
A: This is a compiler bug that we find at least with Sun’s WorkShop Compiler version ‘5.0
00/05/17 FORTRAN 90 2.0 Patch 107356-05’. Upgrade the compiler version (and possi-
bly also the operating system): we find that the code compiles and works with version
‘Sun WorkShop 6 update 2 Fortran 95 6.2 Patch 111690-05 2002/01/17’ under SunOS
version ‘5.8 Generic 108528-11’.
7.2 Compilation 77
7.2.11 After some dirty tricks I got pencil code to compile with MPI, ...
> Before that i installed lam-7.1.4 from source.
Goodness gracious me, you shouldn’t have to compile your own MPI library.
A: Then don’t use the old LAM-MPI. It is long superseded by open-mpi now. Open-mpi
doesn’t need a daemon to be running. I am using the version that ships with Ubuntu
(e.g. 9.04):
frenesi:~> aptitude -w 210 search openmpi | grep ’^i’
i libopenmpi-dev - high performance message passing library -- header files
i A libopenmpi1 - high performance message passing library -- shared library
i openmpi-bin - high performance message passing library -- binaries
i A openmpi-common - high performance message passing library -- common files
i openmpi-doc - high performance message passing library -- man pages
Install that and keep your configuration (Makefile.src and getconf.csh) close to that for
‘frenesi’ or ‘norlx50’. That should work.
7.2.12 Error: Symbol ’mpi comm world’ at (1) has no IMPLICIT type
I installed the pencil code on Ubuntu system and tested "run.csh"
in ...\samples\conv-slab. Here the code worked pretty well.
Nevertheless, running (auto-test), I found there are some errors.
The messages are,
Error: Symbol ’mpi_comm_world’ at (1) has no IMPLICIT type
Fatal Error: Error count reached limit of 25.
make[2]: *** [mpicomm_double.o] Error 1
make[2]: Leaving directory
‘/home/pkiwan/Desktop/pencil-code/samples/2d-tests/selfgravitating-shearwave/src’
make[1]: *** [code] Error 2
make[1]: Leaving directory
‘/home/pkiwan/Desktop/pencil-code/samples/2d-tests/selfgravitating-shearwave/src’
make: *** [default] Error 2
Finally, ### auto-test failed ###
Will it be OK? Or, how can I fix this?
A: Thanks for letting me know about the status, and congratulations on your progress!
Those tests that fail are those that use MPI. If your machine is a dual or multi core
machine, you could run faster by running under MPI. But this is probably not crucial
for you at this point. (I just noticed that there is a ToDo listed in the auto-test command
to implement the option not to run the MPI tests, but this hasn’t been done yet. So I
guess you can start with the science next.
78 THE PENCIL CODE
7.2.13 Error: Can’t open included file ’mpif.h’
It always worked, but now, after some systems upgrade, I get
gfortran -O3 -o mpicomm.o -c mpicomm.f90
Error: Can’t open included file ’mpif.h’
When I say locate mpif.h I only get things like
/scratch/ntest/1.2.7p1-intel/include/mpif.h
But since I use FC=mpif90 I thought I don’t need to worry.
A: Since you use FC=mpif90 there must definitely be something wrong with their setup.
Try mpif90 -showme or mpif90 -show; the ‘-I’ option should say where it looks for ’mpif.h’.
If those directories don’t exist, it’s no wonder that it doesn’t work, and it is time to
complain.
7.3 Pencil check
7.3.1 The pencil check complains for no reason.
A: The pencil check only complains for a reason.
7.3.2 The pencil check reports MISSING PENCILS and quits
A: This could point to a serious problem in the code. Check where the missing pencil
is used in the code. Request the right pencils, likely based on input parameters, by
adapting one or more of the pencil_criteria_MODULE subroutines.
7.3.3 The pencil check reports unnecessary pencils
The pencil check reports possible overcalculation... pencil rho ( 43) is
requested, but does not appear to be required!
A: Such warnings show that your simulation is possibly running too slowly because it is
calculating pencils that are not actually needed. Check in the code where the unneces-
sary pencils are used and adapt one or more of the pencil_criteria_MODULE subroutines
to request pencils only when they are actually needed.
7.3.4 The pencil check reports that most or all pencils are missing
A: This is typically a thing that can happen when testing new code development for the
first time. It is usually an indication that the reference df changes every time you call
pde. Check whether any newly implemented subroutines or functionality has a “mem-