benchmark about Collisions between particles of the same species.

[ercoop]

I set up an electron with a uniform density temperature distribution and used a Coulomb collision module. As the simulation evolved, there were differences in the temperature of the electron in three directions. I currently think that there is an problem with one line of code in BinaryCollision.H


Does anyone know where the problem lies？

[RemiLehe]

Thanks a lot for reporting this!
Are you saying that, when using cell_stop_1 instead of cell_half_1, the above-mentioned problem is solved? (i.e. electrons have the same temperature in all direction)

[ax3l]

Note that the line that @ercoop is referring to was changed locally. Checked in, the original line with cell_half_1 is used.
https://github.com/ECP-WarpX/WarpX/blob/23.12/Source/Particles/Collision/BinaryCollision/BinaryCollision.H#L348

Maybe @Yin-YinjianZhao can comment on this, it looks like he did the local modification?

[RemiLehe]

Btw, @ercoop could you share the input script that you used for this test?

Note also that we do have an automated test that checks for temperature isotropization:
https://github.com/ECP-WarpX/WarpX/blob/development/Examples/Tests/collision/inputs_3d_isotropization
https://github.com/ECP-WarpX/WarpX/blob/development/Examples/Tests/collision/analysis_collision_3d_isotropization.py

Is your test significantly different from that automated test?

[ercoop]

Yes! The problem of temperature anisotropy has been solved, but we found that the rate of temperature relaxation is slightly slower than our theoretical calculation, which is based on smiley. https://smileipic.github.io/Smilei/Understand/collisions.html

Here is my input script:
#################################
####### MY CONSTANTS ######
#################################
my_constants.n0 = 1.116e25 # number of electrons per m^3
my_constants.me = 1 * m_e
my_constants.mi = 1 * m_e
my_constants.de = clight / sqrt(n0q_e^2/epsilon0/me)
my_constants.di = clight / sqrt(n0q_e^2/epsilon0/mi)
my_constants.L0 = 32 * de
my_constants.tmax = 3.0e-9

#################################
####### GENERAL PARAMETERS ######
#################################
max_step = 144
stop_time = tmax
amr.n_cell = 16 16 16
amr.max_grid_size = 8
amr.blocking_factor = 8
#amr.max_grid_size_x = 32
#amr.blocking_factor_x = 16
#amr.max_grid_size_y = 16
#amr.blocking_factor_y = 16
amr.max_level = 0
#amr.restart = "check000005"
geometry.dims = 3
geometry.prob_lo = 0 0 0 # physical domain
geometry.prob_hi = L0 L0 L0
#warpx.fine_tag_lo = 0 0
#warpx.fine_tag_hi = L0 L0

#################################
####### Boundary condition ######
#################################
boundary.field_lo = periodic periodic periodic
boundary.field_hi = periodic periodic periodic
boundary.particle_lo = periodic periodic periodic
boundary.particle_hi = periodic periodic periodic

#################################
############ NUMERICS ###########
#################################
warpx.serialize_initial_conditions = 1
warpx.verbose = 1
warpx.cfl = 0.8
warpx.use_filter = 1
warpx.numprocs = 2 2 2

Order of particle shape factors

algo.particle_shape = 1
algo.maxwell_solver = yee

#################################
############ PLASMA #############
#################################
particles.species_names = electron

electron.charge = -q_e
electron.mass = m_e
electron.injection_style = "NRandomPerCell"
electron.num_particles_per_cell = 1280
electron.profile = constant
electron.density = n0 # number of electrons per m^3
#electron.momentum_distribution_type = "maxwell_boltzmann"
#electron.theta = 0.001
electron.momentum_distribution_type = "gaussian_parse_momentum_function"
electron.momentum_function_ux_m(x,y,z) = 0.0
electron.momentum_function_uy_m(x,y,z) = 0.0
electron.momentum_function_uz_m(x,y,z) = 0.0
electron.momentum_function_ux_th(x,y,z) = 0.0315918518549115
electron.momentum_function_uy_th(x,y,z) = 0.0315918518549115
electron.momentum_function_uz_th(x,y,z) = 0.0315918518549115
#electron.momentum_function_uz_th(x,y,z) = 0.00331633316355356
electron.do_bremsstrahlung_radiation_reaction = 0
electron.do_not_deposit = 1
electron.do_not_gather = 0
electron.do_not_push = 0

#################################
############ FIELDS #############
#################################

External magnetic field

#warpx.B_ext_grid_init_style = parse_B_ext_grid_file
#warpx.B_external_grid_file = turb_helical_512.h5
warpx.B_ext_grid_init_style = parse_B_ext_grid_function
warpx.Bx_external_grid_function(x,y,z) = 0
warpx.By_external_grid_function(x,y,z) = 0
warpx.Bz_external_grid_function(x,y,z) = 0

#################################
############ COLLISION ##########
#################################
collisions.collision_names = collision1
collision1.species = electron electron
collision1.CoulombLog = 50.0
collision1.ndt = 1

#################################
############ OUTPUT #############
#################################

Diagnostics

diagnostics.diags_names = diag1 diag2
diag1.diag_type = Full
diag1.format = plotfile
diag1.write_species = 0
diag1.intervals = 100
diag1.dump_last_timestep = 0
diag1.fields_to_plot = rho_electron

diag2.diag_type = Full
diag2.format = openpmd
diag2.write_species = 0
diag2.intervals = 1
#diag2.coarsening_ratio = 2 2 2
diag2.dump_last_timestep = 0
diag2.fields_to_plot = none
diag2.particle_fields_species = electron
diag2.particle_fields_to_plot = vlx vly vlz
diag2.particle_fields.vlx(x,y,z,ux,uy,uz) = uxux
diag2.particle_fields.vly(x,y,z,ux,uy,uz) = uyuy
diag2.particle_fields.vlz(x,y,z,ux,uy,uz) = uz*uz

check.diag_type = Full
check.format = checkpoint
check.intervals = 50:100:50
check.dump_last_timestep = 0

[pordyna]

Hi @ercoop. i stumpled upon this issue more or less by accident, but those plots defintely look familiar. So we did a simiilar test comparing PIConGPU and smilei and both codes were showing this slower relaxation. Idk which case has been shown in smilei docs and how up to date it is but here https://picongpu.readthedocs.io/en/latest/models/binary_collisions.html#thermalization you can see that the same thing happens in smilei.

[Yin-YinjianZhao]

Which version is that? I guess I was testing the shuffling, but when there is only one species, shuffling the half should be equivalent to shuffling them all.
But anyway I don't see why it would affect the isotropization.

[ercoop]

23.11 and 23.12 both have this problem. I did not carefully study the specific meanings of these two codes, cell_stop_1 and cell_half_1, but after making the modifications, the temperature will no longer exhibit non physical anisotropy. This is also the reason why I asked in the community to clarify why it works.

[JustinRayAngus]

Which version is that? I guess I was testing the shuffling, but when there is only one species, shuffling the half should be equivalent to shuffling them all. But anyway I don't see why it would affect the isotropization.

Shuffling half is not equivalent to shuffling them all. Each particle in a cell must have an equal opportunity of being paired with every other particle in the cell so that each particle can sample the full distribution over many scattering cycles. If you only shuffle half, then you are restricting each particle to only being able to sample half of the distribution function and thus it will not produce the correct scattering physics for anisotropic problems, such as the one reported here.

See PR #5045