STELLOPT Boundary Representations

By Caoxiang Zhu

In fixed-boundary optimizations, STELLOPT varies the boundary coefficients to achieve target properties. There are three boundary representations in STELLOPT: VMEC representation, Hirshman-Breslau representation and Garabedian representation. There is a brief webpage introducing these representations on VMECwiki. In this page, we are going to introduce them with more details.

VMEC representation

General description

The VMEC representation is from the VMEC code. A toroidal surface, $\cal{S}$, is represented in cylinder coordinates $(R, \phi, Z)$ as $$ R = \sum R_{mn} \cos(m \theta -nN\phi) $$ $$ Z = \sum Z_{mn} \sin(m \theta -nN\phi) $$ Here, we assume that the surface has the so-called "stellarator symmetry", i.e. $R(-\theta, -\phi) = R(\theta, \phi)$ and $Z(-\theta, -\phi) = -Z(\theta, \phi)$. $N>0$ is the number of field periods. $\theta$ is a parametric angle variable and is usually called the poloidal angle, while the cylindrical angle $\phi$ is also called the toroidal angle.

The two arrays of VMEC coefficients, $R_{mn}$ and $Z_{mn}$, can now uniquely describe a toroidal surface (of course, if no stellarator symmetry is assumed, we need four arrays, both cosine and sine harmonics for R and Z). The poloidal angle $\theta$ is not constrained, thus there is null space.

Numerical details

In the INDATA namelist (for both VMEC and STELLOPT), the following variables are related to the VMEC representation.

Variable name Type Description
NFP integer Number of toroidal field periods
LASYM logical Switch for stellarator symmetry, T: no; F: yes
MPOL integer Poloidal Mode Number (m)
NTOR integer Toroidal Mode Number (n)
RBC(N,M) real(-NTOR:NTOR, 0:MPOL) Boundary cosine coefficients for R
ZBS(N,M) real(-NTOR:NTOR, 0:MPOL) Boundary sine coefficients for Z
RBS(N,M) real(-NTOR:NTOR, 0:MPOL) Boundary cosine coefficients for R
ZBC(N,M) real(-NTOR:NTOR, 0:MPOL) Boundary sine coefficients for Z

In the OPTIMUM namelist (for STELLOPT optimization), the variables LBOUND_OPT(N,M) are used to control the boundary coefficients. For example, if LBOUND_OPT(-1,1) = T is set, the RBC(-1,1) and ZBS(-1,1) parameters are varied during the optimization (LASYM = F). If LASYM = T in INDATA, then RBS(-1,1) and ZBC(-1,1) are also free variables. The corresponding BOUND variables are RBC_MIN(N,M), RBC_MAX(N,M), etc.

Example

In [23]:
cd /p/stellopt/ANALYSIS/czhu/runs/test
/p/stellopt/ANALYSIS/czhu/runs/test

We use the NCSX boundary as illustration. The crosssections at $\phi=0$ and $\phi=\pi/3$ (half period) are plotted below. The dots on the following figure are equally spaced in poloidal angle $\theta$. There are more points at the position with larger curvature. What should be noted is that the boundary data probably comes from a STELLOPT optimization and the poloidal spectrum has already been condensed.

!----- Grid Parameters -----
  LASYM = F
  NFP = 0003
  MPOL = 0009
  NTOR = 0005
  PHIEDGE =    5.143860000000E-001
!----- Boundary Parameters -----
  RBC(0,0) =   1.3782E+00     ZBS(0,0) =   0.0000E+00
  RBC(1,0) =  -4.1452E-03     ZBS(1,0) =   7.1634E-03
  RBC(2,0) =  -4.1374E-03     ZBS(2,0) =   7.5995E-03
  RBC(3,0) =   3.6116E-03     ZBS(3,0) =  -2.9503E-03
  RBC(4,0) =  -3.3285E-04     ZBS(4,0) =   5.1481E-04
  RBC(-4,1) = -4.5103E-04     ZBS(-4,1) =  7.0486E-05
  RBC(-3,1) =  4.7123E-04     ZBS(-3,1) = -1.5952E-04
  RBC(-2,1) = -2.6156E-03     ZBS(-2,1) =  1.4116E-03
  RBC(-1,1) =  2.0869E-02     ZBS(-1,1) =  9.2873E-03
  RBC(0,1) =   2.7073E-01     ZBS(0,1) =   4.6465E-01
  RBC(1,1) =  -1.3500E-01     ZBS(1,1) =   1.6516E-01
  RBC(2,1) =   4.8519E-03     ZBS(2,1) =  -6.0559E-03
  RBC(3,1) =   4.8906E-04     ZBS(3,1) =  -1.7736E-04
  RBC(4,1) =  -9.4694E-04     ZBS(4,1) =   5.6642E-04
  RBC(-4,2) =  4.4897E-04     ZBS(-4,2) =  1.7521E-04
  RBC(-3,2) = -4.6166E-04     ZBS(-3,2) =  5.9225E-04
  RBC(-2,2) = -1.9988E-04     ZBS(-2,2) =  1.9350E-04
  RBC(-1,2) =  8.9339E-03     ZBS(-1,2) =  1.2754E-02
  RBC(0,2) =   9.1094E-02     ZBS(0,2) =   1.5451E-02
  RBC(1,2) =   9.0684E-02     ZBS(1,2) =   1.1618E-02
  RBC(2,2) =   2.7208E-02     ZBS(2,2) =  -2.7337E-02
  RBC(3,2) =  -4.4665E-03     ZBS(3,2) =   1.3515E-03
  RBC(4,2) =  -1.4812E-04     ZBS(4,2) =   8.6535E-04
  RBC(-4,3) = -2.4030E-05     ZBS(-4,3) =  1.7592E-05
  RBC(-3,3) =  5.9433E-04     ZBS(-3,3) = -5.8654E-04
  RBC(-2,3) =  5.2473E-05     ZBS(-2,3) = -1.0220E-04
  RBC(-1,3) =  3.1952E-03     ZBS(-1,3) = -3.0523E-03
  RBC(0,3) =  -1.4174E-02     ZBS(0,3) =   1.1779E-02
  RBC(1,3) =  -1.7796E-03     ZBS(1,3) =  -1.9260E-03
  RBC(2,3) =  -1.0529E-02     ZBS(2,3) =   1.0663E-02
  RBC(3,3) =  -5.8013E-03     ZBS(3,3) =   5.8095E-03
  RBC(4,3) =   1.4257E-03     ZBS(4,3) =  -1.4443E-03
  RBC(-4,4) =  2.5893E-05     ZBS(-4,4) = -8.7822E-06
  RBC(-3,4) = -3.2665E-04     ZBS(-3,4) =  2.6078E-04
  RBC(-2,4) = -4.8039E-05     ZBS(-2,4) =  2.3453E-05
  RBC(-1,4) = -5.2966E-04     ZBS(-1,4) =  2.9088E-04
  RBC(0,4) =   4.5451E-03     ZBS(0,4) =   1.8260E-04
  RBC(1,4) =  -4.6998E-03     ZBS(1,4) =   9.6416E-03
  RBC(2,4) =   6.4586E-03     ZBS(2,4) =  -3.0495E-03
  RBC(3,4) =   1.7811E-04     ZBS(3,4) =  -5.4173E-04
  RBC(4,4) =   1.7237E-03     ZBS(4,4) =  -1.7870E-03
  RBC(-4,5) = -2.6971E-06     ZBS(-4,5) = -2.6971E-06
  RBC(-3,5) =  7.6520E-05     ZBS(-3,5) =  7.6520E-05
  RBC(-2,5) =  1.0022E-05     ZBS(-2,5) =  1.0022E-05
  RBC(-1,5) =  4.0483E-04     ZBS(-1,5) =  4.0483E-04
  RBC(0,5) =  -1.6817E-03     ZBS(0,5) =  -1.6817E-03
  RBC(1,5) =   9.4841E-06     ZBS(1,5) =   9.4841E-06
  RBC(2,5) =  -1.3733E-03     ZBS(2,5) =  -1.3733E-03
  RBC(3,5) =  -7.5237E-04     ZBS(3,5) =  -7.5237E-04
  RBC(4,5) =   1.8598E-04     ZBS(4,5) =   1.8598E-04
  RBC(-4,6) =  3.4247E-06     ZBS(-4,6) =  3.4247E-06
  RBC(-3,6) = -5.8017E-05     ZBS(-3,6) = -5.8017E-05
  RBC(-2,6) = -7.0609E-06     ZBS(-2,6) = -7.0609E-06
  RBC(-1,6) = -8.1041E-05     ZBS(-1,6) = -8.1041E-05
  RBC(0,6) =   4.3087E-04     ZBS(0,6) =   4.3087E-04
  RBC(1,6) =  -1.4165E-03     ZBS(1,6) =  -1.4165E-03
  RBC(2,6) =   9.3906E-04     ZBS(2,6) =   9.3906E-04
  RBC(3,6) =   7.1096E-05     ZBS(3,6) =   7.1096E-05
  RBC(4,6) =   3.4674E-04     ZBS(4,6) =   3.4674E-04
!-----------------------------------------------------------------------
!          OPTIMIZED QUANTITIES
!-----------------------------------------------------------------------
  LBOUND_OPT( 0:5,0) =  6*T  
  LBOUND_OPT(-5:5,1) = 11*T  
  LBOUND_OPT(-5:5,2) = 11*T  
  LBOUND_OPT(-5:5,3) = 11*T  
  LBOUND_OPT(-5:5,4) = 11*T  
  LBOUND_OPT(-5:5,5) = 11*T  
  LBOUND_OPT(-5:5,6) = 11*T  
  LBOUND_OPT(-5:5,7) = 11*T
In [102]:
plt.figure()
plot_plasma_boundary('ncsx.boundary', zeta=0, marker='o', npoints=64)
plot_plasma_boundary('ncsx.boundary', zeta=np.pi/3, marker='o', color='b', npoints=64)
read 59 Fourier harmonics for plasma boundary and 0 for Bn distribution in ncsx.boundary.
read 59 Fourier harmonics for plasma boundary and 0 for Bn distribution in ncsx.boundary.
Out[102]:
(1.7819387236000002,
 1.0382141884626426,
 0.27003475415163336,
 -0.2700347541516334)

Hirshman-Breslau representation

General description

The Hirshman-Breslau representation is introduced in Hirshman, S. P., & Breslau, J. (1998). Physics of Plasmas, 5(7). In that paper, Hirshman & Breslau introduced two unqiue representations to provide additional contraints on the poloidal angle. The one used in STELLOPT is the so-called "quasi-polar"representation. It uses an effective (quasipolar) "radius" $\rho$ as the only coordinate. $$ \rho = \sum \rho^c_{mn} \cos(m \theta -nN\phi) + \rho^s_{mn} \sin(m \theta -nN\phi)$$ $$ r = \rho \cos(\theta)$$ $$ z = \rho \sin(\theta)$$

The implemention of Hirshman-Breslau representation in STELLOPT is introduced Sec. IV.B in that paper. For a given set of $RBC$ and $ZBS$ (only avalid for stellarator symmetry and only $\rho^s_{mn}$ is used), the covert is carried out as $$ m=0,1 \quad \rho(n,m) = \frac{1}{2} \left [\frac{RBC(n,m+1)+ZBS(n,m+1)}{t_1(m+1)} \right ] $$ $$ m\geq2 \quad \rho(n,m) = \frac{1}{4} \left[\frac{RBC(n,m+1)+ZBS(n,m+1)}{t_1(m+1)} + \frac{RBC(n,m-1)+ZBS(n,m-1)}{t_2(m-1)}\right] $$

where $t_1(m) = (\frac{m-1}{m})^p$, $t_2(m) = (\frac{m+1}{m})^p$. The exponential index $p$ can improve the convergence. $p=0$ is the polar representation and $p=1$ is the equal arc-length representation. Normally, $p=4$ is a good value.

The Hirshman-Breslau representation has a unique poloidal angle.

Numerical details

The LRHO_OPT(N,M) logical array controls the ability to utilize the Hirshman-Breslau representation of the VMEC boundary. For every LRHO_OPT set to true, only one variable is loaded if LASYM=F in the VMEC INDATA namelist. In order to transform to this representation the code must first take the VMEC boundary definition and convert is to Hirshman-Breslau. This is done by the code if any LRHO_OPT is set. STELLOPT will output a conversion accuracy message to the screen when this is done. The resulting spectrum in RHO will have one less poloidal mode than VMEC. Once the RHO(N,M) array is calculated, the modes are loaded into the optimization vector. The major radius is stored in the RHO(0,0) array element so if one wishes this variable to be varied LRHO_OPT(0,0)=T should be set. This representation does not treat the m=0 modes so if the user wishes STELLOPT to vary the m=0 modes the LBOUND_OPT(N,0) modes should be set to true (ignoring the N=0 mode). Setting any of these modes to false is equivalent to setting LFIX_NTOR(N) in the older version of STELLOPT to true. A DRHO_OPT(N,M) may also be set to control the associated auxiliary variable (see LMDIF options).

The variable RHO_EXP is the exponential index $p$ above. The default value is $4$. BOUND_MIN(N,M) and BOUND_MAX(N,M) define the lower and upper bounds for each boundary harmonic. The default behavior in the old version of STELLOPT was to set BOUND_MIN(N,M) = RHO(N,M)*0.3 and BOUND_MAX(N,M) = 2.0*RHO(N,M) with a check to make sure that BOUND_MAX(M,N) > BOUND_MIN(N,M) (if not true, the value were flipped).

Example

Using the same boundary as above, we turn on the option of converting into HB representation by adding LRHO_OPT(N.M) in the input file. In input.ncsx_hb, we now have

!-----------------------------------------------------------------------
!          OPTIMIZED QUANTITIES
!-----------------------------------------------------------------------
  LRHO_OPT(-5:5,0) = 11*T
  LRHO_OPT(-5:5,1) = 11*T
  LRHO_OPT(-5:5,2) = 11*T
  LRHO_OPT(-5:5,3) = 11*T
  LRHO_OPT(-5:5,4) = 11*T
  LRHO_OPT(-5:5,5) = 11*T
  LRHO_OPT(-5:5,6) = 11*T
  LRHO_OPT(-5:5,7) = 11*T
  LBOUND_OPT(1:4,0) = 4*T
  RHO_EXP = 1

Here we set RHO_EXP = 1 and it will use the equal arc-length angle. The converted boundary is shown below. It indicates that

  1. The convert is not performed with very high accurancy, since the original VMEC boundary is obtained with RHO_EXP = 4.
  2. The final angle is indeed a equal arc-length angle.
In [106]:
plt.figure()
plot_plasma_boundary('ncsx.boundary_eqarc', zeta=0, marker='o', npoints=64)
plot_plasma_boundary('ncsx.boundary_eqarc', zeta=np.pi/3, marker='o', color='b', npoints=64)
read 59 Fourier harmonics for plasma boundary and 0 for Bn distribution in ncsx.boundary_eqarc.
read 59 Fourier harmonics for plasma boundary and 0 for Bn distribution in ncsx.boundary_eqarc.
Out[106]:
(1.760784974483333,
 1.0411653636051406,
 0.29378019224460034,
 -0.29378019224460034)

Garabedian representation

General description

The Garabedian representation packs $R$ and $Z$ into a complex number $$ R + iZ = e^{2\pi iu} \sum \Delta_{mn} e^{-2\pi i (mu-nv)} .$$ Here, $R$ and $Z$ are the radial and axial components of cylindrical coordinates, $m$ and $n$ are the poloidal and toroidal mode numbers, and $u$ and $v$ are the normalized poloidal and toroidal angle-like variables, $0 \leq u< 1$ and $0 \leq v< 1$ ($2\pi u = \theta,2\pi v /N = \phi$, $\theta$ and $\phi$ being the poloidal and toroidal angles, respectively).

In this notation, the coefficient $\Delta_{00}$ is a measure of the plasma minor radius and $\Delta_{10}$ the major radius. It's usually imposed the normalization $\Delta_{00}=1.0$ on the small plasma radius. Conformal mapping shows that shape factors $\Delta_{mn}$ with large negative $m$ cause unrealistic cusps to penetrate the plasma, so we choose to leave them out when $m<-1$. The terms with $m =−1$, which define catenoids or crescents, are helpful because they contribute significantly to the magnetic well. Analytic geometry shows that the coefficients $\Delta_{1n}$ specify the helical excursion of the magnetic axis, whereas $\Delta_{2n}$ makes the plasma shape elliptical, $\Delta_{3n}$ makes it triangular, $\Delta_{4n}$ makes it rectangular, and so forth. Experience with the computations establishes that each $\Delta_{mn}$ has a strong influence on the corresponding coefficient Bmn in the magnetic spectrum, which makes it easier to design stellarators with approximate two-dimensional symmetry.

In Garabedian representation, the poloidal angle $u$ (or $\theta$) is not unique. Improved zoning can be obtained by performing an additional substitution $ u = u_1 - \sum Z_{mn} \sin(mu_1-nv) $ on the poloidal angle $u$. The poloidal angle in the Garabedian representation is not unique.

Numerical details

The LDELTAMN_OPT(N,M) logical array controls the ability to utilize the Garabedian representation of the VMEC boundary. This representation utilizes a rho coordinate to define the boundary. Thus for every LDELTAMN_OPT set to true, only one variable is loaded if LASYM=F in the VMEC INDATA namelist. In order to transform to this representation the code must first take the VMEC boundary definition and convert it to that of Garabedian. This is done by the code if any LDELTAMN_OPT is set. STELLOPT will output a conversion accuracy message to the screen when this is done. The resulting spectrum in RHO will have one less poloidal mode than VMEC. Once the DELTAMN(N,M) array is calculated, the modes are loaded into the optimization vector.

DELTA_MIN(N,M) and DELTA_MAX(N,M) define the lower and upper bounds for each boundary harmonic. A DDELTAMN_OPT(N,M) may also be set to control the associated auxiliary variable (see LMDIF options).

To convert VMEC representation to Garabedian representation, STELLOPT will call a subroutine convert_boundary_PG. What should be noted is that STELLOPT will store the majoir radius $RBC(0,0)$ into $\Delta_{00}$, but during the computation $\Delta_{00}=1.0$ is still imposed. Thus, $\Delta_{00}$ is not a free parameter. The dimension of DELTAMN is (-ntord:ntord,-mpol1d:mpol1d), almost twice as RBC, ZBS, and RHOBC.

Example

Using the same boundary as above, we turn on the option of converting into Garabedian representation by adding LDELTA_OPT(N.M) in the input file. In input.ncsx_garabedian, we now have

!-----------------------------------------------------------------------
!          OPTIMIZED QUANTITIES
!-----------------------------------------------------------------------
  LDELTAMN_OPT(-5:5,0) = 11*T
  LDELTAMN_OPT(-5:5,1) = 11*T
  LDELTAMN_OPT(-5:5,2) = 11*T
  LDELTAMN_OPT(-5:5,3) = 11*T
  LDELTAMN_OPT(-5:5,4) = 11*T
  LDELTAMN_OPT(-5:5,5) = 11*T
  LDELTAMN_OPT(-5:5,6) = 11*T
  LDELTAMN_OPT(-5:5,7) = 11*T

Since converting to the Garabedian representation doesn't change the poloidal angle, the converted shape is the same as original.

In [ ]: