Difference between Probe and Field Diagnostics using BTIS3 scheme

[Tissot11]

I am using BTIS3 scheme to study the laser propagation in e-p plasmas. I always preferred in past using Probe because of the data size. However, I had to check the output of Field diagnostic and I noticed that both diagnostics give different answers for the laser field propagation in the plasma. I am using latest commits from GitHub for SMILEI. I paste below the Namelist where I only make change for ellipticity of the laser pulse. I also paste my plotting scripts for Probe and Field diagnostics and corresponding results. You can see for the CP laser pulse, difference is quite huge (factor of 3) and for the LP laser pulse, it is smaller. This is really puzzling...

import math, numpy as np

l0 = 2.0*math.pi              # laser wavelength
t0 = l0                       # optical cicle
Lx = 1024*l0
resx = 128.                    # nb of cells in one laser wavelength
dx = l0/resx                            # space step
dt  = 0.95 * dx                         # timestep (0.95 x CFL)
Tsim = 1024.0*t0                 # duration of the simulation
`

# Plasma density and profile
n0 = 1e-2                     # particle density

#Laser parameters
a0 = 10

start = 0                               # Laser start
fwhm = 5*t0                            # Gaussian time fwhm
duration = 10*t0                        # Laser duration
center = duration*0.5                   # Laser profile center

pusher = "borisBTIS3"

Main(
    geometry = "1Dcartesian",

    interpolation_order = 2 ,

    cell_length = [dx],
    grid_length  = [Lx],

    number_of_patches = [8*1024],

    timestep = dt,
    simulation_time = Tsim,

    EM_boundary_conditions = [ ['silver-muller', 'silver-muller'] ],

    use_BTIS3_interpolation = True,

)

LaserPlanar1D(
    box_side         = "xmin",
    a0               = a0,
    omega            = 1.,
    polarization_phi = 0.,
    ellipticity      = 0.,  #\pm 1 for CP and 0 for LP
    time_envelope    = tgaussian(start=start,duration=duration,fwhm=fwhm,center=center,order=2)
)

Species(
    # name = "electron_" + pusher,
    name = "electron",
    position_initialization = "regular",
    momentum_initialization = "cold",
    particles_per_cell = 128,
    #c_part_max = 1.0,
    mass = 1.0,
    charge = -1.0,
    number_density = trapezoidal(n0, xvacuum=32.0*l0, xplateau=1.0*Lx, xslope1=2.0*l0, xslope2=0.),
    # number_density = densProfile,
    #charge_density = n0,
    mean_velocity = [0., 0.0, 0.0],
    temperature = [0.],
    pusher = pusher,
    boundary_conditions = [["remove"]],
    #is_test = True
)

Species(
    # name = "electron_" + pusher,
    name = "positron",
    position_initialization = "regular",
    momentum_initialization = "cold",
    particles_per_cell = 128,
    #c_part_max = 1.0,
    mass = 1.0,
    charge = 1.0,
    number_density = trapezoidal(n0, xvacuum=32.0*l0, xplateau=1.0*Lx, xslope1=2.0*l0, xslope2=0.),
    #charge_density = n0,
    mean_velocity = [0., 0.0, 0.0],
    temperature = [0.],
    pusher = pusher,
    boundary_conditions = [["remove"]],
    #is_test = True
)

DiagFields(
    every = 100,
    fields = ['Rho_electron', 'Rho_positron', 'Ex', 'Ey','Ez','Bx','By','Bz']
)

DiagProbe(
    every = 100,
    origin = [0.0],
    corners = [ [Lx] ],
    number = [1024],
    fields = ['Rho_electron', 'Rho_positron', 'Ex', 'Ey','Ez','Bx','By','Bz', 'Jx','Jy','Jz']
)

Plotting scripts, first for Probe diagnostic

S = happi.Open(Path, verbose=True)

nCritical = 2*S.namelist.n0

aLaser = S.namelist.a0

# %% Plot laser field using Probe doagnostic

xGridProbe = S.Probe(0, 'Ey').getAxis('axis1')
tGridProbe = S.Probe(0, 'Ey').getTimes()
timeStepsGrid = S.Probe(0, 'Ey').getTimesteps()
dt = S.namelist.Main.timestep

timesteps = timeStepsGrid[1000]

timeStr = np.array2string(np.round(timesteps*dt/S.namelist.t0, 2))

# Renormalizations
timeGridNorm = tGridProbe/S.namelist.t0
xGridNorm = xGridProbe/S.namelist.l0

if laserPol == 'lp':
    DiagEy = S.Probe(0, 'Ey', timesteps=timesteps)
    ProbeDataEy = DiagEy.getData()[0].T
    Eperp = np.sqrt(ProbeDataEy**2)
    yTitle = r'$|E_y|$'

else:
    DiagEy = S.Probe(0, 'Ey', timesteps=timesteps)
    ProbeDataEy = DiagEy.getData()[0].T
    DiagEz = S.Probe(0, 'Ez', timesteps=timesteps)
    ProbeDataEz = DiagEz.getData()[0].T
    Eperp = np.sqrt(ProbeDataEy**2 + ProbeDataEz**2)
    yTitle = r'$E_{\perp}=\sqrt{E_y^2+E_z^2}$'

fig, ax = plt.subplots()

im = ax.plot(xGridNorm, Eperp)

# plt.colorbar(im)
ax.set(xlabel=r'$x/\lambda_L$', ylabel=yTitle,
       title=r'$t/\tau_L$ = ' + timeStr)

and then Field diagnostic

S = happi.Open(Path, verbose=True)

nCritical = 2*S.namelist.n0

aLaser = S.namelist.a0

#%% Plotting the results 

xGridFields = S.Field(0, 'Ey').getAxis('x')
tGridFields = S.Field(0, 'Ey').getTimes()

timeStepsGrid = S.Field(0, 'Ey').getTimesteps()
# Renormalizations
dt = S.namelist.Main.timestep

timesteps = timeStepsGrid[1000]

timeStr = np.array2string(np.round(timesteps*dt/S.namelist.t0,2))

timeGridNorm = tGridFields/S.namelist.t0
xGridNorm = xGridFields/S.namelist.l0
a0Laser = S.namelist.a0

if laserPol == 'lp':
    DiagEy = S.Field(0, 'Ey', timesteps=timesteps)        
    FieldDataEy =  DiagEy.getData()[0].T
    Eperp = np.sqrt(FieldDataEy**2)
    yTitle = r'$|E_y|$'

else:
    DiagEy = S.Field(0, 'Ey', timesteps=timesteps)        
    FieldDataEy =  DiagEy.getData()[0].T        
    DiagEz = S.Field(0, 'Ez', timesteps=timesteps)        
    FieldDataEz =  DiagEz.getData()[0].T        
    Eperp = np.sqrt( FieldDataEy**2 + FieldDataEz**2 )
    yTitle = r'$E_{\perp}=\sqrt{E_y^2+E_z^2}$'

fig1,ax1 = plt.subplots()

im1 = ax1.plot(xGridNorm, Eperp )

ax1.set(xlabel=r'$x/\lambda_L$', ylabel=yTitle,
       title= r'$t/\tau_L$ = ' + timeStr )

And the results for CP laser pulse (first Field and then Probe diagnostics)

and LP laser pulse (first Field and then Probe diagnostics)

[Z10Frank]

Hello, I know it is not efficient in a real use-case, but just to delve deeper in this behaviour: what happens if you use a number of probe points equal to the grid points? This way you should obtain the same number of points along x for the Fields and for the Probes. Do you observe the same difference?

Do you observe the same phenomenon without the BTIS3? Normally it does not affect directly the grid fields, but it can do it indirectly by acting on the macro-particle movement, which in turn affects the fields.

[Tissot11]

I can try to see if choosing same number of grid points for Probe helps. I did do simulations without BTIS3 but did not use Field diagnostic for those simulations. In fact, I always use Probe diagnostic and that's why I am a bit worried...I'll try next week to run a simulation without BTIS3 for the CP laser pulse and use same number of points for Probe diagnostic.

[Z10Frank]

Yes, my guess is that it is something related to the interpolation of the Probes. This would happen with and without B-TIS3.

PS in 1D you should be able to use a dt/dx even closer to 1, to make the B-TIS3 work even better, but I would repeat the tests keeping the dt fixed for the moment for coherence