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("../../../"))
If you have access to a GPU, uncomment the following two lines before any DESC or JAX related imports. You should see about an order of magnitude speed improvement with only these two lines of code!
[ ]:
# from desc import set_device
# set_device("gpu")
As mentioned in DESC Documentation on performance tips, one can use compilation cache directory to reduce the compilation overhead time. Note: One needs to create jax-caches folder manually.
[3]:
# import jax
# jax.config.update("jax_compilation_cache_dir", "../jax-caches")
# jax.config.update("jax_persistent_cache_min_entry_size_bytes", -1)
# jax.config.update("jax_persistent_cache_min_compile_time_secs", 0)
[4]:
import desc.io
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. Note that DESC makes certain assumptions about the VMEC input file when converting, please see the input file docs. 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].
[5]:
%matplotlib inline
from desc.plotting import plot_comparison, plot_section
[6]:
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
FourierRZToroidalSurfaceProfiles, as
Profileobjects (there are several different options, the most common isPowerSeriesProfile)Total flux
Psias a float
[7]:
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,
[8]:
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:
[9]:
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
[10]:
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>`__
[11]:
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.
[12]:
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
[13]:
eq1, info = eq.solve(verbose=3, copy=True)
Building objective: force
Precomputing transforms
Timer: Precomputing transforms = 735 ms
Timer: Objective build = 1.30 sec
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 = 536 ms
Timer: LinearConstraintProjection build = 5.27 sec
Number of parameters: 351
Number of objectives: 2106
Timer: Initializing the optimization = 7.18 sec
Starting optimization
Using method: lsq-exact
Solver options:
------------------------------------------------------------
Maximum Function Evaluations : 501
Maximum Allowed Total Δx Norm : inf
Scaled Termination : True
Trust Region Method : qr
Initial Trust Radius : 8.819e+01
Maximum Trust Radius : inf
Minimum Trust Radius : 2.220e-16
Trust Radius Increase Ratio : 2.000e+00
Trust Radius Decrease Ratio : 2.500e-01
Trust Radius Increase Threshold : 7.500e-01
Trust Radius Decrease Threshold : 2.500e-01
------------------------------------------------------------
Iteration Total nfev Cost Cost reduction Step norm Optimality
0 1 9.095e-01 6.282e-01
1 2 7.594e-02 8.335e-01 2.539e-01 1.139e-01
2 3 5.712e-03 7.023e-02 3.003e-01 2.947e-02
3 4 4.117e-03 1.596e-03 1.775e-01 3.678e-02
4 5 8.399e-04 3.277e-03 1.178e-01 1.255e-02
5 6 1.741e-04 6.658e-04 7.147e-02 4.889e-03
6 7 9.597e-05 7.817e-05 3.119e-02 8.282e-04
7 8 9.373e-05 2.236e-06 1.633e-02 5.819e-04
8 9 9.260e-05 1.131e-06 1.821e-02 2.174e-04
9 10 9.228e-05 3.286e-07 1.102e-02 2.164e-04
Optimization terminated successfully.
`ftol` condition satisfied. (ftol=1.00e-02)
Current function value: 9.228e-05
Total delta_x: 6.849e-01
Iterations: 9
Function evaluations: 10
Jacobian evaluations: 10
Timer: Solution time = 8.27 sec
Timer: Avg time per step = 827 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 (~)
==============================================================================================================
If you would like to store the objective values before and after equilibrium solve, you can access them through info.
[14]:
# info["Objective values"] is a dictionary of the objective values
for key in info["Objective values"].keys():
if isinstance(info["Objective values"][key], dict):
for subkey in info["Objective values"][key].keys():
print(key, subkey, info["Objective values"][key][subkey])
# if there are multiple instances of the same objective type,
# there will be a list of dictionaries, each storing the values of the objective
elif isinstance(info["Objective values"][key], list):
for i in range(len(info["Objective values"][key])):
for subkey in info["Objective values"][key][i].keys():
print(key, subkey, info["Objective values"][key][i][subkey])
Total (sum of squares) f 9.227504236352532e-05
Total (sum of squares) f0 0.9094859837715115
Force error: f_max 59027.56295828799
Force error: f0_max 1584413.5259700886
Force error: f_min 1.1451026388460366
Force error: f0_min 22.41548413517067
Force error: f_mean 2313.6575222125243
Force error: f0_mean 247770.03845029263
Force error: f_max_norm 0.0047472782879893075
Force error: f0_max_norm 0.12742609645513542
Force error: f_min_norm 9.209461855564835e-08
Force error: f0_min_norm 1.8027601990761712e-06
Force error: f_mean_norm 0.00018607537852789608
Force error: f0_mean_norm 0.01992684883129156
R boundary error: f 0.0
R boundary error: f0 0.0
Z boundary error: f 0.0
Z boundary error: f0 0.0
Fixed Psi error: f 0.0
Fixed Psi error: f0 0.0
Fixed pressure profile error: f 0.0
Fixed pressure profile error: f0 0.0
Fixed iota profile error: f 0.0
Fixed iota profile error: f0 0.0
Fixed sheet current error: f 0.0
Fixed sheet current error: f0 0.0
[15]:
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 = 736 ms
Timer: Objective build = 832 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 = 438 ms
Timer: LinearConstraintProjection build = 4.23 sec
Number of parameters: 27
Number of objectives: 98
Timer: Initializing the optimization = 5.57 sec
Starting optimization
Using method: lsq-exact
Solver options:
------------------------------------------------------------
Maximum Function Evaluations : 501
Maximum Allowed Total Δx Norm : inf
Scaled Termination : True
Trust Region Method : qr
Initial Trust Radius : 3.472e+01
Maximum Trust Radius : inf
Minimum Trust Radius : 2.220e-16
Trust Radius Increase Ratio : 2.000e+00
Trust Radius Decrease Ratio : 2.500e-01
Trust Radius Increase Threshold : 7.500e-01
Trust Radius Decrease Threshold : 2.500e-01
------------------------------------------------------------
Iteration Total nfev Cost Cost reduction Step norm Optimality
0 1 1.513e-02 1.414e-01
1 2 5.311e-03 9.815e-03 2.956e-01 5.522e-02
2 3 4.208e-04 4.890e-03 1.389e-01 1.646e-02
3 4 1.747e-05 4.033e-04 2.667e-02 1.320e-03
4 5 1.014e-06 1.646e-05 1.331e-02 4.090e-04
5 6 1.366e-07 8.775e-07 4.450e-03 3.705e-05
6 7 1.026e-07 3.402e-08 4.940e-03 3.497e-05
7 8 9.623e-08 6.334e-09 4.969e-03 5.524e-06
8 9 9.547e-08 7.556e-10 9.957e-04 1.262e-06
Optimization terminated successfully.
`ftol` condition satisfied. (ftol=1.00e-02)
Current function value: 9.547e-08
Total delta_x: 1.608e-01
Iterations: 8
Function evaluations: 9
Jacobian evaluations: 9
Timer: Solution time = 5.54 sec
Timer: Avg time per step = 616 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 = 12.8 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 = 699 ms
Timer: Objective build = 829 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 = 362 ms
Timer: LinearConstraintProjection build = 4.30 sec
Number of parameters: 48
Number of objectives: 162
Timer: Initializing the optimization = 5.56 sec
Starting optimization
Using method: lsq-exact
Solver options:
------------------------------------------------------------
Maximum Function Evaluations : 501
Maximum Allowed Total Δx Norm : inf
Scaled Termination : True
Trust Region Method : qr
Initial Trust Radius : 5.791e+01
Maximum Trust Radius : inf
Minimum Trust Radius : 2.220e-16
Trust Radius Increase Ratio : 2.000e+00
Trust Radius Decrease Ratio : 2.500e-01
Trust Radius Increase Threshold : 7.500e-01
Trust Radius Decrease Threshold : 2.500e-01
------------------------------------------------------------
Iteration Total nfev Cost Cost reduction Step norm Optimality
0 1 9.241e-08 1.213e-04
1 3 2.638e-08 6.603e-08 1.086e-02 5.110e-05
2 4 1.391e-08 1.247e-08 1.881e-02 5.332e-05
3 5 3.581e-09 1.033e-08 1.674e-02 2.827e-05
4 7 1.881e-10 3.393e-09 2.533e-03 8.765e-06
5 9 1.304e-11 1.751e-10 1.234e-03 2.510e-06
6 11 1.393e-12 1.164e-11 6.311e-04 7.111e-07
7 13 2.868e-13 1.106e-12 3.274e-04 2.163e-07
8 15 2.019e-13 8.494e-14 1.646e-04 5.912e-08
9 17 1.970e-13 4.912e-15 8.174e-05 1.472e-08
10 19 1.964e-13 5.380e-16 4.072e-05 3.658e-09
`gtol` condition satisfied. (gtol=1.00e-08)
Current function value: 1.964e-13
Total delta_x: 4.817e-02
Iterations: 10
Function evaluations: 19
Jacobian evaluations: 11
Timer: Solution time = 6.18 sec
Timer: Avg time per step = 562 ms
==============================================================================================================
Start --> End
Total (sum of squares): 9.241e-08 --> 1.964e-13,
Maximum absolute Force error: 1.359e+03 --> 8.679e-01 (N)
Minimum absolute Force error: 3.405e-01 --> 4.216e-04 (N)
Average absolute Force error: 7.427e+01 --> 1.041e-01 (N)
Maximum absolute Force error: 1.093e-04 --> 6.980e-08 (normalized)
Minimum absolute Force error: 2.738e-08 --> 3.391e-11 (normalized)
Average absolute Force error: 5.973e-06 --> 8.372e-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 = 14.7 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.0
Pressure ratio = 0.5
Perturbation Order = 2
Objective: force
Optimizer: lsq-exact
================
Perturbing equilibrium
Building objective: force
Precomputing transforms
Timer: Precomputing transforms = 38.1 ms
Timer: Objective build = 48.8 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 = 104 ms
Perturbing p_l
Factorizing linear constraints
Timer: linear constraint factorize = 76.2 ms
Computing df
Timer: df computation = 4.33 sec
Factoring df
Timer: df/dx factorization = 201 ms
Computing d^2f
Timer: d^2f computation = 3.17 sec
Timer: Objective build = 1.15 ms
||dx||/||x|| = 3.665e-02
Timer: Total perturbation = 12.4 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 = 12.2 ms
Timer: LinearConstraintProjection build = 85.8 ms
Number of parameters: 48
Number of objectives: 162
Timer: Initializing the optimization = 105 ms
Starting optimization
Using method: lsq-exact
Solver options:
------------------------------------------------------------
Maximum Function Evaluations : 501
Maximum Allowed Total Δx Norm : inf
Scaled Termination : True
Trust Region Method : qr
Initial Trust Radius : 8.377e+01
Maximum Trust Radius : inf
Minimum Trust Radius : 2.220e-16
Trust Radius Increase Ratio : 2.000e+00
Trust Radius Decrease Ratio : 2.500e-01
Trust Radius Increase Threshold : 7.500e-01
Trust Radius Decrease Threshold : 2.500e-01
------------------------------------------------------------
Iteration Total nfev Cost Cost reduction Step norm Optimality
0 1 2.645e-01 2.918e-01
1 2 2.038e-02 2.442e-01 3.697e-01 9.280e-02
2 4 5.427e-04 1.984e-02 1.042e-01 1.164e-02
3 6 1.074e-04 4.353e-04 4.440e-02 2.683e-03
4 8 7.560e-05 3.184e-05 2.060e-02 5.207e-04
5 9 5.202e-05 2.358e-05 3.496e-02 1.388e-03
6 11 3.676e-05 1.526e-05 1.823e-02 3.539e-04
7 12 2.857e-05 8.197e-06 3.493e-02 1.570e-03
8 13 2.252e-05 6.041e-06 3.410e-02 1.879e-03
9 14 1.971e-05 2.814e-06 3.041e-02 2.272e-03
10 15 7.530e-06 1.218e-05 8.479e-03 1.650e-04
11 16 7.002e-06 5.281e-07 1.336e-02 5.307e-04
12 17 6.068e-06 9.341e-07 1.291e-02 4.443e-04
13 18 5.311e-06 7.574e-07 1.307e-02 4.346e-04
14 19 4.675e-06 6.359e-07 1.326e-02 4.378e-04
15 20 4.125e-06 5.494e-07 1.346e-02 4.402e-04
16 21 3.643e-06 4.822e-07 1.367e-02 4.413e-04
17 22 3.217e-06 4.257e-07 1.391e-02 4.419e-04
18 23 2.841e-06 3.764e-07 1.416e-02 4.422e-04
19 24 2.508e-06 3.326e-07 1.442e-02 4.424e-04
20 25 2.215e-06 2.935e-07 1.466e-02 4.424e-04
21 26 1.957e-06 2.581e-07 1.489e-02 4.422e-04
22 27 1.731e-06 2.259e-07 1.511e-02 4.419e-04
23 28 1.534e-06 1.967e-07 1.530e-02 4.415e-04
24 29 1.364e-06 1.705e-07 1.549e-02 4.410e-04
25 30 1.217e-06 1.472e-07 1.568e-02 4.403e-04
26 31 1.090e-06 1.266e-07 1.589e-02 4.393e-04
27 32 9.820e-07 1.080e-07 1.613e-02 4.378e-04
28 33 8.914e-07 9.052e-08 1.640e-02 4.360e-04
29 34 8.185e-07 7.297e-08 1.667e-02 4.340e-04
30 35 5.429e-07 2.755e-07 4.354e-03 2.609e-05
31 36 5.284e-07 1.459e-08 8.614e-03 1.073e-04
32 37 5.022e-07 2.616e-08 8.517e-03 1.074e-04
33 38 4.789e-07 2.330e-08 8.484e-03 1.077e-04
34 39 4.579e-07 2.100e-08 8.400e-03 1.079e-04
35 40 4.390e-07 1.890e-08 8.328e-03 1.082e-04
36 41 4.220e-07 1.699e-08 8.255e-03 1.085e-04
37 42 4.068e-07 1.520e-08 8.189e-03 1.087e-04
38 43 3.933e-07 1.352e-08 8.128e-03 1.090e-04
39 44 3.814e-07 1.192e-08 8.069e-03 1.092e-04
40 45 3.710e-07 1.041e-08 8.012e-03 1.094e-04
41 46 3.620e-07 8.975e-09 7.955e-03 1.097e-04
42 47 3.544e-07 7.610e-09 7.897e-03 1.100e-04
43 48 3.481e-07 6.309e-09 7.836e-03 1.103e-04
44 49 3.275e-07 2.056e-08 2.054e-03 6.784e-06
45 50 3.261e-07 1.432e-09 3.913e-03 2.762e-05
Optimization terminated successfully.
`ftol` condition satisfied. (ftol=1.00e-02)
Current function value: 3.261e-07
Total delta_x: 2.792e-01
Iterations: 45
Function evaluations: 50
Jacobian evaluations: 46
Timer: Solution time = 3.61 sec
Timer: Avg time per step = 78.6 ms
==============================================================================================================
Start --> End
Total (sum of squares): 2.645e-01 --> 3.261e-07,
Maximum absolute Force error: 7.314e+05 --> 1.064e+03 (N)
Minimum absolute Force error: 7.338e+02 --> 1.920e-01 (N)
Average absolute Force error: 9.810e+04 --> 1.337e+02 (N)
Maximum absolute Force error: 5.882e-02 --> 8.559e-05 (normalized)
Minimum absolute Force error: 5.902e-05 --> 1.544e-08 (normalized)
Average absolute Force error: 7.890e-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 --> 4.456e-12 (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 = 17.0 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.0
Pressure ratio = 1.0
Perturbation Order = 2
Objective: force
Optimizer: lsq-exact
================
Perturbing equilibrium
Building objective: force
Precomputing transforms
Timer: Precomputing transforms = 40.7 ms
Timer: Objective build = 58.3 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 = 165 ms
Perturbing p_l
Factorizing linear constraints
Timer: linear constraint factorize = 711 ms
Computing df
Timer: df computation = 14.9 ms
Factoring df
Timer: df/dx factorization = 112 ms
Computing d^2f
Timer: d^2f computation = 5.41 ms
Timer: Objective build = 1.13 ms
||dx||/||x|| = 2.431e-02
Timer: Total perturbation = 2.10 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 = 10.2 ms
Timer: LinearConstraintProjection build = 74.7 ms
Number of parameters: 48
Number of objectives: 162
Timer: Initializing the optimization = 92.2 ms
Starting optimization
Using method: lsq-exact
Solver options:
------------------------------------------------------------
Maximum Function Evaluations : 501
Maximum Allowed Total Δx Norm : inf
Scaled Termination : True
Trust Region Method : qr
Initial Trust Radius : 8.167e+01
Maximum Trust Radius : inf
Minimum Trust Radius : 2.220e-16
Trust Radius Increase Ratio : 2.000e+00
Trust Radius Decrease Ratio : 2.500e-01
Trust Radius Increase Threshold : 7.500e-01
Trust Radius Decrease Threshold : 2.500e-01
------------------------------------------------------------
Iteration Total nfev Cost Cost reduction Step norm Optimality
0 1 6.221e-04 1.709e-02
1 2 6.537e-05 5.567e-04 7.710e-02 5.613e-03
2 3 6.773e-06 5.860e-05 3.741e-02 5.270e-04
3 4 5.575e-06 1.198e-06 1.284e-02 2.339e-04
4 6 4.797e-06 7.776e-07 1.011e-02 6.064e-05
5 7 4.734e-06 6.377e-08 1.597e-02 1.933e-04
6 8 4.652e-06 8.119e-08 1.358e-02 2.926e-04
7 9 4.594e-06 5.888e-08 1.470e-02 2.410e-04
8 10 4.420e-06 1.734e-07 1.542e-02 2.692e-04
9 11 4.344e-06 7.653e-08 1.525e-02 2.772e-04
10 12 4.296e-06 4.769e-08 1.457e-02 3.049e-04
11 13 4.280e-06 1.603e-08 1.419e-02 3.091e-04
12 14 4.170e-06 1.097e-07 4.200e-03 1.679e-05
13 15 4.170e-06 4.761e-10 4.206e-03 2.495e-05
Optimization terminated successfully.
`ftol` condition satisfied. (ftol=1.00e-02)
Current function value: 4.170e-06
Total delta_x: 1.934e-01
Iterations: 13
Function evaluations: 15
Jacobian evaluations: 14
Timer: Solution time = 1.04 sec
Timer: Avg time per step = 74.5 ms
==============================================================================================================
Start --> End
Total (sum of squares): 6.221e-04 --> 4.170e-06,
Maximum absolute Force error: 3.328e+04 --> 4.413e+03 (N)
Minimum absolute Force error: 7.372e+00 --> 9.589e+00 (N)
Average absolute Force error: 4.872e+03 --> 4.761e+02 (N)
Maximum absolute Force error: 2.676e-03 --> 3.549e-04 (normalized)
Minimum absolute Force error: 5.929e-07 --> 7.712e-07 (normalized)
Average absolute Force error: 3.918e-04 --> 3.829e-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 = 3.70 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.0
Perturbation Order = 2
Objective: force
Optimizer: lsq-exact
================
Perturbing equilibrium
Building objective: force
Precomputing transforms
Timer: Precomputing transforms = 61.5 ms
Timer: Objective build = 72.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 = 101 ms
Perturbing Rb_lmn, Zb_lmn
Factorizing linear constraints
Timer: linear constraint factorize = 127 ms
Computing df
Timer: df computation = 6.08 sec
Factoring df
Timer: df/dx factorization = 410 ms
Computing d^2f
Timer: d^2f computation = 3.72 sec
Timer: Objective build = 1.14 ms
||dx||/||x|| = 1.066e-02
Timer: Total perturbation = 15.3 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 = 19.0 ms
Timer: LinearConstraintProjection build = 124 ms
Number of parameters: 351
Number of objectives: 2106
Timer: Initializing the optimization = 153 ms
Starting optimization
Using method: lsq-exact
Solver options:
------------------------------------------------------------
Maximum Function Evaluations : 501
Maximum Allowed Total Δx Norm : inf
Scaled Termination : True
Trust Region Method : qr
Initial Trust Radius : 4.875e+01
Maximum Trust Radius : inf
Minimum Trust Radius : 2.220e-16
Trust Radius Increase Ratio : 2.000e+00
Trust Radius Decrease Ratio : 2.500e-01
Trust Radius Increase Threshold : 7.500e-01
Trust Radius Decrease Threshold : 2.500e-01
------------------------------------------------------------
Iteration Total nfev Cost Cost reduction Step norm Optimality
0 1 3.206e-04 7.784e-03
1 2 6.178e-05 2.589e-04 4.781e-02 5.459e-03
2 3 3.613e-05 2.565e-05 5.057e-02 3.128e-03
3 4 5.111e-06 3.102e-05 1.515e-02 2.105e-04
4 6 4.974e-06 1.366e-07 3.942e-03 1.495e-05
5 8 4.973e-06 1.889e-09 1.737e-03 3.898e-06
Optimization terminated successfully.
`ftol` condition satisfied. (ftol=1.00e-02)
Current function value: 4.973e-06
Total delta_x: 1.973e-02
Iterations: 5
Function evaluations: 8
Jacobian evaluations: 6
Timer: Solution time = 753 ms
Timer: Avg time per step = 125 ms
==============================================================================================================
Start --> End
Total (sum of squares): 3.206e-04 --> 4.973e-06,
Maximum absolute Force error: 9.703e+04 --> 5.674e+03 (N)
Minimum absolute Force error: 5.729e-01 --> 4.809e-01 (N)
Average absolute Force error: 3.571e+03 --> 5.292e+02 (N)
Maximum absolute Force error: 7.804e-03 --> 4.563e-04 (normalized)
Minimum absolute Force error: 4.607e-08 --> 3.868e-08 (normalized)
Average absolute Force error: 2.872e-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 = 16.9 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.0
Perturbation Order = 2
Objective: force
Optimizer: lsq-exact
================
Perturbing equilibrium
Building objective: force
Precomputing transforms
Timer: Precomputing transforms = 65.1 ms
Timer: Objective build = 81.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 = 177 ms
Perturbing Rb_lmn, Zb_lmn
Factorizing linear constraints
Timer: linear constraint factorize = 87.3 ms
Computing df
Timer: df computation = 30.0 ms
Factoring df
Timer: df/dx factorization = 52.8 ms
Computing d^2f
Timer: d^2f computation = 6.03 ms
Timer: Objective build = 1.21 ms
||dx||/||x|| = 1.034e-02
Timer: Total perturbation = 1.07 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 = 15.7 ms
Timer: LinearConstraintProjection build = 83.2 ms
Number of parameters: 351
Number of objectives: 2106
Timer: Initializing the optimization = 106 ms
Starting optimization
Using method: lsq-exact
Solver options:
------------------------------------------------------------
Maximum Function Evaluations : 501
Maximum Allowed Total Δx Norm : inf
Scaled Termination : True
Trust Region Method : qr
Initial Trust Radius : 5.949e+01
Maximum Trust Radius : inf
Minimum Trust Radius : 2.220e-16
Trust Radius Increase Ratio : 2.000e+00
Trust Radius Decrease Ratio : 2.500e-01
Trust Radius Increase Threshold : 7.500e-01
Trust Radius Decrease Threshold : 2.500e-01
------------------------------------------------------------
Iteration Total nfev Cost Cost reduction Step norm Optimality
0 1 2.538e-04 4.595e-03
1 2 5.615e-05 1.977e-04 4.992e-02 4.991e-03
2 4 7.819e-06 4.833e-05 1.554e-02 2.576e-04
3 6 7.402e-06 4.170e-07 7.480e-03 9.095e-05
4 8 7.337e-06 6.506e-08 4.040e-03 2.335e-05
Optimization terminated successfully.
`ftol` condition satisfied. (ftol=1.00e-02)
Current function value: 7.337e-06
Total delta_x: 5.199e-02
Iterations: 4
Function evaluations: 8
Jacobian evaluations: 5
Timer: Solution time = 678 ms
Timer: Avg time per step = 135 ms
==============================================================================================================
Start --> End
Total (sum of squares): 2.538e-04 --> 7.337e-06,
Maximum absolute Force error: 6.469e+04 --> 8.267e+03 (N)
Minimum absolute Force error: 2.999e+00 --> 1.432e-01 (N)
Average absolute Force error: 3.330e+03 --> 6.332e+02 (N)
Maximum absolute Force error: 5.203e-03 --> 6.649e-04 (normalized)
Minimum absolute Force error: 2.412e-07 --> 1.152e-08 (normalized)
Average absolute Force error: 2.678e-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 = 2.40 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.0
Perturbation Order = 2
Objective: force
Optimizer: lsq-exact
================
Perturbing equilibrium
Building objective: force
Precomputing transforms
Timer: Precomputing transforms = 63.9 ms
Timer: Objective build = 78.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 = 154 ms
Perturbing Rb_lmn, Zb_lmn
Factorizing linear constraints
Timer: linear constraint factorize = 81.8 ms
Computing df
Timer: df computation = 16.4 ms
Factoring df
Timer: df/dx factorization = 57.9 ms
Computing d^2f
Timer: d^2f computation = 6.22 ms
Timer: Objective build = 1.18 ms
||dx||/||x|| = 1.051e-02
Timer: Total perturbation = 546 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 = 17.4 ms
Timer: LinearConstraintProjection build = 93.6 ms
Number of parameters: 351
Number of objectives: 2106
Timer: Initializing the optimization = 118 ms
Starting optimization
Using method: lsq-exact
Solver options:
------------------------------------------------------------
Maximum Function Evaluations : 501
Maximum Allowed Total Δx Norm : inf
Scaled Termination : True
Trust Region Method : qr
Initial Trust Radius : 7.738e+01
Maximum Trust Radius : inf
Minimum Trust Radius : 2.220e-16
Trust Radius Increase Ratio : 2.000e+00
Trust Radius Decrease Ratio : 2.500e-01
Trust Radius Increase Threshold : 7.500e-01
Trust Radius Decrease Threshold : 2.500e-01
------------------------------------------------------------
Iteration Total nfev Cost Cost reduction Step norm Optimality
0 1 3.039e-04 5.906e-03
1 2 5.936e-05 2.445e-04 5.153e-02 4.126e-03
2 4 1.496e-05 4.440e-05 1.460e-02 6.771e-04
3 6 1.408e-05 8.796e-07 6.569e-03 1.952e-04
4 8 1.396e-05 1.159e-07 3.320e-03 5.484e-05
Optimization terminated successfully.
`ftol` condition satisfied. (ftol=1.00e-02)
Current function value: 1.396e-05
Total delta_x: 6.362e-02
Iterations: 4
Function evaluations: 8
Jacobian evaluations: 5
Timer: Solution time = 523 ms
Timer: Avg time per step = 104 ms
==============================================================================================================
Start --> End
Total (sum of squares): 3.039e-04 --> 1.396e-05,
Maximum absolute Force error: 7.907e+04 --> 1.600e+04 (N)
Minimum absolute Force error: 6.983e-01 --> 5.589e-01 (N)
Average absolute Force error: 3.768e+03 --> 9.117e+02 (N)
Maximum absolute Force error: 6.359e-03 --> 1.287e-03 (normalized)
Minimum absolute Force error: 5.616e-08 --> 4.495e-08 (normalized)
Average absolute Force error: 3.030e-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 = 1.73 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.0
Perturbation Order = 2
Objective: force
Optimizer: lsq-exact
================
Perturbing equilibrium
Building objective: force
Precomputing transforms
Timer: Precomputing transforms = 64.7 ms
Timer: Objective build = 81.7 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 = 169 ms
Perturbing Rb_lmn, Zb_lmn
Factorizing linear constraints
Timer: linear constraint factorize = 101 ms
Computing df
Timer: df computation = 17.3 ms
Factoring df
Timer: df/dx factorization = 57.9 ms
Computing d^2f
Timer: d^2f computation = 6.74 ms
Timer: Objective build = 1.12 ms
||dx||/||x|| = 1.056e-02
Timer: Total perturbation = 562 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 = 14.9 ms
Timer: LinearConstraintProjection build = 90.4 ms
Number of parameters: 351
Number of objectives: 2106
Timer: Initializing the optimization = 112 ms
Starting optimization
Using method: lsq-exact
Solver options:
------------------------------------------------------------
Maximum Function Evaluations : 501
Maximum Allowed Total Δx Norm : inf
Scaled Termination : True
Trust Region Method : qr
Initial Trust Radius : 1.030e+02
Maximum Trust Radius : inf
Minimum Trust Radius : 2.220e-16
Trust Radius Increase Ratio : 2.000e+00
Trust Radius Decrease Ratio : 2.500e-01
Trust Radius Increase Threshold : 7.500e-01
Trust Radius Decrease Threshold : 2.500e-01
------------------------------------------------------------
Iteration Total nfev Cost Cost reduction Step norm Optimality
0 1 7.340e-04 1.586e-02
1 2 8.730e-05 6.467e-04 7.815e-02 3.569e-03
2 3 6.511e-05 2.219e-05 2.973e-02 1.934e-03
3 4 5.596e-05 9.146e-06 2.849e-02 4.838e-04
4 5 5.331e-05 2.651e-06 2.139e-02 3.374e-04
5 6 5.276e-05 5.437e-07 1.210e-02 1.320e-04
6 7 5.243e-05 3.362e-07 1.378e-02 6.876e-05
Optimization terminated successfully.
`ftol` condition satisfied. (ftol=1.00e-02)
Current function value: 5.243e-05
Total delta_x: 1.335e-01
Iterations: 6
Function evaluations: 7
Jacobian evaluations: 7
Timer: Solution time = 456 ms
Timer: Avg time per step = 65.2 ms
==============================================================================================================
Start --> End
Total (sum of squares): 7.340e-04 --> 5.243e-05,
Maximum absolute Force error: 1.918e+05 --> 3.269e+04 (N)
Minimum absolute Force error: 8.094e-01 --> 1.567e+00 (N)
Average absolute Force error: 5.547e+03 --> 1.811e+03 (N)
Maximum absolute Force error: 1.542e-02 --> 2.629e-03 (normalized)
Minimum absolute Force error: 6.510e-08 --> 1.260e-07 (normalized)
Average absolute Force error: 4.461e-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 = 1.66 sec
====================
Done
Timer: Total time = 1.26 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
[16]:
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.
[17]:
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_automatic: 7.3832e-03