Ex, Ey, Bx, By

Theory

We want to solve

\[ \begin{align}\begin{aligned}\Delta_\perp E_x &= \frac{\partial \rho}{\partial x} - \frac{\partial j_x}{\partial \xi}\\\Delta_\perp E_y &= \frac{\partial \rho}{\partial y} - \frac{\partial j_y}{\partial \xi}\\\Delta_\perp B_x &= \frac{\partial j_y}{\partial \xi} - \frac{\partial j_z}{\partial y}\\\Delta_\perp B_y &= \frac{\partial j_z}{\partial x} - \frac{\partial j_x}{\partial \xi}\end{aligned}\end{align} \]

with mixed boundary conditions.

Todo

DOCS: specify the boundary conditions.

Unfortunately, what we actually solve is no less than three steps away from these.

Helmholtz equations

The harsh reality of numerical stability forces us to solve these Helmholtz equations instead:

\[ \begin{align}\begin{aligned}\Delta_\perp E_x - E_x &= \frac{\partial \rho}{\partial x} - \frac{\partial j_x}{\partial \xi} - E_x\\\Delta_\perp E_y - E_y &= \frac{\partial \rho}{\partial y} - \frac{\partial j_y}{\partial \xi} - E_y\\\Delta_\perp B_x - B_x &= \frac{\partial j_y}{\partial \xi} - \frac{\partial j_z}{\partial y} - B_x\\\Delta_\perp B_y - B_y &= \frac{\partial j_z}{\partial x} - \frac{\partial j_x}{\partial \xi} - B_y\end{aligned}\end{align} \]

Note

The behaviour is actually configurable with config.field_solver_subtraction_trick (what a mouthful). 0 or False corresponds to Laplace equation, and while any floating-point values or whole matrices of them should be accepted, it’s recommended to simply use 1 or True instead.

config_example.field_solver_subtraction_trick = 1

0 for Laplace eqn., Helmholtz otherwise

Method

The algorithm can be succinctly written as iMIX2D(mixed_matrix * MIX2D(RHS)), where MIX2D and iMIX2D are Type-1 Forward and Inverse Discrete Trasforms, Sine in one direction and Cosine in the other. RHS is the right-hand side of the equiation above, and dirichlet_matrix is a ‘magical’ matrix that does all the work.

lcode.mixed_matrix(grid_steps, grid_step_size, subtraction_trick)

Calculate a magical matrix that solves the Helmholtz or Laplace equation (subtraction_trick=True and subtraction_trick=False correspondingly) if you elementwise-multiply the RHS by it “in DST-DCT-transformed-space”. See Samarskiy-Nikolaev, p. 189 and around.

Todo

DOCS: expand with method description (Kargapolov, Shalimova)

lcode.calculate_Ex_Ey_Bx_By(config, Ex_avg, Ey_avg, Bx_avg, By_avg, beam_ro, ro, jx, jy, jz, jx_prev, jy_prev)

Calculate transverse fields as iDST-DCT(mixed_matrix * DST-DCT(RHS.T)).T, with and without transposition depending on the field component.

Note that some outer cells do not participate in the calculations, and the result is simply padded with zeroes in the end. We don’t define separate functions for separate boundary condition types and simply transpose the input and output data.

DST-DCT Hybrid

lcode.mix2d(a)

Calculate a DST-DCT-hybrid transform (DST in first direction, DCT in second one), jury-rigged from padded rFFT (anti-symmetrically in first direction, symmetrically in second direction).

As cupy currently ships no readily available function for calculating even 1D DST/DCT on the GPU, we, once again, roll out our own FFT-based implementation.

We don’t need a separate function for the inverse transform, as it matches the forward one up to the normalization multiplier, which is taken into account in mixed_matrix().

Variant B

But wait, the complications don’t stop here.

While we do have a succesfully implemented \((\Delta_\perp - 1)\) inverse operator, there’s still an open question of supplying an unknown value to the RHS.

The naive version (internally known as ‘Variant B’) is to pass the best known substitute to date, i.e. previous layer fields at the predictor phase and the averaged fields at the corrector phase. \(\xi\)-derivatives are taken at half-steps, transverse derivatives are averaged at half-steps or taken from the previous layer if not available.

\[ \begin{align}\begin{aligned}(\Delta_\perp - 1) E_x^{next} &= \frac{\partial \rho^{prev}}{\partial x} - (\frac{\partial j_x}{\partial \xi})^{\mathrm{halfstep}} - E_x^{avg}\\(\Delta_\perp - 1) E_y^{next} &= \frac{\partial \rho^{prev}}{\partial y} - (\frac{\partial j_y}{\partial \xi})^{\mathrm{halfstep}} - E_y^{avg}\\(\Delta_\perp - 1) B_x^{next} &= (\frac{\partial j_y}{\partial \xi})^{\mathrm{halfstep}} - \frac{\partial j_z^{prev}}{\partial y} - B_x^{avg}\\(\Delta_\perp - 1) B_y^{next} &= \frac{\partial j_z^{prev}}{\partial x} - (\frac{\partial j_x}{\partial \xi})^{\mathrm{halfstep}} - B_y^{avg}\end{aligned}\end{align} \]

Variant A

The more correct version (known as ‘Variant A’) mutates the equations once again to take everything at half-steps:

\[ \begin{align}\begin{aligned}(\Delta_\perp - 1) E_x^{\mathrm{halfstep}} &= \frac{\partial \rho^{\mathrm{avg}}}{\partial x} - (\frac{\partial j_x}{\partial \xi})^{\mathrm{halfstep}} - E_x^{\mathrm{avg}}\\(\Delta_\perp - 1) E_y^{\mathrm{halfstep}} &= \frac{\partial \rho^{\mathrm{avg}}}{\partial y} - (\frac{\partial j_y}{\partial \xi})^{\mathrm{halfstep}} - E_y^{\mathrm{avg}}\\(\Delta_\perp - 1) B_x^{\mathrm{halfstep}} &= (\frac{\partial j_y}{\partial \xi})^{\mathrm{halfstep}} - \frac{\partial j_z^{\mathrm{avg}}}{\partial y} - B_x^{\mathrm{avg}}\\(\Delta_\perp - 1) B_y^{\mathrm{halfstep}} &= \frac{\partial j_z^{\mathrm{avg}}}{\partial x} - (\frac{\partial j_x}{\partial \xi})^{\mathrm{halfstep}} - B_y^{\mathrm{avg}}\end{aligned}\end{align} \]

and calculates the fields at next step in the following fashion: \(E_x^{next} = 2 E_x^{avg} - E_x^{prev}\), e.t.c.

Solving these is equivalent to solving Variant B equations with averaged fields, \(\rho\) and \(j_z\) and applying the above transformation to the result. See lcode.step() for the wrapping code that does that.

config_example.field_solver_variant_A = True

Use Variant A or Variant B for Ex, Ey, Bx, By