<jupyter_start><jupyter_text>STELLOPT hands-on instructionSTELLOPT is written in Fortran, a historic language commonly used in scientific computing.Unlike Python, you need to compile the code first and then execute it.The source code of STELLOPT and other distributed parts can be found at https://github.com/PrincetonUniversity/STELLOPT.For convenience, we have installed STELLOPT and related software on the remote server.You can directly use the executable by typing `xstelloptv2` etc.STELLOPT uses the message passing interface (MPI) to accelerate the computation.You can specify the number of CPUs by `mpiexec -n 4 xstelloptv2`.More detailed documentation can be found on the wiki pages for STELLOPT (https://princetonuniversity.github.io/STELLOPT/).Although not required for this tutorial, you can learn more about  - Fortran: [Beginner Fortran 90 tutorial](https://people.ucsc.edu/~dlee79/2019/fall/am129_209/_downloads/19b82f919ae387db3090603199f2bd45/section1.pdf)  - Linux: [Introduction to Linux for HPC](https://portal.tacc.utexas.edu/c/document_library/get_file?uuid=a0c33bf1-b2a4-48b4-a23d-fc8a79c887ec&groupId=13601)  - MPI: [MPI turtorial](https://mpitutorial.com/tutorials/)  In this instruction, we are going to use my personal python package [coilpy](https://github.com/zhucaoxiang/CoilPy) for post-processing. There are two more comprehensive packages using GUI that are publicly available, [matlabVMEC](https://github.com/lazersos/matlabVMEC), [pySTEL](https://github.com/PrincetonUniversity/STELLOPT/tree/develop/pySTEL).  VMEC  [VMEC](https://princetonuniversity.github.io/STELLOPT/VMEC) uses a variational method to find the minimum of energy functional.As an equilibrium calculation code, VMEC takes three essential inputs:  1. plasma boundary   2. pressure profile  3. rotational transform or toroidal current density profileLike most Fortran codes, VMEC uses an input namelist to read these information (as well as others).  VMEC input namelistLet's open an example of VMEC input file `input.QAS`. VMEC will read all the variables in the namelist `INDATA`.<jupyter_code>less input.QAS<jupyter_output><empty_output><jupyter_text>Plasma boundaryVMEC uses Fourier series to represent the plasma boundary in cylinder coordinates $(R, \phi, Z)$,$$ 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. The related variables in the input file are| Variable | Type | Size | Description ||---|---|---|---|| nfp | integer | 1 | Number of toroidal field periods | | lasym | logical | 1 | Non-stellarator symmetric configuration (False) | | mpol | integer | 1 | Poloidal Mode Number (m) | | ntor | integer | 1 | Toroidal Mode Number (n) | | rbc | real(12,100) | -61:61,0:60 | Boundary cosine coefficients for $R= \cos(m\theta-n\zeta)$ | | rbs | real(12,100) | -61:61,0:60 | Boundary sine coefficients for R | | zbc | real(12,100) | -61:61,0:60 | Boundary cosine coefficients for Z | | zbs | real(12,100) | -61:61,0:60 | Boundary sine coefficients for Z|<jupyter_code># plot the boundary shape
import matplotlib.pyplot as plt
import numpy as np
from coilpy import FourSurf
qas = FourSurf.read_vmec_input('input.QAS')
qas.plot(zeta=0, label=r'$\phi=0$')
qas.plot(zeta=np.pi/6, label=r'$\phi=\pi/6$')
qas.plot(zeta=np.pi/3, label=r'$\phi=\pi/3$')
plt.legend()
# plot it as 3D surfa
qas.plot3d();
# adjust the shape
print(qas.xn[0], qas.xm[0], qas.rbc[0])
qas.plot(zeta=0, label=r'$r_{00} = 1.41$')
qas.rbc[0] = 1.6
qas.plot(zeta=0, label=r'$r_{00} = 1.6$')
plt.legend()<jupyter_output><empty_output><jupyter_text>Pressure profileThe pressure profile in VMEC is specified as a function of radial flux space coordinates (s).Where for the default case is the normalized toroidal flux (normalized to PHIEDGE). The AM parameter determines the polynomial coefficients (0..10) used to calculate pressure profile: $$ p=\sum_{n=0}^{10} am(n) * s^n . $$ The PRES_SCALE value is a scale factor applied to the profile allowing it to be scaled up and down.<jupyter_code># plot pressure profile
am = np.array([6.85517649352426E+04, -5.12027745123057E+03, -3.61510451745464E+04, -4.74263014113066E+05,
  1.78878195473870E+06, -3.21513828868170E+06,  2.69041023837233E+06, -8.17049854168367E+05,
  0.00000000000000E+00,  0.00000000000000E+00,  0.00000000000000E+00])
pressure = np.polynomial.polynomial.Polynomial(am)
flux = np.linspace(0,1,100)
plt.figure()
plt.plot(flux, pressure(flux))
plt.xlabel('normalized flux')
plt.ylabel('pressure [Pa]')<jupyter_output><empty_output><jupyter_text>Rotational transform or toroidal current density profileThe VMEC code provides the user an option to specify either a rotational transform radial profile or a toroidal current density profile. The NCURR parameter determine which form of the profile to use (0: Rotational Transform, 1: Toroidal Current Density). The AI parameter specifies the polynomial coefficients (0..10) used to calculate the rotational transform profile (NCURR=0) $$ \iota=\sum_{n=0}^{10} ai(n) * s^n . $$ The AC_FORM parameter determines the form of the current profiles used (NCURR=1). For AC_FORM=0 the toroidal current profile is power series in s defined by the AC parameter $$ j=\sum_{n=0}^{10} ac(n) * s^n . $$<jupyter_code># plot current profile
ac = np.array([8.18395699999999E+03,  1.43603560000000E+06, -1.07407140000000E+07,  7.44389200000000E+07,
 -3.22215650000000E+08,  8.81050800000000E+08, -1.49389660000000E+09,  1.52746800000000E+09,
 -8.67901590000000E+08,  2.10351200000000E+08,  0.00000000000000E+00])
current = np.polynomial.polynomial.Polynomial(ac)
flux = np.linspace(0,1,100)
plt.figure()
plt.plot(flux, current(flux))
plt.xlabel('normalized flux')
plt.ylabel('current [A]')<jupyter_output><empty_output><jupyter_text>Other input variablesThere are also other input variables required for a VMEC run.More details can be found in [VMEC input namelist](https://princetonuniversity.github.io/STELLOPT/VMEC%20Input%20Namelist%20(v8.47)), [Tutorial: Prepare VMEC Input Namelist](https://princetonuniversity.github.io/STELLOPT/Tutorial%20VMEC%20Input%20Namelist).  Run VMECNow run VMEC using```mpiexec -n 4 xvmec2000 input.QAS 2>&1 | tee log.vmec```On the screen, you should be able to see the output information.<jupyter_code>less log.vmec<jupyter_output><empty_output><jupyter_text>Check VMEC resultsYou can now check the VMEC outputs.When executed normally, VMEC will produce several files.  - `jxbout.QAS' contains values for various quantities on a grid throughout the simulation domain.  - 'mercier.QAS' contains radial profiles (radial index in VMEC is denoted by the variable `s`) of various quantities.   - 'threed1.QAS' can be considered an expanded log file where various quantities are calculated which were not output to the screen. This file is fairly self explanatory.   - 'wout_QAS.nc' is a [NetCDF](https://www.unidata.ucar.edu/software/netcdf/) file contains main data. It contains the Fourier Coefficients for the magnetic field along with various quantities.<jupyter_code># Read the netcdf file
import xarray
wout = 'wout_QAS.nc'
vmec = xarray.open_dataset(wout)
print(vmec)
# plot iota
plt.figure()
plt.plot(np.linspace(0, 1, vmec['ns'].values), vmec['iotaf'])
# plot |B|
from coilpy import VMECout
qas = VMECout(wout)
plt.figure()
ns = 10 # surface label [0,127]
plt.imshow(qas.data['b'][10], origin='lower', extent=[0, 2*np.pi/3, 0, 2*np.pi], aspect='auto')
plt.colorbar()
# plot flux surface
zeta = 0
plt.figure()
for isurf in qas.surface[::8]:
    isurf.plot(zeta=zeta)<jupyter_output><empty_output><jupyter_text>STELLOPT  STELLOPT input namelistLike VMEC, [STELLOPT](https://princetonuniversity.github.io/STELLOPT/STELLOPT) reads input variables from the Fortran namelist [`&OPTIMUM`](https://princetonuniversity.github.io/STELLOPT/STELLOPT%20Input%20Namelist).The two namelists have to placed into the same file.STELLOPT has a long list of input variables.Here we are going to only introduce some of them. Runtime controlThe following parameters control how the code runs. Depending on the type of optimization the parameters may change functionality.```&OPTIMUM!------------------------------------------------------------------------!       Optimizer Run Control Parameters!------------------------------------------------------------------------  NFUNC_MAX    = 5000                       ! Maximum number of function evaluations  EQUIL_TYPE   = 'VMEC2000'                 ! Equilibrium Code VMEC2000  OPT_TYPE     = 'LMDIF'                    ! Optimization Type (LMDIF),GADE,MAP,PSO  NOPTIMIZERS  = -1                         ! Number of optimizers  FTOL         = 1.0E-4                     ! Absolute Tolerance  XTOL         = 1.0E-4                     ! Relative Tolerance  GTOL         = 1.0E-30                    ! Orthogonality Tolerance  EPSFCN       = 1.0E-4                     ! Finite Difference Stepsize  FACTOR       = 100.0                      ! Initial step scaling  MODE         = 1                          ! Mode Parameter  LKEEP_MINS   = T                          ! Logical to keep minimum states``` VariablesSpecify the variables you want to vary.By default, all the variables will be fixed.The variables are named in the style of `LVAR_OPT` (`VAR` is variable name).If the optimizer supports, you can specify the upper and lower bounds for the variable, `VAR_MIN`, `VAR_MAX`, with a scaling value `DVAR`.```!------------------------------------------------------------------------!       Optimized Quantities!------------------------------------------------------------------------  LPHIEDGE_OPT  = T    PHIEDGE_MIN = 0.01    PHIEDGE_MAX = 1.00    DPHIEDGE    = 1.0  LCURTOR_OPT   = T  LPSCALE_OPT   = T  LBCRIT_OPT    = T  LEXTCUR_OPT   = T T T T  LAPHI_OPT     = T T T T  LAM_OPT       = T T T T  LAC_OPT       = F T F T  LAI_OPT       = T F T F  LBOUND_OPT(1,0) = T      ! Optimize RBC,ZBS Arrays  LRHO_OPT(1,0) = T        ! Optimize RHOMN (Hirshman/Breslau representation)  LDELTAMN_OPT(1,0) = T    ! Optimize DELTAMN (Garabedian representation) ```The three different boundary representations will be introduced in Project_02 in the computer labs tomorrow.  TargetsThere are various chi-squared functionals and each is a unique set of parameters. In general for each target value there is a sigma associated with it. Here they are some examples.```!------------------------------------------------------------------------!       Equilibrium / Geometry Optimization Targets!------------------------------------------------------------------------  TARGET_PHIEDGE = 2.5    SIGMA_PHIEDGE = 0.025   ! Enclosed Toroidal Flux [Wb]  TARGET_CURTOR  = 1.0E6  SIGMA_CURTOR  = 1.0E3   ! Total Toroidal Current [A]  TARGET_CURVATURE = 1.0E-3  SIGMA_CURVATURE = 1.0 ! Flux surface curvature  TARGET_RBTOR   = 7.2    SIGMA_RBTOR   = 0.01    ! R*Btor [T-m]  TARGET_R0      = 3.6    SIGMA_R0      = 0.01    ! Magnetic Axis R (phi=0) [m]  TARGET_Z0      = 0.5    SIGMA_Z0      = 0.01    ! Magnetic Axis Z (phi=0) [m]  TARGET_VOLUME  = 1.0E2  SIGMA_VOLUME  = 1.0     ! Plasma Volume [m^-3]  TARGET_BETA    = 0.02   SIGMA_BETA    = 0.0001  ! Average Plasma Beta  TARGET_WP      = 1.0E3  SIGMA_WP      = 1.0     ! Stored Energy [J]  TARGET_ASPECT  = 5.0    SIGMA_ASPECT  = 0.5     ! Aspect Ratio (R/a)!------------------------------------------------------------------------!       Boozer Coordinate Helicity!         Note that helicity targeting is by surface.  Axis (01) is ignored.!         (X,0): Quasi-Axisymetry!         (0,X): Quasi-Poloidal Symmetry!         (L,K): Quasi-Helical Symmetry (m *K + n*L)!------------------------------------------------------------------------  HELICITY = (1,0)  TARGET_HELICITY(02:128) = 127*0.0  SIGMA_HELICITY(02:128) = 127*0.01!------------------------------------------------------------------------!       Ballooning Stability (as calculated by COBRA_VMEC)!         Note that ballooning stability is by surface.  Axis (01) is ignored.!         THETA, ZETA: Ballooning angle perturbations!------------------------------------------------------------------------  BALLOON_THETA = 0.0  3.14  5.50  BALLOON_ZETA    = 0.0 3.14  5.50  TARGET_BALLOON(02:128) = 127*0.0   SIGMA_BALLOON(02:128) = 127*0.2!------------------------------------------------------------------------!       Neoclassical Transport Calculation (as calculated by NEO)!         Note that neoclassical transport is by surface. Axis (01) is ignored.!------------------------------------------------------------------------  TARGET_NEO(02:128) = 127*0.0  SIGMA_NEO(02:128) = 127*0.1```  One iteration evaluationLet's first just do a "one-iteration" optimization, by using `OPT_TYPE = 'ONE_ITER'`.We have pasted the above namelist and append it at the end of `input.QAS`.Now, we can evaluate the properties of the equilibrium.```mpiexec -n 4 xstelloptv2 input.QAS 2>&1 | tee log.stellopt```<jupyter_code>less log.stellopt<jupyter_output><empty_output><jupyter_text>Check STELLOPT resultsOnce a run has been completed, STELLOPT will write a ASCII file `stellopt.QAS` containing the main data.STELLOPT will also save all the intermediate outputs from sub-codes, like `wout_QAS.00000.nc`.<jupyter_code># Read stellopt.QAS file
from coilpy import STELLout
stell_qas = STELLout('stellopt.QAS')

# plot the convergence
stell_qas.plot()
# plot Ballooning stability
stell_qas.plot_balloon()
# Read BOOZ_XFORM output
from coilpy import BOOZ_XFORM
booz = BOOZ_XFORM('boozmn_QAS.00000.nc')

# plot non-axisymmetric terms
booz.plot()
# plot |B| contour on one surface
booz.plot2d(ns=10, contour=True)
# plot epsilon_eff
plt.figure()
plt.semilogy(stell_qas['NEO_K'][0], stell_qas['NEO_equil'][0])<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.16246215943491424, 'fraction_numerical': 0.058094277065013696, 'mean_word_length': 3.057017543859649, 'pattern_counts': {'":': 0, '<': 33, '<?xml version=': 0, '>': 36, 'https://': 14, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 14, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '156167', 'license': 'MIT', 'n_tokens_mistral': 5025, 'n_tokens_neox': 4350, 'n_words': 1347, 'path': '00_stellopt_instructions/00_stellopt_instructions.ipynb', 'repo_name': 'zhucaoxiang/stellopt_labs', 'size': '13874'}