Control variate for hybrid5D

Solves the following issue(s):

Closes #255 (closed) and #264 (closed)

Core changes:

• fix particle_to_grid.py

• make all the moments of background kinetic distribution particles.f_backgr.n, particles.f_backgr.vth, ... are always callable

Currently it is float when our background kinetic distribution is uniform without perturbation

Model-specific changes:

• renew&fix control variate moment calculation at the model LinearMHDVlasovCC and LinearMHDVlasovPC

• implement control variate moment calculation at the model LinearMHDDriftkineticCC

Maxwellain distribution for Particles5D

n(\mathbf x) = \iint \mathcal{M}(\mathbf x, v_\parallel, v_\perp) d v_\parallel d v_\perp = \iint \frac{1}{\sqrt{2\pi} v_{th,\parallel}} \frac{1}{v_{th, \perp}²} v_\perp \text{exp} \left[ -\frac{(v_\parallel - u_\parallel)²}{2 v_{th \parallel}²} -\frac{(v_\perp - u_\perp)²}{2 v_{th \perp}²}\right] d v_\parallel d v_\perp

Considering the Jacobian of the velocity transformation (v_\perp, \theta) \rightarrow (\mu, \theta) \left| \begin{pmatrix} \frac{v_\perp}{B} & 0 \\ 0 & 1 \end{pmatrix} \right| = \frac{v_\perp}{B}

n(\mathbf x) = \iint \mathcal{M}'(\mathbf x, v_\parallel, \mu) d v_\parallel d \mu = \iint \frac{1}{\sqrt{2\pi} v_{th,\parallel}} \frac{1}{v_{th, \perp}²} B \text{exp} \left[ -\frac{(v_\parallel - u_\parallel)²}{2 v_{th \parallel}²} -\frac{(v_\perp - u_\perp)²}{2 v_{th \perp}²}\right] d v_\parallel d \mu

DriftKinetic moment calculation with Maxwellian distribution

1. 1st order in \mu (assume u_\perp = 0)

\iint \mathcal{M}'\mu \mathbf{b}_0(\mathbf{x}) d v_\parallel d \mu = n(\mathbf x) \int^\infty_0 \frac{1}{v_{th, \perp}²} \frac{v^2_\perp}{2} \text{exp} \left[-\frac{v_\perp²}{2 v²_{th \perp}}\right] d \mu = n(\mathbf x) \int^\infty_0 \frac{1}{v_{th, \perp}²} \frac{v_\perp³}{2B(\mathbf x)} \text{exp} \left[-\frac{v_\perp²}{2 v_{th \perp}²}\right] d v_\perp = n(\mathbf x) \frac{v_{th, \perp}²}{B}

Gaussian integral \int^\infty_0 x^{2n+1} \text{exp}\left[-\frac{x²}{a²}\right] dx = \frac{n!}{2}a^{2n+2} is used.

2. 2nd order in v_\parallel

\iint \mathcal{M}'\frac{1}{B^*_\parallel(\mathbf{x}, v_\parallel)} v_\parallel² d v_\parallel d \mu = n(\mathbf x) \int^\infty_{-\infty} \frac{1}{\sqrt{2\pi}} \frac{1}{B^*_\parallel(\mathbf x, v_\parallel)} v_\parallel² \text{exp} \left[ -\frac{(v_\parallel - u_\parallel)²}{2 v_{th \parallel}²} \right] d v_\parallel

approximate B^*_\parallel(\mathbf x, v_\parallel) \approx B^*_\parallel(\mathbf x) = B(\mathbf x) + \epsilon u_\parallel (\nabla \times \mathbf b_0 (\mathbf x))

= n(\mathbf x) \frac{v_{th \parallel}² + u_\parallel²}{B^*_\parallel (\mathbf x)}

Gaussian integral \int^\infty_0 x^{2n} \text{exp}\left[-\frac{x²}{a²}\right] dx = \sqrt{\pi} \frac{a^{2n+1} (2n-1)!!}{2^{n+1}} is used.

• fix wrong implementation at the coupling step CurrentCOupling5DGradB

So far the step is solved with our standard implicit Crank-Nicolson with Schur-solver method. However, it was found that the step cannot be solved with the Schur-solver strategy, since the PDEs are not a linear.

\begin{bmatrix}
\dot{\mathbf u} \\
\dot{\mathbf H} \\
\end{bmatrix} = 
\begin{bmatrix}
0 & (\mathbb{M}^{2,n})^{-1} \mathbb{L}² \frac{1}{\bar{\sqrt{g}}} \frac{1}{\bar B^{*0}_\parallel}\bar{B}^\times_f \bar{G}^{-1} \bar{b}^\times_0 \bar{G}^{-1} \\
-\bar{G}^{-1} \bar{b}^\times_0 \bar{G}^{-1}  \bar{B}^\times_f \frac{1}{\bar B^{*0}_\parallel} \frac{1}{\bar{\sqrt{g}}} (\mathbb{L}²)^\top (\mathbb{M}^{2,n})^{-1} & 0
\end{bmatrix}
\begin{bmatrix}
\mathbb{M}^{2,n} \mathbf{u}
\\
\underbrace{\frac{A_h}{A_b} \bar M \bar W \bar{\nabla B_\parallel}}_{\text{not a function of   } \mathbf H}
\end{bmatrix}

So the step was replaced with the explicit solver.

Documentation changes:

None