LWFA simulation with PSATD solver produces asymmetry

[prkkumar]

LWFA simulation with PSATD solver produces asymmetry. The following input script results in

#################################
############# PSATD ##############
##################################
#warpx.numprocs = 1 1 12
algo.maxwell_solver = psatd
warpx.do_nodal = 1
psatd.use_default_v_galilean = 1
psatd.nox = 12
psatd.noy = 12
psatd.noz = 12
psatd.nx_guard = 12
psatd.ny_guard = 12
psatd.nz_guard = 12
particles.use_fdtd_nci_corr = 0
algo.current_deposition = direct
warpx.do_single_precision_comms = 0

amr.blocking_factor = 64
amr.max_grid_size = 64

#################################
########## MESH PATCH ###########
#################################
amr.max_level = 0

#################################
######### BOX PARAMETERS ########
#################################
my_constants.nx = 128
my_constants.nz = 3008
my_constants.xmax = 50e-6
my_constants.zmin = -150e-6
my_constants.zmax = 1e-6
warpx.zmax_plasma_to_compute_max_step = 0.05
amr.n_cell = nx nx nz
geometry.dims = 3
boundary.field_lo = periodic periodic damped
boundary.field_hi = periodic periodic damped
# physical domain when running in the lab frame
geometry.prob_lo = -xmax -xmax zmin
geometry.prob_hi = xmax   xmax zmax

#################################
########## DIAGNOSTICS ##########
#################################
diagnostics.diags_names = diag diag_lab
diag.intervals = 500
diag.format = plotfile
diag.fields_to_plot = Ex Ey Ez Bx By Bz jx jy jz rho
diag.diag_type = Full

diag_lab.diag_type = BackTransformed
diag_lab.do_back_transformed_fields = 1
diag_lab.num_snapshots_lab = 20
diag_lab.dz_snapshots_lab = 0.0025
diag_lab.fields_to_plot = Ex Ey Ez Bx By Bz jx jy jz rho
diag_lab.format = plotfile
diag_lab.diag_lo = -xmax -2*xmax/nx zmin
diag_lab.diag_hi = xmax   2*xmax/nx   0

#################################
############ NUMERICS ###########
#################################
warpx.verbose = 1
algo.charge_deposition = standard
algo.particle_pusher = vay
algo.particle_shape = 3
warpx.use_filter = 1
warpx.filter_npass_each_dir = 1 1 1
warpx.cfl = 1.

#################################
####### Moving Window ########
#################################
warpx.do_moving_window = 1
warpx.moving_window_dir = z
warpx.moving_window_v = 1.0

#################################
####### BOOST PARAMETERS ########
#################################
warpx.gamma_boost = 30.0
warpx.boost_direction = z

#################################
############ PLASMA #############
#################################
particles.species_names = ions1 electrons1 beam
particles.rigid_injected_species = beam

electrons1.charge = -q_e
electrons1.mass = m_e
electrons1.injection_style = NUniformPerCell
electrons1.num_particles_per_cell_each_dim = 1 1 1
electrons1.momentum_distribution_type = "gaussian"
electrons1.xmin = -46.e-6
electrons1.xmax =  46.e-6
electrons1.ymin = -46.e-6
electrons1.ymax =  46.e-6
electrons1.zmin = 0.0
electrons1.zmax = 1.
electrons1.profile                 = "predefined"
electrons1.predefined_profile_name = "parabolic_channel"
#         predefined_profile_params = z_start   ramp_up   plateau   ramp_down   rc       n0
electrons1.predefined_profile_params = 0.0 1.e-9 .05 1.e-9 20.e-6 3.4e23
electrons1.do_continuous_injection = 1
ions1.charge = q_e
ions1.mass = 1.6726219236900001e-29
ions1.injection_style = NUniformPerCell
ions1.num_particles_per_cell_each_dim = 1 1 1
ions1.momentum_distribution_type = "gaussian"
ions1.xmin = -46.e-6
ions1.xmax =  46.e-6
ions1.ymin = -46.e-6
ions1.ymax =  46.e-6
ions1.zmin = 0.0
ions1.zmax = 1.
ions1.profile                 = "predefined"
ions1.predefined_profile_name = "parabolic_channel"
#    predefined_profile_params = z_start   ramp_up   plateau   ramp_down   rc       n0
ions1.predefined_profile_params = 0.0 1.e-9 .05 1.e-9 20.e-6 3.4e23
ions1.do_continuous_injection = 1

beam.charge = -q_e
beam.mass = m_e
beam.injection_style = "gaussian_beam"
beam.x_rms = 2.0e-6
beam.y_rms = 2.0e-6
beam.z_rms = 5.0e-6
beam.x_m = 0.
beam.y_m = 0.
beam.z_m = -80.e-6
beam.npart = 50000
beam.q_tot = -3.5e-10
beam.momentum_distribution_type = "gaussian"
beam.ux_m = 0.
beam.uy_m = 0.
beam.uz_m = 1.e9
beam.ux_th = 0.
beam.uy_th = 0.
beam.uz_th = 0.
beam.zinject_plane = 0.0
beam.rigid_advance = true

#################################
############# LASER #############
#################################
lasers.names        = laser1
laser1.profile      = Gaussian
laser1.position = 0. 0. -1.e-9  # This point is on the laser plane
laser1.direction    = 0. 0. 1.      # The plane normal direction
laser1.polarization = 0. 1. 0.      # The main polarization vector
laser1.e_max        = 6.82274e12       # Maximum amplitude of the laser field (in V/m)
laser1.profile_waist = 20.e-6       # The waist of the laser (in meters)
laser1.profile_duration = 7.33841e-14   # The duration of the laser (in seconds)
laser1.profile_t_peak = 1.46764864e-13 # The time at which the laser reaches its peak (in seconds)
laser1.profile_focal_distance = 0.
laser1.wavelength = 0.8e-6         # The wavelength of the laser (in meters)
laser1.do_continuous_injection = 1

[EZoni]

Thank you for reporting this issue, Prabhat. In the past we observed asymmetries when some of the subdomains have odd numbers of cells (odd as opposed to even). Would you be able to double check if that is the case here, by looking at the grids summary in the standard output of the simulation?

[prkkumar]

I used even number of grids by specifying

amr.blocking_factor = 64
amr.max_grid_size = 64

in the input file. Grid summary from the output file reads:

Grids Summary:
  Level 0   188 grids  49283072 cells  100 % of domain
            smallest grid: 64 x 64 x 64  biggest grid: 64 x 64 x 64

[prkkumar]

Update : Using only one box in the transverse direction by specifying

amr.blocking_factor = 128
amr.max_grid_size = 128

I do not see the asymmetry:

[EZoni]

Ok, thank you. It would be interesting to see if we are using too few guard cells transversally. I can try to use the Python script in Tools/DevUtils/Stencil.py to get a rough estimate quickly. Could you write down for me the values of dx, dy, dz, and dt?

[prkkumar]

dt = 2.605969494e-15 ; dx = 7.8125e-07 ; dy = 7.8125e-07 ; dz = 2.948399295e-06

[EZoni]

Ok, thank you. One thing you could try is using 20 guard cells in x and y, instead of 12. This is one estimate I'm getting from the script. The asymmetry issue seems quite significant, so I'm not sure if this will make any difference, but maybe it is worth trying. This will work only provided that each subdomain is larger than 20 cells in x and y, otherwise the code will abort saying that the number of guard cells exceeds the number of valid cells.

[prkkumar]

I tried 20 guard cells instead of 12 and the results are perfectly symmetric - identical to the image I posted using single box in x. Thank you so much! I will familiarize myself with Python script in Tools/DevUtils/Stencil.py and use that to specify number of guard cells when using PSATD solver. This issue may be closed.