Zigzag structure of gamma-z phase space

[mulingLHY]

After simulating with hollow plasma, I observed a zigzag structure in the gamma-z phase space. I suspect this might be due to the internal push of particles without proper interpolation of Ez? Or should I enable any additional options to address this issue?

Here are the simulation parameters:

amr.n_cell = 255 255 255
amr.max_level = 0

# n0 = 2.83184e15 cm^-3
my_constants.kp_inv = 100e-6
my_constants.kp = 1. / kp_inv
my_constants.wp = clight * kp
my_constants.ne = wp^2 * m_e * epsilon0 / q_e^2

max_step = 500
hipace.dt = 4/wp
random_seed = 12345

hipace.verbose = 2

boundary.field = Dirichlet
boundary.particle = Absorbing
geometry.prob_lo     = -6.*kp_inv -6.*kp_inv -4.*kp_inv  # physical domain
geometry.prob_hi     =  6.*kp_inv  6.*kp_inv  4.*kp_inv


beams.names = beam
beam.injection_type = fixed_weight_pdf
beam.num_particles = 500000
beam.pdf = "(z<=300e-6)*(z>20e-6)*(1+z/3000e-6) + 1.1*(z>300e-6)*(z<=320e-6)*(320e-6-z)/20e-6 + (z>0)*(z<=20e-6)*z/20e-6"
beam.total_charge = 500e-12
beam.u_mean = 0. 0. 2935.42
beam.u_std = 0.02 0.02 0.2935
beam.position_mean = 0. 0.
beam.position_std = 50.e-6 50.e-6
beam.radius = 300.e-6


plasmas.names = plasma
# n = 6e16 cm^-3
plasma.density(x,y,z) = "r=sqrt(x*x+y*y); 21.19*ne* if(r>280e-6 and r<480e-6, 1., 6.-abs(r-380e-6)/20e-6)"
plasmas.min_density = 0
plasma.ppc = 2 2
plasma.u_mean = 0.0 0.0 0.0
plasma.element = electron


hipace.output_input = 1
diagnostic.diag_type = xyz
diagnostic.field_data = ExmBy By Ez rho_plasma
diagnostic.output_period = 4
hipace.file_prefix = diags/hollow
hipace.deposit_rho = 1

Could anyone provide insights into this issue or confirm if the lack of Ez interpolation could be the cause?

Thank you!

[mulingLHY]

Sorry, because I'm not familiar with GitHub, I submitted two useless issues before, and I don't know if I can delete them.

[mulingLHY]

Here is the post-processing code in Jupyter Notebook.

import h5py
import matplotlib.pyplot as plt
import numpy as np
from scipy.constants import c,e

step = 500
with h5py.File('./diags/hollow/openpmd_%06d.h5'%step, 'r') as f:
    # display_h5_struct(f)

    plt.figure(figsize=(11,5))
    plt.subplot(121)
    plasma = f['/data/%d/meshes/rho_plasma'%step][:]
    plt.imshow(plasma[:,:,plasma.shape[2]//2].T, extent=[-400,400,-600,600], aspect='auto')
    plt.xlabel('z (um)')
    plt.ylabel('x (um)')
    # plt.colorbar()

    plt.sca(plt.gca().twinx())

    Ez = f['/data/%d/meshes/Ez'%step][:]
    
    Ex = (f['/data/%d/meshes/ExmBy'%step][:] + f['/data/%d/meshes/By'%step][:])[100,Ez.shape[1]//2,:]
    Ez = Ez[:,Ez.shape[1]//2, Ez.shape[2]//2]
    plt.plot(np.linspace(-400,400,len(Ex)), Ez/1e9, linewidth=1.0)
    plt.ylabel('E (GV/m)')

    plt.subplot(122)
    z = f['/data/%d/particles/beam/position/z'%step][:]
    x = f['/data/%d/particles/beam/position/x'%step][:]
    px = f['/data/%d/particles/beam/momentum/x'%step][:]/c
    gamma = np.sqrt(f['/data/%d/particles/beam/momentum/z'%step][:]**2 \
                     + f['/data/%d/particles/beam/momentum/x'%step][:]**2 \
                        + f['/data/%d/particles/beam/momentum/y'%step][:]**2)/c
    weighting = f['/data/%d/particles/beam/weighting'%step][:]
    hist,edges = np.histogram(z, bins=200)

    plt.plot((edges[:-1]+edges[1:])/2*1e6, hist*weighting[0]*e/(edges[1]-edges[0])*c, linewidth=1.0, label='current')
    plt.ylabel('current (A)')
    plt.xlabel('z (um)')

    plt.sca(plt.gca().twinx())
    nslice = 400
    sort_idx = np.argsort(z)[:(z.shape[0]//nslice)*nslice]
    z, x, px, gamma = z[sort_idx].reshape(-1,nslice), x[sort_idx].reshape(-1,nslice), \
                        px[sort_idx].reshape(-1,nslice), \
                            gamma[sort_idx].reshape(-1,nslice)

    plt.plot(np.mean(z, axis=1)*1e6, np.mean(gamma, axis=1), linewidth=1.0, color='red', label='gamma')
    
    # plt.plot(np.mean(z, axis=1)*1e6, np.sqrt(np.mean(x**2, axis=1)*np.mean(px**2, axis=1)-np.mean(x*px, axis=1)**2)*1e6, linewidth=1.0, label='x')
    # plt.ylim(0, 2.0)

    cent:slice = slice(120,-120)
    print("beta", np.mean(x[cent]*x[cent])/(1e-6/np.mean(gamma[cent])))
    print("alpha", -np.mean(x[cent]*px[cent]/(1e-6)))
    
    plt.legend()

    plt.tight_layout()
    plt.show()

[MaxThevenet]

Hi @mulingLHY,
Thanks for opening an issue! Just to make sure I understand, is your question about the little steps than one can see on the red curve showing the Lorentz factor in your right plot? If so, that comes from the interpolation indeed: the order of interpolation in the longitudinal direction is 0 (it is 2 by default in the transverse directions x and y). Roughly speaking, all particles within the same longitudinal slice gather from this slice only, regardless of their longitudinal position within the cell. This typically gives rise to this behavior. In most cases, this does not really matter. If it matters for your case, you could consider increasing the number of cells longitudinally, from amr.n_cell = 255 255 255 to amr.n_cell = 255 255 500.

[mulingLHY]

Thanks for answering @MaxThevenet ,It indeed resolved my doubts.
Is this determined by hipace.depos_order_z and hipace.depos_order_xy in the input file? It seems that there is a plan to add depos_order_z based on the document.

hipace.depos_order_xy (int) optional (default 2)
Transverse particle shape order. Currently, 0,1,2,3 are implemented.

hipace.depos_order_z (int) optional (default 0)
Longitudinal particle shape order. Currently, only 0 is implemented.

[MaxThevenet]

Sounds good.
No real plan to go to higher order, no. Due to the algorithm, this would be a significant endeavor, and for most of our current applications we don't expect a significant gain, so really not high priority. The plots show these little steps, but as long as the beam cover many cells, which is basically always the case, the physics is well captured for our applications (plasma acceleration). Nevertheless, let us know if you think your application would benefit from this.

[mulingLHY]

@MaxThevenet Sorry for bothering again. Recently, after applying R56 to the beam from HiPACE++, I found that these little steps lead to micro-bunching in the beam.

This micro-bunching has a significant impact on the behavior of the beam. Therefore, I reviewed some of the source code and I'm puzzled as to why the configuration parameter depos_order_z is used for both field gather and current deposition.

hipace/src/Hipace.H

Lines 70 to 75 in 0ab9f68

	/** Order of the field gather and current deposition shape factor in the transverse directions
	*/
	inline static int m_depos_order_xy = 2;
	/** Order of the field gather and current deposition shape factor in the longitudinal direction
	*/
	inline static int m_depos_order_z = 0;

Simply adding linear interpolation for Ez during the field gather can smooth the output gamma-z phase space, but I am not sure if there are issues with physical self-consistency.

[MaxThevenet]

Dear @mulingLHY,
Indeed your example seems to be affected by the interpolation order in z. But as you can imagine, things are not that simple, otherwise we would have implemented z interpolation already. Long story short, we compute the plasma response in a loop over longitudinal cells, starting from the head of the beam, such that at any given point we only have the data of the current cell, hence the 0th order interpolation. More details at https://www.sciencedirect.com/science/article/pii/S0010465522001400. Workarounds are possible, but non-trivial, we will discuss it internally. For you, the simplest solution is to increase the longitudinal resolution, which should fix the issue. Could you do some convergence tests in longitudinal cell size, and let us know whether the problem can be fixed in this way?

[mulingLHY]

Thanks @MaxThevenet,
Increasing the longitudinal resolution will only lead to smaller micro-bunching, and smaller micro-bunching has a greater impact on the free-electron laser (FEL) which I am concerned with.

Regarding the acquisition of fields data from previous cell, by examining the source code, I discovered that the fields can be 'shifted' to WhichSlice::Previous in m_fields.ShiftSlices (within Hipace::SolveOneSlice), so that they can be read in the next slice. Ignoring the potential increase in memory usage for m_fields, I have made a simple linear interpolation for Ez, which cannot be pulled because of my unfamiliarity with C++ and the entire HiPACE++ project. I will organize my rough modifications for your review and correction.

[MaxThevenet]

Great to see that you're looking into the code!
You are correct, the information of the previous slice is stored, which could be used for a non-centred linear interpolation, which should be easier to implement. Whether this will improve the accuracy enough for your use cases is not clear, due to the order of interpolation, but certainly worth a try. If more accuracy is needed, there are also several option. Any chance you could use https://hipace.readthedocs.io/en/latest/run/chat.html for easier discussion? Otherwise we can simply proceed here.