Picongpu: Jumps in output data of plugins (in moving window simulations)

[steindev]

I recognized in two different kind of simulations with two different plugins jumps in output data of longitudinal (y) quantities. I suspect these are related to calculations of output data within the moving window.

The two cases are:

Resetting of electron bunch longitudinal position in phase space output of LWFA sims

• picongpu version: 0.4.0rc2
• plugin: e_phaseSpace y-py

The electron position in y-direction (laser propagation and moving window direction) is repeatedly reset at a period of about twice the simulation box size in y-direction. Here, electron position resets take place every 3k to 4k (mostly) time steps. The simulation box is 2304 cells in y direction, with 8 gpus along this direction. The ratio of y cell size to time step is DY / ( cspeed * DT ) = 1.5275. That is, it seems like a reset occurs when the moving window moved once over the simulation box, i.e. every n_reset = 2304 * DY / ( cspeed * DT ) = 2304 * 1.5275 = 3519 time steps. This value fits to the observed reset period.

Increase in longitudinal E- and B-field energy in field energy plugin

• picongpu version: 0.3.0
• plugin: fields_energy

The longitudinal field energy jumps with a period of 20kDT. This time the simulation volume is 8192 cells in y-direction distributed over 16 gpus. The ratio of y cell size to time step is DY / ( cspeed * DT ) = 0.99976.

Zoom in to Ey data. Note that the sampling rate is high compared to period length of jumps.

As this second appearance of jumps does not really match the first one I will start a primitive bunch acceleration with piconpu-0.4.0 and see how the fields develop and if this behavior is reproduced.

Still it is funny that these jumps appear at a time the window needs to pass about twice the simbox.

[psychocoderHPC]

@steindev I checked the phase space plugin and it looks like sliding window is not supported from the plugin.

[ax3l]

Afair, the plugin dumps the whole box but also adds the meta information to slice out the window accordingly. Our Python plugin should do this reliably, so users don't need to fiddle with the HDF5 attributes and the cutting: https://picongpu.readthedocs.io/en/0.4.0-rc2/usage/plugins/phaseSpace.html#analysis

[PrometheusPi]

offline discussion with @psychocoderHPC:
The field energy plugin provides an output every 1k time step. The sliding occurs every 8192 (#cells) / 16 (# y-GPUs) = 512 ~0.5k time steps. Assuming the plugin wrote out the field energy just right before a GPU slide. The total field energy is thus at a (local) maximum. For the next plugin output, the energy is slightly lower (because 100ß%512 until the plugin output reaches a point, right before sliding, where the energy is at a local minimum. The next output contains a rapid increase in energy.

The effect seen is thus a aliasing effect due to the low sampling rate of the field energy plugin compared to the sliding rate.

A visualization of this effect is shown below.
The field energy is given by E(t) = 1.0 + sin(t/10000) + 0.5 t%512.
The energy thus fluctuates with the sliding every 512 time steps by 0.5 arb.u..

The result is a rapid increase and a slow decrease as seen by @steindev:

[PrometheusPi]

The effect can be resolved if the field energy plugin output occurs more frequently (at least every ~200 time steps (~512/2 = Nyquist sampling rate)) or can be avoided if the plugin only works on the moving window.

[PrometheusPi]

Or with an energy increase:

N = 100000
delta_plugin = 1000
delta_slide = 512

t = np.arange(N)

E = 1.0 + np.exp(t/50000.) + 0.5 * (t%delta_slide)/delta_slide

[ax3l]

Yep, many plugins need that filter still: #89 and it is now easy to implement with our filters: #2131 #2582 (note: not as a user option; I mean default-filtered)

Nevertheless, the phase space plugin has this implemented in simple post-processing. @steindev Can you please show the Python snippet used for PS visualization? I want to see if there is a bug hidden or if this is just not processing the data correctly.

[PrometheusPi]

@steindev or for comparison a "forward (green) and backward (orange) in time propagating" aliasing:

[steindev]

@ax3l: I used my own script for visualization of the phase space data:

tstep = 130000
ps = "ypy"

## Read data
f = h5py.File(path + "phaseSpace/PhaseSpace_e_all_%s_%i.h5"%(ps,tstep), "r")

# Unit data
Dr_SI = f["/data/%i/%s"%(tstep,ps)].attrs['dr_unit'] # position unit
Dp_SI = f["/data/%i/%s"%(tstep,ps)].attrs['p_unit'] / m_e / cspeed # momentum unit in gamma*beta

# Grid data 
p_min = f["/data/%i/%s"%(tstep,ps)].attrs['p_min'] * Dp_SI 
p_max = f["/data/%i/%s"%(tstep,ps)].attrs['p_max'] * Dp_SI 
N_r, N_p = f["/data/%i/%s"%(tstep,ps)].shape

# Data
rp = f["/data/%i/%s"%(tstep,ps)][...]

f.close()

## Plot data
figure(figsize=(8,6))
imshow(abs(rp).T, cmap = mpl.cm.cubehelix_r#, vmax = 2, vmin = 0
       , origin='lower', interpolation='none', aspect='auto'
       , extent=(0, N_r, p_min, p_max)
       , norm = matplotlib.colors.LogNorm())

# Plot setup
colorbar()
grid()
xticks(arange(0,N_r,200))
xlabel(r"position [cells]")
ylabel(r"momentum [$\gamma\beta_{komp}$]")
title(r"%s space @ %ikDt"%(ps,tstep*1.e-3))

show()

That is, I plot the complete data set stored in the phase space hdf5 file as is shown in the phaseSpace plugin docs. If there is a script for analysis provided by picongpu, I would wish there is a link or comment about it on this page. Further, the documentation page could list under known limitations the observed behavior and/or link to #89.

[ax3l]

Ah, that's the reason it jumps. Older versions (dev between 0.3.X and 0.4 - lol) of PIConGPU had an explanation on how to manually crop & visualize the phase-space output. Check this out: https://github.com/ComputationalRadiationPhysics/picongpu/blob/60506f84c0268ec8cd2fb92d1ce896fca02116c5/docs/source/usage/plugins/phaseSpace.rst#spatial-offset

But you should not do this anymore. Use the new python lib that you can add to your PYTHONPATH and which is in the latest manual. It will remove these nasty details from your sight :)

documentation page could list under known limitations the observed behavior

No, it now documents how to read the data. That's not a limitation for users: you are currently hacking the recommended workflow ;-) I don't recommend to anyone to parse HDF5 files by hand, it's verbose and error-prune.

[ax3l]

If you are in the mood to be extra-lazy (tm), which I can only recommend to be, then you could also replace the documented

from picongpu.plugins.phase_space import PhaseSpace

with an even lazier

from picongpu.plugins.plot_mpl.phase_space_visualizer import Visualizer

That object has the same arguments as the data reader API, but generates you immediately an MPL plot :D
https://github.com/ComputationalRadiationPhysics/picongpu/blob/2f265eb9715a4b4b4cbe82b4d09ab214773824e8/lib/python/picongpu/plugins/plot_mpl/phase_space_visualizer.py#L161-L164

I am recommending you the undocumented feature, since you seem to skip the recommended, documented ones ;-p

[steindev]

My apologies. I remember seeing the script which loads the python lib. I simply did not know where to find it... :-(
Simply loading the hdf5 data as I have done before seemed faster than searching for the lib. I think by doing so I already qualified for lazy (tm), at least.

[ax3l]

Nothing to be sorry for, we just moved all the plugin docs from the wiki to the mainline docs during the last year. They were hard to find. Everything is more shiny now! :sparkles:

But really, switch to the new Python imports. It's just one line in your picongpu.profile (or postprocessing.profile or whatever you use for jupyter) but will declutter your analysis notebooks so you can focus on the relevant things :)

[ax3l]

(Feel free to comment further if I missed something. The other issue is a duplicate of the energy plugin not yet implementing a moving window filter, cross-linked above.)