Warwick-Plasma / epoch

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

Relativistic electrons #531

Closed tobiassh0 closed 12 months ago

tobiassh0 commented 1 year ago

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

Status-Mirror commented 1 year ago

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

tobiassh0 commented 1 year ago

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

Status-Mirror commented 1 year ago

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