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

Validating collisions #463

Open bzdjordje opened 1 year ago

bzdjordje commented 1 year ago

I have been running high density cases, 100-1000*nc, and decided to do some sanity checks for collisions. I noticed that despite being exactly energy conserving the Sentoku model seems to result in growing field energy in the EPOCH implementation. Likewise, both Nanbu and Sentolu have almost twice the total particle energy versus the collisionless case example. I initialized the plasma at Te=1000 eV and Ti=100 eV to see how thermal equilibration is achieved. Collisionless temperatures remain constant, Nanbu seems to follow the simple theoretical model per the NRL formulary although not exactly, not sure if my Lambda=4.3 is correct, and then Sentoku is very off.

image

image

image

I noticed a similar discussion in #94 but I am not sure if this strong discrepancy was observed.

Are there any caveats for when field ionization is also considered?

I noticed in the Epoch implementation for collisions there are threshold values for the Coulomb log, should these be modified?

calc_coulomb_log = 0.0_num
DO j = 1-ng, ny+ng
DO i = 1-ng, nx+ng
  local_ekbar1 = MAX(ekbar1(i,j), 100.0_num * q0)
  local_temp2 = MAX(temp2(i,j), 100.0_num)
  IF (dens1(i,j) <= 1.0_num .OR. dens2(i,j) <= 1.0_num) THEN
    calc_coulomb_log(i,j) = 1.0_num
  ELSE
    bmax = SQRT(epsilon0 * q0 * local_temp2 / (ABS(q2) * q0 * dens2(i,j)))
    b0 = ABS(q1 * q2) / (8.0_num * pi * epsilon0 * local_ekbar1)
    gamm = (local_ekbar1 / (m1 * cc)) + 1.0_num
    dB = 2.0_num * pi * h_bar / (SQRT(gamm**2 - 1.0_num) * m1 * c)
    bmin = MAX(b0, dB)
    calc_coulomb_log(i,j) = MAX(1.0_num, LOG(bmax / bmin))
  END IF
END DO
END DO

Attached is a simple input deck I used.

Thank you very much, Blagoje

input.deck.txt

Status-Mirror commented 1 year ago

Hi Blagoje,

I believe the Sentoku method is only included for backwards compatability - we do not recommend the use of this today.

I also think the threshold values prevent the Coulomb logarithm from drifting towards unphysical values. For example, if $b{max} < b{min}$, then the cumulative elastic scattering model may not be entirely appropriate.

In your input deck, you set the initial electron and ion temperatures to 1000 and 100 respectively. EPOCH uses SI units, so these are in [K], while your figure axis uses [eV]. Could this explain the NRL difference?

The main issue I see here is the potential numerical heating with the particle energy. Recently a paper was published warning of statistical issues with the Nanbu model by Higginson. We will eventually implement this fix, which may address the numerical heating, but I can't give you a time-scale of when to expect this.

Hopefully some of this was helpful, Stuart

bzdjordje commented 1 year ago

That makes sense, thank you!

bzdjordje commented 1 year ago

Hi Stuart, I am still having some issue validating basic NP collisions. I tried fixing the coulomb_log in the input deck compared to a theoretical calculation as well as tuning the basic thermal equilibration calculation to match the initial conditions of the epoch run based on what I see for calc_coulomb_log() in collisions.F90. There is still essentially the same discrepancy as before. However, I noticed that we can output the coulomb log as a function of grid point, but when I do that I just see zeros, whether in toto or per species.

print(data['Derived/Coulomb Logarithm'].data) print(data['Derived/Coulomb Logarithm/Electron'].data) print(data['Derived/Coulomb Logarithm/Hydrogen'].data)

[[0. 0. 0.] [0. 0. 0.] [0. 0. 0.]] [[0. 0. 0.] [0. 0. 0.] [0. 0. 0.]] [[0. 0. 0.] [0. 0. 0.] [0. 0. 0.]]

Not sure if this is just an error with the diagnostic or something else.

Thank you, Blagoje