SmileiPIC / Smilei

Particle-in-cell code for plasma simulation
https://smileipic.github.io/Smilei
334 stars 119 forks source link

Laser pulse time center is changing when changing the incidence_angle along y propagation #646

Closed Sanjeev273 closed 1 year ago

Sanjeev273 commented 1 year ago

@mccoys ## Description

I am running a 2D simulation with four counter Propagating laser starting from xmin, xmax, ymin, ymax. Two of laser at an incidence_angle of 1.5 degree ( one from ymin and one xmax). I have find the issue, when I am changing the incidence_angle, Pulse time centre is changing for the laser starting from the ymin or ymax boundary. See the attached video here. you can see in video, three lasers from xmin, xmax and ymax in beginging of video at correct time that was introduced by me in the script. but Laser from ymin boundary at an incidence_angle of 1.5 degree is starting after 2ps approx. In principle the time center souldn't change. The cell length and the grid length are same in both the x and y direction. I have tested it with the 1 particle per cell. the issue isame.

A script from the smilei benchmark ( tst2d_04_laser_wake.py) also tested with the 1.5 degree incidence_angle and 0 degree incidence_angle. I have find the same issue. I think it is bug.

https://github.com/SmileiPIC/Smilei/assets/67734147/b52cb4e2-6df5-4805-80b7-cd78a12abe3b


dx = 0.125
dt = 0.08
nx = 896
Lx = nx * dx
npatch_x = 128
laser_fwhm = 19.80

Main(
    geometry = "2Dcartesian",

    interpolation_order = 2,

    timestep = dt,
    simulation_time = int(2*Lx/dt)*dt,

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

    number_of_patches = [npatch_x, npatch_x],

    cluster_width = nx/npatch_x,

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

    solve_poisson = False,
    print_every = 100,

)

MovingWindow(
    time_start = Main.grid_length[0]*0.98,
    velocity_x = 0.9997
)

LoadBalancing(
    initial_balance = False,
    every = 20,
    cell_load = 1.,
    frozen_particle_load = 0.1
)

Species(
    name = "electron",
    position_initialization = "regular",
    momentum_initialization = "maxwell-juettner",
    particles_per_cell = 16,
    c_part_max = 1.0,
    mass = 1.0,
    charge = -1.0,
    charge_density = 0.000494,
    mean_velocity = [0.0, 0.0, 0.0],
    temperature = [0.000001],
    pusher = "boris",
    time_frozen = 0.0,
    boundary_conditions = [
        ["remove", "remove"],
        ["remove", "remove"],
    ],
)

LaserGaussian2D(
    box_side         = "ymin",
    a0              = 2.,
    focus           = [Main.grid_length[0]/2., 0.],
    waist           = 26.16,
    #incidence_angle  = 1.5*math.pi/180.,
    time_envelope   = tgaussian(center=2**0.5*laser_fwhm, fwhm=laser_fwhm)
)

Checkpoints(
    dump_step = 0,
    dump_minutes = 0.0,
    exit_after_dump = False,
)

list_fields = ['Ex','Ey','Rho','Jx']

DiagFields(
    every = 100,
    fields = list_fields
)

DiagProbe(
    every = 10,
    origin = [0., Main.grid_length[1]/2.],
    corners = [
        [Main.grid_length[0], Main.grid_length[1]/2.],
    ],
    number = [nx],
    fields = ['Ex','Ey','Rho','Jx', 'Rho_electron']
)

DiagScalar(
    every = 10,
    vars=[
        'Uelm','Ukin_electron',
        'ExMax','ExMaxCell','EyMax','EyMaxCell','RhoMin','RhoMinCell',
        'Ukin_bnd','Uelm_bnd','Ukin_out_mvw','Ukin_inj_mvw','Uelm_out_mvw','Uelm_inj_mvw'
    ]
)

DiagParticleBinning(
    deposited_quantity = "weight_charge",
    every = 50,
    species = ["electron"],
    axes = [
        ["moving_x", 0, Lx, 300],
        ["px", -1, 4., 100]
    ]
)

DiagPerformances(
    every = 100,
)

Steps to reproduce the problem

If relevant, provide a step-by-step guide

Parameters

Sanjeev273 commented 1 year ago

I run the test smulation. Temporal center is decreases with increase in incidence angle. with incidence_angle 1.5: temporal center is 2ps; incidence_angle 10: temporal center is 500fs; with incidence angle 30 degree: temporal center is 300fs. Presentation2

I have repeated the same process for the Laser starting from xmin boundary having same grid length, cell size and times. But a negligible chang in temporal center is observed. That is true in principle.

Hi @Z10Frank, Please, I rquest you to look into this matter. A job script is here:


import numpy as np
from numpy import vectorize, real
import cmath
import math
l0 = 2.0*math.pi         # laser wavelength [in code units]
t0 = l0                  # optical cycle
Lsimx = 42.*l0          # length of the simulation box in x
Lsimy = 42.*l0           # length of the simulation box in y
Tsim =  200.*t0            # duration of the simulation (= 28714.*2.67e-15 sec)
resx = 16.               # nb of cells in one laser wavelength along x-axis
resy = 16.
rest = resx/0.7         # nb of timesteps in one optical cycle 
dt = t0/rest             #0.19
nx = 42.*32.             # number of cells along x-axis
ny = 42.*32.
npatch_x = 16.         # number of patches along x-axis
npatch_y = 16.           # number of patches along y-axis
# Plasma density
n0 = 0.040121511

lambda0_laser     = 0.8e-6              # laser central frequency, m
c_over_omega0     = lambda0_laser/2./math.pi # converts from c/omega0 to um

# Laser parameters

######Pump Lasers parameter##############
######Spatial profile######
waistinI_pump =  12.5*l0       #(50 um) 
FWHMinI_pump  = waistinI_pump*(2.0*math.sqrt(math.log(2.0)))
FWHMinE_pump  = waistinI_pump*math.sqrt(2.0)
waistinE_pump = FWHMinI_pump/math.sqrt(2.0*math.log(2.0))
#sigma = FWHMinI/2.*math.sqrt(math.log(2.0))

#########Temporal profile##################
FWHMtinI_pump = 37.45*t0                  #(= 100 fs)
FWHMtinE_pump = FWHMtinI_pump*math.sqrt(2.0)
sigmatinI = FWHMtinI_pump/2.*math.sqrt(math.log(2.0)) 
sigmatinE = FWHMtinE_pump/2.*math.sqrt(math.log(2.0)) 
omega = 1.
a0 = 0.022        # Laser amplitude
center_laser= 3.*sigmatinE
focus_laser = 1.*FWHMinE_pump

############ Probe Laser parameters############

###########Spatial profile for probe ##########
waistinI_probe =  2.5*l0       #(10 um) 
FWHMinI_probe  = waistinI_probe*(2.0*math.sqrt(math.log(2.0)))
waistinE_probe = FWHMinI_probe/math.sqrt(2.0*math.log(2.0))
#sigma = FWHMinI/2.*math.sqrt(math.log(2.0))

###########Temporal profile for probe ##########
FWHMtinI_probe = 18.7266*t0                  #(= 50 fs)
FWHMtinE_probe = FWHMtinI_probe*math.sqrt(2.0)
sigmatinI_probe = FWHMtinI_probe/2.*math.sqrt(math.log(2.0)) 
sigmatinE_probe = FWHMtinE_probe/2.*math.sqrt(math.log(2.0)) 
#tau = FWHMtinI/2.*math.sqrt(math.log(2.0))
probe_delay = 10.*t0              #( 800 fs)
#center_probe= probe_delay + 3.*sigmatinE_probe

Main(
    geometry = "2Dcartesian",

    interpolation_order = 2,

    timestep = dt,
    simulation_time = Tsim, #int(1*Lx/dt)*dt,

    cell_length = [l0/resx, l0/resy],
    grid_length = [ Lsimx, Lsimy ],

    number_of_patches = [npatch_x, npatch_y],

    #clrw = nx/npatch_x,

    reference_angular_frequency_SI = 2.0*math.pi*3e8/0.8e-6,
    EM_boundary_conditions = [
        ["PML","PML"],
        ["PML","PML"],
    ],
    number_of_pml_cells =[[20, 20],[20, 20]], 
    solve_poisson = False,

    solve_relativistic_poisson = True,

    print_every = 100,

    random_seed = smilei_mpi_rank
)

#MovingWindow(
#    time_start =250.,                      # time_start = ( (0.9*LX / c0) + 5*tau_FWHM) / ref_t,     ref_t = 1/omega_L
#    velocity_x = 0.9999999
#)

LoadBalancing(
    initial_balance = False,
    every = 20,
    cell_load = 1.,
    frozen_particle_load = 0.1
)

### Transverse plasma-channel

Plasma_xmax = 20.*l0          # half length of plasma box in x direction
#bunch_sigma_r = 1.75*l0   
#delta_n= 0.25*n0
Radius_plasma = 20.*l0         # half length of plasma box in y direction 
def nplasma(x,y):
    longitudinal_profile = constant(n0) # gaussian(1., xvacuum=0., xlength =None, xfwhm=None, xcenter=None, xorder=2)  #constant(n0)
    profile_r = 0.
    if ((y-Main.grid_length[1]/2.)**2<Radius_plasma**2) and ((x-Main.grid_length[0]/2.)**2<Plasma_xmax**2):
        profile_r = 1.
    return profile_r*longitudinal_profile(x,y)

Te_keV = 0.005              # electron temperature in keV
Te  = Te_keV/511.           # Te normalised in mec^2 (code units)
Tp  = Te/10.

Species(
    name = "electron",
    position_initialization = "regular",
    momentum_initialization = "maxwell-juettner",
    particles_per_cell = 1,
    c_part_max = 1.0,
    mass = 1.0,
    atomic_number = 0., 
    charge = -1.0,
    charge_density = nplasma, #polygonal(xpoints=[focus_laser+20,focus_laser+30.,15000.,15100.],xvalues=[0.,n0,n0,0.]),
    mean_velocity = [0.0, 0.0, 0.0],
    temperature = [Te],
    pusher = "boris",
    time_frozen = 0.,
    boundary_conditions = [
       ["remove", "remove"],
       ["remove", "remove"],
    ],
)

Species(
    name = "proton",
    position_initialization = "regular",
    momentum_initialization = "maxwell-juettner",
    particles_per_cell = 1,
    c_part_max = 1.0,
    mass = 1836.,
    atomic_number = 0.,
    charge = 1.0,
    charge_density = nplasma, # polygonal(xpoints=[focus_laser+20,focus_laser+30.,15000.,15100.],xvalues=[0.,n0,n0,0.]),
    mean_velocity = [0.0, 0.0, 0.0],
    temperature = [Tp],
    pusher = "boris",
    time_frozen = 0.0,
    boundary_conditions = [
       ["remove", "remove"],
       ["remove", "remove"],
    ],
)

# We build a gaussian laser from scratch instead of using LaserGaussian2D
# The goal is to test the space_time_profile attribute
focus = [Main.grid_length[0]/2., Main.grid_length[1]/2.]

time_envelope = tgaussian(fwhm=FWHMtinE_pump, center=center_laser)   #(center=center_laser, fwhm=laser_fwhm)

ellipticity = 0
### Laser1 from ymin boundary #######
LaserGaussian2D(
    box_side         = "xmin",
    a0               = a0,
    omega            = 1.,
    focus            = focus,
    waist            = waistinE_pump,
    polarization_phi = 0.,  #90.*math.pi/180.,
    ellipticity      = 0.,
    time_envelope    = time_envelope
    )

###Laser2 from xmax boundary #######
LaserGaussian2D(
    box_side         = "xmax",
    a0               = a0,
    omega            = 1.,
    focus            = focus,
    waist            = waistinE_pump,
    incidence_angle  = 1.5*math.pi/180.,
    polarization_phi = 0.,  #90.*math.pi/180.,
    ellipticity      = 0.,
    time_envelope    = time_envelope
    )

### Laser3 from ymin boundary #######
LaserGaussian2D(
    box_side         = "ymin",
    a0               = a0,
    omega            = 1.,
    focus            = focus,
    waist            = waistinE_pump,
    incidence_angle  = 1.5*math.pi/180.,
    polarization_phi = 0.,  #90.*math.pi/180.,
    ellipticity      = 0.,
    time_envelope    = time_envelope
    )

###Laser4 from xmax boundary #######
LaserGaussian2D(
    box_side         = "ymax",
    a0               = a0,
    omega            = 1.,
    focus            = focus,
    waist            = waistinE_pump,
    polarization_phi = 0.,  #90.*math.pi/180.,
    ellipticity      = 0.,
    time_envelope    = time_envelope
    )
Checkpoints(
    dump_step = 0,
    dump_minutes = 0.0,
    exit_after_dump = False,
)

list_fields = ["Ex", "Ey", "Bx", "By", "Rho_electron", "Rho", "Rho_proton"] #all fields

DiagFields(
    every = 100.,
    fields = list_fields
)

DiagProbe(
        every = 100.,
        origin = [0., Main.grid_length[1]/2.],
        corners = [
            [Main.grid_length[0], Main.grid_length[1]/2.]
        ],
        number = [nx],
        fields = list_fields
)

DiagProbe(
        every = 100.,
        origin = [Main.grid_length[0]/2. - math.tan(30.*math.pi/180.)*Main.grid_length[1]/2., 0.],
        corners = [
            [Main.grid_length[0]/2.+ math.tan(30.*math.pi/180.)*Main.grid_length[1]/2., Main.grid_length[1]]
        ],
        number = [nx],
        fields = list_fields
)
DiagScalar(every = 10, vars=['Ukin','Ukin_electron'])

DiagParticleBinning(
    deposited_quantity = "weight",
    every = 10.,
  time_average = 1,
    species = ["electron", "proton"],
    axes = [
        ["x", 0., Main.grid_length[0], 500]
    ]
)

DiagPerformances (
         every = 100

)
mccoys commented 1 year ago

I don't think this really is a bug, as lasers do propagate correctly. It is expected that the delay changes with angle of incidence because lasers need more or less time to reach the center.

However, I agree that this large change vs. angle is not nice to handle, and the sharp change at angle=0 is not a good behavior. I will thus try to find a better way to define the laser delay depending on the angle.

mccoys commented 1 year ago

@MickaelGrechX @beck-llr @Z10Frank I have found a way to make the delay of LaserGaussian2D consistent from one boundary to the other, as requested in this issue. The result is that lasers with large incidence, or more specifically lasers coming from the ymin/ymax boundaries, will have a different delay after these changes in the code. I believe this is not a big problem for users as very few use this feature. For this reason I propose to integrate this change directly in develop without worrying about backwards compatibility. Is that ok with you guys?

beck-llr commented 1 year ago

Make sure that a correction is not applied just to a single geometry and not to the others. It is really a pain when behaviours are not the same in all geometries. Apart from that no problem for me.

mccoys commented 1 year ago

I did it for 2D and 3D. I don't think there is an angle of incidence in AM

Z10Frank commented 1 year ago

No problem for me. Indeed there is no angle of incidence in AM.

mccoys commented 1 year ago

This issue has been addressed in the develop branch. The delay now changes more smoothly with respect to the incidence angle.