Other Laser-Plasma example decks

Now that you have a basic understanding of how the input decks work, you should be able to work through the remaining example decks by referring to the User manual for a description of any new parameters used. Several experienced users of the code will be available throughout the duration of the workshop, so if you want help with anything please don’t hesitate to ask. The decks are:

• 01-1d_laser.deck Described in notes above (here)
• 02-2d_laser.deck Described in notes above (here)
• 03-1d_two_stream.deck Described in notes above (here)
• 04-1d_two_stream_io.deck This is the same as the previous deck but with the addition of more sophisticated output diagnostics
• 05-2d_moving_window.deck This deck contains an example of firing a laser into a plasma and then using the moving window facility to track the wave front as it moves beyond the edge of the original domain.
• 06-2d_ramp.deck This deck contains an example of firing a laser at a plasma with a ramped density profile.
• 07-1d_heating.deck This deck contains a setup for investigating the anomalous heating of a plasma that occurs for purely resolved systems.

Other things to try

1. Landau damping predicts collisionless damping of electrostatic waves. Setup a 1D problem with an electrostatic wave and check for a range of wavelengths. Points to note:
   1. Does the answer depend on whether the initial condition is a travelling wave or standing wave? How are these setup?
   2. Look for trapping in the Langmuir wave
   3. Check the damping rate against published formulae. Try for a range of $k\lambda_D$ as the most commonly reported formulae assume $kλ_D≪1$
   4. The answer is more accurate, assuming you have enough grid-points and particles to get a good answer, if you ignore the first maxima or two – why?
2. A more realistic instability than the two cold beams tested above is the bump-on-tail instability. Setup a 1D bump-on-tail distribution and check that the simple formula for the growth-rates is correctly reproduced. The main problem with the initial conditions is how to setup a suitable initial distribution.
3. Try setting up the initial conditions for a problem of direct relevance to your research. This may be too computationally demanding to run on the workshop computers but it is a good exercise as you can get some help on trickier input decks and diagnostic planning than the simple exercises so far.
4. Check that EPOCH works as expected on your host institution computer. If not we may be able to help before you leave.

Copies of the decks

The decks can be downloaded here and viewed or copied from here:

begin:control
  nx = 400

  # Size of domain
  x_min = 0
  x_max = 5.0e5

  # Final time of simulation
  t_end = 1.5e-1

  stdout_frequency = 400
end:control


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


begin:constant
  drift_p = 2.5e-24
  temp = 273
  dens = 10
end:constant


begin:species
  # Rightwards travelling electrons
  name = Right
  charge = -1
  mass = 1.0
  temp = temp
  drift_x = drift_p
  number_density = dens
  npart = 4 * nx
end:species


begin:species
  # Leftwards travelling electrons
  name = Left
  charge = -1
  mass = 1.0
  temp = temp
  drift_x = -drift_p
  number_density = dens
  npart = 4 * nx
end:species


begin:output
  name = normal

  # Number of timesteps between output dumps
  dt_snapshot = 1.5e-3

  # Properties at particle positions
  particles = always
  px = always

  # Properties on grid
  grid = always
  ex = always
  ey = always
  ez = always
  bx = always
  by = always
  bz = always
  jx = always
  ekbar = always
  mass_density = never + species
  charge_density = always
  number_density = always + species
  temperature = always + species

  distribution_functions = always
end:output


begin:output
  name = restart

  # Number of timesteps between output dumps
  dt_snapshot = 0.15

  restartable = T
end:output


begin:dist_fn
  name = x_px
  ndims = 2

  direction1 = dir_x
  direction2 = dir_px

  # Range is ignored for spatial coordinates
  range1 = (1, 1)
  range2 = (-5e-24, 5e-24)

  # Resolution is ignored for spatial coordinates
  resolution1 = 1
  resolution2 = 200

  include_species:Left
  include_species:Right
end:dist_fn

begin:constant
  x0 = 20 * micron
  lambda = 10 * micron

  t_laser = 120 * femto
  sigma_t = t_laser / 2 / sqrt(loge(2))

  w0_laser = 30 * micron
  sigma_w0 = w0_laser / 2 / sqrt(loge(2))

  den_peak = 5.0e19 * 1.0e6
  win_start = 340 * femto
end:constant


begin:control
  nx = 1550 / 8
  ny = 600 / 8
  npart = (60e6) / 8

  # Size of domain
  x_min = 0
  x_max = 155 * micron
  y_min = -30 * micron
  y_max = -y_min

  # Final time of simulation
  t_end = 1600 * femto

  stdout_frequency = 1
  print_eta_string = T
end:control


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


begin:species
  name = electron
  charge = -1.0
  mass = 1.0
  number_density = if((x lt x0), 0.0, den_peak)
  frac = 0.5
end:species


begin:species
  name = proton
  charge = 1.0
  mass = 1836.2
  number_density = number_density(electron)
  frac = 0.5
end:species


begin:output
  name = normal
  dt_snapshot = 50 * femto

  grid = always
  ex = always
  ey = always
  ez = always
  bx = always
  by = always
  bz = always
  jx = always
  jy = always
  ekbar = always
  mass_density = never + species
  charge_density = always
  number_density = always + species
  temperature = always + species
end:output


begin:output
  name = large
  dt_snapshot = 500 * femto

  particles = always
  particle_weight = always
end:output


begin:laser
  boundary = x_min
  intensity_w_cm2 = 1.9e18
  lambda = lambda
  t_profile = gauss(time, 2*sigma_t, sigma_t)
  profile = gauss(y, 0, sigma_w0)
end:laser


begin:window
  move_window = T
  window_v_x = c * 0.87
  window_start_time = win_start
  bc_x_min_after_move = simple_outflow
  bc_x_max_after_move = simple_outflow
end:window

begin:constant
  # Particles per cell
  part = 32

  las_lambda = 1 * micron
  las_omega = 2.0 * pi * c / las_lambda
  las_time = 2.0 * pi / las_omega

  n_crit = critical(las_omega)

  max_dens = 0.8 * n_crit

  scale_x = 20 * micron

  las_scale_y = 8 * micron

  xmin = -4 * micron

  # Gaussian Beam stuff
  w0 = las_scale_y
  rayleigh_range = pi * w0^2 / las_lambda
  wz = w0 * sqrt(1 + (x_start / rayleigh_range)^2)
  radius_of_curvature = x_start * (1.0 + (rayleigh_range / x_start)^2)
end:constant


begin:control
  nx = 1024 / 4
  ny = 512 / 4

  # Final time of simulation
  t_end = 0.4 * pico

  # Size of domain
  x_min = xmin
  x_end = scale_x + 20 * micron

  y_min = -20 * micron
  y_max = -y_min

  stdout_frequency = 10
end:control


begin:laser
  boundary = x_min
  intensity_w_cm2 = 1.0e16
  omega = las_omega
  t_profile = if (time lt 2*las_time, gauss(time, 2*las_time, 2*las_time), 1)
  profile = (1.0 + 0.05 * sin(32.0*pi*y/lengthy)) * gauss(y, 0, las_scale_y)
end:laser


begin:boundaries
  bc_x_min = simple_laser
  bc_x_max = simple_outflow
  bc_y_min = periodic
  bc_y_max = periodic
end:boundaries


begin:species
  # Electron
  name = electron
  charge = -1.0
  mass = 1.0
  npart = nx * ny * part

  number_density = max_dens * (exp(x/scale_x) - 1) / (exp(1) - 1)
  number_density = if(x lt 0, 0.0, number_density(electron))
  number_density = if(number_density(electron) gt max_dens, max_dens, \
                      number_density(electron))
  number_density = if(x gt 75*micron, 0.0, number_density(electron))
  #number_density = number_density(electron) \
                    * (0.8 + 0.2 * gauss(y, 0, 0.5*las_scale_y))
  number_density_min = 0.0001 * n_crit
  number_density_max = n_crit

  temp_ev = 10^3
end:species


begin:species
  # Protons
  name = proton
  charge = 1.0
  mass = 1836.2
  npart = nx * ny * part

  number_density = number_density(electron)
  number_density_min = 0.0001 * n_crit
  number_density_max = 1.2 * n_crit

  temp_ev = 40
end:species


begin:output
  name = normal

  # Number of timesteps between output dumps
  dt_snapshot = 5 * femto

  # Properties at particle positions
  particles = always
  px = always
  particle_weight = always

  # Properties on grid
  grid = always
  ex = always
  ey = always
  ez = always
  bx = always
  by = always
  bz = always
  jx = always
  jy = always
  jz = always
  ekbar = always + species
  mass_density = never + species
  charge_density = always # + average + snapshot
  number_density = always + species
  temperature = never + species

  # Extended io
  distribution_functions = always
end:output


begin:dist_fn
  name = en
  ndims = 1

  direction1 = dir_en

  range1 = (0, 15*kev)

  resolution1 = 5000

  include_species:electron
end:dist_fn


begin:dist_fn
  name = x_en
  ndims = 2

  direction1 = dir_x
  direction2 = dir_en

  # Range is ignored for spatial coordinates
  #range1 = (1, 1)
  range2 = (0, 15*kev)

  # Resolution is ignored for spatial coordinates
  #resolution1 = 1
  resolution2 = 1500

  include_species:electron
end:dist_fn


begin:dist_fn
  name = x_px
  ndims = 2

  direction1 = dir_x
  direction2 = dir_px

  # Range is ignored for spatial coordinates
  #range1 = (1, 1)
  range2 = (-5e-23, 5e-23)

  # Resolution is ignored for spatial coordinates
  #resolution1 = 1
  resolution2 = 1500

  include_species:electron
end:dist_fn


begin:probe
  name = electron_probe

  point = (0.5 * (x_max + x_min), y_min)
  normal = (1, 0)

  include_species:electron
  include_species:proton
end:probe

begin:constant
  dl = 74.33942 * micron
end:constant


begin:control
  nx = 10

  # Size of domain
  x_min = 0
  x_max = 14000 * dl

  # Final time of simulation
  t_end = 1.5e-2

  stdout_frequency = 10000
  print_eta_string = T
end:control


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


begin:species
  name = electron
  charge = -1
  mass = 1.0
  temp_x_ev = 1
  number_density = 1e16
  npart = nx * 5
end:species


begin:output
  name = normal

  # Number of timesteps between output dumps
  dt_snapshot = 1.5e-3

  # Properties on grid
  grid = always
  ekbar = always
  temperature = always
end:output


begin:output
  name = large

  # Number of timesteps between output dumps
  dt_snapshot = 75e-3

  # Properties at particle positions
  particles = always
  px = always
  py = always
  pz = always
end:output

Remote Visualisation with VisIt

If the local workstation you are using isn’t big enough for your test problems you may also use a your host institutes HPC cluster.

Remote Visualisation with VisIt

Most large simulations are carried out on a remotely located machine. Often this machine is located many miles away, perhaps even in a different country. Viewing data on remote systems can be awkward and poor network speeds can often make it nearly impossible. The VisIt visualisation tool solves this problem by using a client-server model. The program which reads, processes and renders the data is completely separated from the program which displays the results on the screen. It is therefore possible to run VisIt on your local machine and look at data located on a different machine. The method of setting this up varies depending on the configuration of the remote machine so we will not go into details here. However, the desktop machines have been setup to be able to view data located on remote clusters so you can try it out.

In the VisIt control window, click the “Open” button which launches a file browser window. The first entry is called “Host” and contains a drop-down list of all configure remote machines.

If you want to know more about how to set up remote visualisation in VisIt, you can ask one of the Warwick staff members.

When viewing data across a slow network connection, there is one more useful thing to know. VisIt has two methods of drawing plots generated on a remote machine. The first method is to construct the polygons used in drawing the plot on the remote machine and send them across the network. The local machine then turns these into a plot image. This makes manipulating the figure very fast (zooming, rotating, etc), since all the polygons that generate the image are on the local machine. However, if there are a lot of polygons then they can be slow to transfer across the network. They can also use up a lot of memory. For these cases, the alternative is to render the image on the remote machine and just transfer the image across the network. The downside of this approach is that whenever you manipulate the plot, it must be re-drawn on the remote machine and then transferred across the network again. The options controlling this behaviour are to be found under “Options->Rendering” in the “Advanced” tab. The feature is called “scalable rendering”.

Collisions in EPOCH

EPOCH now contains a collision routine based on the technique outlined in Sentoku & Kemp¹

Collisions are enabled using the output block named collisions which accepts the following three parameters.

• use_collisions – This is a logical flag which determines whether or not to call the collision routine. If omitted, the default is “true” if any of the frequency factors are non-zero (see below) and “false” otherwise.

• coulomb_log – This may either be set to a real value, specifying the Coulomb logarithm to use when scattering the particles or to the special value “auto”. If “auto” is used then the routine will calculate a value based on the properties of the two species being scattered. If omitted, the default value is “auto”.

• collide – This sets up a symmetric square matrix of size nspecies*nspecies containing the collision frequency factors to use between particle species. The element (s1,s2) gives the frequency factor used when colliding species s1 with species s2. If the factor is less than zero, no collisions are performed. If it is equal to one, collisions are performed normally. For any value between zero and one, the collisions are performed using a frequency multiplied by the given factor. If “collide” has a value of “all” then all elements of the matrix are set to one. If it has a value of “none” then all elements are set to minus one. If the syntax “species1 species2 ” is used, then the (species1,species2) element of the matrix is set to the factor “”. This may either be a real number, or the special value “on” or “off”. The “collide” parameter may be used multiple times. The default value is “all” (ie. all elements of the matrix are set to one).

For example:

begin:collisions
  use_collisions = T
  coulomb_log = auto
  collide = all
  collide = spec1 spec2 off
  collide = spec2 spec3 0.1
end:collisions

With this block, collisions are turned on and the Coulomb logarithm is automatically calculated. All values of the frequency array are set to one except (spec1,spec2) is set to minus one (and also (spec2,spec1)) and (spec2,spec3) is set to 0.1

Ionisation in EPOCH

EPOCH includes field ionization which can be activated by defining “field_ionisation = T” in the control block along with ionisation energies and an electron for the ionising species in one of the species blocks. This is done via the species block in the “ionisation_energies” and “electron_species” parameter respectively. “ionisation_energies” should be given as a list in joules, and “electron_species” should be the name of the species to be used as the electron species. For example, ionising carbon species might appear in the input deck as:

begin:species
  charge = 0.0
  mass = 1837.2
  name = carbon
  ionisation_energies = \
    (11.26*ev, 24.38*ev, 47.89*ev, 64.49*ev, 392.1*ev, 490.0*ev)
  electron_species = electron
  number_density = den_gas
end:species

begin:species
  charge = -1.0
  mass = 1.0
  name = electron
  number_density = 0.0
end:species

It is possible to define different electron species for each ionisation level, which is particularly useful in monitoring specific ionisation levels. If we wished to monitor the fourth ionisation level of carbon in the above example, the above example might appear:

begin:species
  charge = 0.0
  mass = 1837.2
  name = carbon
  ionisation_energies = \
    (11.26*ev, 24.38*ev, 47.89*ev, 64.49*ev, 392.1*ev, 490.0*ev)
  electron_species = (electron, electron, electron, fourth, electron, electron)
  number_density = den_gas
end:species

begin:species
  charge = -1.0
  mass = 1.0
  name = electron
  number_density = 0.0
end:species

begin:species
  charge = -1.0
  mass = 1.0
  name = fourth
  number_density = 0.0
end:species

Field ionisation consists of three distinct regimes; multiphoton in which ionisation is best described as absorption of multiple photons, tunneling in which deformation of the atomic coulomb potential is the dominant factor, and barrier suppression ionisation in which the electric field is strong enough for an electron to escape classically. It is possible to turn off multiphoton or barrier suppression ionisation through the input deck by adding “use_multiphoton=F” and/or “use_bsi=F” to the control block.

QED Effects in EPOCH

EPOCH has recently been extended to include some quantum electrodynamic effects that are important for high intensity (>) lasers. The two processes that are included are

• Gamma ray production by QED corrected synchrotron emission (Also called magnetic bremsstrahlung or nonlinear Compton scattering).
• Electron positron pair production by the Breit-Wheeler process from these gamma ray photons.

For more information on the theory see Duclous et al. ²

Simulating the QED effects increases EPOCH’s memory requirements and so the code has to be compiled with the correct compilation options to turn the module on. To turn the module on, open “Makefile” in an editor and find the commented out line #DEFINES += $(D)PHOTONS. Uncomment this line, then type “make clean” and then “make” (remember to include the COMPILER= if you haven’t specified the environment variable) to rebuild the code with QED support.

Once the code is built with QED support, actually turning on QED for a specific simulation requires the addition of a new block into the input deck. This block is simply called qed and starts with the usual “begin:qed” and “end:qed” markers of the other blocks. The parameters which can go into the block are:

• use_qed - Turns QED on or off. If you don’t want QED effects at all then compile the code without the “-DPHOTONS” lines in the makefile.
• qed_start_time - Specifies the time after which QED effects should be turned on. For example you can turn off the routines until a laser has crossed the vacuum region in front of the target.
• produce_photons - Specifies whether you’re interested in the photons generated by synchrotron emission. If this is F then the radiation reaction force is calculated but the properties of the emitted photons are not tracked.
• photon_energy_min - Minimum energy of produced photons. Radiation reaction is calculated for photons of all energies, but photons with energy below this cutoff are not tracked.
• photon_dynamics - If F then photons are generated, but their motion through the domain is not simulated and they stay where they were generated. Photon motion is often less interesting than photon generation unless you want to simulate pair production. In these cases set this to F.
• produce_pairs - Whether or not to simulate the process of pair generation from gamma ray photons. Both produce_photons and photon_dynamics must be T for this to work.
• qed_table_location - EPOCH’s QED routines use lookup tables to calculate gamma ray emission and pair production. If you want to use tables in a different location from the default put the location in this parameter.

QED also requires that the code now know which species are electrons, positrons and photons. Rather than try to do this automatically the user has to specify the type of a species. This is done by using a single “identify” tag in a species block. To specify an electron the block in the deck would look like

begin:species
  name = electron
  frac = 0.5
  number_density = 7.7e29
  identify:electron
end:species

Once the identity of a species is set then the code automatically assigns mass and charge states for the species. At present, the user cannot override these. Possible identities are

• electron : A normal electron species. All species of electrons in the simulation must be identified in this way or they will not generate photons.
• positron : A normal positron species. All species of positron in the simulation must be identified in this way or they will not generate photons.
• photon : A normal photon species. One species of this type is needed for photon production to work. If multiple species are present then generated photons will appear in the first species of this type.
• bw_electron : The electron species for pair production. If a species of this type exists then electrons from the pair production module will be created in this species. If no species of this type is specified then pair electrons will be generated in the first electron species.
• bw_positron : The positron species for pair production. If a species of this type exists then positrons from the pair production module will be created in this species. If no species of this type is specified then pair positrons will be generated in the first positron species.

A species should be identified only once, so a “bw_electron” species does not need to also be identified as an “electron” species. If the code is running with “produce_photons=T” then a photon species must be created by user and identified. If the code is running with “produce_pairs=T” then the code must specify at least one electron (or bw_electron) species and one positron (or bw_positron) species. The code will fail to run if the needed species are not specified.

Other Useful Info

Bug reports, feature requests and questions

All questions and requests after the workshop should be posted on the GitHub EPOCH project web page

The VisIt programme

The VisIt programme is free. It can be downloaded from https://wci.llnl.gov/simulation/computer-codes/visit/ There are many pre-compiled binaries so this ought to be easy. If you have any problems post a question on the GitHub EPOCH project.

GDL not IDL

If you don’t have IDL, or don’t want to pay for it!, then the free GDL is available from http://gnudatalanguage.sourceforge.net/

Updating EPOCH

To update to the latest version of EPOCH simple cd into your Epoch directory and enter ‘git pull’. This will work fine provided you haven’t edited any of the Fortran source code. If you have edited the source code then you need to learn git.

Getting Old Copies of EPOCH

You can also checkout an old version of EPOCH, you may want to get the version used 18 months ago to reproduce some previous simulations exactly for example. In this case it is best to checkout a new branch in the EPOCH repository. If you wanted the version from 10 February 2010 for example you would first enter git log --before=2010-02-11 This will give you the log of commits in reverse order, starting on the 11th of February. Identify the commit you want and copy the commit hash (the long string of numbers and letters following the word “commit”). To checkout a copy of this version of the code, type git checkout -b old-code After this your repository will reflect the state of the code at that point in time. To get back to the current version, just type git checkout master

References