Skip to content

COIL - BUG FIX - Use the exact field of each straight coil segment (PrincetonUniversity/GPEC#361) #9

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
krystophny merged 1 commit into from
Nov 2, 2026
Merged

COIL - BUG FIX - Use the exact field of each straight coil segment (PrincetonUniversity/GPEC#361) #9

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
85 changes: 57 additions & 28 deletions coil/field.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,8 @@
c-----------------------------------------------------------------------
c 0. field_mod.
c 1. field_bs_psi.
c 2. field_bs_rz.
c 2. field_bs_rzphi.
c 3. field_bs_segment.
c-----------------------------------------------------------------------
c subprogram 0. field_mod.
c module declarations.
Expand Down Expand Up @@ -137,12 +138,9 @@ SUBROUTINE field_bs_psi(psi,bmn,wegt,op_start,op_stop, op_verbose)
DO j=1,coil(i)%ncoil
DO k=1,coil(i)%s
DO iseg=1,nseg
xa(j,k,iseg)=(coil(i)%x(j,k,iseg+1)+
$ coil(i)%x(j,k,iseg))/2.0
ya(j,k,iseg)=(coil(i)%y(j,k,iseg+1)+
$ coil(i)%y(j,k,iseg))/2.0
za(j,k,iseg)=(coil(i)%z(j,k,iseg+1)+
$ coil(i)%z(j,k,iseg))/2.0
xa(j,k,iseg)=coil(i)%x(j,k,iseg)
ya(j,k,iseg)=coil(i)%y(j,k,iseg)
za(j,k,iseg)=coil(i)%z(j,k,iseg)

dx(j,k,iseg)=coil(i)%x(j,k,iseg+1)-coil(i)%x(j,k,iseg)
dy(j,k,iseg)=coil(i)%y(j,k,iseg+1)-coil(i)%y(j,k,iseg)
Expand Down Expand Up @@ -202,13 +200,12 @@ SUBROUTINE field_bs_psi(psi,bmn,wegt,op_start,op_stop, op_verbose)
ry=yobs(itheta,izeta)-ya(j,k,iseg)
rz=zobs(itheta,izeta)-za(j,k,iseg)

dl=SQRT(rx**2+ry**2+rz**2)
dbx=dbx+coil(i)%cur(j)*coil(i)%nw*1e-7/
$ dl**3*(dy(j,k,iseg)*rz-dz(j,k,iseg)*ry)
dby=dby+coil(i)%cur(j)*coil(i)%nw*1e-7/
$ dl**3*(dz(j,k,iseg)*rx-dx(j,k,iseg)*rz)
dbz=dbz+coil(i)%cur(j)*coil(i)%nw*1e-7/
$ dl**3*(dx(j,k,iseg)*ry-dy(j,k,iseg)*rx)
dl=coil(i)%cur(j)*coil(i)%nw*1e-7*
$ field_bs_segment(rx,ry,rz,dx(j,k,iseg),
$ dy(j,k,iseg),dz(j,k,iseg),0.0_r8)
dbx=dbx+dl*(dy(j,k,iseg)*rz-dz(j,k,iseg)*ry)
dby=dby+dl*(dz(j,k,iseg)*rx-dx(j,k,iseg)*rz)
dbz=dbz+dl*(dx(j,k,iseg)*ry-dy(j,k,iseg)*rx)

ENDDO
ENDDO
Expand Down Expand Up @@ -356,12 +353,9 @@ SUBROUTINE field_bs_rzphi(nr,nz,np,gdr,gdz,cbr,cbz,cbp,op_verbose)
DO j=1,coil(i)%ncoil
DO k=1,coil(i)%s
DO iseg=1,nseg
xa(j,k,iseg)=(coil(i)%x(j,k,iseg+1)+
$ coil(i)%x(j,k,iseg))/2.0
ya(j,k,iseg)=(coil(i)%y(j,k,iseg+1)+
$ coil(i)%y(j,k,iseg))/2.0
za(j,k,iseg)=(coil(i)%z(j,k,iseg+1)+
$ coil(i)%z(j,k,iseg))/2.0
xa(j,k,iseg)=coil(i)%x(j,k,iseg)
ya(j,k,iseg)=coil(i)%y(j,k,iseg)
za(j,k,iseg)=coil(i)%z(j,k,iseg)

dx(j,k,iseg)=coil(i)%x(j,k,iseg+1)-coil(i)%x(j,k,iseg)
dy(j,k,iseg)=coil(i)%y(j,k,iseg+1)-coil(i)%y(j,k,iseg)
Expand Down Expand Up @@ -399,14 +393,12 @@ SUBROUTINE field_bs_rzphi(nr,nz,np,gdr,gdz,cbr,cbz,cbp,op_verbose)
ry=yobs(ir,iz,ip)-ya(j,k,iseg)
rz=zobs(ir,iz,ip)-za(j,k,iseg)

dl=MAXVAL((/SQRT(rx**2+ry**2+rz**2),dlmin/))

dbx=dbx+coil(i)%cur(j)*coil(i)%nw*1e-7/
$ dl**3*(dy(j,k,iseg)*rz-dz(j,k,iseg)*ry)
dby=dby+coil(i)%cur(j)*coil(i)%nw*1e-7/
$ dl**3*(dz(j,k,iseg)*rx-dx(j,k,iseg)*rz)
dbz=dbz+coil(i)%cur(j)*coil(i)%nw*1e-7/
$ dl**3*(dx(j,k,iseg)*ry-dy(j,k,iseg)*rx)
dl=coil(i)%cur(j)*coil(i)%nw*1e-7*
$ field_bs_segment(rx,ry,rz,dx(j,k,iseg),
$ dy(j,k,iseg),dz(j,k,iseg),dlmin)
dbx=dbx+dl*(dy(j,k,iseg)*rz-dz(j,k,iseg)*ry)
dby=dby+dl*(dz(j,k,iseg)*rx-dx(j,k,iseg)*rz)
dbz=dbz+dl*(dx(j,k,iseg)*ry-dy(j,k,iseg)*rx)

ENDDO
ENDDO
Expand Down Expand Up @@ -459,6 +451,43 @@ SUBROUTINE field_bs_rzphi(nr,nz,np,gdr,gdz,cbr,cbz,cbp,op_verbose)
c-----------------------------------------------------------------------
RETURN
END SUBROUTINE field_bs_rzphi
c-----------------------------------------------------------------------
c subprogram 3. field_bs_segment.
c
c *DESCRIPTION:
c Exact Biot-Savart field of a straight filament segment
c (Hanson & Hirshman, Phys. Plasmas 9, 4410 (2002)).
c The field of a current cur along the segment from xs to
c xs+(dx,dy,dz) at the point xs+(rx,ry,rz) is
c mu0*cur/(4*pi)*field_bs_segment*((dx,dy,dz) x (rx,ry,rz)).
c Far from the segment the factor tends to 1/|r|**3 of the
c midpoint rule. A positive rmin regularizes the field like a
c wire of radius rmin, with |B| bounded by mu0*cur/(2*pi*rmin).
c
c *ARGUMENTS:
c rx,ry,rz : real
c Observation point relative to the segment start.
c dx,dy,dz : real
c Segment vector from its start to its end.
c rmin : real
c Regularization radius, 0 for none.
c-----------------------------------------------------------------------
PURE REAL(r8) FUNCTION field_bs_segment(rx,ry,rz,dx,dy,dz,rmin)

REAL(r8), INTENT(IN) :: rx,ry,rz,dx,dy,dz,rmin
REAL(r8) :: ri,rf,sl,eps

sl=SQRT(dx**2+dy**2+dz**2)
ri=SQRT(rx**2+ry**2+rz**2)
rf=SQRT((rx-dx)**2+(ry-dy)**2+(rz-dz)**2)
IF(rmin>0)THEN
ri=MAX(ri,rmin)
rf=MAX(rf,rmin)
ENDIF
eps=ri+rf-sl
IF(rmin>0)eps=MAX(eps,rmin**2*(ri+rf)/(2*ri*rf))
field_bs_segment=2*(ri+rf)/(ri*rf*eps*(ri+rf+sl))
END FUNCTION field_bs_segment

END MODULE field_mod

Expand Down
123 changes: 123 additions & 0 deletions regression/test/coil_segment_field.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,123 @@
! Biot-Savart field of polygonal coils (coil/field.F).
! 1. field_bs_segment against the textbook closed form of a straight segment,
! B = mu0 I/(4 pi rho) (cos(theta1) - cos(theta2)), from near to far points.
! 2. field_bs_rzphi for a rectangular coil passed as its 4 sides and with each
! side split into 16 collinear pieces. The conductor is the same, so the
! n = 1 field must agree (the former midpoint rule differs by 10 %).
! 3. With rmin > 0, points within rmin of the filament give a finite field
! bounded by that of a wire of radius rmin.
! Build after `make` in install/ (gfortran):
! gfortran -fopenmp -I../coil -I../equil ../regression/test/coil_segment_field.f90 \
! -L../deps/lib -lcoil -lequil -lopenblas -lnetcdff -lnetcdf -o coil_segment_field
program coil_segment_field
use coil_mod, only: coil, coil_num, cro, cnn, r8
use field_mod, only: field_bs_rzphi, field_bs_segment
implicit none
integer, parameter :: nr = 4, nz = 4, np = 128
real(r8), parameter :: a(3) = [0.1_r8, -0.2_r8, -0.5_r8]
real(r8), parameter :: d(3) = [0.2_r8, 0.1_r8, 1.0_r8]
real(r8) :: x(3), b(3), bref(3), err, rmin
real(r8) :: gdr(0:nr, 0:nz), gdz(0:nr, 0:nz)
complex(r8), dimension(0:nr, 0:nz) :: br1, bz1, bp1, br16, bz16, bp16
integer :: i, ir, iz

err = 0
do i = 1, 8
x = a + [0.02_r8*3**(i - 1), 0.01_r8, 0.3_r8*i - 0.9_r8]
b = field_bs_segment(x(1) - a(1), x(2) - a(2), x(3) - a(3), &
d(1), d(2), d(3), 0.0_r8)*cross(d, x - a)
bref = segment_closed_form(a, a + d, x)
err = max(err, norm2(b - bref)/norm2(bref))
end do
print '(a,es10.2)', "segment vs closed form, max rel. error:", err
if (err > 1.0e-12_r8) error stop "field_bs_segment differs from the closed form"

do iz = 0, nz
do ir = 0, nr
gdr(ir, iz) = 0.95_r8 + 0.05_r8*ir
gdz(ir, iz) = -0.3_r8 + 0.15_r8*iz
end do
end do
cro = 1.0_r8
cnn = 1
coil_num = 1
allocate(coil(1))
call rectangle_field(1, br1, bz1, bp1)
call rectangle_field(16, br16, bz16, bp16)
err = max(maxval(abs(br1 - br16)), maxval(abs(bz1 - bz16)), &
maxval(abs(bp1 - bp16))) &
/max(maxval(abs(br16)), maxval(abs(bz16)), maxval(abs(bp16)))
print '(a,es10.2)', "rzphi 4 vs 64 pieces, max rel. difference:", err
if (err > 1.0e-12_r8) error stop "field_bs_rzphi depends on the polyline split"

rmin = 1.0e-3_r8
err = 0
do i = 0, 2
x = a + 0.3_r8*d + 0.5_r8*i*rmin*cross(d, [1.0_r8, 0.0_r8, 0.0_r8]) &
/norm2(cross(d, [1.0_r8, 0.0_r8, 0.0_r8]))
b = field_bs_segment(x(1) - a(1), x(2) - a(2), x(3) - a(3), &
d(1), d(2), d(3), rmin)*cross(d, x - a)
if (.not. norm2(b) <= 2/rmin*(1 + 1.0e-6_r8)) &
error stop "regularized field is not bounded by the rmin wire"
err = max(err, norm2(b)*rmin/2)
end do
print '(a,f8.5)', "near-filament max |B|/(mu0 I/(2 pi rmin)):", err

contains

pure function cross(u, v) result(w)
real(r8), intent(in) :: u(3), v(3)
real(r8) :: w(3)
w = [u(2)*v(3) - u(3)*v(2), u(3)*v(1) - u(1)*v(3), u(1)*v(2) - u(2)*v(1)]
end function cross

! Field per mu0 I/(4 pi) from the angles subtended at x by the segment p--q.
pure function segment_closed_form(p, q, x) result(bv)
real(r8), intent(in) :: p(3), q(3), x(3)
real(r8) :: bv(3), e(3), rho(3), r, z1, z2
e = (q - p)/norm2(q - p)
rho = (x - p) - dot_product(x - p, e)*e
r = norm2(rho)
z1 = -dot_product(x - p, e)
z2 = z1 + norm2(q - p)
bv = (z2/hypot(z2, r) - z1/hypot(z1, r))/r*cross(e, rho/r)
end function segment_closed_form

! Rectangle 0.6 m x 0.4 m in the plane x = 1.25 m, I = 1 kA, each side in m pieces.
subroutine rectangle_field(m, cbr, cbz, cbp)
integer, intent(in) :: m
complex(r8), dimension(0:nr, 0:nz), intent(out) :: cbr, cbz, cbp
real(r8), parameter :: cy(5) = [-0.3_r8, 0.3_r8, 0.3_r8, -0.3_r8, -0.3_r8]
real(r8), parameter :: cz(5) = [-0.2_r8, -0.2_r8, 0.2_r8, 0.2_r8, -0.2_r8]
integer :: iside, ipiece, ipt
real(r8) :: t

if (associated(coil(1)%x)) &
deallocate(coil(1)%x, coil(1)%y, coil(1)%z, coil(1)%cur)
coil(1)%coil_name = "rectangle"
coil(1)%ncoil = 1
coil(1)%s = 1
coil(1)%nsec = 4*m + 1
coil(1)%nw = 1.0_r8
allocate(coil(1)%x(1, 1, 4*m + 1), coil(1)%y(1, 1, 4*m + 1), &
coil(1)%z(1, 1, 4*m + 1), coil(1)%cur(1))
coil(1)%cur(1) = 1.0e3_r8
coil(1)%x = 1.25_r8
ipt = 0
do iside = 1, 4
do ipiece = 0, m - 1
t = real(ipiece, r8)/m
ipt = ipt + 1
coil(1)%y(1, 1, ipt) = cy(iside) + t*(cy(iside + 1) - cy(iside))
coil(1)%z(1, 1, ipt) = cz(iside) + t*(cz(iside + 1) - cz(iside))
end do
end do
coil(1)%y(1, 1, ipt + 1) = cy(5)
coil(1)%z(1, 1, ipt + 1) = cz(5)
cbr = 0
cbz = 0
cbp = 0
call field_bs_rzphi(nr, nz, np, gdr, gdz, cbr, cbz, cbp)
end subroutine rectangle_field

end program coil_segment_field