Numerics

This section describes how the equations are discretized and evaluated.

Grid Conventions

The surface is parameterized by toroidal angle phi and poloidal angle theta. SIMSOPT and the reference BIEST code use normalized angles of period 1, not 2*pi. The grids are uniform:

  • theta in [0, 1) with ntheta points.

  • phi in [0, 1/nfp) with nphi points for a field period.

For stellarator symmetry (half_period=True), the toroidal grid is shifted by half a grid point. This is essential for spectral accuracy with uniform weights.

Spectral Resampling

The core surface operators are Fourier based:

  • Upsample: zero-pad Fourier coefficients.

  • Resample: upsample then decimate.

  • Grad2D: spectral differentiation in both angles.

  • RotateToroidal: phase shift in Fourier space.

The JAX implementation uses unitary FFTs to match the normalization in SCTL. The reference FFT wrapper scales by 1/sqrt(N) on both forward and inverse transforms.

Surface Differentiation and Normals

Surface derivatives are computed spectrally. For a Fourier mode exp(2*pi*i*(m*theta + n*phi)):

  • partial_theta multiplies the coefficient by -2*pi*i*m.

  • partial_phi multiplies the coefficient by -2*pi*i*n.

The negative sign matches the FFT sign convention used in BIEST’s Grad2D implementation.

Toroidal frequencies use signed indices m with the Nyquist index treated as positive, matching the reference implementation:

m = t for t <= Nt/2 and m = t - Nt otherwise.

Given X_theta and X_phi, the unit normal and area element are:

\[\begin{split}n = \\frac{X_{\\theta} \\times X_{\\phi}} {\\lVert X_{\\theta} \\times X_{\\phi} \\rVert}, \\quad dA = \\frac{\\lVert X_{\\theta} \\times X_{\\phi} \\rVert}{N},\end{split}\]

with N = Nt * Np grid points. The orientation is chosen so that the normal component corresponding to the maximum coordinate among x,y,z is positive, matching BIEST’s convention.

Field-Period Target Selection

The BIEST field-period operator evaluates on a subset of the quadrature grid. Given:

  • quad_nt / quad_np: quadrature grid sizes for the full surface.

  • trg_nt / trg_np: target grid sizes for one field period.

  • nfp: number of field periods.

The target indices are selected by uniform strides:

\[\begin{split}\\text{skip}_t = \\frac{\\text{quad\\_nt}}{n_{fp}\\,\\text{trg\\_nt}}, \\quad \\text{skip}_p = \\frac{\\text{quad\\_np}}{\\text{trg\\_np}}\end{split}\]

and (i,j) maps to the quadrature index (t,p) = (i * skip_t, j * skip_p) with i \\in [0, trg_nt). This picks points from the first field period only.

Singular Quadrature

The boundary integrals involve kernels singular at r = r'. BIEST uses a high-order scheme:

  1. Partition of unity (POU) on a local patch around each target point.

  2. The non-singular part is integrated by the trapezoidal rule.

  3. The singular part is evaluated in local polar coordinates with a high-order radial quadrature.

The scheme is described in detail in Malhotra et al. (2020), and the implementation in this repository mirrors the reference code in biest/singular_correction.hpp.

Implementation Notes

The JAX port implements the partition-of-unity correction for Laplace FxdU (Hedgehog order 1) and Laplace Fxd2U (Hedgehog order > 1). The patch size is chosen using the same thresholding rules as BIEST:

\[\begin{split}\\text{PDIM} = \\lfloor 1.6\\,\\text{digits}\\,\\text{cond} \\rfloor\end{split}\]

and then rounded up to the nearest supported template value {6, 8, 12, 16, ..., 64}. The base polar quadrature order is RAD_DIM_0 = \\lfloor 1.6\\,\\text{PDIM} \\rfloor.

For Laplace Fxd2U, BIEST uses Hedgehog quadrature with a tripled radial resolution while keeping the angular resolution tied to the base order:

\[\begin{split}\\text{RAD\\_DIM} = 3\\,\\text{RAD\\_DIM}_0, \\quad \\text{ANG\\_DIM} = 2\\,\\text{RAD\\_DIM}_0\end{split}\]

The JAX implementation mirrors this convention and uses the same hedgehog extrapolation weights at nodes 1..16 when HedgehogOrder=8.

For quadrature selection, the JAX port also implements the singular-corrected Laplace DxU (double-layer) operator with Hedgehog order 1. This matches the BIEST self-test that verifies the expected 1/2 jump for constant density and drives adaptive quadrature refinement.

Adaptive Quadrature Resolution

The field-period operator chooses a quadrature resolution quad_Nt quad_Np based on:

  • Surface anisotropy estimates.

  • A double-layer self-test that checks the expected 1/2 jump for constant density.

This adaptive selection is replicated in the JAX implementation to ensure parity.