Warwick-Plasma / epoch

Particle-in-cell code for plasma physics simulations
https://epochpic.github.io
GNU General Public License v3.0
186 stars 59 forks source link

Injectors cannot inject photons #502

Closed Status-Mirror closed 1 year ago

Status-Mirror commented 1 year ago

When using standard injectors to inject photons (not file injectors), the positions of photons are set to Inf, and we output the terminal message: "Note: The following floating-point exceptions are signalling: IEEE_DIVIDE_BY_ZERO".

This error was achieved using a code compiled with -DPHOTONS, and the input deck:

begin:control
  nx = 600
  t_end = 10 * femto
  x_min = -20 * micron
  x_max = 100 * micron
end:control

begin:boundaries
  bc_x_min = periodic
  bc_x_max = periodic
end:boundaries

begin:species
  name = Photon
  npart_per_cell = 0
  number_density = if((x gt 40*micron)and(x lt 1000*micron),1e30,0.0)
  identify:photon
end:species

begin:qed
  photon_dynamics = T
end:qed

begin:injector
  boundary = x_min
  species = Photon
  number_density = 1.0e20
  temp_x = 0
  drift_x = 1.0e-20
  npart_per_cell = 10
end:injector

begin:output
  dt_snapshot = 1.0 * femto
  particle_grid = always
end:output
Status-Mirror commented 1 year ago

Problems have also been reported using another input deck, in EPOCH1D. The user introduces 3 species: Photon (identify:photon), TestPhoton (mass = 1.0e-41kg), and TestElectron (identify:electron). All are modelled as two bunches - the first bunch is initialised inside the simulation window by x_min with drift, and the second comes from x_min injectors at a later time. TestElectron and TestPhoton behave as expected, while both Photon bunches fail:

(Note: the following input deck has been edited for brevity - some default and unused paramters have been removed).

begin:control
  nx = 550
  npart = 1e4
  nsteps = -1
  t_end = 70e-15
  x_min = -1e-6
  x_max = 10e-6
  dt_multiplier = 0.8
end:control

begin:qed
  use_qed = T
  qed_start_time = 0
  produce_photons = T
  photon_energy_min = 0
  produce_pairs = F
  photon_dynamics = T
end:qed

begin:boundaries
  bc_x_min = simple_outflow
  bc_x_max = simple_laser
end:boundaries

begin:species
  name = Photon
  dump = T
  mass = 0
  charge = 0
  npart = 1e4
  density = if((x gt -1e-6) and (x lt 0e-6),1e4*1e6,0)
  drift_x = 1e9 * qe/c
  identify:photon
end:species

begin:injector
  boundary = x_min
  species = Photon
  t_start = 33.33e-15
  t_end = 36.67e-15
  npart_per_cell = 1e4/50
  density = 1e4*1e6
  drift_x = 1e9*qe/c
  temp = 0
end:injector

begin:species
  name = TestPhoton
  dump = T
  mass = 1e-10
  charge = 0
  npart = 1e4
  density = if((x gt -1e-6) and (x lt 0e-6),1e4*1e6,0)
  drift_x = 1e9 * qe/c
end:species

begin:injector
  boundary = x_min
  species = TestPhoton
  t_start = 33.333e-15
  t_end = 36.666e-15
  npart_per_cell = 1e4/50
  density = 1e4*1e6
  drift_x = 1e9*qe/c
  temp = 0
end:injector

begin:species
  name = TestElectron
  dump = T
  mass = 1
  charge = -1
  npart = 1e4
  density = if((x gt -1e-6) and (x lt 0e-6),1e4*1e6,0)
  drift_x = 1e9 * qe/c
  identify:electron
end:species

begin:injector
  boundary = x_min
  species = TestElectron
  t_start = 33.333e-15
  t_end = 36.666e-15
  npart_per_cell = 1e4/50
  density = 1e4*1e6
  drift_x = 1e9*qe/c
  temp = 0
end:injector

begin:output
  use_offset_grid = F
  dt_snapshot = 1.0e-15
  full_dump_every = 1
  restart_dump_every = -1
  force_final_to_be_restartable = F
  particles = full
  px = full
  py = full
  pz = full
  particle_weight = full
  grid = always
  ey = always
  ez = always
  ekbar = always + species
  number_density = always + species
  distribution_functions = always
end:output

begin:dist_fn
  name = momentum
  ndims = 3
  dumpmask = always
  direction1 = dir_px
  direction2 = dir_py
  direction3 = dir_pz
  range1 = (-6e-19, 6e-19)
  range2 = (-6e-19, 6e-19)
  range3 = (-6e-19, 6e-19)
  resolution1 = 200
  resolution2 = 200
  resolution3 = 200
  include_species:Photon
  include_species:Electron
end:dist_fn
Status-Mirror commented 1 year ago

This issue has been traced back to the push_photons subroutine in particles.F90. Unlike the push_particles routine, push_photons uses the additional particle variable particle_energy, which is only defined when the code is compiled with -DPHOTONS or -DBREMSSTRAHLUNG. This variable is set upon creation of a photon using the physics packages, but not when photons are either injected or initialised. In these cases particle_energy=0, which leads to infinite photon displacement, which immediately pushes photons outside the simulation window.

In this fix, we properly set the particle_energy variable for both injectors, and in the auto_load subroutine. Pull request #510 has been submitted with this fix.