Simple LUME-Genesis example for Genesis 1.3 Version 2 (Genesis2)¶

In [ ]:

# Useful for debugging
%load_ext autoreload
%autoreload 2

# Useful for debugging %load_ext autoreload %autoreload 2

In [ ]:

from genesis import Genesis2

# Make genesis object with some input file
G = Genesis2("data/basic/genesis.in", verbose=True)

from genesis import Genesis2 # Make genesis object with some input file G = Genesis2("data/basic/genesis.in", verbose=True)

In [ ]:

# Turn on field output
G["idmpfld"] = 1
# Turn on particle output
G["idmppar"] = 1
G["npart"] = 2048

# Turn on history
G["ippart"] = 10
G["ipradi"] = 0

# Turn on field output G["idmpfld"] = 1 # Turn on particle output G["idmppar"] = 1 G["npart"] = 2048 # Turn on history G["ippart"] = 10 G["ipradi"] = 0

In [ ]:

# Run genesis with default lattice
G.run()

# Run genesis with default lattice G.run()

In [ ]:

G.output["data"].keys()

G.output["data"].keys()

In [ ]:

G.write_wavefront()

G.write_wavefront()

In [ ]:

G.archive()

G.archive()

In [ ]:

# This is the working path
G.path

# This is the working path G.path

In [ ]:

# Contents of this
!ls -ahl {G.path}

# Contents of this !ls -ahl {G.path}

In [ ]:

# Get a hash of the input
G.fingerprint()

# Get a hash of the input G.fingerprint()

In [ ]:

# Lattice parameters
G.lattice["param"]

# Lattice parameters G.lattice["param"]

In [ ]:

# Get a list of z from the output
G.output.keys()

# Get a list of z from the output G.output.keys()

In [ ]:

# These are the available data
G.output["data"].keys()

# These are the available data G.output["data"].keys()

In [ ]:

zlist = G.output["data"]["z"]
zlist.shape

zlist = G.output["data"]["z"] zlist.shape

In [ ]:

# Get power. This is a 2d array of: slice, z
power = G.output["data"]["power"]
power.shape

# Get power. This is a 2d array of: slice, z power = G.output["data"]["power"] power.shape

In [ ]:

# Simpler plot

import matplotlib.pyplot as plt

%config InlineBackend.figure_format = 'retina'

plt.xlabel("Z (m)")
plt.ylabel("Power (GW)")

plt.plot(zlist, power[0] / 1e9)

# Simpler plot import matplotlib.pyplot as plt %config InlineBackend.figure_format = 'retina' plt.xlabel("Z (m)") plt.ylabel("Power (GW)") plt.plot(zlist, power[0] / 1e9)

Archiving¶

In [ ]:

afile = G.archive()

afile = G.archive()

In [ ]:

G2 = Genesis2()
G2.load_archive(afile)

G2 = Genesis2() G2.load_archive(afile)

In [ ]:

# Check that fingetprints (hash on input) are the same
G.fingerprint() == G2.fingerprint()

# Check that fingetprints (hash on input) are the same G.fingerprint() == G2.fingerprint()

In [ ]:

import numpy as np


# Check that all output data are the same
for k in G.output["data"]:
    print(k, np.all(G.output["data"][k] == G2.output["data"][k]))

import numpy as np # Check that all output data are the same for k in G.output["data"]: print(k, np.all(G.output["data"][k] == G2.output["data"][k]))

Wavefront (dfl) in openPMD-wavefront format¶

This will write the loaded dfl data to a proper openPMD-wavefront file

In [ ]:

wfile = G.write_wavefront()

wfile = G.write_wavefront()

In [ ]:

# Read back
import h5py

with h5py.File(wfile, "r") as h5:
    component = h5["data/000000/meshes/electricField/x"]
    E = component[:]
    attrs = dict(component.attrs)

# This is the conversion factor to the proper SI unit for electric field
attrs

# Read back import h5py with h5py.File(wfile, "r") as h5: component = h5["data/000000/meshes/electricField/x"] E = component[:] attrs = dict(component.attrs) # This is the conversion factor to the proper SI unit for electric field attrs

In [ ]:

# Check the shapes
E.shape

# Check the shapes E.shape

In [ ]:

# Plot the phase
%matplotlib inline
import matplotlib.pyplot as plt

plt.rcParams["figure.figsize"] = [8, 8]


ndat = np.angle(E[:, :, 0])  # Complex angle
plt.imshow(ndat, extent=[1e3 * G["dgrid"] * i for i in [-1, 1, -1, 1]])
plt.xlabel("x (mm)")
plt.ylabel("y (mm)")
plt.show()

# Plot the phase %matplotlib inline import matplotlib.pyplot as plt plt.rcParams["figure.figsize"] = [8, 8] ndat = np.angle(E[:, :, 0]) # Complex angle plt.imshow(ndat, extent=[1e3 * G["dgrid"] * i for i in [-1, 1, -1, 1]]) plt.xlabel("x (mm)") plt.ylabel("y (mm)") plt.show()

Raw particle history¶

In [ ]:

import numpy as np
import os

import numpy as np import os

In [ ]:

# gamma, phase, x, y, px/mc, py/mc
b = G.output["data"]["par"]

# gamma, phase, x, y, px/mc, py/mc b = G.output["data"]["par"]

In [ ]:

p = b[-1]
plt.hist2d(
    p[1] % (2 * np.pi), 0.511 * p[0], bins=[100, 100], cmap=plt.get_cmap("plasma")
)
plt.xlabel("phase")
plt.ylabel("energy (MeV)")
# plt.savefig('frame.png')
plt.show()

p = b[-1] plt.hist2d( p[1] % (2 * np.pi), 0.511 * p[0], bins=[100, 100], cmap=plt.get_cmap("plasma") ) plt.xlabel("phase") plt.ylabel("energy (MeV)") # plt.savefig('frame.png') plt.show()

In [ ]:

def frame(
    i,
):
    p = b[i]
    plt.hist2d(
        p[1] % (2 * np.pi), 0.511 * p[0], bins=[200, 200], cmap=plt.get_cmap("plasma")
    )
    plt.xlabel("phase")
    plt.ylabel("energy (MeV)")


# plt.savefig('frame_'+str(i)+'.png')

frame(100)

def frame( i, ): p = b[i] plt.hist2d( p[1] % (2 * np.pi), 0.511 * p[0], bins=[200, 200], cmap=plt.get_cmap("plasma") ) plt.xlabel("phase") plt.ylabel("energy (MeV)") # plt.savefig('frame_'+str(i)+'.png') frame(100)

In [ ]:

# Make frames
# for i in range(nbunch):
#    frame(i);

# Make frames # for i in range(nbunch): # frame(i);

In [ ]:

# Make movie
#!ffmpeg -framerate 10 -i frame_%d.png -c:v libx264 -c:a libfdk_aac output.mp4

# Make movie #!ffmpeg -framerate 10 -i frame_%d.png -c:v libx264 -c:a libfdk_aac output.mp4

Final particles¶

If a final particle (.dpa file) data was made, this will process the internal arrays into a high-level ParticleGroup object.

In [ ]:

P = G.final_particles()
P

P = G.final_particles() P

In [ ]:

P.plot("t", "energy")

P.plot("t", "energy")

Cleanup¶

In [ ]:

os.remove(afile)
os.remove(wfile)

os.remove(afile) os.remove(wfile)