Skip to content

ForcingTerms - FEATURE - Use the exact straight-segment field in the coil Biot-Savart sum #590

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
matt-pharr merged 1 commit into from
Oct 26, 2026
Merged

ForcingTerms - FEATURE - Use the exact straight-segment field in the coil Biot-Savart sum #590

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
170 changes: 95 additions & 75 deletions src/ForcingTerms/BiotSavart.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,86 +3,115 @@

Biot-Savart field calculation for 3D coil geometries.

The fundamental formula is (SI units):
Each conductor is a polyline of straight filament segments, and each segment contributes
its exact field (Hanson & Hirshman 2002, Phys. Plasmas 9, 4410, Eq. 8; SI units):

dB⃗ = (μ₀/4π) × I × nw × (dl⃗ × r̂) / |r|²
B⃗ = (μ₀/4π) × I × nw × 2(Rᵢ + R_f) / (Rᵢ R_f ((Rᵢ + R_f)² - L²)) × (dl⃗ × R⃗ᵢ)

where μ₀/4π = 1e-7 [T⋅m/A], I is the current [A], nw is the winding multiplier,
dl⃗ is the current element direction (segment vector), and r⃗ is the displacement
from the segment midpoint to the observation point.

Matches the Fortran `field_bs_psi` discretization in `field.F`.
dl⃗ is the segment vector of length L, and R⃗ᵢ, R⃗_f are the displacements from the
segment start and end to the observation point, with lengths Rᵢ and R_f.
"""

const MU0_OVER_4PI = 1e-7 # T⋅m/A

# Minimum squared distance [m²] to skip in single-point evaluation (r ~ 1 μm); avoids
# singularity when an observation point coincides with a conductor segment midpoint.
# Minimum squared distance [m²] from a segment (r ~ 1 μm) below which its contribution is
# skipped; avoids the singularity when an observation point lies on a conductor.
const BIOT_SAVART_MIN_DIST_SQ = 1e-12

"""
_segment_field(rx, ry, rz, dlx, dly, dlz, len) -> (bx, by, bz)

Exact field of a straight unit-current filament segment of vector `dl⃗` and length `len`,
without the μ₀/4π prefactor, at displacement `r⃗` from the segment start
(Hanson & Hirshman 2002, Eq. 8). With ε = Rᵢ + R_f - L, the effective squared distance
ρ² = 2RᵢR_f ε / (Rᵢ + R_f) reduces to the perpendicular distance squared near the
segment, so the contribution is zero for ρ² below `BIOT_SAVART_MIN_DIST_SQ`.
"""
@inline function _segment_field(rx::Float64, ry::Float64, rz::Float64,
dlx::Float64, dly::Float64, dlz::Float64, len::Float64)
ri = sqrt(rx * rx + ry * ry + rz * rz)
fx = rx - dlx
fy = ry - dly
fz = rz - dlz
rf = sqrt(fx * fx + fy * fy + fz * fz)
sum_r = ri + rf
rho2_num = 2 * ri * rf * (sum_r - len)
# Guard: zero contribution if observer is on conductor (avoids Inf/NaN)
on_wire = rho2_num <= BIOT_SAVART_MIN_DIST_SQ * sum_r
factor = ifelse(on_wire, 0.0, 4 * sum_r / (ifelse(on_wire, 1.0, rho2_num) * (sum_r + len)))
return (dly * rz - dlz * ry) * factor, (dlz * rx - dlx * rz) * factor, (dlx * ry - dly * rx) * factor
end

"""
PrecomputedSegments

Pre-computed segment vectors and midpoints for one conductor strand, stored
Pre-computed segment vectors, lengths and start points for one conductor strand, stored
in contiguous arrays for cache-friendly access in the inner Biot-Savart loop.
"""
struct PrecomputedSegments
dlx::Vector{Float64}
dly::Vector{Float64}
dlz::Vector{Float64}
midx::Vector{Float64}
midy::Vector{Float64}
midz::Vector{Float64}
len::Vector{Float64}
x0::Vector{Float64}
y0::Vector{Float64}
z0::Vector{Float64}
prefactor::Float64 # MU0_OVER_4PI * current * nw
nseg::Int
end

"""
precompute_segments(cs, j, k, current_nw) → PrecomputedSegments

Build contiguous segment-vector and midpoint arrays for strand `(j,k)` of coil set `cs`.
Build contiguous segment-vector, length and start-point arrays for strand `(j,k)` of coil set `cs`.
"""
function precompute_segments(cs::CoilSet, j::Int, k::Int, current_nw::Float64)
nsec = cs.nsec
nseg = nsec - 1
dlx = Vector{Float64}(undef, nseg)
dly = Vector{Float64}(undef, nseg)
dlz = Vector{Float64}(undef, nseg)
midx = Vector{Float64}(undef, nseg)
midy = Vector{Float64}(undef, nseg)
midz = Vector{Float64}(undef, nseg)
dlx = Vector{Float64}(undef, nseg)
dly = Vector{Float64}(undef, nseg)
dlz = Vector{Float64}(undef, nseg)
len = Vector{Float64}(undef, nseg)
x0 = Vector{Float64}(undef, nseg)
y0 = Vector{Float64}(undef, nseg)
z0 = Vector{Float64}(undef, nseg)

@inbounds for l in 1:nseg
x1 = cs.x[j, k, l]; x2 = cs.x[j, k, l+1]
y1 = cs.y[j, k, l]; y2 = cs.y[j, k, l+1]
z1 = cs.z[j, k, l]; z2 = cs.z[j, k, l+1]
dlx[l] = x2 - x1
dly[l] = y2 - y1
dlz[l] = z2 - z1
midx[l] = (x1 + x2) * 0.5
midy[l] = (y1 + y2) * 0.5
midz[l] = (z1 + z2) * 0.5
x1 = cs.x[j, k, l]
x2 = cs.x[j, k, l+1]
y1 = cs.y[j, k, l]
y2 = cs.y[j, k, l+1]
z1 = cs.z[j, k, l]
z2 = cs.z[j, k, l+1]
dlx[l] = x2 - x1
dly[l] = y2 - y1
dlz[l] = z2 - z1
len[l] = sqrt(dlx[l]^2 + dly[l]^2 + dlz[l]^2)
x0[l] = x1
y0[l] = y1
z0[l] = z1
end

return PrecomputedSegments(dlx, dly, dlz, midx, midy, midz,
MU0_OVER_4PI * current_nw, nseg)
return PrecomputedSegments(dlx, dly, dlz, len, x0, y0, z0,
MU0_OVER_4PI * current_nw, nseg)
end

"""
accumulate_strand_field!(Bx, By, Bz, i, obs_x, obs_y, obs_z, xs, ys, zs, current_nw)

Accumulate the Biot-Savart contribution from one conductor strand onto observation point `i`.

Uses midpoint-rule integration: each segment from point `l` to `l+1` contributes:
dB = μ₀/(4π) × current_nw × (dl⃗ × r⃗) / |r|³
where `r⃗` is the displacement from the segment midpoint to the observer.
Each straight segment from point `l` to `l+1` contributes its exact filament field,
see `_segment_field`.

## Arguments
- `Bx`, `By`, `Bz`: pre-allocated output arrays (Cartesian field, Tesla); updated in-place
- `i`: observation point index
- `obs_x/y/z`: observer Cartesian coordinates [m]
- `xs`, `ys`, `zs`: strand point coordinates [nsec], Cartesian [m]
- `current_nw`: effective current = I × nw [A⋅turns]

- `Bx`, `By`, `Bz`: pre-allocated output arrays (Cartesian field, Tesla); updated in-place
- `i`: observation point index
- `obs_x/y/z`: observer Cartesian coordinates [m]
- `xs`, `ys`, `zs`: strand point coordinates [nsec], Cartesian [m]
- `current_nw`: effective current = I × nw [A⋅turns]
"""
function accumulate_strand_field!(
Bx::AbstractVector{Float64},
Expand All @@ -98,33 +127,29 @@ function accumulate_strand_field!(
nsec = length(xs)
prefactor = MU0_OVER_4PI * current_nw

@inbounds for l in 1:(nsec - 1)
@inbounds for l in 1:(nsec-1)
# Segment vector dl⃗
dlx = xs[l+1] - xs[l]
dly = ys[l+1] - ys[l]
dlz = zs[l+1] - zs[l]

# Displacement r⃗ from segment midpoint to observer
rx = obs_x - (xs[l] + xs[l+1]) * 0.5
ry = obs_y - (ys[l] + ys[l+1]) * 0.5
rz = obs_z - (zs[l] + zs[l+1]) * 0.5
# Displacement r⃗ from segment start to observer
rx = obs_x - xs[l]
ry = obs_y - ys[l]
rz = obs_z - zs[l]

r2 = rx*rx + ry*ry + rz*rz
r2 < BIOT_SAVART_MIN_DIST_SQ && continue # avoid singularity (observer on conductor)
inv_r3 = prefactor / (r2 * sqrt(r2))

# dB = prefactor × (dl⃗ × r⃗) / |r|³ [cross product]
Bx[i] += (dly*rz - dlz*ry) * inv_r3
By[i] += (dlz*rx - dlx*rz) * inv_r3
Bz[i] += (dlx*ry - dly*rx) * inv_r3
bx, by, bz = _segment_field(rx, ry, rz, dlx, dly, dlz, sqrt(dlx * dlx + dly * dly + dlz * dlz))
Bx[i] += bx * prefactor
By[i] += by * prefactor
Bz[i] += bz * prefactor
end
end

"""
_accumulate_precomputed(ox, oy, oz, seg) -> (dbx, dby, dbz)

Inner kernel: accumulate Biot-Savart field from pre-computed segments.
Returns `(dbx, dby, dbz)` increments. Matches Fortran `field_bs_psi` midpoint rule.
Returns `(dbx, dby, dbz)` increments using the exact segment field.

Operates on `PrecomputedSegments` with contiguous arrays for maximum
cache efficiency.
Expand All @@ -139,22 +164,14 @@ cache efficiency.
pf = seg.prefactor

@inbounds @simd for l in 1:seg.nseg
rx = ox - seg.midx[l]
ry = oy - seg.midy[l]
rz = oz - seg.midz[l]

r2 = rx*rx + ry*ry + rz*rz
# Guard: zero contribution if observer is on conductor (avoids Inf/NaN)
safe_r2 = ifelse(r2 < BIOT_SAVART_MIN_DIST_SQ, one(r2), r2)
scale = ifelse(r2 < BIOT_SAVART_MIN_DIST_SQ, zero(pf), pf)
inv_r3 = scale / (safe_r2 * sqrt(safe_r2))

bx_acc += (seg.dly[l]*rz - seg.dlz[l]*ry) * inv_r3
by_acc += (seg.dlz[l]*rx - seg.dlx[l]*rz) * inv_r3
bz_acc += (seg.dlx[l]*ry - seg.dly[l]*rx) * inv_r3
bx, by, bz = _segment_field(ox - seg.x0[l], oy - seg.y0[l], oz - seg.z0[l],
seg.dlx[l], seg.dly[l], seg.dlz[l], seg.len[l])
bx_acc += bx
by_acc += by
bz_acc += bz
end

return bx_acc, by_acc, bz_acc
return bx_acc * pf, by_acc * pf, bz_acc * pf
end

"""
Expand All @@ -168,13 +185,14 @@ returned field components are also cylindrical.
Uses `Threads.@threads` over observation points -- each thread writes to a unique
index so no synchronization is required.

Performance: segment vectors and midpoints are pre-computed once outside the
Performance: segment vectors, lengths and start points are pre-computed once outside the
observation-point loop, then accessed from contiguous arrays in the inner kernel.

## Arguments
- `B_R`, `B_phi`, `B_Z`: output field arrays [Tesla], length `nobs`; overwritten
- `obs_R`, `obs_phi`, `obs_Z`: observation point cylindrical coordinates [m, rad, m]
- `coil_sets`: conductor geometry and current data

- `B_R`, `B_phi`, `B_Z`: output field arrays [Tesla], length `nobs`; overwritten
- `obs_R`, `obs_phi`, `obs_Z`: observation point cylindrical coordinates [m, rad, m]
- `coil_sets`: conductor geometry and current data
"""
function compute_biot_savart_boundary!(
B_R::AbstractVector{Float64},
Expand Down Expand Up @@ -204,8 +222,8 @@ function compute_biot_savart_boundary!(
# Pre-compute sincos for cylindrical ↔ Cartesian conversion (avoids redundant trig)
cos_phi = Vector{Float64}(undef, nobs)
sin_phi = Vector{Float64}(undef, nobs)
obs_x = Vector{Float64}(undef, nobs)
obs_y = Vector{Float64}(undef, nobs)
obs_x = Vector{Float64}(undef, nobs)
obs_y = Vector{Float64}(undef, nobs)

@inbounds for i in 1:nobs
sp, cp = sincos(obs_phi[i])
Expand All @@ -218,7 +236,9 @@ function compute_biot_savart_boundary!(
# Each thread writes exclusively to its own index i — no synchronization needed
Threads.@threads for i in 1:nobs
@inbounds begin
ox = obs_x[i]; oy = obs_y[i]; oz = obs_Z[i]
ox = obs_x[i]
oy = obs_y[i]
oz = obs_Z[i]
bx_total = 0.0
by_total = 0.0
bz_total = 0.0
Expand All @@ -233,9 +253,9 @@ function compute_biot_savart_boundary!(
# Convert Cartesian → cylindrical in-place
cp = cos_phi[i]
sp = sin_phi[i]
B_R[i] = bx_total * cp + by_total * sp
B_R[i] = bx_total * cp + by_total * sp
B_phi[i] = -bx_total * sp + by_total * cp
B_Z[i] = bz_total
B_Z[i] = bz_total
end
end
end
Expand Down
54 changes: 54 additions & 0 deletions test/runtests_coils.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -122,6 +122,60 @@ end
@test abs(By[1]) < 1e-14
end

@testset "BiotSavart: straight segment field is exact" begin
# Segment from (0,0,-0.5) to (0,0,0.5) m, unit current. Textbook finite-wire field at
# distance d, height z: B_φ = μ₀I/(4πd) × ((L/2 - z)/√(d² + (L/2 - z)²) + (L/2 + z)/√(d² + (L/2 + z)²)).
xs = [0.0, 0.0]
ys = [0.0, 0.0]
zs = [-0.5, 0.5]
for d in (0.05, 0.25, 0.5, 2.0), z in (0.0, 0.2, 0.9)
Bx = zeros(1)
By = zeros(1)
Bz = zeros(1)
ForcingTerms.accumulate_strand_field!(Bx, By, Bz, 1, d, 0.0, z, xs, ys, zs, 1.0)
a = 0.5 - z
b = 0.5 + z
B_exact = 1e-7 / d * (a / sqrt(d^2 + a^2) + b / sqrt(d^2 + b^2))
@test By[1] ≈ B_exact rtol = 1e-12
@test abs(Bx[1]) < 1e-14 * B_exact && abs(Bz[1]) < 1e-14 * B_exact
end
@inferred ForcingTerms._segment_field(1.0, 0.2, 0.3, 0.0, 0.0, 1.0, 1.0)

# Observer on the conductor: the segment contribution is skipped, never Inf/NaN
Bx = zeros(1)
By = zeros(1)
Bz = zeros(1)
ForcingTerms.accumulate_strand_field!(Bx, By, Bz, 1, 0.0, 0.0, 0.1, xs, ys, zs, 1.0)
@test Bx[1] == 0 && By[1] == 0 && Bz[1] == 0
end

@testset "BiotSavart: field independent of polyline subdivision" begin
# 0.6 m × 0.4 m rectangle in the plane x = 1.25 m with each side split into m collinear
# pieces; the conductor is the same for every m, so the field must be too.
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 = [cy[s] + p / m * (cy[s+1] - cy[s]) for s in 1:4 for p in 0:(m-1)]
zs = [cz[s] + p / m * (cz[s+1] - cz[s]) for s in 1:4 for p in 0:(m-1)]
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
obs_R = [r for r in 0.95:0.05:1.15 for z in -0.3:0.15:0.3]
obs_Z = [z for r in 0.95:0.05:1.15 for z in -0.3:0.15:0.3]
obs_phi = fill(0.1, length(obs_R))
B = map((1, 16)) do m
B_R = zeros(length(obs_R))
B_phi = zeros(length(obs_R))
B_Z = zeros(length(obs_R))
ForcingTerms.compute_biot_savart_boundary!(B_R, B_phi, B_Z, obs_R, obs_phi, obs_Z, [rectangle(m)])
hcat(B_R, B_phi, B_Z)
end
@test maximum(abs, B[1] .- B[2]) < 1e-12 * maximum(abs, B[2])
end

@testset "BiotSavart: compute_biot_savart_boundary! threading" begin
# Same circular loop; verify threaded result matches single-obs result
a = 0.8
Expand Down
Loading