Enable TAE-EP simulations (fix several bugs)

[bkna0327]

Core changes:
Allow TAE-EP simulations with devel.

• Renew CanonicalMaxwellian.
• Pass n_cols_diag and n_cols_aux to the class Particles.

Model-specific changes:

Bug fix in the model LinearMHDDriftkineticCC

[spossann]

Looks good! Have written some suggestions.

[spossann]

A type KineticBackground should be evaluated at (eta1, eta2, ...) etc. thus on the unit cube, not at psic. Wouldn't it be more consistent to have a function eta_to_psic and use it in the evaluations?

[bkna0327]

If we do that, we need to recalculate $\psi_c$ for all particles every time, whereas it is currently calculated only once with the kernel and stored in the marker array. Just so you know, psic_to_eta is not $\psi_c \rightarrow \eta$ but $\psi_c \rightarrow r_c$ where $r_c = \sqrt{\frac{\psi_c - \psi_\text{axis}}{\psi_\text{edge} - \psi_\text{axis}}}$.

[bkna0327]

May be we should rename it to psic_to_rc.

[spossann]

We have two options here:

1. We keep the hierarchy class CanonicalMaxwellian(KineticBackground):, then CanonicalMaxwellian needs to respect what the base class demands (see docstring).
2. Do not inherit, thus make a new class class CanonicalMaxwellian:; then we loose some consistency with the rest of the code and have to adapt a bit.

I prefer option 1. to have a unique class for backgrounds, which can then be added/subtracted etc. If I understand correctly you say we do not need to re-evaluate the canonical Maxwellian, and just evaluate it at the beginning and then store the value in the marker array? This could be done also in eta-space and the problem is solved, right?

[bkna0327]

I mean, in the case of control_variate=True scheme, we need to evaluate the canonical Maxwellian whenever we push particles. So far the approach was that first we evaluate constants of motion of all particles $\psi_c$ and $\varepsilon$ with the pyccelized kernels and store them in the marker array and then use the stored values when we evaluate the canonical Maxwellian. With your suggested approach, if I understand correctly, we should implement the function which converts an array of etas into an array of $\varepsilon$ and $\psi_c$ and keep etas as an input argument right?

[spossann]

I mean, in the case of control_variate=True scheme, we need to evaluate the canonical Maxwellian whenever we push particles. So far the approach was that first we evaluate constants of motion of all particles ψ c and ε with the pyccelized kernels and store them in the marker array and then use the stored values when we evaluate the canonical Maxwellian. With your suggested approach, if I understand correctly, we should implement the function which converts an array of etas into an array of ε and ψ c and keep etas as an input argument right?

What I don't understand is the following: When the value of the canoncial Maxwellian never changes for a particle, because the arguments are constants of motion, why not store the value of the Canonical Maxwellian in the marker array?

On the other hand, if the value does change, and we need to recompute anyway, then why have an intermediate step for storing the changed "constants of motion", and not evaluate the Maxwellian directly? In the eval we can compute the constants, as you suggested.

[spossann]

Here it is not clear whether one evaluates at eta or at psic. I would opt for eta for consistency.

[bkna0327]

Currently, args includes the canonical toroidal momentum psic and the energy $\varepsilon$.

[spossann]

Eval at eta here?

[bkna0327]

In fact, we use $r_c$ for the evaluation.

[spossann]

Please adapt the name to rc to be clear.

[bkna0327]

But the input argument is still psic. In the function, we calculate r_c from the psic, then evaluate at the r_c.

[spossann]

Ok, so in general there is nothing wrong with having this helper function that evaluate at psic, but these should be called inside the methods demanded by the base class.

[bkna0327]

Hi @spossann, are you considering any further modifications for this MR?

[spossann]

Hi @spossann, are you considering any further modifications for this MR?

I still want to get some consistency for the arguments of the CanonicalMaxwellian.

[bkna0327]

That is possible, but then we need to pass $v_\parallel$ and $\mu$ in addition to $\eta_1, \eta_2, \eta_3$ to calculate $\psi_c$ internally. Would it be fine?

[spossann]

That is possible, but then we need to pass v ∥ and μ in addition to η 1 , η 2 , η 3 to calculate ψ c internally. Would it be fine?

Yes I think it is fine.

[bkna0327]

After implementing this, I realized that the performance degradation might be more significant than I initially expected. The current fix introduces multiple large 2D array computations (roughly number_of_particle_sizes × 6). Moreover, these calculations must be done every time the particle positions or velocities are updated when using the control variate method.

[max-models]

After implementing this, I realized that the performance degradation might be more significant than I initially expected. The current fix introduces multiple large 2D array computations (roughly number_of_particle_sizes × 6). Moreover, these calculations must be done every time the particle positions or velocities are updated when using the control variate method.

We should definitely finish the performance testing to make sure this doesn't happen 😅

[bkna0327]

After implementing this, I realized that the performance degradation might be more significant than I initially expected. The current fix introduces multiple large 2D array computations (roughly number_of_particle_sizes × 6). Moreover, these calculations must be done every time the particle positions or velocities are updated when using the control variate method.

This comment refers specifically to the most recent commit that I made in this MR, not to the overall approach in general ;)