Basic Equilibrium
DESC is a 3D MHD equilibrium and optimization code suite, which solves the ideal MHD equilibrium equations to find stellarator equilibria.
Like VMEC, DESC requires 4 main inputs to define the equilibrium problem:
Pressure Profile
Rotational Transform or Toroidal Current Profile
Last Closed Flux Surface Boundary Shape
Total toroidal magnetic flux enclosed by the LCFS
DESC can be run both with a text input file from the command line (e.g. python -m desc INPUT_FILE) or through python scripts (which offers more options than the text file for equilibrium solving and optimization). This tutorial notebook will focus on basic functionality for solving an equilibrium.
Installing DESC
The installation instructions are located in the installation documentation page.
Once these instructions are followed and DESC is installed correctly, you should be able to import a module from desc such as desc.io without any errors.
[1]:
import sys
import os
sys.path.insert(0, os.path.abspath("."))
sys.path.append(os.path.abspath("../../../"))
[2]:
import desc.io
DESC version 0.12.3+56.ge61501fcc.dirty,using JAX backend, jax version=0.4.31, jaxlib version=0.4.30, dtype=float64
Using device: CPU, with 3.48 GB available memory
Running DESC from a VMEC input
Perhaps the simplest way to run DESC if you already have a VMEC input file for your equilibrium is to simply pass that to DESC on the command line.
In this directory is an input.HELIOTRON VMEC file which contains a 19 field-period, finite beta HELIOTRON fixed boundary stellarator. To convert the input, we simply run DESC from the command line with
python -m desc -vv input.HELIOTRON -o input.HELIOTRON_output.h5.
You can leave out python -m if you’ve added DESC to your python path. The -vv flag is for double verbosity. The -o flag specifies the filename of the output HDF5 file (XXX_output.h5 by default). See the command line docs for more information.
The conversion creates a new file with _desc appended to the file name. The above command automatically converts the VMEC input file to a DESC input file and begins solving it. There are many DESC solver options without VMEC analogs, such as the multi-grid continuation method steps, so DESC will automatically choose a conservative continuation method when ran in this way, which is generally recommended for most equilibria to ensure robust convergence.
Generally, the most important parameters to tweak are related to the spectral resolution and the continuation method. DESC leverages a multi-grid continuation method to allow for robust convergence of highly shaped equilibria. The following parameters control this method:
Spectral Resolution
L_rad(int): Maximum radial mode number for the Fourier-Zernike basis.default =
M_polifspectral_indexing = ANSIdefault =
2*M_polifspectral_indexing = FringeFor more information see Basis functions and collocation nodes.
M_pol(int): Maximum poloidal mode number for the Fourier-Zernike basis. Required input.N_tor(int): Maximum toroidal mode number for the Fourier-Zernike basis.default = 0
L_grid(int): Radial resolution of nodes in collocation grid.default =
M_gridifspectral_indexing = ANSIdefault =
2*M_gridifspectral_indexing = Fringe
M_grid(int): Poloidal resolution of nodes in collocation grid.default =
2*M_pol
N_grid(int): Toroidal resolution of nodes in collocation grid.default =
2*N_tor
When M_grid = M_pol, the number of collocation nodes in each toroidal cross-section is equal to the number of Zernike polynomials in a FourierZernike basis set of the same resolution L_rad = M_pol. When N_grid = N_tor, the number of nodes with unique toroidal angles is equal to the number of terms in the toroidal Fourier series. Convergence is typically superior when the number of nodes exceeds the number of spectral coefficients, so by default the collocation grids are oversampled
with respect to the spectral resolutions.
Loading the results
DESC provides utility functions to load and compare equilibria. These will be covered in detail later in the tutorial. For now, notice that DESC saves solutions as EquilibriaFamily objects. An EquilibriaFamily is essentially a list of equilibria, which can be indexed to retrieve individual equilibria. Higher indices hold equilibria solved later in the continuation process. You can retrieve the final state of a DESC equilibrium solve with eq = desc.io.load("XXX_output.h5")[-1].
[3]:
%matplotlib inline
from desc.plotting import plot_comparison, plot_section
[4]:
eq_fam = desc.io.load("input.HELIOTRON_output.h5")
print("Number of equilibria in the EquilibriaFamily:", len(eq_fam))
fig, ax = plot_comparison(
eqs=[eq_fam[1], eq_fam[3], eq_fam[-1]],
labels=[
"Axisymmetric w/o pressure",
"Axisymmetric w/ pressure",
"Nonaxisymmetric w/ pressure",
],
)
Number of equilibria in the EquilibriaFamily: 8
Creating an Equilibrium from scratch
The recommended way to work with DESC is through python scripts, in which we construct and solve an Equilibrium object.
We initialize the Equilibrium with the desired resolution, and the inputs specified above: - Boundary shape, as a FourierRZToroidalSurface - Profiles, as Profile objects (there are several different options, the most common is PowerSeriesProfile) - Total flux Psi as a float
[5]:
from desc.equilibrium import Equilibrium
from desc.geometry import FourierRZToroidalSurface
from desc.profiles import PowerSeriesProfile
When you already have the VMEC or DESC input file, you can instantiate an Equilibrium directly from the input file, using Equilibrium.from_input_file,
[6]:
eq = Equilibrium.from_input_file("input.HELIOTRON")
Converting VMEC input to DESC input
Generated DESC input file input.HELIOTRON_desc:
Or, you may extract just the boundary from the input file with the FourierRZToroidalSurface.from_input_file:
[7]:
surf1 = FourierRZToroidalSurface.from_input_file("input.HELIOTRON")
Converting VMEC input to DESC input
Generated DESC input file input.HELIOTRON_desc:
If you do not have an input file to work from, you can see the below steps to create an Equilibrium from scratch.
When starting from scratch, you can construct the surface by specifying the Fourier coefficients and mode numbers manually.
The boundary is represented by a double Fourier series for R and Z in terms of a poloidal angle \(\theta\) and the geometric toroidal angle \(\zeta\). We specify the mode numbers for R and Z as 2D arrays of [m,n] pairs, and the coefficients as a 1D array.
In DESC the double Fourier series for R and Z are defined in a slightly different manner than VMEC:
\[\begin{split}\begin{align}
R^b(\theta,\zeta) &= \sum^M \sum^N R^b_{mn} \mathcal{G}^{m}_{n}(\theta,\zeta)\\
Z^b(\theta,\zeta) &= \sum^M \sum^N Z^b_{mn} \mathcal{G}^{m}_{n}(\theta,\zeta)
\end{align}\end{split}\]
where
\[\begin{split}\mathcal{G}^{m}_{n}(\theta,\zeta) = \begin{cases}
\cos(|m|\theta)\cos(|n|N_{FP}\zeta) &\text{for }m\ge0, n\ge0 \\
\cos(|m|\theta)\sin(|n|N_{FP}\zeta) &\text{for }m\ge0, n<0 \\
\sin(|m|\theta)\cos(|n|N_{FP}\zeta) &\text{for }m<0, n\ge0 \\
\sin(|m|\theta)\sin(|n|N_{FP}\zeta) &\text{for }m<0, n<0.
\end{cases}\end{split}\]
Note: in DESC, radial modes are indexed by l, poloidal modes by m, and toroidal modes by n
[8]:
surf2 = FourierRZToroidalSurface(
R_lmn=[10.0, -1.0, -0.3, 0.3],
modes_R=[
(0, 0),
(1, 0),
(1, 1),
(-1, -1),
], # (m,n) pairs corresponding to R_mn on previous line
Z_lmn=[1, -0.3, -0.3],
modes_Z=[(-1, 0), (-1, 1), (1, -1)],
NFP=19,
)
For profiles, we will use the standard PowerSeriesProfile which represents a function as a power series in the radial coordinate \(\rho\) which is the square root of the normalized toroidal flux: \(\rho = \sqrt{\Psi/\Psi_b}\) (Note this is different from VMEC which uses the normalized toroidal flux without the square root). We could also use splines (SplineProfile), or Zernike polynomials (FourierZernikeProfile), for more options, see the `desc.profiles
module <https://desc-docs.readthedocs.io/en/stable/api.html#profiles>`__
[9]:
pressure = PowerSeriesProfile(
[1.8e4, 0, -3.6e4, 0, 1.8e4]
) # coefficients in ascending powers of rho
iota = PowerSeriesProfile([1, 0, 1.5]) # 1 + 1.5 r^2
Finally, we create the equilibrium using the objects we just made, leaving the other parameters at their defaults. The number of field periods by default is inferred from the given surface.
[10]:
eq = Equilibrium(
L=8, # radial resolution
M=8, # poloidal resolution
N=3, # toroidal resolution
surface=surf1,
pressure=pressure,
iota=iota,
Psi=1.0, # total flux, in Webers
)
We now have an Equilibrium object, but it is generally not actually in equilibrium, as the default initial guess is just to scale the outermost flux surface.
To find the actual solution, we must solve the equilibrium. There are two primary ways to do this in DESC:
eq.solveis a method of theEquilibriumclass that is best used when starting close to the correct solution, or for refining a solution after optimization, for example. Generally not recommended for “cold start” solves.desc.continuation.solve_continuation_automaticuses a multi-grid continuation method to arrive at a particular desired equilibrium. It is generally very robust, and is the recommended method.
Here, we will use both methods
[11]:
eq1, info = eq.solve(verbose=3, copy=True)
Building objective: force
Precomputing transforms
Timer: Precomputing transforms = 533 ms
Timer: Objective build = 822 ms
Building objective: lcfs R
Building objective: lcfs Z
Building objective: fixed Psi
Building objective: fixed pressure
Building objective: fixed iota
Building objective: fixed sheet current
Building objective: self_consistency R
Building objective: self_consistency Z
Building objective: lambda gauge
Building objective: axis R self consistency
Building objective: axis Z self consistency
Timer: Objective build = 591 ms
Timer: Linear constraint projection build = 3.00 sec
Number of parameters: 351
Number of objectives: 2106
Timer: Initializing the optimization = 4.43 sec
Starting optimization
Using method: lsq-exact
Iteration Total nfev Cost Cost reduction Step norm Optimality
0 1 9.095e-01 3.829e+01
1 2 7.594e-02 8.335e-01 2.539e-01 1.197e+01
2 3 5.712e-03 7.023e-02 3.003e-01 1.557e+00
3 4 4.117e-03 1.596e-03 1.775e-01 5.823e+00
4 5 8.399e-04 3.277e-03 1.178e-01 1.623e+00
5 6 1.741e-04 6.658e-04 7.147e-02 6.174e-01
6 7 9.597e-05 7.817e-05 3.119e-02 1.320e-01
7 8 9.373e-05 2.236e-06 1.633e-02 8.389e-02
8 9 9.260e-05 1.131e-06 1.821e-02 3.462e-02
9 10 9.228e-05 3.286e-07 1.102e-02 3.377e-02
Optimization terminated successfully.
`ftol` condition satisfied.
Current function value: 9.228e-05
Total delta_x: 6.849e-01
Iterations: 9
Function evaluations: 10
Jacobian evaluations: 10
Timer: Solution time = 5.83 sec
Timer: Avg time per step = 583 ms
==============================================================================================================
Start --> End
Total (sum of squares): 9.095e-01 --> 9.228e-05,
Maximum absolute Force error: 1.584e+06 --> 5.903e+04 (N)
Minimum absolute Force error: 2.242e+01 --> 1.145e+00 (N)
Average absolute Force error: 2.478e+05 --> 2.314e+03 (N)
Maximum absolute Force error: 1.274e-01 --> 4.747e-03 (normalized)
Minimum absolute Force error: 1.803e-06 --> 9.209e-08 (normalized)
Average absolute Force error: 1.993e-02 --> 1.861e-04 (normalized)
R boundary error: 0.000e+00 --> 0.000e+00 (m)
Z boundary error: 0.000e+00 --> 0.000e+00 (m)
Fixed Psi error: 0.000e+00 --> 0.000e+00 (Wb)
Fixed pressure profile error: 0.000e+00 --> 0.000e+00 (Pa)
Fixed iota profile error: 0.000e+00 --> 0.000e+00 (dimensionless)
Fixed sheet current error: 0.000e+00 --> 0.000e+00 (~)
==============================================================================================================
[12]:
from desc.continuation import solve_continuation_automatic
eqf = solve_continuation_automatic(eq.copy(), verbose=3)
================
Step 1
================
Spectral indexing: ansi
Spectral resolution (L,M,N)=(6,6,0)
Node resolution (L,M,N)=(12,12,0)
Boundary ratio = 0
Pressure ratio = 0
Perturbation Order = 2
Objective: force
Optimizer: lsq-exact
================
Building objective: force
Precomputing transforms
Timer: Precomputing transforms = 697 ms
Timer: Objective build = 750 ms
Building objective: lcfs R
Building objective: lcfs Z
Building objective: fixed Psi
Building objective: fixed pressure
Building objective: fixed iota
Building objective: fixed sheet current
Building objective: self_consistency R
Building objective: self_consistency Z
Building objective: lambda gauge
Building objective: axis R self consistency
Building objective: axis Z self consistency
Timer: Objective build = 392 ms
Timer: Linear constraint projection build = 2.04 sec
Number of parameters: 27
Number of objectives: 98
Timer: Initializing the optimization = 3.21 sec
Starting optimization
Using method: lsq-exact
Iteration Total nfev Cost Cost reduction Step norm Optimality
0 1 1.513e-02 1.947e+00
1 2 5.311e-03 9.815e-03 2.956e-01 5.699e-01
2 3 4.208e-04 4.890e-03 1.389e-01 1.699e-01
3 4 1.747e-05 4.033e-04 2.667e-02 2.859e-02
4 5 1.014e-06 1.646e-05 1.331e-02 4.222e-03
5 6 1.366e-07 8.775e-07 4.450e-03 5.429e-04
6 7 1.026e-07 3.402e-08 4.940e-03 9.463e-04
7 8 9.623e-08 6.334e-09 4.969e-03 1.430e-04
8 9 9.547e-08 7.556e-10 9.957e-04 1.413e-05
Optimization terminated successfully.
`ftol` condition satisfied.
Current function value: 9.547e-08
Total delta_x: 1.608e-01
Iterations: 8
Function evaluations: 9
Jacobian evaluations: 9
Timer: Solution time = 2.97 sec
Timer: Avg time per step = 330 ms
==============================================================================================================
Start --> End
Total (sum of squares): 1.513e-02 --> 9.547e-08,
Maximum absolute Force error: 2.067e+05 --> 9.478e+02 (N)
Minimum absolute Force error: 5.178e+01 --> 8.292e-01 (N)
Average absolute Force error: 3.277e+04 --> 7.425e+01 (N)
Maximum absolute Force error: 1.663e-02 --> 7.623e-05 (normalized)
Minimum absolute Force error: 4.164e-06 --> 6.669e-08 (normalized)
Average absolute Force error: 2.636e-03 --> 5.972e-06 (normalized)
R boundary error: 0.000e+00 --> 0.000e+00 (m)
Z boundary error: 0.000e+00 --> 0.000e+00 (m)
Fixed Psi error: 0.000e+00 --> 0.000e+00 (Wb)
Fixed pressure profile error: 0.000e+00 --> 0.000e+00 (Pa)
Fixed iota profile error: 0.000e+00 --> 0.000e+00 (dimensionless)
Fixed sheet current error: 0.000e+00 --> 0.000e+00 (~)
==============================================================================================================
Timer: Iteration 1 total = 7.25 sec
================
Step 2
================
Spectral indexing: ansi
Spectral resolution (L,M,N)=(8,8,0)
Node resolution (L,M,N)=(16,16,0)
Boundary ratio = 0
Pressure ratio = 0
Perturbation Order = 2
Objective: force
Optimizer: lsq-exact
================
Building objective: force
Precomputing transforms
Timer: Precomputing transforms = 676 ms
Timer: Objective build = 729 ms
Building objective: lcfs R
Building objective: lcfs Z
Building objective: fixed Psi
Building objective: fixed pressure
Building objective: fixed iota
Building objective: fixed sheet current
Building objective: self_consistency R
Building objective: self_consistency Z
Building objective: lambda gauge
Building objective: axis R self consistency
Building objective: axis Z self consistency
Timer: Objective build = 442 ms
Timer: Linear constraint projection build = 1.91 sec
Number of parameters: 48
Number of objectives: 162
Timer: Initializing the optimization = 3.11 sec
Starting optimization
Using method: lsq-exact
Iteration Total nfev Cost Cost reduction Step norm Optimality
0 1 9.241e-08 2.600e-03
1 3 2.638e-08 6.603e-08 1.086e-02 1.154e-03
2 4 1.398e-08 1.240e-08 1.879e-02 1.407e-03
3 5 3.708e-09 1.027e-08 1.679e-02 1.396e-03
4 7 1.869e-10 3.521e-09 2.535e-03 5.663e-04
5 9 1.299e-11 1.740e-10 1.235e-03 1.540e-04
6 11 1.396e-12 1.159e-11 6.316e-04 4.710e-05
7 13 2.847e-13 1.112e-12 3.275e-04 1.608e-05
8 15 2.015e-13 8.317e-14 1.647e-04 3.824e-06
9 17 1.966e-13 4.902e-15 8.176e-05 9.518e-07
10 19 1.961e-13 5.381e-16 4.073e-05 2.366e-07
Optimization terminated successfully.
`ftol` condition satisfied.
Current function value: 1.961e-13
Total delta_x: 4.820e-02
Iterations: 10
Function evaluations: 19
Jacobian evaluations: 11
Timer: Solution time = 2.76 sec
Timer: Avg time per step = 250 ms
==============================================================================================================
Start --> End
Total (sum of squares): 9.241e-08 --> 1.961e-13,
Maximum absolute Force error: 1.359e+03 --> 8.676e-01 (N)
Minimum absolute Force error: 3.405e-01 --> 9.701e-05 (N)
Average absolute Force error: 7.427e+01 --> 1.040e-01 (N)
Maximum absolute Force error: 1.093e-04 --> 6.977e-08 (normalized)
Minimum absolute Force error: 2.738e-08 --> 7.802e-12 (normalized)
Average absolute Force error: 5.973e-06 --> 8.366e-09 (normalized)
R boundary error: 0.000e+00 --> 0.000e+00 (m)
Z boundary error: 0.000e+00 --> 0.000e+00 (m)
Fixed Psi error: 0.000e+00 --> 0.000e+00 (Wb)
Fixed pressure profile error: 0.000e+00 --> 0.000e+00 (Pa)
Fixed iota profile error: 0.000e+00 --> 0.000e+00 (dimensionless)
Fixed sheet current error: 0.000e+00 --> 0.000e+00 (~)
==============================================================================================================
Timer: Iteration 2 total = 7.72 sec
================
Step 3
================
Spectral indexing: ansi
Spectral resolution (L,M,N)=(8,8,0)
Node resolution (L,M,N)=(16,16,0)
Boundary ratio = 0
Pressure ratio = 0.5
Perturbation Order = 2
Objective: force
Optimizer: lsq-exact
================
Perturbing equilibrium
Building objective: force
Precomputing transforms
Timer: Precomputing transforms = 30.2 ms
Timer: Objective build = 38.1 ms
Building objective: lcfs R
Building objective: lcfs Z
Building objective: fixed Psi
Building objective: fixed pressure
Building objective: fixed iota
Building objective: fixed sheet current
Building objective: self_consistency R
Building objective: self_consistency Z
Building objective: lambda gauge
Building objective: axis R self consistency
Building objective: axis Z self consistency
Timer: Objective build = 63.2 ms
Perturbing p_l
Factorizing linear constraints
Timer: linear constraint factorize = 23.3 ms
Computing df
Timer: df computation = 2.88 sec
Factoring df
Timer: df/dx factorization = 21.3 ms
Computing d^2f
Timer: d^2f computation = 2.78 sec
Timer: Objective build = 960 us
||dx||/||x|| = 3.670e-02
Timer: Total perturbation = 7.11 sec
Building objective: self_consistency R
Building objective: self_consistency Z
Building objective: lambda gauge
Building objective: axis R self consistency
Building objective: axis Z self consistency
Timer: Objective build = 8.44 ms
Timer: Linear constraint projection build = 22.7 ms
Number of parameters: 48
Number of objectives: 162
Timer: Initializing the optimization = 36.2 ms
Starting optimization
Using method: lsq-exact
Iteration Total nfev Cost Cost reduction Step norm Optimality
0 1 2.654e-01 1.197e+01
1 2 2.078e-02 2.446e-01 3.706e-01 4.065e+00
2 4 5.501e-04 2.023e-02 1.044e-01 6.551e-01
3 6 1.087e-04 4.414e-04 4.467e-02 3.784e-01
4 8 7.638e-05 3.231e-05 2.072e-02 7.234e-02
5 9 5.247e-05 2.391e-05 3.517e-02 1.443e-01
6 11 3.705e-05 1.543e-05 1.827e-02 3.400e-02
7 12 2.867e-05 8.377e-06 3.512e-02 1.177e-01
8 13 2.265e-05 6.015e-06 3.423e-02 1.380e-01
9 14 1.989e-05 2.764e-06 3.048e-02 1.633e-01
10 15 7.558e-06 1.233e-05 8.505e-03 1.059e-02
11 16 7.037e-06 5.213e-07 1.335e-02 4.062e-02
12 17 6.098e-06 9.386e-07 1.288e-02 3.656e-02
13 18 5.338e-06 7.606e-07 1.304e-02 3.273e-02
14 19 4.699e-06 6.390e-07 1.323e-02 3.296e-02
15 20 4.147e-06 5.522e-07 1.343e-02 3.314e-02
16 21 3.662e-06 4.847e-07 1.365e-02 3.323e-02
17 22 3.234e-06 4.279e-07 1.389e-02 3.328e-02
18 23 2.856e-06 3.784e-07 1.414e-02 3.331e-02
19 24 2.521e-06 3.344e-07 1.440e-02 3.333e-02
20 25 2.226e-06 2.951e-07 1.465e-02 3.333e-02
21 26 1.967e-06 2.596e-07 1.488e-02 3.332e-02
22 27 1.739e-06 2.272e-07 1.510e-02 3.330e-02
23 28 1.541e-06 1.979e-07 1.529e-02 3.327e-02
24 29 1.370e-06 1.715e-07 1.548e-02 3.323e-02
25 30 1.222e-06 1.481e-07 1.568e-02 3.317e-02
26 31 1.094e-06 1.274e-07 1.589e-02 3.310e-02
27 32 9.857e-07 1.088e-07 1.613e-02 3.298e-02
28 33 8.944e-07 9.124e-08 1.641e-02 3.284e-02
29 34 8.207e-07 7.367e-08 1.668e-02 3.269e-02
30 35 5.445e-07 2.763e-07 4.357e-03 1.964e-03
31 36 5.297e-07 1.474e-08 8.624e-03 8.084e-03
32 37 5.034e-07 2.632e-08 8.529e-03 8.089e-03
33 38 4.800e-07 2.344e-08 8.497e-03 8.107e-03
34 39 4.589e-07 2.112e-08 8.414e-03 8.129e-03
35 40 4.398e-07 1.901e-08 8.342e-03 8.148e-03
36 41 4.228e-07 1.708e-08 8.268e-03 8.168e-03
37 42 4.075e-07 1.529e-08 8.202e-03 8.186e-03
38 43 3.939e-07 1.360e-08 8.140e-03 8.205e-03
39 44 3.819e-07 1.199e-08 8.081e-03 8.224e-03
40 45 3.714e-07 1.048e-08 8.024e-03 8.243e-03
41 46 3.624e-07 9.037e-09 7.967e-03 8.263e-03
42 47 3.547e-07 7.667e-09 7.908e-03 8.285e-03
43 48 3.483e-07 6.360e-09 7.847e-03 8.309e-03
44 49 3.277e-07 2.065e-08 2.057e-03 5.107e-04
45 50 3.262e-07 1.451e-09 3.919e-03 2.080e-03
Optimization terminated successfully.
`ftol` condition satisfied.
Current function value: 3.262e-07
Total delta_x: 2.786e-01
Iterations: 45
Function evaluations: 50
Jacobian evaluations: 46
Timer: Solution time = 2.28 sec
Timer: Avg time per step = 49.6 ms
==============================================================================================================
Start --> End
Total (sum of squares): 2.654e-01 --> 3.262e-07,
Maximum absolute Force error: 7.301e+05 --> 1.065e+03 (N)
Minimum absolute Force error: 6.385e+02 --> 3.135e-01 (N)
Average absolute Force error: 9.822e+04 --> 1.337e+02 (N)
Maximum absolute Force error: 5.872e-02 --> 8.564e-05 (normalized)
Minimum absolute Force error: 5.135e-05 --> 2.521e-08 (normalized)
Average absolute Force error: 7.899e-03 --> 1.075e-05 (normalized)
R boundary error: 0.000e+00 --> 0.000e+00 (m)
Z boundary error: 0.000e+00 --> 0.000e+00 (m)
Fixed Psi error: 0.000e+00 --> 0.000e+00 (Wb)
Fixed pressure profile error: 0.000e+00 --> 0.000e+00 (Pa)
Fixed iota profile error: 0.000e+00 --> 0.000e+00 (dimensionless)
Fixed sheet current error: 0.000e+00 --> 0.000e+00 (~)
==============================================================================================================
Timer: Iteration 3 total = 10.4 sec
================
Step 4
================
Spectral indexing: ansi
Spectral resolution (L,M,N)=(8,8,0)
Node resolution (L,M,N)=(16,16,0)
Boundary ratio = 0
Pressure ratio = 1.0
Perturbation Order = 2
Objective: force
Optimizer: lsq-exact
================
Perturbing equilibrium
Building objective: force
Precomputing transforms
Timer: Precomputing transforms = 52.5 ms
Timer: Objective build = 82.4 ms
Building objective: lcfs R
Building objective: lcfs Z
Building objective: fixed Psi
Building objective: fixed pressure
Building objective: fixed iota
Building objective: fixed sheet current
Building objective: self_consistency R
Building objective: self_consistency Z
Building objective: lambda gauge
Building objective: axis R self consistency
Building objective: axis Z self consistency
Timer: Objective build = 69.2 ms
Perturbing p_l
Factorizing linear constraints
Timer: linear constraint factorize = 482 ms
Computing df
Timer: df computation = 5.67 ms
Factoring df
Timer: df/dx factorization = 646 us
Computing d^2f
Timer: d^2f computation = 1.46 ms
Timer: Objective build = 810 us
||dx||/||x|| = 2.432e-02
Timer: Total perturbation = 578 ms
Building objective: self_consistency R
Building objective: self_consistency Z
Building objective: lambda gauge
Building objective: axis R self consistency
Building objective: axis Z self consistency
Timer: Objective build = 7.30 ms
Timer: Linear constraint projection build = 24.9 ms
Number of parameters: 48
Number of objectives: 162
Timer: Initializing the optimization = 37.2 ms
Starting optimization
Using method: lsq-exact
Iteration Total nfev Cost Cost reduction Step norm Optimality
0 1 6.139e-04 9.516e-01
1 2 6.380e-05 5.501e-04 7.726e-02 2.453e-01
2 3 6.783e-06 5.702e-05 3.730e-02 5.517e-02
3 4 5.567e-06 1.216e-06 1.276e-02 2.213e-02
4 6 4.798e-06 7.693e-07 1.009e-02 5.617e-03
5 7 4.732e-06 6.511e-08 1.596e-02 1.489e-02
6 8 4.652e-06 8.079e-08 1.353e-02 2.256e-02
7 9 4.593e-06 5.863e-08 1.466e-02 2.231e-02
8 10 4.420e-06 1.729e-07 1.537e-02 2.077e-02
9 11 4.343e-06 7.704e-08 1.522e-02 2.135e-02
10 12 4.295e-06 4.846e-08 1.453e-02 2.350e-02
11 13 4.278e-06 1.699e-08 1.415e-02 2.378e-02
12 14 4.170e-06 1.073e-07 4.172e-03 1.303e-03
13 15 4.170e-06 2.917e-10 4.594e-03 2.356e-03
14 16 4.169e-06 1.215e-09 7.249e-04 3.884e-05
Optimization terminated successfully.
`ftol` condition satisfied.
Current function value: 4.169e-06
Total delta_x: 1.932e-01
Iterations: 14
Function evaluations: 16
Jacobian evaluations: 15
Timer: Solution time = 180 ms
Timer: Avg time per step = 12.0 ms
==============================================================================================================
Start --> End
Total (sum of squares): 6.139e-04 --> 4.169e-06,
Maximum absolute Force error: 3.300e+04 --> 4.335e+03 (N)
Minimum absolute Force error: 7.970e+00 --> 1.156e+01 (N)
Average absolute Force error: 4.844e+03 --> 4.768e+02 (N)
Maximum absolute Force error: 2.654e-03 --> 3.486e-04 (normalized)
Minimum absolute Force error: 6.410e-07 --> 9.296e-07 (normalized)
Average absolute Force error: 3.895e-04 --> 3.835e-05 (normalized)
R boundary error: 0.000e+00 --> 0.000e+00 (m)
Z boundary error: 0.000e+00 --> 0.000e+00 (m)
Fixed Psi error: 0.000e+00 --> 0.000e+00 (Wb)
Fixed pressure profile error: 0.000e+00 --> 0.000e+00 (Pa)
Fixed iota profile error: 0.000e+00 --> 0.000e+00 (dimensionless)
Fixed sheet current error: 0.000e+00 --> 0.000e+00 (~)
==============================================================================================================
Timer: Iteration 4 total = 1.13 sec
================
Step 5
================
Spectral indexing: ansi
Spectral resolution (L,M,N)=(8,8,3)
Node resolution (L,M,N)=(16,16,6)
Boundary ratio = 0.25
Pressure ratio = 1
Perturbation Order = 2
Objective: force
Optimizer: lsq-exact
================
Perturbing equilibrium
Building objective: force
Precomputing transforms
Timer: Precomputing transforms = 35.3 ms
Timer: Objective build = 44.2 ms
Building objective: lcfs R
Building objective: lcfs Z
Building objective: fixed Psi
Building objective: fixed pressure
Building objective: fixed iota
Building objective: fixed sheet current
Building objective: self_consistency R
Building objective: self_consistency Z
Building objective: lambda gauge
Building objective: axis R self consistency
Building objective: axis Z self consistency
Timer: Objective build = 76.3 ms
Perturbing Rb_lmn, Zb_lmn
Factorizing linear constraints
Timer: linear constraint factorize = 802 ms
Computing df
Timer: df computation = 4.94 sec
Factoring df
Timer: df/dx factorization = 85.9 ms
Computing d^2f
Timer: d^2f computation = 3.94 sec
Timer: Objective build = 901 us
||dx||/||x|| = 1.065e-02
Timer: Total perturbation = 12.1 sec
Building objective: self_consistency R
Building objective: self_consistency Z
Building objective: lambda gauge
Building objective: axis R self consistency
Building objective: axis Z self consistency
Timer: Objective build = 23.0 ms
Timer: Linear constraint projection build = 571 ms
Number of parameters: 351
Number of objectives: 2106
Timer: Initializing the optimization = 633 ms
Starting optimization
Using method: lsq-exact
Iteration Total nfev Cost Cost reduction Step norm Optimality
0 1 3.216e-04 5.160e-01
1 2 6.162e-05 2.600e-04 4.750e-02 4.500e-01
2 3 3.637e-05 2.525e-05 5.032e-02 2.124e-01
3 4 5.114e-06 3.126e-05 1.529e-02 1.830e-02
4 6 4.974e-06 1.400e-07 3.913e-03 1.264e-03
5 8 4.972e-06 1.824e-09 1.721e-03 3.344e-04
Optimization terminated successfully.
`ftol` condition satisfied.
Current function value: 4.972e-06
Total delta_x: 2.028e-02
Iterations: 5
Function evaluations: 8
Jacobian evaluations: 6
Timer: Solution time = 3.83 sec
Timer: Avg time per step = 639 ms
==============================================================================================================
Start --> End
Total (sum of squares): 3.216e-04 --> 4.972e-06,
Maximum absolute Force error: 9.822e+04 --> 5.672e+03 (N)
Minimum absolute Force error: 1.369e+00 --> 7.219e-02 (N)
Average absolute Force error: 3.584e+03 --> 5.292e+02 (N)
Maximum absolute Force error: 7.899e-03 --> 4.562e-04 (normalized)
Minimum absolute Force error: 1.101e-07 --> 5.806e-09 (normalized)
Average absolute Force error: 2.882e-04 --> 4.256e-05 (normalized)
R boundary error: 0.000e+00 --> 0.000e+00 (m)
Z boundary error: 0.000e+00 --> 0.000e+00 (m)
Fixed Psi error: 0.000e+00 --> 0.000e+00 (Wb)
Fixed pressure profile error: 0.000e+00 --> 0.000e+00 (Pa)
Fixed iota profile error: 0.000e+00 --> 0.000e+00 (dimensionless)
Fixed sheet current error: 0.000e+00 --> 0.000e+00 (~)
==============================================================================================================
Timer: Iteration 5 total = 18.6 sec
================
Step 6
================
Spectral indexing: ansi
Spectral resolution (L,M,N)=(8,8,3)
Node resolution (L,M,N)=(16,16,6)
Boundary ratio = 0.5
Pressure ratio = 1
Perturbation Order = 2
Objective: force
Optimizer: lsq-exact
================
Perturbing equilibrium
Building objective: force
Precomputing transforms
Timer: Precomputing transforms = 41.0 ms
Timer: Objective build = 51.6 ms
Building objective: lcfs R
Building objective: lcfs Z
Building objective: fixed Psi
Building objective: fixed pressure
Building objective: fixed iota
Building objective: fixed sheet current
Building objective: self_consistency R
Building objective: self_consistency Z
Building objective: lambda gauge
Building objective: axis R self consistency
Building objective: axis Z self consistency
Timer: Objective build = 70.0 ms
Perturbing Rb_lmn, Zb_lmn
Factorizing linear constraints
Timer: linear constraint factorize = 887 ms
Computing df
Timer: df computation = 232 ms
Factoring df
Timer: df/dx factorization = 59.8 ms
Computing d^2f
Timer: d^2f computation = 4.39 ms
Timer: Objective build = 920 us
||dx||/||x|| = 1.034e-02
Timer: Total perturbation = 2.36 sec
Building objective: self_consistency R
Building objective: self_consistency Z
Building objective: lambda gauge
Building objective: axis R self consistency
Building objective: axis Z self consistency
Timer: Objective build = 36.6 ms
Timer: Linear constraint projection build = 967 ms
Number of parameters: 351
Number of objectives: 2106
Timer: Initializing the optimization = 1.05 sec
Starting optimization
Using method: lsq-exact
Iteration Total nfev Cost Cost reduction Step norm Optimality
0 1 2.534e-04 4.496e-01
1 2 5.619e-05 1.972e-04 4.996e-02 4.518e-01
2 4 7.822e-06 4.837e-05 1.555e-02 2.485e-02
3 6 7.402e-06 4.200e-07 7.498e-03 8.701e-03
4 8 7.337e-06 6.505e-08 4.047e-03 2.247e-03
Optimization terminated successfully.
`ftol` condition satisfied.
Current function value: 7.337e-06
Total delta_x: 5.202e-02
Iterations: 4
Function evaluations: 8
Jacobian evaluations: 5
Timer: Solution time = 1.61 sec
Timer: Avg time per step = 322 ms
==============================================================================================================
Start --> End
Total (sum of squares): 2.534e-04 --> 7.337e-06,
Maximum absolute Force error: 6.527e+04 --> 8.269e+03 (N)
Minimum absolute Force error: 9.393e-01 --> 1.625e-01 (N)
Average absolute Force error: 3.331e+03 --> 6.333e+02 (N)
Maximum absolute Force error: 5.250e-03 --> 6.650e-04 (normalized)
Minimum absolute Force error: 7.554e-08 --> 1.307e-08 (normalized)
Average absolute Force error: 2.679e-04 --> 5.093e-05 (normalized)
R boundary error: 0.000e+00 --> 0.000e+00 (m)
Z boundary error: 0.000e+00 --> 0.000e+00 (m)
Fixed Psi error: 0.000e+00 --> 0.000e+00 (Wb)
Fixed pressure profile error: 0.000e+00 --> 0.000e+00 (Pa)
Fixed iota profile error: 0.000e+00 --> 0.000e+00 (dimensionless)
Fixed sheet current error: 0.000e+00 --> 0.000e+00 (~)
==============================================================================================================
Timer: Iteration 6 total = 5.34 sec
================
Step 7
================
Spectral indexing: ansi
Spectral resolution (L,M,N)=(8,8,3)
Node resolution (L,M,N)=(16,16,6)
Boundary ratio = 0.75
Pressure ratio = 1
Perturbation Order = 2
Objective: force
Optimizer: lsq-exact
================
Perturbing equilibrium
Building objective: force
Precomputing transforms
Timer: Precomputing transforms = 41.4 ms
Timer: Objective build = 74.9 ms
Building objective: lcfs R
Building objective: lcfs Z
Building objective: fixed Psi
Building objective: fixed pressure
Building objective: fixed iota
Building objective: fixed sheet current
Building objective: self_consistency R
Building objective: self_consistency Z
Building objective: lambda gauge
Building objective: axis R self consistency
Building objective: axis Z self consistency
Timer: Objective build = 80.1 ms
Perturbing Rb_lmn, Zb_lmn
Factorizing linear constraints
Timer: linear constraint factorize = 712 ms
Computing df
Timer: df computation = 234 ms
Factoring df
Timer: df/dx factorization = 61.0 ms
Computing d^2f
Timer: d^2f computation = 3.79 ms
Timer: Objective build = 811 us
||dx||/||x|| = 1.051e-02
Timer: Total perturbation = 2.11 sec
Building objective: self_consistency R
Building objective: self_consistency Z
Building objective: lambda gauge
Building objective: axis R self consistency
Building objective: axis Z self consistency
Timer: Objective build = 30.5 ms
Timer: Linear constraint projection build = 552 ms
Number of parameters: 351
Number of objectives: 2106
Timer: Initializing the optimization = 635 ms
Starting optimization
Using method: lsq-exact
Iteration Total nfev Cost Cost reduction Step norm Optimality
0 1 3.043e-04 4.975e-01
1 2 5.966e-05 2.447e-04 5.163e-02 5.137e-01
2 4 1.496e-05 4.470e-05 1.463e-02 8.413e-02
3 6 1.408e-05 8.821e-07 6.584e-03 2.403e-02
4 8 1.396e-05 1.164e-07 3.327e-03 6.503e-03
Optimization terminated successfully.
`ftol` condition satisfied.
Current function value: 1.396e-05
Total delta_x: 6.371e-02
Iterations: 4
Function evaluations: 8
Jacobian evaluations: 5
Timer: Solution time = 1.56 sec
Timer: Avg time per step = 313 ms
==============================================================================================================
Start --> End
Total (sum of squares): 3.043e-04 --> 1.396e-05,
Maximum absolute Force error: 7.930e+04 --> 1.600e+04 (N)
Minimum absolute Force error: 7.137e-01 --> 6.245e-01 (N)
Average absolute Force error: 3.772e+03 --> 9.117e+02 (N)
Maximum absolute Force error: 6.378e-03 --> 1.287e-03 (normalized)
Minimum absolute Force error: 5.740e-08 --> 5.022e-08 (normalized)
Average absolute Force error: 3.034e-04 --> 7.332e-05 (normalized)
R boundary error: 0.000e+00 --> 0.000e+00 (m)
Z boundary error: 0.000e+00 --> 0.000e+00 (m)
Fixed Psi error: 0.000e+00 --> 0.000e+00 (Wb)
Fixed pressure profile error: 0.000e+00 --> 0.000e+00 (Pa)
Fixed iota profile error: 0.000e+00 --> 0.000e+00 (dimensionless)
Fixed sheet current error: 0.000e+00 --> 0.000e+00 (~)
==============================================================================================================
Timer: Iteration 7 total = 4.76 sec
================
Step 8
================
Spectral indexing: ansi
Spectral resolution (L,M,N)=(8,8,3)
Node resolution (L,M,N)=(16,16,6)
Boundary ratio = 1.0
Pressure ratio = 1
Perturbation Order = 2
Objective: force
Optimizer: lsq-exact
================
Perturbing equilibrium
Building objective: force
Precomputing transforms
Timer: Precomputing transforms = 42.6 ms
Timer: Objective build = 54.4 ms
Building objective: lcfs R
Building objective: lcfs Z
Building objective: fixed Psi
Building objective: fixed pressure
Building objective: fixed iota
Building objective: fixed sheet current
Building objective: self_consistency R
Building objective: self_consistency Z
Building objective: lambda gauge
Building objective: axis R self consistency
Building objective: axis Z self consistency
Timer: Objective build = 76.5 ms
Perturbing Rb_lmn, Zb_lmn
Factorizing linear constraints
Timer: linear constraint factorize = 676 ms
Computing df
Timer: df computation = 300 ms
Factoring df
Timer: df/dx factorization = 62.3 ms
Computing d^2f
Timer: d^2f computation = 4.50 ms
Timer: Objective build = 873 us
||dx||/||x|| = 1.056e-02
Timer: Total perturbation = 1.76 sec
Building objective: self_consistency R
Building objective: self_consistency Z
Building objective: lambda gauge
Building objective: axis R self consistency
Building objective: axis Z self consistency
Timer: Objective build = 156 ms
Timer: Linear constraint projection build = 686 ms
Number of parameters: 351
Number of objectives: 2106
Timer: Initializing the optimization = 884 ms
Starting optimization
Using method: lsq-exact
Iteration Total nfev Cost Cost reduction Step norm Optimality
0 1 7.358e-04 2.552e+00
1 2 8.750e-05 6.483e-04 7.822e-02 2.978e-01
2 3 6.510e-05 2.240e-05 2.974e-02 2.970e-01
3 4 5.595e-05 9.155e-06 2.846e-02 4.422e-02
4 5 5.330e-05 2.643e-06 2.140e-02 2.283e-02
5 6 5.276e-05 5.418e-07 1.211e-02 1.807e-02
6 7 5.243e-05 3.353e-07 1.377e-02 5.833e-03
Optimization terminated successfully.
`ftol` condition satisfied.
Current function value: 5.243e-05
Total delta_x: 1.335e-01
Iterations: 6
Function evaluations: 7
Jacobian evaluations: 7
Timer: Solution time = 953 ms
Timer: Avg time per step = 136 ms
==============================================================================================================
Start --> End
Total (sum of squares): 7.358e-04 --> 5.243e-05,
Maximum absolute Force error: 1.919e+05 --> 3.270e+04 (N)
Minimum absolute Force error: 1.282e+00 --> 1.546e+00 (N)
Average absolute Force error: 5.554e+03 --> 1.811e+03 (N)
Maximum absolute Force error: 1.544e-02 --> 2.630e-03 (normalized)
Minimum absolute Force error: 1.031e-07 --> 1.243e-07 (normalized)
Average absolute Force error: 4.467e-04 --> 1.456e-04 (normalized)
R boundary error: 0.000e+00 --> 0.000e+00 (m)
Z boundary error: 0.000e+00 --> 0.000e+00 (m)
Fixed Psi error: 0.000e+00 --> 0.000e+00 (Wb)
Fixed pressure profile error: 0.000e+00 --> 0.000e+00 (Pa)
Fixed iota profile error: 0.000e+00 --> 0.000e+00 (dimensionless)
Fixed sheet current error: 0.000e+00 --> 0.000e+00 (~)
==============================================================================================================
Timer: Iteration 8 total = 3.96 sec
====================
Done
Timer: Total time = 1.04 min
====================
solve_continuation_automatic starts with a low resolution vacuum axisymmetric solution, and proceeds to increase the pressure, boundary shaping, and resolution until the final desired configuration is reached. It returns not just the final equilibrium, but each step along the way, as an EquilibriaFamily.
Finally, we can look at the differences between the two methods, and the initial guess
[13]:
plot_comparison(
[eq, eq1, eqf[-1]], labels=["Initial", "eq.solve", "solve_continuation_automatic"]
);
If we compute the normalized force balance error for each case, we see that using the continuation method gives ~20% lower error, indicating a better solution. For more complex equilibria this difference will often be much larger, which is why the continuation method is usually recommended.
[ ]:
f1 = (
eq1.compute("<|F|>_vol")["<|F|>_vol"]
/ eq1.compute("<|grad(|B|^2)|/2mu0>_vol")["<|grad(|B|^2)|/2mu0>_vol"]
)
f2 = (
eqf[-1].compute("<|F|>_vol")["<|F|>_vol"]
/ eqf[-1].compute("<|grad(|B|^2)|/2mu0>_vol")["<|grad(|B|^2)|/2mu0>_vol"]
)
print(f"Force error after eq.solve(): {f1:.4e}")
print(f"Force error after solve_continuation_automatic: {f2:.4e}")
Force error after eq.solve(): 9.3711e-03
Force error after solve_continuation_autmatic: 7.3831e-03