[krystophny]

Fixes #360.

field_bs_psi and field_bs_rzphi evaluated each straight polyline segment once, at its midpoint (dx x r/|r|^3), which has an O((L/d)^2) error for a segment of length L at distance d. Each segment now contributes its exact closed-form field (Hanson and Hirshman, Phys. Plasmas 9, 4410 (2002)) through a new pure function field_bs_segment, so the field of the polygon in the coil file is exact at any distance. Interfaces, units, coil file format and outputs are unchanged. The loops keep their structure, with xa now holding the segment start instead of the midpoint.

In field_bs_rzphi the former distance floor dlmin = 2e-3*cro is kept as the radius of a regularizing wire: within dlmin of a filament the field is bounded by mu0 I/(2 pi dlmin) and goes to zero on the filament instead of becoming singular. field_bs_psi had no floor and still has none.

Test

regression/test/coil_segment_field.f90 (build command in its header, no input files). It checks

1. field_bs_segment against the textbook closed form mu0 I/(4 pi rho) (cos theta1 - cos theta2) at eight points from near the segment to far beyond its ends: max relative error 5.9e-15 (tolerance 1e-12);
2. field_bs_rzphi for a rectangular coil given as its four sides and with each side split into 16 pieces, which is the same conductor: max relative n = 1 difference 2.3e-16 (tolerance 1e-12). On develop the same comparison differs by 0.11 (reproducer in COIL - FEATURE - Exact straight-segment field instead of the midpoint rule per coil segment #360);
3. that points within rmin of the filament give a finite field below the wire bound.

It also passes with the debug flags of the gfortran CI job (-O0 -fcheck=all -fsanitize=bounds -finit-real=snan ...) with no new warnings in coil/field.F.

Before and after

	develop 5499124	this branch
reproducer of #360, 4 sides vs 1024 pieces	1.1e-1	6.8e-16
MAST-U correction coils (our case), resonant vacuum flux / independent Biot-Savart, q = 2-7	1.130-1.137	0.9975-1.0005
same, coil files densified 16x	0.998-1.001	identical to the undensified run to 1.4e-12
regression_solovev_ideal_example (finely discretized coil)		outputs change by at most 1.2e-4 relative

The remaining 0.25 % in the MAST-U row is the difference between the codes after the coil field, the same as with densified coils on develop.

Cost: one segment evaluation is about twice as expensive (two square roots instead of one, 4.9 s vs 9.8 s single-threaded for 2.5e8 segment-point pairs in field_bs_rzphi). Reaching 0.3 % accuracy with the midpoint rule needed 16x more segments in the case above.

Built with the native Make (make, gfortran 16, -O3 -fopenmp); regression_solovev_ideal_example runs to normal termination. The CI builds with gfortran and ifort pass. Not run: the kinetic and resistive Soloviev examples, which use the same coil routines.

Chris&AI

[matt-pharr]

@krystophny I would generally classify this as a feature add, not a bugfix, as the previous definition of Biot-Savart fields is still correct and does converge with sufficient resolution :)

This is handy, though, and partially resolves a backburner item I have wanted to look at for a long time (that being investigating the convergence of vacuum source fields).

I would be interested in investigation or analysis of how this affects convergence for round coils; EFC coils are frequently straight or slightly round, but we also use this Biot-Savart evaluator with TF/PF/CS coils and feeds. Does this new formula impact the convergence of nearby fields when a coil is circular?

[krystophny]

Agreed, this is an accuracy improvement rather than a bug fix, since the midpoint rule converges with more segments. I will retitle #360 and this pull request as a feature, and the same for the Julia port (#589, #590).

For a round coil both versions converge as N⁻² once the segments are short compared with the distance d from the conductor, because the remaining error is mostly the polygon replacing the circle. The two differ close to the coil, where the segment length is comparable to d: the midpoint rule (each segment as one current element at its midpoint) then overshoots, while the exact field of each straight segment does not. The plot compares both with the analytic field of a thin circular loop of radius a, as the largest relative error of |B| over azimuth in the coil plane, at distance d outside the wire:

• At d = a the exact formula is 10 to 30 % more accurate, so for TF, PF and CS coils seen from the plasma little changes.
• At d = 0.1 a it reaches 1 % at 64 segments, where the midpoint rule is still at 3 %.
• At d = 0.02 a, for example near a feed, the midpoint rule is 40 % off at 128 segments and the exact formula about 1 %.

The script and data are in this branch.

As a rule of thumb, the midpoint rule is fine once the segment length is well below d. For a straight coil side the exact formula has no discretisation error at all, whatever the number of segments.

Chris&AI

[matt-pharr]

This is a useful change and already reviewed in Julia. The formulae here match appropriately and are accurate for all applications of this subroutine that I am aware of. Approve