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

[Tissot11]

I can confirm that if I choose as many points as in the Field diagnostic, the laser field (for circular polarization) between Probe and Field looks identical. So the problem is indeed with interpolation in Probe diagnostic. I hope this does not pose a problem for simulations without the BTIS-3 scheme. Please let me know what is your finding.

[Z10Frank]

Thank you, to know this you can repeat the test with and without the B-TIS3 scheme, although I doubt that the result will change.

Remember that the B field in Fields is defined with a dt/2 shift in time, while the B you see in the Probes are those interpolated for the macro-particle pusher, so it is the field actually acting on the macro-particles.
You can also add e.g. BzBTIS3 in the Probe diagnostic and Bz_mBTIS3 in the Field diagnostic.

Probes are not just the field at a certain grid point, they interpolate the field as if they were macro-particles, thus using a shape function with finite extent which is larger than one grid point. For this reason they may damp high-frequency oscillations seen in the Field diagnostic.

[Tissot11]

Ok. Indeed, Probe diagnostic does damp high-frequency components. But the results should not vary that much between the two. For me the biggest reason to use Probe to have light file size and Probe diagnostic used to quite fast compared to Field. Does it mean I should abandon the Probe diagnostic?

[Z10Frank]

If there are high frequencies in the Field, the interpolation made by Probes will inevitably reduce the oscillation amplitudes. If the high-frequency components are significant, I am not surprised that a low-pass filter like an interpolation yields significantly different results.

To reduce the file size, is there a reason for not dumping a subset of the grid in the Field diagnostic?
https://smileipic.github.io/Smilei/Use/namelist.html#subgrid
This way you would have exactly the fields on the grid (without interpolation of the Probes), but if your sampling step is larger than the period of those high-frequency, high-amplitude oscillations then you will not see them in the outputs.

If you think that these oscillations are only numerical noise, you may try to further increase the number of macro-particles.

[Tissot11]

Field diagnostic used to be really slow in past. That's why I started using Probes. I guess now the Field diagnostic is also fast enough that I could use it instead of the Probe, with this subset feature.

I mean these high-frequency oscillations are only important for a certain class of laser-plasma interactions. As long as the diagnostic is able to resolve laser frequency/oscillations, it is fine. In my case, I had an extremely high resolution lambda_L/128 and in Probe diagnostic, I was essentially saving lambda_L/8, which is good enough for resolving laser oscillations and that's why I find this difference a bit troubling. I also had 128 particles per cell. It is not only the laser field, plasma density oscillations also show markedly different results on plasma density compression between two diagnostics but I could ignore this difference due to number of points chosen for different diagnostics and this difference was not as stark as for the laser field amplitude for circular polarisation.

[Z10Frank]

I'll close the issue then, feel free to reopen it or to open a new one if you need.

[Tissot11]

So you suggest that despite high resolution chosen for the Probe diagnostic, a factor of 3 lower laser amplitude compared to the Field diagnostic is justified? If this is true then Probe diagnostic is not useful for the purpose it was created in SMILEI. Also if it is indeed an interpolation issue then what is your confidence level for using the subset feature in the Field diagnostic? In that case too, there is an interpolation involved.

[Z10Frank]

The behaviour of Probes is to sample fields as if they were macro-particles at the requested positions, i.e. with an interpolation which uses a finite extent shape function. By mathematical definition an interpolation of a highly varying field can yield a value significantly different from the field at one point, so I don't see a contradiction mathematically, independently of what could have been their initial purpose.

All diagnostics have a defined behaviour, which may or may not be suited for a given case. If you think that it may not be suited for what you want for this particular case, you are free to choose to use them or not.
You may try increasing the number of macro-particles, maybe you need more than 128 if that is a numerical noise.

The subgrid feature of the Fields does not use an interpolation as Probes. I cannot give an ad-hoc confidence level, you are the best placed to test it on your specific case.