Bz

Equations

We want to solve

\[\Delta_\perp B_z = \frac{\partial j_x}{\partial y} - \frac{\partial j_y}{\partial x}\]

with Neumann boundary conditions (derivative = 0).

Method

The algorithm can be succinctly written as iDCT2D(neumann_matrix * DCT2D(RHS)), where DCT2D and iDCT2D are Type-1 Forward and Inverse Discrete Sine 2D Trasforms respectively, RHS is the right-hand side of the equiation above, and neumann_matrix is a ‘magical’ matrix that does all the work.

lcode.neumann_matrix(grid_steps, grid_step_size)

Calculate a magical matrix that solves the Laplace equation if you elementwise-multiply the RHS by it “in DST-space”. See Samarskiy-Nikolaev, p. 187.

In addition to the magic values, it also hosts the DCT normalization multiplier.

Todo

DOCS: expand with method description (Kargapolov, Shalimova)

lcode.calculate_Bz(config, jx, jy)

Calculate Bz as iDCT2D(dirichlet_matrix * DCT2D(djx/dy - djy/dx)).

Note that this time the outer cells do not participate in the calculations, so the RHS derivatives are padded with zeroes in the beginning.

DCT2D

lcode.dct2d(a)

Calculate DCT-Type1-2D, jury-rigged from symmetrically-padded rFFT.

As cupy currently ships no readily available function for calculating the DCT2D on the GPU, we roll out our own FFT-based implementation.

We don’t need to make a separate iDCT2D function as (for Type-1) it matches DCT2D up to the normalization multiplier, which is taken into account in neumann_matrix().