Get em back

[Antoinehoff]

This PR is related to gkeyllorg/gkeyll#867 and implements the Maxima code dedicated to generating electromagnetic (EM) kernels.

apar contribution:

The apar contribution is related to the effective magnetic field

$$B^*_h = B_{0h} + \frac{m v_{\parallel h}}{q} \nabla \times \hat{ b} + \nabla \times \left(A_{\parallel h} \hat{ b}\right).$$

Volume kernels

In the volume kernels, the contribution of the $\nabla \times \left(A_{\parallel h} \hat{ b}\right)$ term is added to the characteristics as

  /* Expand Apar.*/
  Apar_e : doExpand1(apar,bC),

  /* Expand dBperp/Bmag. */
  dBperpoverB_x : (rdy2*diff(Apar_e*b_z_e,y) - rdz2*diff(Apar_e*b_y_e,z))*jacobTotInv_e,
  dBperpoverB_y : (rdz2*diff(Apar_e*b_x_e,z) - rdx2*diff(Apar_e*b_z_e,x))*jacobTotInv_e,
  dBperpoverB_z : (rdx2*diff(Apar_e*b_y_e,x) - rdy2*diff(Apar_e*b_x_e,y))*jacobTotInv_e,
  dBperpoverB_list : [dBperpoverB_x, dBperpoverB_y, dBperpoverB_z],

Flux kernels

In the flux kernels, the contribution is added following

$$\nabla \times \left(A_{\parallel h} \hat{ b}\right) = A_{\parallel h} \nabla \times \hat{ b} + \nabla A_{\parallel h} \times \hat{ b},$$

which will be then used to compute a dot product with $\nabla H$.
The implementation of the first term, $A_{\parallel h} \nabla \times \hat{ b}$, is straightforward since we can just reuse the dualcurlbhat_quad variable to have $\nabla \times \hat{ b}$ and multiply it by Apar.
For the second term, $\nabla A_{\parallel h} \times \hat{ b}$, we evaluate the gradient of Apar as it is done for the Hamiltonian. Then we need to compute the covariant components

$$e^m \cdot \nabla A_{\parallel h} \times \hat{ b} = \nabla A_{\parallel h} \cdot \hat{ b} \times e^m = \nabla A_{\parallel h} \cdot \frac{1}{\sqrt{g_{33}}} e_3 \times e^m = \nabla A_{\parallel h} \cdot \frac{g_{3i}}{\sqrt{g_{33}}} e^i\times e^m.$$

Writing the gradient of $A_{\parallel h}$ in terms of its covariant components, $\nabla A_{\parallel h} = \partial_i A_{\parallel h} e^i$, and using $e^j \times e^k = e_i/J_c$, we can compute the covariant components of the second term as,

$$e^1 \cdot \nabla A_{\parallel h} \times \hat{ b} = \frac{1}{J_c \sqrt{g_{33}}} \left(g_{33} \partial_{x_2} A_{\parallel h} - g_{32} \partial_{x_3} A_{\parallel h}\right),$$ $$e^2 \cdot \nabla A_{\parallel h} \times \hat{ b} = \frac{1}{J_c \sqrt{g_{33}}} \left(g_{31} \partial_{x_3} A_{\parallel h} - g_{33} \partial_{x_1} A_{\parallel h}\right),$$ $$e^3 \cdot \nabla A_{\parallel h} \times \hat{ b} = \frac{1}{J_c \sqrt{g_{33}}} \left(g_{32} \partial_{x_1} A_{\parallel h} - g_{31} \partial_{x_2} A_{\parallel h}\right).$$

We can now compute the dot product $\nabla H \cdot (\nabla \times (A_\parallel b))$, which yields for cdim=3

/* EM term curl(Apar*b) = Apar * curl(b) grad(Apar) x b */
for k : 1 thru cdim do (
  /* Apar * curl(b) */
  printf(fh, "  alpha_quad += -(~a)*m_bmag_inv * (~a)*dualcurlbhat_quad[~a]; ~%", dH_dz_nodes[k][j1index], apar_nodes[i1index], k-1)
),
/* grad(Apar) x b */
printf(fh, "  alpha_quad += -(~a)*m_bmag_inv * 1/mag_e_3/Jc_quad * (g_33 * (~a) - g_23 * (~a)); ~%", dH_dz_nodes[1][j1index], dA_dx_nodes[2][i1index], dA_dx_nodes[3][i1index]),
printf(fh, "  alpha_quad += -(~a)*m_bmag_inv * 1/mag_e_3/Jc_quad * (g_13 * (~a) - g_33 * (~a)); ~%", dH_dz_nodes[2][j1index], dA_dx_nodes[3][i1index], dA_dx_nodes[1][i1index]),
printf(fh, "  alpha_quad += -(~a)*m_bmag_inv * 1/mag_e_3/Jc_quad * (g_23 * (~a) - g_13 * (~a)); ~%", dH_dz_nodes[3][j1index], dA_dx_nodes[1][i1index], dA_dx_nodes[2][i1index])

apardot contribution:

The apardot contribution is related to the parallel electric field term \dot A_\parallel in the characteristics and is quite straightforward as it does not require any differential operator or geometric terms.