NaN in output spectrum by using the trajectroy of FBPIC code

[QJohn2017]

Hi, I am using the output trajectories of the FBPIC code to calculate the radiation by the CLARA2 code, but the output spectrum has many NaN values. At first, I thought that this NaN appeared was due to the poor output accuracy of the particle trajectories in FBPIC code, So I Interpolate the particle trajectories to increase the accuracy (dt, dz, et al). But when I used the Interpolated particle trajectories to calculate the radiation spectrum, the NaN still appeared. So, does anyone know how this mistake came about? And how to solve this issue? Thank you very much!

[PrometheusPi]

Sorry for the late reply @QJohn2017, I am rarely checking ClaRa2 since I am not actively developing it. If you want to get a quick reply by addressing me directly via @PrometheusPi (but I will try to adjust my mail notification for this repo).

NaNs occur if:

• the floating point precision of the trajectories is not fine enough to
  • result in particle velocity v < speed of light c
  • the distance between two consecutive time steps is larger c * time step
• relativistic dynamics is not correctly handled (might be just rounding in the output for PIC codes) thus triggering the above errors
• the interpolation rate (required for FFT) is too low

• particles have a too high velocity so that $1 - \frac{v}{c}$ is rounded to zero.

  Depending on the reason, the NaNs occur only in specific directions or all directions.

  For what frequencies and directions do you observe NaNs? If you want, you can share a sample trajectories that produces NaNs with my and I will have a look.

[QJohn2017]

@PrometheusPi Thank you very much for replying to this message. Because the NaN always appears in calculations using this code, I have given up this issue. I will familiarize myself with the problem again in a few days and send you my examples, thanks again.

[PrometheusPi]

Thanks a lot for your quick reply @QJohn2017 and that you will look into it. And once again, I am very sorry for my late reply to your question.

[QJohn2017]

Hi, @PrometheusPi , one of my particle trace examples was in the attached file, which contains three particle trajectories. The order in each file was [x, y, z, betax, betay, betaz, t], and z is the propagation axis. These files are uniformly interpolated by time from particle tracks tracked by FBPIC code, in order to increase time accuracy. By using these trace files, a large amount of NaN appears in the output radiation results by the CLARA2 code. Please check these results in your code, thank you very much! trace5.zip

[PrometheusPi]

@QJohn2017 Thank you for providing the trajectories. I will have a look at it today.

[PrometheusPi]

@QJohn2017 I had a lock at the trajectories you provided. For example trace_17.txt describes an particle moving faster than the speed of light, for some time steps.

Below you see the squared normalized velocity beta^2 plotted over time.

I generated this analysis plot using the following python code:

import numpy as np
import matplotlib.pyplot as plt

# load data
data_17 = np.loadtxt("trace_17.dat")

# compute beta^2
beta_squared = data_17[:, 3]**2 + data_17[:, 4]**2 + data_17[:, 5]**2
# check where beta^2 ≥ 1
wrong_filter = np.greater_equal(beta_squared, 1.0)

# plot 
plt.figure(facecolor="white")
plt.title("trace_17.dat", fontsize=18)

# entire trace
plt.plot(data_17[:, -1]/1e-12, beta_squared, ) 

# only faster than c part
plt.plot(data_17[wrong_filter, -1]/1e-12, beta_squared[wrong_filter], 
         "o", color="red", ms=1, alpha=0.5, label=r"$\vec{\beta}^2 \geq 1$")

# labels etc.
plt.xlabel(r"$t \, \mathrm{[ps]}$", fontsize=18)
plt.xticks(fontsize=14)
plt.ylabel(r"$\vec{\beta}^2$", fontsize=18)
plt.yticks(fontsize=14)
plt.ylim(0.997, 1.0005)
plt.xlim(150, 205)
plt.legend(fontsize=18, loc=4)

plt.tight_layout()
plt.show()

Based on this finding, it is very likely that you divide by zero when beta^2 ≥ 1 at some point and your observer vec{n} is close to parallel to the velocity (pointing in forward direction) thus that the denominator of (eq. 2.10 in or see reduced equation below) becomes zero and produces a NaN.

The faster-than-light particles in your trace might not necessarily be a problem with FBPIC but might originate from how the interpolation is done to make the traces equidistant in time. When doing a linear interpolation, physics constrains might not necessarily be fulfilled.

I hope this helped you.

[PrometheusPi]

PS: I would be interested whether the faster than light particles are comping (as I suspect) from the interpolation or from FBPIC. In the latter case not using boosted frame might solve the problem - (if you are using boosted frame simulations with FBPIC).

If it is FBPIC, please inform @RemiLehe via an issue in the FBPIC repo.

[PrometheusPi]

PPS: I correct myself. Looking at trace_50.dat, there are long stretches of faster than light velocity.

I doubt that this is just an interpolation error when preparing the FBPIC data for ClaRa2.

Highly interesting problem on the trajectory side.

Note to myself: adding a beta^2 ≥ 1 in ClaRa2 might allow to debug faulty input faster in the future.

[QJohn2017]

@PrometheusPi , Thank you very much for showing my issue, it is true that I used the boosted frame in FBPIC simulations, and the original particle trajectories were traced in boosted frame simulation. Because the simulation length is very long (~6 cm). To save calculation time, only boosted frame can be used. After obtaining the particle trace files, I used the time uniform interpolation method (e.g. dt is the same) to interpolate the trajectory, because the FBPIC could not track the particles in a very small time interval (e.g. the computational accuracy in FBPIC, or I don't know if I didn't find the right way to use particle tracking method). The original particle trajectory has only just over 100 data points. So, for this situation of mine, is there any way to solve it? Thank you again!

[PrometheusPi]

@QJohn2017 You are welcome. Sorry for just showing the origin of the problem and not having a solution.

I think for now, the only way to do your analysis is to either:

• exclude the entire trajectory from analysis if it contains a beta^2≥1 point

• or split the beta^2≥1 trajectories into several shorter trajectories that do not contain wrong points (with some extra correct time steps cut out around the beta^2≥1 regions to be save)

  The first is the safest option. The latter might provide some more details. But since quite often the relevant radiation for e.g. light sources comes from the fastest electrons, excluding either entire trajectories or only the "fast" parts will alter the spectra in the high frequency region.

  The only "good" approach would be to debug the beta^2≥1 problem.

  Was the interpolation ("time uniform interpolation method") done by FBPIC or code provided by FBPIC?

[QJohn2017]

@PrometheusPi Thank you for giving me the solution suggestions. I have a question for the second method you provided, if I do some extra correct time steps cut out around the beta^2≥1 regions, the time step should be different from the "good" trajectory part, is this allowed by CLARA2 code (that is the time step is different in the whole particle trajectory file). secondly, the interpolation ("time uniform interpolation method") was done by myself in the MATLAB software through the interp1 function : t_fit = linspce(min(t_old), max(t_old), 50000); x_fit = interp1(t, x, t_fit, 'pchip'); y_fit = interp1(t, y, t_fit, 'pchip'); z_fit = interp1(t, z, t_fit, 'pchip'); ux_fit = interp1(t, ux, t_fit, 'pchip'); uy_fit = interp1(t, uy, t_fit, 'pchip'); uz_fit = interp1(t, uz, t_fit, 'pchip');

[PrometheusPi]

@QJohn2017 Sorry for my imprecise wording. With "cut out around the beta^2≥1 regions" I meant:

• take your wrong trajectory trace_1.dat
• for each (large enough) region i of N correct data regions in the trajectory:
  • save only the correct data as trace_1_i.dat
• remove trace_1.dat from analysis with ClaRa2 and take instead trace_1_i.dat till trace_1_N.dat (I used Fortran indexing here.)

Thanks for the information on the interpolation method. Can you check whether FBPIC had already some ux^2 + uy^2 + uz^2 ≥ 1 data?

[QJohn2017]

@PrometheusPi The following figure is the original trajectory data of trace_17: The second figure above is the z-x trajectories of trace_17. It can be seen that by the interpolation, the trajectory fits well, but the velocity figure is different compared to what you plotted above. Besides, in the original output data, there is no beta^2≥1 appeared. This shows that this NaN originated from the interpolation error. Finally, the attached file was the original output data of the FBPIC code, which corresponding to the interpolated data that I provided above. trace6.zip ,

[PrometheusPi]

Thanks a lot for clarifying this.

Thus I would correct my recommendation: I would suggest correcting/improving the interpolation method. A straight forward way to do this could be to convert beta_x, beta_y and beta_z in momenta values. Interpolate on those values and then convert back to normalized velocities. That should prevent beta^2≥1 - but it might come at the cost of numerical inaccuracies when doing the back and forth sqrt for the Lorentz factor.

[QJohn2017]

@PrometheusPi Thank you for your patience in explaining my question to me. Do you mean that, if there is 5 correct data regions in 'trace_1.dat', then I save only the correct data as trace_1_0.dat to trace_1_4.dat, and take these 5 files into the CLARA2 code to calculate the radiation?

[PrometheusPi]

@PrometheusPi Thank you for your patience in explaining my question to me. Do you mean that, if there is 5 correct data regions in 'trace_1.dat', then I save only the correct data as trace_1_0.dat to trace_1_4.dat, and take these 5 files into the CLARA2 code to calculate the radiation?

Yes, exactly. Thus you will have 4 smaller traces with non-overlapping times. But they ae equidistant in time. (There will be some artifacts from stopping the Fourier transform at high acceleration values though - but no NaNs)

[QJohn2017]

@PrometheusPi , Okay, thank you very much! The interpolation method is very good, I've been ignoring this before. In fact, only position and momenta values are the correct parameters, not the normalized value. Thank you again, and I will try the interpolation again by your method. Thank you very much for helping me so much!

[PrometheusPi]

@QJohn2017 You are welcome.

[PrometheusPi]

@QJohn2017 Did implementing a momentum-based interpolation help?

[QJohn2017]

@PrometheusPi, sorry for the late reply, I have done a momentum-based interpolation and the NaN disappears, thank you very much! But the output radiation results showed only ~30eV radiated photon energy, which seems to be inconsistent with the theoretical prediction.

According to the Betatron radiation, the critical radiation energy is about ~60 keV, so the simulation results is far lower than the analysis result. The attached file was the momentum-based interpolated trajectories. trace7.zip