Improvements on the adaptive time step calculation

[MaxThevenet]

This PR proposes 2 improvements of the adaptive time step algorithm:

1. prediction feature: In the current algorithm, the beam properties at time iteration it are used to predict the time step at the next iteration performed by the current rank, which is it+nranks. This can be a problem for largely parallel simulations, as the beam energy and local density can be quite off. Here, we use current Lorentz factor and accelerating field to estimate the beam position and energy at it+nranks, and use it to calculate the adaptive time step. This should make both the beam energy and local density evaluations more accurate. Addresses Improve adaptive time step #883.
2. phase feature: The temporal resolution to resolve the effect of density ramps depends on the local density and its gradient, not only the local density. At the beginning of a time step, the dt is corrected if needed to resolve the gradient. The new dt is chosen such that the resulting phase advance equals the default one in a uniform plasma (1/nt_per_betatron), to a certain tolerance level.

The figure below shows the physical time and time step for a z-dependent density using the dev version (old), the new code with both features above turned off (off, right on top of dev), only the phase correction on, only the prediction method on, both on, and both on with a 3x larger tolerance on the phase advance error (3x tol). As expected, the phase feature changes the behavior in the upramp. With prediction, the time step is modified in the density bumps, rather than after them.

This input file can be useful to test further improvements on the time step. Few python lines for the analysis are below.

import scipy.constants as scc
import numpy as np
import matplotlib.pyplot as plt
from openpmd_viewer.addons import LpaDiagnostics

ts0 = LpaDiagnostics('test_adaptive/diags0/hdf5/')
ts1 = LpaDiagnostics('test_adaptive/diags/hdf5/')

n0 = 2.e22
lramp = 0.005
z0 = 0.01
db = 0.0004
dn = 0.6*n0
def dens(z):
    return np.where(z>=0, np.where(z>=lramp,n0,n0*.1+(n0*.9)*z/lramp) ,n0*.1) - dn*np.exp(-(z-z0)**2/db**2)

plt.figure(figsize=(8,8))
plt.subplot(311)
plt.plot(scc.c*ts0.t[:-1], '+-', label='old')
plt.plot(scc.c*ts1.t[:-1], '+-', label='new')
plt.xlabel('time step')
plt.ylabel('z (m)')
plt.legend()
plt.subplot(312)
plt.plot(scc.c*ts0.t[:-1], scc.c*(ts0.t[1:]-ts0.t[:-1]), '+-')
plt.plot(scc.c*ts1.t[:-1], scc.c*(ts1.t[1:]-ts1.t[:-1]), '+-')
plt.xlim(0,.02)
plt.xlabel('z (m)')
plt.ylabel('c*dt (m)')
plt.subplot(313)
z = np.linspace(0,0.02,10000)
plt.plot(z, dens(z), 'k-')
plt.xlim(0,.02)
plt.xlabel('z (m)')
plt.ylabel('density (m-3)')
plt.tight_layout()

[MaxThevenet]

CI is fixed.

[SeverinDiederichs]

Looks great! See few minor comments

[MaxThevenet]

Thanks! I included your comments, and incorporated a few additional changes:

• Replace MPI Wait/Notify with Bcast
• Clean up behavior in the first few time steps

[SeverinDiederichs]

Thanks for the effort! Looking forward to the convergence scans and new defaults..