This file is part of postpic.

postpic is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.

postpic is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.

You should have received a copy of the GNU General Public License along with postpic. If not, see http://www.gnu.org/licenses/.

Copyright Stephan Kuschel 2016-2018

First, download the example data

In [1]:

def download(url, file):
    import urllib3
    import shutil
    import os
    if os.path.isfile(file):
        return
    urllib3.disable_warnings()
    http = urllib3.PoolManager()
    print('downloading {:} ...'.format(file))
    with http.request('GET', url, preload_content=False) as r, open(file, 'wb') as out_file:
        shutil.copyfileobj(r, out_file)

# download the example data
import os
if not os.path.exists('_epochexample'):
    os.mkdir('_epochexample')
    download('https://gist.github.com/skuschel/9c78c47130a181a1e2bec13900ab4293/raw/'
           + '6588b96090caa9032a6a9ed8d785f2034f446ad2/epoch_exampledata.tar.xz',
            '_epochexample/epoch_exampledata.tar.xz')

import tarfile
tar = tarfile.open('_epochexample/epoch_exampledata.tar.xz')
tar.extractall('_epochexample/')
print('done.')

downloading _epochexample/epoch_exampledata.tar.xz ...
done.

In [2]:

%pylab inline
import postpic as pp
pp.__version__

Populating the interactive namespace from numpy and matplotlib

Out[2]:

'v0.3.1+318.ga7cb0c4'

In [3]:

pp.chooseCode('epoch')
# open a specific dump
dr = pp.readDump('_epochexample/epoch_exampledata/0420.sdf')
print(dr)
# the dumpreader knwos all the information of the simulation
print('The simulations was running on {} spatial dimensions.'.format(dr.simdimensions()))

<Sdfreader at "_epochexample/epoch_exampledata/0420.sdf">
The simulations was running on 2 spatial dimensions.

In [4]:

# if needed, postpic can be bypassed and the underlying datastructure (sdf in this case)
# can be accessed directly via keys
print(dr.keys())
dr['Header']

['Header', 'Electric Field/Ex', 'Electric Field/Ey', 'Magnetic Field/Bz', 'Particles/Weight/electron', 'Particles/Weight/electront', 'Particles/Px/electron', 'Particles/Px/electront', 'Particles/Py/electron', 'Particles/Py/electront', 'Particles/Pz/electron', 'Particles/Pz/electront', 'Particles/ID/electron', 'Particles/ID/electront', 'Grid/Particles/electron', 'Grid/Particles/electront', 'Derived/Number_Density', 'Grid/Grid', 'Grid/Grid_mid', 'CPUs/Original rank', 'CPU split/electron', 'CPU split/electront', 'CPUs/Current rank']

Out[4]:

{'code_io_version': 1,
 'code_name': 'Epoch2d',
 'file_revision': 4,
 'file_version': 1,
 'jobid1': 1520354594,
 'jobid2': 572,
 'other_domains': False,
 'restart_flag': False,
 'station_file': False,
 'step': 1489,
 'time': 5.001526760099565e-13}

The Field data¶

Field data is data, which is already on a grid, examples are Electric Field, Magnetic Field, Number Density, A velocity field of the streaming particles... Such data can be accessed by the datareader directly. It is represented in postpic by a postpic.Field object, which is a numpy array holding the data plus information about the axis.

In [5]:

def addcolorbar(ax, im):
    from mpl_toolkits.axes_grid1 import make_axes_locatable
    divider = make_axes_locatable(ax)
    cax = divider.append_axes("right", size="5%", pad=0.05)
    cbar = colorbar(im, cax=cax)
    return cbar

In [6]:

ez = dr.Ex()
print(ez)
fig, ax = subplots()
im = ax.imshow(ez.T, origin='lower', extent=ez.extent*1e6)
addcolorbar(ax, im)
ax.set_xlabel('y [$\mu m$]');
ax.set_ylabel('z [$\mu m$]');

<Field "Ex" (240, 150)>

Particle data¶

Particle data is represented in postpic by a postpic.MultiSpecies object. It contains a number of particles and it offers a very easy to use interface to access the properties or derived quantities or even use derived quantities to create a postpic.Field object, which can finally be plotted.

In [7]:

# the multispecies Object is used to access particle data
print(dr.listSpecies())
ms = pp.MultiSpecies(dr, 'electron')
# now, ms is representing the species "electrons"
print(ms)

['electron', 'electront']
<MultiSpecies including all "electron" (36248)>

In [8]:

# we can access the properties for each individual particle, like the x coordinate
x = ms('x')
print(len(x))

In [9]:

# or do something more complicated such as:
pr = ms('sqrt(px**2 + py**2)')
# actually ridiculous things will work:
pr = ms('x + x**2 + (gamma**2 - 2) - sin(px/pz)')
# you can look at  for a list of predefined values
pp.particle_scalars

Out[9]:

Ekin = gamma_m1 * mass * c**2
Ekin_MeV = Ekin / elementary_charge / 1e6
Ekin_MeV_amu = Ekin / elementary_charge / 1e6 / mass_u
Ekin_MeV_qm = Ekin / elementary_charge / 1e6 / mass_u * charge_e
Ekin_keV = Ekin / elementary_charge / 1e3
Ekin_keV_amu = Ekin / elementary_charge / 1e3 / mass_u
Ekin_keV_qm = Ekin / elementary_charge / 1e3 / mass_u * charge_e
Eruhe = mass * c**2
_np2 = (px**2 + py**2 + pz**2)/(mass * c)**2
angle_xaxis = arctan2(sqrt(py**2 + pz**2), px)
angle_xy = arctan2(py, px)
angle_xz = arctan2(pz, px)
angle_yaxis = arctan2(sqrt(pz**2 + px**2), py)
angle_yx = arctan2(px, py)
angle_yz = arctan2(pz, py)
angle_zaxis = arctan2(sqrt(px**2 + py**2), pz)
angle_zx = arctan2(px, pz)
angle_zy = arctan2(py, pz)
beta = sqrt(gamma**2 - 1) / gamma
betax = vx / c
betay = vy / c
betaz = vz / c
charge = charge
charge_e = charge / elementary_charge
gamma = _np2 / (sqrt(1 + _np2) + 1) + 1
gamma_m = gamma * mass
gamma_m1 = _np2 / (sqrt(1 + _np2) + 1)
id = id
m = mass
m_u = mass / atomic_mass
mass = mass
mass_u = mass / atomic_mass
p = sqrt(px**2 + py**2 + pz**2)
px = px
py = pz
pz = py
q = charge
q_e = charge / elementary_charge
r_xy = sqrt(x**2 + y**2)
r_xyz = sqrt(x**2 + y**2 + z**2)
r_yz = sqrt(y**2 + z**2)
r_zx = sqrt(z**2 + x**2)
t = time
time = time
v = beta * c
vx = px / (gamma * mass)
vy = py / (gamma * mass)
vz = pz / (gamma * mass)
w = weight
weight = weight
x = x
x_um = x * 1e6
y = y
y_um = y * 1e6
z = z
z_um = y * 1e6
--> 56 known particle scalars.

In [10]:

# we can use the particle properties to create a Field for plotting.
# Particle Shapes as in the Simulation will be included
# calculate the number density
nd = ms.createField('x', 'y', bins=[240, 150], shape=3)
# note, that particle weights are included automatically
# and plot
fig, ax = subplots()
im = ax.imshow(nd.T/1e6, origin='lower', extent=nd.extent*1e6)
addcolorbar(ax, im)
ax.set_xlabel('x [$\mu m$]');
ax.set_ylabel('y [$\mu m$]');

In [11]:

# create a Field object qd holding the charge density
# note, that particle weights are included automatically
qd = ms.createField('x', 'y', weights='charge', bins=[240, 150], shape=3)
# and plot
fig, ax = subplots()
im = ax.imshow(qd.T, origin='lower', extent=qd.extent*1e6)
addcolorbar(ax, im)
ax.set_xlabel('x [$\mu m$]');
ax.set_ylabel('y [$\mu m$]');

In [12]:

# average kinetic energy on grid
ekin = ms.createField('x', 'y', weights='Ekin_MeV', bins=[240, 150])
ekinavg = ekin/nd
# and plot
fig, ax = subplots()
im = ax.imshow(ekinavg.T, origin='lower', extent=ekinavg.extent*1e6, vmax=20)
addcolorbar(ax, im)
ax.set_title('Avg kin. Energy on Grid [MeV]')
ax.set_xlabel('x [$\mu m$]');
ax.set_ylabel('y [$\mu m$]');

/glusterhome/stephan/.local/lib/python3.6/site-packages/postpic/datahandling.py:687: RuntimeWarning: invalid value encountered in true_divide
  result = getattr(ufunc, method)(*inputs, **kwargs)

In [13]:

f = ms.createField('x', 'p', bins=[240, 150])
# and plot
fig, ax = subplots()
im = ax.imshow(f.T, origin='lower', extent=f.extent, aspect='auto', vmax=0.5e41)
addcolorbar(ax, im)
ax.set_xlabel('x [m]');
ax.set_ylabel('p');

In [14]:

f = ms.createField('x', 'gamma', bins=[240, 150])
# and plot
fig, ax = subplots()
im = ax.imshow(f.T, origin='lower', extent=f.extent, aspect='auto', vmax=1e19)
addcolorbar(ax, im)
ax.set_xlabel('z [$\mu m$]');
ax.set_ylabel('$\gamma$');

In [15]:

f = ms.createField('x', 'beta', bins=[240, 150])
# and plot
fig, ax = subplots()
im = ax.imshow(f.T, origin='lower', extent=f.extent, aspect='auto', vmax=5e21)
addcolorbar(ax, im)
ax.set_xlabel('z [$\mu m$]');
ax.set_ylabel(r'$\beta$');

The Simulationreader¶

The Simulationreader represents an entire simulation, i.e. a series of dumps.

In [16]:

sr = pp.readSim('_epochexample/epoch_exampledata/normal.visit')
print('There are {:} dumps in this simulationreader object.'.format(len(sr)))
print(sr)

There are 5 dumps in this simulationreader object.
<Visitreader at "_epochexample/epoch_exampledata/normal.visit" (5 dumps)>

In [17]:

for dr in sr:
    print('Simulation time of current dump t = {:.2e} s'.format(dr.time()))

Simulation time of current dump t = 1.68e-16 s
Simulation time of current dump t = 1.25e-13 s
Simulation time of current dump t = 2.50e-13 s
Simulation time of current dump t = 3.75e-13 s
Simulation time of current dump t = 5.00e-13 s

In [18]:

fig, axs = subplots(1,5)
for dr, ax in zip(sr, axs):
    f = dr.Ex()
    im = ax.imshow(f.T, origin='lower', extent=f.extent*1e6)
    #addcolorbar(ax, im)
    ax.set_xlabel('x [$\mu m$]');
    ax.set_ylabel('y [$\mu m$]');
    ax.set_title('{:.2f} fs'.format(dr.time()*1e15))
fig.set_size_inches(13,3)
fig.tight_layout()

In [ ]: