search

google-deepmind / torax

TORAX: Tokamak transport simulation in JAX
https://torax.readthedocs.io
Other
334 stars 29 forks source link

[bug] `NaN` in `d_face_el` and `v_face_el` when using Bohm-gyroBohm #403

Open theo-brown opened 2 days ago

theo-brown commented 2 days ago

Behaviour: simulation exits due to nan in core profiles when run with the BgB transport model. Diagnosis: In some scenarios, chi terms are zero on axis. This causes the D_face on axis to be nan (as it divides by chi), and hence the V_face to be nan (as it is proportional to D_face).

theo-brown commented 2 days ago

A simple solution would be to add a small term to the chis. I have done this manually and the simulations run; however, I'm going to do some investigating to see how other codes (RAPTOR, JETTO) handle the same problem.

jcitrin commented 2 days ago

Yes, please compare to RAPTOR and JETTO and we can decide on the fix.

But for D it seems pretty clear that the zero-gradient boundary condition means D=0 on-axis: chi_e=chi_i=0 on axis (see eq 8+9 in Tholerus 2024), and thus in eq.10, we have 0^2/0, i.e. 0. So can define d_face_el as something like:

d_face_el = jnp.concatenate((jnp.zeros(1), weighting[1:] * chi_e[1:] * chi_i[1:] / (chi_e[1:] + chi_i[1:]))))

Same approach for V. It scales with S^2/(Volume*Volume') , and since (roughly) S= pi r^2 and Volume = 2pi^2r^2R , V scales with r and is zero on axis.

NOTE: there seems to be a bug in v_face_el. It's missing a geo.rho_b term, when comparing to equation 12 in Tholerus 2024. Can be fixed in the same opportunity as fixing the divs by zero.

jcitrin commented 2 days ago

If you'd like to send the PR go ahead, otherwise can do it on our side. Would be good to add a BgB case with particle transport to the test suite.

theo-brown commented 2 days ago

RAPTOR adds its inner transport patch, a neoclassical term, and a sawtooth term before D and v are calculated (line 276 in chi_BgB.m).

I'm not proficient in FORTRAN, but from my initial scan of the JETTO repo I think it sets D=0 where chi=0 (line 224 in callcgm.f).

Not sure whether there are any edge cases where you might get chi = 0 elsewhere in the profile; if there are, it would be better to do a jnp.where to catch these points too.

theo-brown commented 2 days ago

Thanks for catching the missing rho_b!

jcitrin commented 2 days ago

To catch any other edge cases, then it's sufficient to add a small term to the denominator of the D formula. i.e. add the eps from constants to chii+chie. That will set both D and V to zero safely when chi=0 also for rhon>0 cases