Notes on benchmarking several coil design codes (NESCOIL, REGCOIL, COILOPT++, FOCUS)¶
By Caoxiang Zhu
Generally, coil design codes optimize external currents to reproduce the target magnetic field. Currently, the four popular codes, NESCOIL, REGCOIL, COILOPT++ and FOCUS, are all trying to fit a given Bn distribution on a given toroidal boundary, despite the difference on solving the current potential or coil filaments and the optimization algorithms they used. Thus, the information of plasma currents is needed as a input and this is provided by the code BNORM. These five codes are developed independently and each code has its own numerical conventions. We checked the details in each code and benchmark the codes with one particular case.
For your convenience, you could just have a quick view at the Summary.
Table of Contents¶
BNORM ¶
BNORM code calculates the Fourier coefficients of the magnetic field (from plasma currents) normal to the VMEC surface. The only input file is VMEC output file, wout_xxx.nc. According to Merkel's note on BNORM, it actually calculates the exterior normal field on the boundary produced by surface currents.
Plasma boundary representation ¶
Although BNORM reads VMEC outputs, it uses (mu+nv) convention for representing the surface, where u and v are normalized poloidal and toroidal angle. Therefore, it flips the sign of n when reading VMEC coefficients. Here is the snippet of code from bn_read_vmecf90.f.
do mn = 1, mnmax
ixm(mn) = nint(xm(mn))
ixn(mn) =-nint(xn(mn))/nfp !!Flip sign: NESCOIL convention
m = ixm(mn)
n = ixn(mn)
bsubu(m,n) = 1.5_dp*bsubumnc(mn,ns) - 0.5_dp*bsubumnc(mn,ns-1)
bsubv(m,n) = 1.5_dp*bsubvmnc(mn,ns) - 0.5_dp*bsubvmnc(mn,ns-1)
cl(m,n) = 1.5_dp*lmns(mn,ns) - 0.5_dp*lmns(mn,ns-1)
cr(m,n) = rmnc(mn,ns)
cz(m,n) = zmns(mn,ns)
end do
Normal vector of plasma boundary ¶
To calculate the normal magnetic field on plasma boundary, it is essential to compute the normal vector of the surface. In BNORM, the normal vector $\bf N = X_u \times X_v$, where $\bf X_u$ and $\bf X_v$ are the derivatives with respect to angles. This part can be found in bn_surface.f.
do i = 1 , nuv
snx = yu(i)*zv(i)-zu(i)*yv(i)
sny = zu(i)*xv(i)-xu(i)*zv(i)
snz = xu(i)*yv(i)-yu(i)*xv(i)
sqf(i) = sqrt(snx*snx+sny*sny+snz*snz)
guu(i) = xu(i)*xu(i)+yu(i)*yu(i)+zu(i)*zu(i)
guv(i) = xu(i)*xv(i)+yu(i)*yv(i)+zu(i)*zv(i)
gvv(i) = xv(i)*xv(i)+yv(i)*yv(i)+zv(i)*zv(i)
dju(i) = bu(i)*guv(i)- bv(i)*guu(i)
djv(i) = bu(i)*gvv(i)- bv(i)*guv(i)
enddo
Bn convention ¶
BNORM takes a Fourier transformation over the normal field and calculates the Fourier coefficients. It represents the normal field as $B_n = \sum_m \sum_n bns_{m,n} \sin(mu+nv)$. This part can be found in bn_fouri.f
do m = 0,mf
faz = 2*fnuv
if(m.eq.0) faz = fnuv
do n = -nf,nf
bnfou(m,n) = 0
bnfou_c(m,n) = 0
do i = 1,nuv
cofp = conu(i,m)*conv(i,n)-sinu(i,m)*sinv(i,n) !cos(u+v)
sifp = sinu(i,m)*conv(i,n)+conu(i,m)*sinv(i,n) !sin(u+v)
bnfou(m,n) = bnfou(m,n) +bn(i)/sqf(i)*sifp*faz
! bnfou_c(m,n)=bnfou_c(m,n) +bn(i)/sqf(i)*cofp*faz
enddo
enddo
enddo
Normalization factor ¶
BNORM results are normalized so that the poloidal current per period, $I_p$, is one. This is found in bnormal.f.
!----------------------------------------------------------------------
c 15.08.98
c purpose:
c
c compuce
c /
c 1 | [ j', x-x'] . n
c B.n = ---- | ----------------- df , x,x' on the boundary
c 4 pi | |x-x'|^(3/2)
c /
c
c
c j = B_subu . xv - B_subv . xu
c
c the result ist normalized so that the toroidal integral
c
c /
c |
c I_p = | B. ds per field period is I_p =1.
c |
c /
c ----------------------------------------------------------------------And also in bn_parallel.f and bn_read_vmecf90.f.
m = ixm(mn)
n = ixn(mn)
bsubu(m,n) = 1.5_dp*bsubumnc(mn,ns) - 0.5_dp*bsubumnc(mn,ns-1)
bsubv(m,n) = 1.5_dp*bsubvmnc(mn,ns) - 0.5_dp*bsubvmnc(mn,ns-1)
curpol = alp*bsubv(0,0) ! alp=pi2/np
do m= 0,md
do n=-nd,nd
bsubu(m,n) = pi2*bsubu(m,n)/curpol
bsubv(m,n) = alp*bsubv(m,n)/curpol
bsubus(m,n) = pi2*bsubus(m,n)/curpol
bsubvs(m,n) = alp*bsubvs(m,n)/curpol
enddo
enddo
Comments ¶
In the above code comments, BNORM is calculating the interior normal field. In the code documentation (link), eq. 1 is indicating interior while eq. 12 is exterior. Based on FOCUS test, BNORM is calculating the interior normal field.
NESCOIL ¶
NESCOIL solves current potential on a given current-carrying-surface to produce the target magnetic field.
Plasma boundary representation ¶
NESCOIL reads boundary coefficients from nescin.xxx file, which is generally prepared by BNORM. It uses the plus convention for representing the boundary. This could be found in the original paper and in the source code surface_plas.f.
do 20 m = 0,ms1
do 20 n = 0,ns1
cm = m*pi2
cn = n*pi2
i = 0
do 20 kv = 1 , nv1/2 + 1
do 20 ku = 1 , nu1
i = i + 1
cofp = comu1(ku,m)*conv1(kv,n)-simu1(ku,m)*sinv1(kv,n)
cofm = comu1(ku,m)*conv1(kv,n)+simu1(ku,m)*sinv1(kv,n)
sifp = simu1(ku,m)*conv1(kv,n)+comu1(ku,m)*sinv1(kv,n)
sifm = simu1(ku,m)*conv1(kv,n)-comu1(ku,m)*sinv1(kv,n)
r1(i) = r1(i) + cr1(m,n)*cofp + cr1(m,-n)*cofm
z1(i) = z1(i) + cz1(m,n)*sifp + cz1(m,-n)*sifm
r1u(i) = r1u(i) - cm *( cr1(m,n)*sifp + cr1(m,-n)*sifm )
y1v(i) = y1v(i) - cn *( cr1(m,n)*sifp - cr1(m,-n)*sifm )
z1u(i) = z1u(i) + cm *( cz1(m,n)*cofp + cz1(m,-n)*cofm )
z1v(i) = z1v(i) + cn *( cz1(m,n)*cofp - cz1(m,-n)*cofm )
r1uu(i) = r1uu(i) -cm*cm*( cr1(m,n)*cofp + cr1(m,-n)*cofm )
z1uu(i) = z1uu(i) -cm*cm*( cz1(m,n)*sifp + cz1(m,-n)*sifm )
20 continue
Normal vector of plasma boundary ¶
In NESCOIL, the normal vector $\bf N = X_u \times X_v$, where $\bf X_u$ and $\bf X_v$ are the derivatives with respect to angles. This part can be found in surface_plas.f.
snx1(i) = y1u(i)*z1v(i) - z1u(i)*y1v(i)
sny1(i) = z1u(i)*x1v(i) - x1u(i)*z1v(i)
snz1(i) = x1u(i)*y1v(i) - y1u(i)*x1v(i)
! snx1, sny1, snz1 are later normalized to |N|
Bn convention ¶
NESCOIL reads bnorm.xxx file to get the Bn information on the boundary. As BNORM, it uses plus sign to compute the Bn on grids from Fourier coefficients. However, it has a factor, $-4\pi$, multiplied to Bn. It comes from $\mu_0/4\pi$, but $\mu_0$ has already been taken out in VMEC. This can be found in bnfld.f.
! Turn bnorm_mn into bn_ext(nuvh1) on plasma surface
bnorm_mn(:ms,0) = .5*bnorm_mn(:ms,0)
do m = 0, ms
do n = 0, ns
! cm = m*pi2
! cn = n*pi2
i = 0
do kv = 1 , nv1/2 + 1
do ku = 1 , nu1
i = i + 1
sifp = simu(ku,m)*conv(kv,n)+comu(ku,m)*sinv(kv,n)
sifm = simu(ku,m)*conv(kv,n)-comu(ku,m)*sinv(kv,n)
bn_ext(i) = bn_ext(i) + bnorm_mn(m,n)*sifp
1 + bnorm_mn(m,-n)*sifm
end do
end do
end do
end do
!
if (iloop .le. 0) deallocate (bnorm_mn)
fact = -2*pi2
bn_ext(:nuvh1) = fact*bn_ext(:nuvh1)
Normalization ¶
Although NESCOIL reads curpol in nescin.xxx, it doesn't use it and the parameters that determines the poloidal and toroidal current are cup and cup. Since BNORM has already normalized the Bn coefficients to let cup=1.0, there is no extra normalizations in NESCOIL. Of course, you could still provide cutomized cup and cut.
Total magnetic field ¶
After calculating the magnetic field from plasma, NESCOIL will add it the normal field produced by the total current potential, just as shown in accuracy.f.
call besurfcur(spx, spy, spz)
berr = snx1(i)*bx+sny1(i)*by+snz1(i)*bz + bn_ext(i)
Comments ¶
- It's strange that in the code NESCOIL uses opposite normal vector as it desicribes in the paper.
- bn_ext is multiplied by a factor of $-4\pi$ and this changes the sign of bn_ext. This is understandable if BNORM is calculating exterior normal component and NESCOIL uses interior normal vector. As to $4\pi$, I think it's from the constant $\mu_0 / 4\pi$ in Biot-Savart Law, while $\mu_0$ is taken out.
REGCOIL ¶
REGCOIL solves current potential on a given current-carrying-surface to produce the target magnetic field, regularized to the total current density.
Plasma boundary representation ¶
After reading boundary coefficients, REGCOIL uses the minus convention for representing the boundary, the same as VMEC. This could be found in the original paper and in the source code regcoil_init_plasma_mod.f90.
do izeta = 1,nzetal_plasma
angle2 = zetal_plasma(izeta)
sinangle2 = sin(angle2)
cosangle2 = cos(angle2)
dsinangle2dzeta = cosangle2
dcosangle2dzeta = -sinangle2
do itheta = 1,ntheta_plasma
do imn = 1,mnmax_plasma
angle = xm_plasma(imn)*theta_plasma(itheta) - xn_plasma(imn)*zetal_plasma(izeta) ! minus convention
sinangle = sin(angle)
cosangle = cos(angle)
dsinangledtheta = cosangle*xm_plasma(imn)
dcosangledtheta = -sinangle*xm_plasma(imn)
dsinangledzeta = -cosangle*xn_plasma(imn)
dcosangledzeta = sinangle*xn_plasma(imn)
r_plasma(1,itheta,izeta) = r_plasma(1,itheta,izeta) + rmnc_plasma(imn) * cosangle * cosangle2
r_plasma(2,itheta,izeta) = r_plasma(2,itheta,izeta) + rmnc_plasma(imn) * cosangle * sinangle2
r_plasma(3,itheta,izeta) = r_plasma(3,itheta,izeta) + zmns_plasma(imn) * sinangle
...
end do
end do
end do
Normal vector of plasma boundary ¶
In REGCOIL, the normal vector is computed as $\bf N = X_{\zeta} \times X_{\theta}$, where $\bf X_{\theta}$ and $\bf X_{\zeta}$ are the derivatives with respect to angles. This part can be found in regcoil_init_plasma_mod.f90.
! Evaluate cross product
normal_plasma(1,:,:) = drdzeta_plasma(2,:,:) * drdtheta_plasma(3,:,:) - drdtheta_plasma(2,:,:) * drdzeta_plasma(3,:,:)
normal_plasma(2,:,:) = drdzeta_plasma(3,:,:) * drdtheta_plasma(1,:,:) - drdtheta_plasma(3,:,:) * drdzeta_plasma(1,:,:)
normal_plasma(3,:,:) = drdzeta_plasma(1,:,:) * drdtheta_plasma(2,:,:) - drdtheta_plasma(1,:,:) * drdzeta_plasma(2,:,:)
Bn convention ¶
REGCOIL reads bnorm.xxx file to get the Bn information on the boundary. As BNORM, it uses plus sign to compute the Bn on grids from Fourier coefficients. It will multiply a factor of curpol to Bn. This can be found in regcoil_init_bnorm.f90.
num_modes_added = 0
do
read(iunit,*,iostat = i) mm, nn, bf
if (i .ne. 0) exit
num_modes_added = num_modes_added+1
!print *," Adding a mode with m=",mm,", n=",nn
do izeta = 1,nzeta_plasma
do itheta = 1,ntheta_plasma
Bnormal_from_plasma_current(itheta,izeta) = Bnormal_from_plasma_current(itheta,izeta) + &
bf*sin(mm*theta_plasma(itheta) + nn*nfp*zeta_plasma(izeta))
! To see that it should be (mu+nv) rather than (mu-nv) in the above line, you can examine
! BNORM/Sources/bn_fouri.f (where the arrays in the bnorm files are computed)
! or
! either NESCOIL/Sources/bnfld.f (where bnorm files are read)
end do
end do
end do
close(iunit)
! BNORM scales B_n by curpol=(2*pi/nfp)*bsubv(m=0,n=0)
! where bsubv is the extrapolation to the last full mesh point of
! bsubvmnc. Let's undo this scaling now.
Bnormal_from_plasma_current = Bnormal_from_plasma_current * curpol
Normalization ¶
In REGCOIL , the Bn coefficients are multiplied by curpol, which is obtained from reading VMEC outputs.
Total magnetic field ¶
The net coil currents, i.e. secular poloidal and toroidal currents from external coils, are also read in the unit of Amper. The produced magnetic field is calculated and stored in Bnormal_from_net_coil_currents. The total magnetic field is the sum of each field, as shown in regcoil_diagnostics.f90.
Bnormal_total(:,:,ilambda) = (reshape(matmul(g,solution),(/ ntheta_plasma, nzeta_plasma /)) / norm_normal_plasma) &
+ Bnormal_from_plasma_current + Bnormal_from_net_coil_currents
Comments ¶
There are two main differences between NESCOIL and REGCOIL that might cause problems.
- The normal vectors are in opposite directions.
- NESCOIL doesn't use SI units ($\mu_0$ is out), while REGCOIL does.
COILOPT++ ¶
COILOPT++ code calculates coil filaments on a given winding surface using B-spline representation.
Plasma boundary representation ¶
After reading boundary coefficients, COILOPT++ uses the plus convention for representing both the plasma boundary and the winding surface. Thus, it flips the sign of n. This could be found in the source code toroidalSurface.C.
// In readVMECcdf.C
// Rectify toroidal mode numbers (VMEC series is in terms of -phi, not +v).
for (int imode=0; imode<nmodes; imode++) {
assert(nnum[imode] % nfper == 0);
nnum[imode] /= -nfper;
}
/////////////////////////////////////
// Parametric surface evaluator (cylindrical)
// void vmecSurface::uv2Rz(const real u, const real v, real& R, real& z) const
// {
// real arg;
//
// R = z = 0.0;
// for (int mode=0; mode<nmodes; mode++) {
// arg = twopi*(mnum[mode]*u + nnum[mode]*v);
// R += rmn[mode] * cos(arg);
// z += zmn[mode] * sin(arg);
// }
// }
Normal vector of plasma boundary ¶
In COILOPT++, the normal vector $\bf N = X_v \times X_u$, where $\bf X_u$ and $\bf X_v$ are the derivatives with respect to angles. This part can be found in toroidalSurface.C.
////////////////////////////////////////////////////////////////////////
// Parametric surface evaluator with derivatives (Cartesian)
vector vmecSurface::uv2xyzprime(const real u, const real v,
vector prime[5])
const
{
real R[6], z[6], phi, cp, sp;
uv2Rzprime(u, v, R, z); //1st get derivs of cylindrical coords
phi = dphidv * v;
cp = cos(phi); sp = sin(phi);
vector x0(R[0]*cp, R[0]*sp, z[0]); // x, y, z
prime[0] = vector(R[1]*cp, R[1]*sp, z[1]); // d/du
prime[1] = vector(R[2]*cp - dphidv*x0.y,
R[2]*sp + dphidv*x0.x, z[2]); // d/dv
prime[2] = vector(R[3]*cp, R[3]*sp, z[3]); // d2/du2
prime[3] = vector(R[4]*cp - dphidv*prime[0].y,
R[4]*sp + dphidv*prime[0].x, z[4]); // d2/dudv
prime[4] = vector(R[5]*cp - dphidv*(2.0*R[2]*sp + dphidv*x0.x),
R[5]*sp + dphidv*(2.0*R[2]*cp - dphidv*x0.y),
z[5]); //d2/dv2
return x0;
}
////////////////////////////////////////////////////////////////
// Normal calculator
vector vmecSurface::unitnormal(const real u, const real v) const
{
vector prime[5];
uv2xyzprime(u, v, prime);
vector normal = prime[1] * prime[0];
real len = normal.mag();
assert(len > 0.0);
return normal/len;
}
Bn convention ¶
COILOPT++ reads bnorm.xxx file to get the Bn information on the boundary. As BNORM, it uses plus sign to compute the Bn on grids from Fourier coefficients. It will multiply a factor of curpol to Bn. This can be found in plasma.C.
// Read data, scale to VMEC poloidal current value.
rewind(fp);
for (int imode=0; imode<numBcoefs; imode++) {
fscanf(fp, "%d %d %le", mnum+imode, nnum+imode, bnormcoef+imode);
bnormcoef[imode] *= curpol_vmec;
}
Normalization factor ¶
In COILOPT++ , the Bn coefficients are multiplied by curpol, which is obtained from reading VMEC outputs.
Total magnetic field ¶
COILOPT++ calculates the total magnetic field and add the normal field to Bn_plasma, as shown in plasma.C.
vector Bcoils(0.0), sp(surfpos[impt]);
// Contributions from varying-geometry coils
for (int iws=0; iws<nws; iws++)
for (int iloop=0; iloop<nloops[iws]; iloop++)
Bcoils += loop[iws][iloop]->Bfield_T(sp);
// Contributions from fixed-geometry coils
for (int ifix=0; ifix<nfixed; ifix++)
Bcoils += fixedcurrent[ifix][0] * B_fixed[ifix][impt];
// Contribution from background TF
if (settings.bgtfR0B0 != 0.0)
Bcoils += ((settings.bgtfR0B0/(sp.x*sp.x + sp.y*sp.y)) *
vector(-sp.y, sp.x, 0.0));
// Contribution from background VF
Bcoils.z += bgBz;
real dBoB = bnorm[impt] + (Bcoils % nhat[impt]);
Comments ¶
COILOPT++ uses NESCOIL convention on plasma boundary, but flips the direction of normal vector.
FOCUS ¶
FOCUS nonlinearly optimizes coil filaments using 3D representation.
Plasma boundary representation ¶
After reading boundary coefficients, FOCUS uses the minus convention for representing the boundary, the same as VMEC. This could be found in the original paper and in the source code rdsurf.h.
! The center point value was used to discretize grid;
do ii = 0, Nteta-1; teta = ( ii + half ) * pi2 / Nteta
do jj = 0, Nzeta-1; zeta = ( jj + half ) * pi2 / ( Nzeta*Nfp )
RR(0:2) = zero ; ZZ(0:2) = zero
do imn = 1, Nfou ; arg = bim(imn) * teta - bin(imn) * zeta ! mu-nv
RR(0) = RR(0) + Rbc(imn) * cos(arg) + Rbs(imn) * sin(arg)
ZZ(0) = ZZ(0) + Zbc(imn) * cos(arg) + Zbs(imn) * sin(arg)
RR(1) = RR(1) + ( - Rbc(imn) * sin(arg) + Rbs(imn) * cos(arg) ) * bim(imn)
ZZ(1) = ZZ(1) + ( - Zbc(imn) * sin(arg) + Zbs(imn) * cos(arg) ) * bim(imn)
RR(2) = RR(2) - ( - Rbc(imn) * sin(arg) + Rbs(imn) * cos(arg) ) * bin(imn)
ZZ(2) = ZZ(2) - ( - Zbc(imn) * sin(arg) + Zbs(imn) * cos(arg) ) * bin(imn)
...
end do
end do
end do
Normal vector of plasma boundary ¶
In FOCUS, the normal vector is computed as $\bf N = X_{\zeta} \times X_{\theta}$, where $\bf X_{\theta}$ and $\bf X_{\zeta}$ are the derivatives with respect to angles. This part can be found in rdsurf.h.
xx(1:3) = (/ RR(0) * czeta, RR(0) * szeta, ZZ(0) /) ! x, y, z
xt(1:3) = (/ RR(1) * czeta, RR(1) * szeta, ZZ(1) /) ! dx/dtheta
xz(1:3) = (/ RR(2) * czeta, RR(2) * szeta, ZZ(2) /) + (/ - RR(0) * szeta, RR(0) * czeta, zero /) ! dx/dzeta
ds(1:3) = -(/ xt(2) * xz(3) - xt(3) * xz(2), & ! minus sign for theta counterclockwise direction;
xt(3) * xz(1) - xt(1) * xz(3), &
xt(1) * xz(2) - xt(2) * xz(1) /)
dd = sqrt( sum( ds(1:3)*ds(1:3) ) )
! x, y, z coordinates for the surface;
surf(1)%xx(ii,jj) = xx(1)
surf(1)%yy(ii,jj) = xx(2)
surf(1)%zz(ii,jj) = xx(3)
! dx/dt, dy/dt, dz/dt (dt for d theta)
surf(1)%xt(ii,jj) = xt(1)
surf(1)%yt(ii,jj) = xt(2)
surf(1)%zt(ii,jj) = xt(3)
! surface normal vectors and ds for the jacobian;
surf(1)%nx(ii,jj) = ds(1) / dd
surf(1)%ny(ii,jj) = ds(2) / dd
surf(1)%nz(ii,jj) = ds(3) / dd
surf(1)%ds(ii,jj) = dd
Bn convention ¶
FOCUS reads Bn coefficients to get the Bn information on the boundary. Unlike BNORM, it uses minus sign to compute the Bn on grids from Fourier coefficients. It will multiply a factor of curpol to Bn. This can be found in vmec2focus.py and rdsurf.h.
...
curpol = 2.0*np.pi / Nfp * wout['rbtor'].values # not rbtor0, it's rbtor
...
for imn in range(Nbnf):
fofile.write("{:d} \t {:d} \t {:23.15E} \t {:23.15E} \n".format(-bn[imn], bm[imn], 0.0, bns[imn]))
# negative n, because FOCUS and BNORM are taking different rep.
...
!calculate target Bn with input harmonics; 05 Jan 17;
if(NBnf > 0) then
do jj = 0, Nzeta-1 ; zeta = ( jj + half ) * pi2 / (Nzeta*Nfp)
do ii = 0, Nteta-1 ; teta = ( ii + half ) * pi2 / Nteta
do imn = 1, NBnf
arg = Bnim(imn) * teta - Bnin(imn) * zeta
surf(1)%pb(ii,jj) = surf(1)%pb(ii,jj) + Bnc(imn)*cos(arg) + Bns(imn)*sin(arg)
enddo
enddo
enddo
endif
Normalization ¶
Inside FOCUS, there is no normalization on Bn. But you have to pay attention to providing the correct Bn.
Total magnetic field ¶
FOCUS reads all the other magnetic field except the coils to be solved and treats them as the target boundary condition. Therefore, the total magnetic field is obtained by subtraction, as shown in bnormal.h.
lbn(iteta, jzeta) = lbx(iteta, jzeta)*surf(1)%nx(iteta, jzeta) &
& + lby(iteta, jzeta)*surf(1)%ny(iteta, jzeta) &
& + lbz(iteta, jzeta)*surf(1)%nz(iteta, jzeta) &
& - surf(1)%pb(iteta, jzeta)
Comments ¶
FOCUS relies on the input Bn information and when discretizing the plasma surface there is a half grid shift.
Summary ¶
The following table briefly summaries each code's numerical conventions.
| Code | Plasma surface | Normal vector | Bn Fourier rep. | Internal normalization | Total Bn |
|---|---|---|---|---|---|
| BNORM | (mu+nv) | $Bn_{plasma}$ (interior) | (mu+nv) | 1/curpol | NA |
| NESCOIL | (mu+nv) | $\bf X_u \times X_v$ (interior) | (mu+nv) | $-4 \pi$ | add |
| REGCOIL | $(m\theta - n\zeta)$ | $\bf X_{\zeta} \times X_{\theta}$ (exterior) | $(m\theta + n\zeta)$ | curpol | add |
| COILOPT++ | (mu+nv) | $\bf X_v \times X_u$ (exterior) | (mu+nv) | curpol | add |
| FOCUS | $(m\theta - n\zeta)$ | $\bf X_{\zeta} \times X_{\theta}$ (exterior) | $(m\theta - n\zeta)$ | 1 | subtract |
The conclusions are
- The total Bn, which is the normal component of the total field, should be zero in all coil design codes.
- BNORM calculates the exterior normal components from the plasma currents. Although it calculates the interior normal vector, it only uses the length of normal vector for later computations.
- NESCOIL, REGCOIL and COILOPT++ have different interpretations on plasma boundary, normal vector and normalizations. They all assume BNORM is calculating the exterior normal field, while my observation is not.
- NESCOIL is not using the SI unit. The current potential solution need to be multiplied with a factor of $I_p/nfp$, where $I_p$ is the actual number of net poloidal current.
- FOCUS relies on user to provide correct input Bn. If using the output of BNORM, user should flip the sign of toroidal modes (n) and multiply a factor of (curpol) to the Fourier coefficients.
Appendix A: Numerical verifications among NESCOIL, REGCOIL and FOCUS ¶
Equiblirium ¶
We choose a simple axisymmetric configuration with finite beta, Solov'ev equilibrium, to test the three codes. The case is used for benchmarking VMEC. It's used in the original VMEC paper.
Here is the input.SOLOVEV file for VMEC fixed-boundary run.
&INDATA
! SOLOV'EV Equilibrium from Hirshman Whitson 1983
MGRID_FILE = 'none'
LFREEB = F
DELT = 9.00E-01
TCON0 = 2.00E+00
LASYM = F
NFP = 1
NCURR = 0
MPOL = 12
NTOR = 0
NZETA = 200
NS_ARRAY = 16 32 64 128
NITER_ARRAY = 1000 2000 4000 20000
FTOL_ARRAY = 1E-30 1E-30 1E-30 1E-12
NITER = 5000
NSTEP = 250
NVACSKIP = 6
GAMMA = 0.000000E+00
PHIEDGE = 1.00000000000000E+00
BLOAT = 1.000000E+00
CURTOR = 0.00000000000000E+00
SPRES_PED = 1.00000000000000E+00
PMASS_TYPE = "power_series"
AM = 0.125 -0.125
AI = 1.0
RAXIS = 3.999
ZAXIS = 0.00000000000000E+00
RBC( 0,0) = 3.999 ZBS( 0,0) = 0.000
RBC( 0,1) = 1.026 ZBS( 0,1) = 1.580
RBC( 0,2) = -0.068 ZBS( 0,2) = 0.010
/
&ENDThe iota profile is flat, $\iota = 1.0$ and pressure profile is $p = (1-s^2)/8$. The soloved equilibrium is shown in the below, with $B_0 = 0.20334$T, $I_{toroidal}=-0.44$MA and $\beta_{total} = 3.87\times 10^{-6}$.
for i in np.arange(0,129,16)-1:
plot_vmec_surface('wout_SOLOVEV.nc',i,color='b', lbl=None)
plot_vmec_surface('wout_SOLOVEV.nc',color='b',lbl='VMEC fixed boundary');
plt.legend()
BNORM results ¶
BNORM is ran to calculate the normal component of plasma currents. The resolution is $Mpol=24$ and $Ntor=14$. We used BNORM to produce a uniformly expanded surface with a separation of 0.5m. Reading the Fourier coefficients and plot out the distribution of Bn_plasma.
plt.figure()
plt.contourf(tbn,zbn,bnorm, 18)
plt.xlabel('zeta', fontsize=16)
plt.ylabel('theta', fontsize=16)
plt.title('Bnormal from BNORM code')
plt.colorbar()
fig = plt.figure()
plot_plasma_boundary('plasma.boundary', color='r', lbl='plasma boundary');
plot_winding_surface('test.winding', color='b', lbl='winding surface');
plt.legend()
NESCOIL results ¶
Run NESCOIL with the input file nescin.SOLOVEV prepared by BNORM. Some resolution parameters are as below.
------ Spatial dimensions ----
nu, nv, nu1, nv1, npol, ntor, lasym_bn
128 128 128 128 64 10 F
------ Fourier Dimensions ----
mf, nf, md, nd (max in surf and bnorm files)
8 8 24 20Here is part of the output.
Complexity = 1.13817350E+00
Time in Solver: 19.8 sec
----- Calling Surfcur_Diag -----
J Surf max, min, ave = 6.67491028E-02 3.14638915E-02 4.39897927E-02
Curv_R of J min, max = 1.12452252E+00 4.22091801E+00
Time in Surfcur_Diag: 0.149 sec
----- Calling Accuracy -----
Berr ave, max, var = 3.65157304E-05 7.60790295E-05 2.03785308E-05
Bmod ave, rms = 2.58046901E-01 3.68682330E-01
ACCURACY took 10.6 sec
ONE NESCOIL RUN took 30.6 secNESCOIL only writes the single-valued current potential. We can read in the Fourier coefficients and plot out the distribution of $\Phi_{SV}$. Of course, by adding the secular terms, $G$ & $I$, we could obtain the total current potential.
plt.figure()
plt.subplot(121)
plt.contourf(zeta,theta,cp, 18)
plt.xlabel('zeta', fontsize=16)
plt.ylabel('theta', fontsize=16)
plt.title('Single-valued current potential by NESCOIL')
plt.colorbar()
######
plt.subplot(122)
plt.contourf(zeta,theta,tcp, 18)
plt.xlabel('zeta', fontsize=16)
plt.ylabel('theta', fontsize=16)
plt.title('Total current potential by NESCOIL')
plt.colorbar()
REGCOIL results ¶
Run REGCOIL with the following input file.
®coil_nml
nlambda = 1
lambda_min = 1e-15
lambda_max = 1e-14
ntheta_plasma = 128
ntheta_coil = 128
nzeta_plasma = 128
nzeta_coil = 128
mpol_potential = 8
ntor_potential = 8
geometry_option_plasma = 2
wout_filename = 'wout_SOLOVEV.nc'
geometry_option_coil=2
nescin_filename = 'test.winding'
separation = 0.5
load_bnorm=.t.
bnorm_filename='bnorm.SOLOVEV'
net_poloidal_current_Amperes = 0.0
net_toroidal_current_Amperes = 0.0
symmetry_option = 1
save_level = 1
/It is running without any regularization on current density and the same resolution with NESCOIL. Here is one part of the screen output.
Solving system for lambda= 0.000E+00 ( 1 of 1)
Additions: 0.00000000 sec.
DSYSV: 1.00000005E-03 sec.
Diagnostics: 2.19999999E-02 sec.
chi2_B: 8.831E-09, chi2_K: 6.755E+12, chi2_Laplace_Beltrami: 1.684E+11
max(B_n): 1.130E-05, max(K): 2.413E+05, rms K: 1.537E+05
REGCOIL complete. Total time= 19.229000091552734 sec.Plot out the results.
reg_scp = plot_regcoil_potential('regcoil_out.test.nc',0,num_contours=18, total=False);
plt.title('Single-valued current potential by REGCOIL')
reg_tcp = plot_regcoil_potential('regcoil_out.test.nc',0,num_contours=18, total=True);
plt.title('Total current potential by REGCOIL')
Benchmark between NESCOIL and REGCOIL ¶
We can now taking the conclusion we derive above, i.e. normalize REGCOIL results with a factor of $I_p/nfp = 3.6152$MA in this case. The numbers are consistent, to some precision.
plt.figure()
plt.subplot(121)
plt.contourf(zeta,theta,reg_scp/fac-cp, 18)
plt.xlabel('zeta', fontsize=16)
plt.ylabel('theta', fontsize=16)
plt.title('Difference between NESCOIL and norm. REGCOIL (Phi_SV)')
plt.colorbar()
plt.subplot(122)
plt.contourf(zeta,theta,cp, 18)
plt.xlabel('zeta', fontsize=16)
plt.ylabel('theta', fontsize=16)
plt.title('Original NESCOIL results (Phi_SV)')
plt.colorbar()
plt.figure()
plt.subplot(121)
plt.contourf(zeta,theta,reg_tcp/fac-tcp, 18)
plt.xlabel('zeta', fontsize=16)
plt.ylabel('theta', fontsize=16)
plt.title('Difference between NESCOIL and norm. REGCOIL results (Phi_total)')
plt.colorbar()
plt.subplot(122)
plt.contourf(zeta,theta,tcp, 18)
plt.xlabel('zeta', fontsize=16)
plt.ylabel('theta', fontsize=16)
plt.title('Original NESCOIL results (Phi_total)')
plt.colorbar()
Free-boundary VMEC results ¶
We can now cut coils from the contours of total current potential. We assume there are 32 coils in total. The amount of current flowing in each coil is $I_p/N = 0.113$MA. Here is a picture of final coils (using cutCoilsFromREGCOIL, num_of_contours=32, no theta shift).
## colors on plasma is the Bn produced by plasma currents
fig = mlab.figure(size=[800,800])
plot_coils(read_coils('coils.test_regcoil'));
plot_plasma_boundary('plasma.boundary', 'bnormal');
mlab.colorbar()
fig
A mgrid file is generated by this modular coils and then used in VMEC for free-boundary running. Somehow, the sign of phiedge in VMEC has to be flipped. Except that, the free-boundary VMEC results are consistent with fixed-boundary results, with surprisingly good precision.
for i in np.arange(0,129,16)-1:
plot_vmec_surface('wout_SOLOVEV.nc',i,color='b', lbl=None)
plot_vmec_surface('wout_regcoil_fb.nc',i,color='r', style='--',lbl=None)
plot_vmec_surface('wout_SOLOVEV.nc',color='b',lbl='fixed boundary');
plot_vmec_surface('wout_regcoil_fb.nc', color='r', style='--', lbl='REGCOIL free boundary');
plt.legend()
When using FOCUS to read the modular coils and diagnose the magnetic field, the surface integral of residual normal field is 1.15695E-03.
half = 0.5
theta_focus = (np.arange(128) + half) * (2*np.pi)/128
zeta_focus = (np.arange(256) + half) * (2*np.pi)/256
plt.figure()
plt.subplot(121)
plt.contourf(zeta_focus, theta_focus, np.transpose(regcoil.Bn + regcoil.plas_Bn), 18)
plt.xlabel('zeta', fontsize=16)
plt.ylabel('theta', fontsize=16)
plt.title('Bn from discretized modular coils (calculated by FOCUS)')
plt.colorbar()
plt.subplot(122)
plt.contourf(zeta_regcoil, theta_regcoil, np.transpose(reg_Bn), 18)
plt.xlabel('zeta', fontsize=16)
plt.ylabel('theta', fontsize=16)
plt.title('Bn_plasma (readed from REGCOIL)')
plt.colorbar()
f = netcdf.netcdf_file('regcoil_out.test.nc','r',mmap=False)
reg_Bn = f.variables['Bnormal_from_plasma_current'][()]
theta_regcoil = f.variables['theta_plasma'][()]
zeta_regcoil = f.variables['zeta_plasma'][()]
plt.figure()
plt.subplot(121)
plt.contourf(zeta_regcoil, theta_regcoil, np.transpose(reg_Bn), 18)
plt.xlabel('zeta', fontsize=16)
plt.ylabel('theta', fontsize=16)
plt.title('Bn_plasma (readed from REGCOIL)')
plt.colorbar()
plt.subplot(122)
plt.contourf(zeta_focus, theta_focus, np.transpose(regcoil.plas_Bn), 18)
plt.xlabel('zeta', fontsize=16)
plt.ylabel('theta', fontsize=16)
plt.title('Bn_plasma (input for FOCUS)')
plt.colorbar()
FOCUS results ¶
When using FOCUS, we initialized 32 TF coils plus 6 PF coils. With BNORM Fourier coefficients, we flipped the sign of toroidal modes and multiply the coefficients with a factor of ($- curpol$), as we described in the above. FOCUS got a result with relatively low Bn error, $f_B = \int_S \frac{1}{2} ({\bf B} \cdot {\bf n})^2 = 1.02 \times 10^{-7}$.
from mayavi import mlab
mlab.init_notebook()
%matplotlib notebook
%run ~/python/scripts/coilpy.py
import seaborn as sns
sns.set(style="white", context="paper")
cd /p/focus/users/czhu/SPEC/test/SOLOVEV
vmec2focus('wout_solovev.nc', focus_file='solovev.boundary_noflip', bnorm_file='bnorm.solovev', flipsign=False)
vmec2focus('wout_solovev.nc', focus_file='solovev.boundary_flip', bnorm_file='bnorm.solovev', flipsign=True)
cd /p/focus/users/czhu/NCSX/c09r00/
vmec2focus('wout_c09r00_fb.nc', focus_file='c09r00.boundary_noflip', bnorm_file='bnorm.c09r00_fb', flipsign=False)
vmec2focus('wout_c09r00_fb.nc', focus_file='c09r00.boundary_flip', bnorm_file='bnorm.c09r00_fb', flipsign=True)
3.61519089469759/32
modular = hdf5('focus_solovev_20200110_modular.h5')
fig = mlab.figure(size=[800,800])
plot_coils(read_coils('solovev_20200110_modular.coils'), currents=False);
plot_hdf5_surface(modular, 'Bmod');
mlab.colorbar()
fig
for i in np.arange(0,129,16)-1:
plot_vmec_surface('wout_solovev.nc',i,color='b', lbl=None)
plot_vmec_surface('wout_solovev_20200110.nc',i,color='r', style='--',lbl=None)
plot_vmec_surface('wout_solovev.nc',color='b',lbl='fixed boundary');
plot_vmec_surface('wout_solovev_20200110.nc', color='r', style='--', lbl='FOCUS free boundary');
plt.legend(fontsize=15)
Appendix B: A simple verification of FOCUS ¶
This is a small demo showing that FOCUS is calculating magnetic field correctly. There is a circular plasma boundary with $R_p=4.0$ m, $r_p=1.0$ m. Assuming there are 128 circular TF coils ($r_c = 2.0$ m) with $I_c = 0.1$MA. The toroidal magnetic field can be computed by Ampere's Law, $B_t = \frac{\mu_0 I}{2\pi r}$. In this case, $I = 128*Ic=1.28\times 10^7$A and $B_t = 2.56/r$.
The toroidal flux should be $\Psi_t = \int_ 3^5 2*\sqrt {1 - (x - 4)^2} \frac{2.56}{x}\, dx = -2.043057082564284$ Wb. It is negative because the coil current flows in the counter-clockwise direction, which means the total current is flowing in negative Z direction and creates a toroidal field in the clock-wise direction.
The magnetic field calculated by FOCUS shows good consistency with the analytical values. The calculated toroidal flux is $-2.043584789045241$ Wb, which is also close to the analytical value up to numerical discretization erros.
test = hdf5('../focus_test.h5')
fig = mlab.figure(size=[800,800])
plot_coils(read_coils('../test.coils'));
plot_hdf5_surface(test, 'Bmod')
mlab.colorbar()
fig
plt.figure()
Bt = np.sqrt(test.Bx*test.Bx + test.By*test.By)
r = np.sqrt(test.xsurf*test.xsurf + test.ysurf*test.ysurf)
plt.plot(r[0,:], Bt[0,:], label='FOCUS value')
plt.plot(r[0,:], 2.56/r[0,:], marker='s', color='r', linestyle='none', label='Analytical value')
plt.ylabel('Bt [T]', fontsize=16)
plt.xlabel('R [m]', fontsize=16)
plt.legend()
