From e52d54d67e5f62744f132f3e8a8fb5d28eda2195 Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Sun, 25 Oct 2026 19:49:28 +0100 Subject: [PATCH] COIL - BUG FIX - Use the exact field of each straight coil segment field_bs_psi and field_bs_rzphi evaluated the Biot-Savart integrand of each polyline segment once, at its midpoint (dl x r/|r|^3). The error is of order (L/d)^2 for a segment of length L seen from distance d, so coils with long straight sides gave fields that were wrong by about 10 % on nearby surfaces. Each segment now contributes its closed-form field (Hanson and Hirshman, Phys. Plasmas 9, 4410 (2002)), which is exact for the polygon in the coil file. In field_bs_rzphi the former distance floor dlmin becomes the radius of a regularizing wire. regression/test/coil_segment_field.f90 checks the segment field against the textbook closed form and that field_bs_rzphi does not depend on how the straight sides of a rectangular coil are split. (cherry picked from commit 5212ddc541ca5e6d4eb263ca73b2e1cc0b37d0c4) --- coil/field.F | 85 +++++++++++------ regression/test/coil_segment_field.f90 | 123 +++++++++++++++++++++++++ 2 files changed, 180 insertions(+), 28 deletions(-) create mode 100644 regression/test/coil_segment_field.f90 diff --git a/coil/field.F b/coil/field.F index 02733b760..8d3dd853b 100644 --- a/coil/field.F +++ b/coil/field.F @@ -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. @@ -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) @@ -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 @@ -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) @@ -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 @@ -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 diff --git a/regression/test/coil_segment_field.f90 b/regression/test/coil_segment_field.f90 new file mode 100644 index 000000000..60baaa04a --- /dev/null +++ b/regression/test/coil_segment_field.f90 @@ -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