ForcingTerms - FEATURE - Exact straight-segment field instead of a midpoint current element per polyline segment #589

src/ForcingTerms/BiotSavart.jl replaces every straight segment of a coil polyline by a single current element at its midpoint, both in accumulate_strand_field! (lines 101-120) and in the threaded kernel _accumulate_precomputed used by compute_biot_savart_boundary! (lines 140-155):

r = x_obs - (x[l] + x[l+1])/2,   dl = x[l+1] - x[l]
B += μ₀ I nw/(4π) * dl × r / |r|^3

This is a port of the same discretisation in the Fortran coil/field.F (field_bs_psi, field_bs_rzphi), reported in GPEC#360 and fixed in GPEC#361. It is the one-point midpoint quadrature of the Biot-Savart line integral along the segment, so its relative error is of order (L/d)^2 for a segment of length L seen from a distance d. For a single 1 m segment the midpoint |B| is 21 % too large at d/L = 0.5, 9.5 % at d/L = 1 and 2.9 % at d/L = 2. Coil files with long straight sides, such as window-frame or correction coils given by their corners, therefore produce a wrong field near the plasma boundary, although the polygon in the coil file has a closed-form field per segment.

Reproducer

A 0.6 m × 0.4 m rectangle in the plane x = 1.25 m, I = 1 kA, field on a grid R = 0.95-1.15 m, Z = -0.3-0.3 m, φ = -0.5-0.5 rad (225 points). Each side is passed as m collinear pieces, so the conductor is the same for every m and the field must not depend on m.

using GeneralizedPerturbedEquilibrium
using GeneralizedPerturbedEquilibrium.ForcingTerms
using Printf

function rectangle(m)
    cy = [-0.3, 0.3, 0.3, -0.3, -0.3]
    cz = [-0.2, -0.2, 0.2, 0.2, -0.2]
    ys = Float64[]; zs = Float64[]
    for side in 1:4, p in 0:m-1
        t = p / m
        push!(ys, cy[side] + t * (cy[side+1] - cy[side]))
        push!(zs, cz[side] + t * (cz[side+1] - cz[side]))
    end
    push!(ys, cy[5]); push!(zs, cz[5])
    n = length(ys)
    return ForcingTerms.CoilSet("rectangle", 1, 1, 1.0, n, reshape(fill(1.25, n), 1, 1, :),
        reshape(ys, 1, 1, :), reshape(zs, 1, 1, :), [1.0e3])
end

R = Float64[]; P = Float64[]; Z = Float64[]
for r in 0.95:0.05:1.15, z in -0.3:0.15:0.3, phi in range(-0.5, 0.5; length=9)
    push!(R, r); push!(Z, z); push!(P, phi)
end
field(m) = (BR = zeros(length(R)); BP = zero(BR); BZ = zero(BR);
            ForcingTerms.compute_biot_savart_boundary!(BR, BP, BZ, R, P, Z, [rectangle(m)]); hcat(BR, BP, BZ))
B0 = field(256); scale = maximum(abs, B0)
println("  m  max|B_m - B_256| / max|B_256|")
for m in (1, 2, 4, 8, 16)
    @printf("%4d %12.3e\n", m, maximum(abs, field(m) .- B0) / scale)
end

Output on develop 81fb41d (Julia 1.13.2):

  m  max|B_m - B_256| / max|B_256|
   1    1.437e+00
   2    3.396e-01
   4    4.230e-02
   8    2.288e-03
  16    5.690e-04

With the corners alone (m = 1) the field is wrong by more than 100 % of its maximum at the grid points closest to the coil, and the error falls roughly as 1/m^2, as expected for the midpoint rule. The existing circular-loop test does not see this, because its 500 segments are short compared with the distance to the observer.

Proposed fix

Use the exact field of a straight filament segment (Hanson and Hirshman, Phys. Plasmas 9, 4410 (2002)),

B = μ₀ I/(4π) * 2 (Ri + Rf) / (Ri Rf ((Ri + Rf)^2 - L^2)) * (dl × R_i)

with R_i and R_f the vectors from the segment start and end to the observer, as in the Fortran fix. A pull request follows.

Chris&AI