pencil-code Wiki

Brought to you by: tobson

PythonForPencil

Authors: Anonymous

Python with the Pencil Code

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, I strongly recommend to also install IPython, which is a more convenient console for python. You will also need the numpy, matplotlib and tk library.

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 you ~/.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()

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 command 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 shabe

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.

ts = pc.read_ts()

Most commands take some arguments. For most of them there is a default value, e.g.

pc.read_ts(filename='time_series.dat', datadir='data', double=0, print_std=0, quiet=0, plot_data=True, comment_char='#')

You can change the values by simply typing e.g.

pc.read_ts(plot_data = False)

Reading and Plotting Time Series

Reading the time series file is very easy. Simply type

ts = pc.read_ts()

and python plots parts of the time series and 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:

plt.plot(ts.t, ts.emag)

The standard plots are not perfect and need a little polishing. See further down about making pretty plots. You can save the plot into a file using the GUI or with

plt.savefig('plot.eps')

Reading and Plotting VAR files and slice files

Read var files:

var = pc.read_var()

Read slice files:

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:

plt.imshow(zip(*var.bb[0,128,:,:]), origin = 'lower', extent = [-4,4,-4,4], interpolation = 'nearest', cmap = 'hot')

For a complete list of arguments of plt.imshow refer to its documentation.

For a more interactive function plot use:

pc.animate_interactive(slices, t)

Be aware: arrays from the reading routines are ordered f[nvar,mz,my,mx], i.e. reversed to IDL. This affects reading var files and slice files.

Examples

Standard plots with any plotting library are not the prettiest ones. The same is true for Matplotlib. Here are a few pretty examples of plots where the default style is changed. You can add your commands into a script e.g. plot_results.py and execute it in IPython with execfile('plot_results.py').

Simple plot:

import pencil as pc
import numpy as np
import pylab as plt

# read the time_series.dat
ts = pc.read_ts(quiet = True, plot_data = False)

# prepare the plot
# set the size and margins
width = 8
height = 6
plt.rc("figure.subplot", left=0.2)
plt.rc("figure.subplot", right=0.95)
plt.rc("figure.subplot", bottom=0.15)
plt.rc("figure.subplot", top=0.90)
figure = plt.figure(figsize = (width, height))
axes = plt.subplot(111)

# make the actual plot
plt.semilogy(ts.t, ts.brms/ts.brms[0], linestyle = '-', linewidth = 2, color = 'black', label = r'$\langle\bar{B}\rangle/\langle\bar{B}\rangle(0)$')
plt.semilogy(ts.t, ts.jrms/ts.jrms[0], linestyle = '--', linewidth = 2, color = 'blue', label = r'$\langle\bar{J}\rangle/\langle\bar{J}\rangle(0)$')
plt.semilogy(ts.t, ts.jmax/ts.jmax[0], linestyle = ':', linewidth = 2, color = 'red', label = r'$J_{\rm max}/J_{\rm max}(0)$')

plt.xlabel(r'$t$', fontsize = 25)
plt.ylabel(r'$\langle\bar{B}\rangle, \langle\bar{J}\rangle, J_{\rm max}$', fontsize = 25)
plt.title('various quantities', fontsize = 25, family = 'serif')

# legend
plt.legend(loc = 1, shadow = False, fancybox = False, numpoints = 1)
leg = plt.gca().get_legend()
# change the font size of the legend
ltext = leg.get_texts() # all the text.Text instance in the legend
for k in range(len(ltext)):
    legLine = ltext[k]
    legLine.set_fontsize(25)
frame = leg.get_frame()
frame.set_facecolor('1.0')
leg.draw_frame(False)

# make plot pretty
plt.xticks(fontsize=20, family = 'serif')
plt.yticks(fontsize=20, family = 'serif')
axes.tick_params(axis = 'both', which = 'major', length = 8)
axes.tick_params(axis = 'both', which = 'minor', length = 4)

# create an offset between the xylabels and the axes
for label in axes.xaxis.get_ticklabels():
    label.set_position((0,-0.03))
for label in axes.yaxis.get_ticklabels():
    label.set_position((-0.03,0))

Simple 2d plot:

import pencil as pc
import numpy as np
import pylab as plt

# read the slices
slices, t = pc.read_slices(field = 'bb1', extension = 'xy')

# read the grid size
grid = pc.read_grid()
x0 = grid.x[3]
x1 = grid.x[-4]
y0 = grid.y[3]
y1 = grid.y[-4]

# prepare the plot
# set the size and margins
width = 8
height = 6
plt.rc("figure.subplot", left=0.15)
plt.rc("figure.subplot", right=0.95)
plt.rc("figure.subplot", bottom=0.15)
plt.rc("figure.subplot", top=0.95)
figure = plt.figure(figsize = (width, height))
axes = plt.subplot(111)

# make the actual plot
plt.imshow(zip(*slices[0,:,:]), origin = 'lower', interpolation = 'nearest', cmap = 'hot', extent = [x0, x1, y0, y1])
plt.xlabel(r'$x$', fontsize = 25)
plt.ylabel(r'$y$', fontsize = 25)

# colorbar
cb = plt.colorbar()
cb.set_label(r'$B_{x}(x,y,z=0)$', fontsize = 25)
cbytick_obj = plt.getp(cb.ax.axes, 'yticklabels')
plt.setp(cbytick_obj, fontsize = 15, family = 'serif')

# make plot pretty
plt.xticks(fontsize=20, family = 'serif')
plt.yticks(fontsize=20, family = 'serif')
axes.tick_params(axis = 'both', which = 'major', length = 8)
axes.tick_params(axis = 'both', which = 'minor', length = 4)

# create an offset between the xylabels and the axes
for label in axes.xaxis.get_ticklabels():
    label.set_position((0,-0.03))
for label in axes.yaxis.get_ticklabels():
    label.set_position((-0.03,0))

Troubleshooting

I'm an a cluster and the library LIBNAME could not be loaded.
Typically system administrators don't install all the software you need. Just contact the person in charge and ask for installing it.

I'm getting complaints about a 'tk' library.
Try launchin python with

ipython --pylab='qt'

If this doesn't work or you have only access to the python console try in Python:

plt.switch_backend('qt')

or any other backend like qtk. If you are still out of luck you can still save the plot into a file with

plt.savefig('plot.eps')

There is nothing displayed when I try plotting.
Try:

plt.show()
plt.draw()