Warwick-Plasma/epoch

Relativistic electrons

Closed this issue · 3 comments

Hi,

I'm trying to run the attached input deck but EPOCH is unable to create a single output file, running on Avon (32 nodes, 48 ntasks/node and 3700 mem/cpu).

I imagine it is struggling because of the relativistic ring-beam electrons, or perhaps just the scale of the initialisation. The CFL condition is maintained in this case.

Thanks

Hey Tobias,

You are right - dist_fn sampling is the issue here. The sampling is performed using the sample_from_deck_expression subroutine in src/user_interaction/particle_temperature.F90. We don't have developer documentation for this, so I looked into it with a few simple tests.

The dist_fn sampler works as follows: a set of px, py and pz values are randomly picked within the ranges you specify (missing a range forces that component to 0). These momenta are substituted into your dist_fn expression, which returns a number that EPOCH saves to the probability variable. A random number between 0 and 1 is sampled, and if the random number is less than probability, then the px, py and pz values are saved. For example, if dist_fn returned 0.2, then we would only have a 20% chance of setting the particle momentum to this. Upon failure, a new set of px, py and pz values are sampled.

This method means the dist_fn in EPOCH has a particular importance, which isn't mentioned in the user documentation.
For optimal performance, you should normalise dist_fn such that the maximum value is 1. If the maximum was, say, 0.1 - then your generated distribution function would still be correct, but the code would arbitrarily throw out 9 out of every 10 guesses and take 10x longer.

For your case, you use something equivalent to:

  dist_fn = exp(-0.5 * ((sqrt(px^2 + py^2) + 15.2) / 1.7)^2) * exp(-0.5 * ((pz - 7.8) / 0.17)^2)
  dist_fn_px_range = (-34.2, 34.2)
  dist_fn_py_range = (-34.2, 34.2)
  dist_fn_pz_range = (-2.49, 18.0)

This has a maximum value of 4.3e-18, when px=py=0. This means a px, py set can only be picked with a probability of, at best, 4.3e-18. Hence, you will never initialise this deck in a reasonable time. I suspect your main problem here is the dist_fn term: exp(-0.5 * ((sqrt(px^2 + py^2) - pring) / p_ring_th)^2), where pring is negative (which you may not be expecting). The pring term is defined as:

pring_beam = sqrt(2 * m3 * qe * ring_beam_energy)  # Equals: 1.708498854424136e-22
pitch_angle = 153.0
pring = pring_beam * cos(pitch_angle*pi/180)       # Equals: -1.522283625860258e-22

Your dist_fn would have the correct normalisation if you used:

  dist_fn = exp(-0.5 * ((sqrt(px^2 + py^2) - abs(pring)) / p_ring_th)^2) \
          * exp(-0.5 * ((pz - abs(pbeam)) / p_beam_th)^2)

I will keep this issue open until I update the EPOCH documentation to make it easier for people to figure out what this subroutine does.

Cheers,
Stuart

Hi Stuart,

Thanks for this, input deck now runs. What I'm getting now is a smaller simulation that completes at ~500 files rather than the 12000 files that I specify. I understand that the time-step is merely a "minimum", as per EPOCH documentation, so I tried reducing it by an order of magnitude or so and it still only outputs ~500 files.

Presume this is because of the relativistic electrons? What else can I do to help this run with increased time-step granularity?

Thanks,
Tobias

Hey Tobias,

In your first attached input deck, my calculations suggest you're using:

nx = 28394
t_end = 112.81e-12
x_min = 0.0
x_max = 2.0

From this, we can deduce your cell-size is $70.44 \mu m$, and for the default field order you use, dt = 1.0 * dx / c, or 0.235 ps. This should yield the roughly 480 time-steps you observe.

I'm not sure how you went about reducing the time-step, but I would advise trying the dt_multiplier in the control block. If you want to boost the number of steps from 500 to 12000, then you need 24x as many steps, so set dt_multiplier = 1/24.

Cheers,
Stuart