diff --git a/src/ForcingTerms/BiotSavart.jl b/src/ForcingTerms/BiotSavart.jl index c03ab07c65..bfbb245509 100644 --- a/src/ForcingTerms/BiotSavart.jl +++ b/src/ForcingTerms/BiotSavart.jl @@ -3,36 +3,60 @@ 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 @@ -40,32 +64,37 @@ 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 """ @@ -73,16 +102,16 @@ end 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}, @@ -98,25 +127,21 @@ 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 @@ -124,7 +149,7 @@ 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. @@ -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 """ @@ -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}, @@ -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]) @@ -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 @@ -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 diff --git a/test/runtests_coils.jl b/test/runtests_coils.jl index e01d57f46f..4bb102afcf 100644 --- a/test/runtests_coils.jl +++ b/test/runtests_coils.jl @@ -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