Free boundary equilibrium
In this example we’ll walk through solving a free boundary problem for a Solovev tokamak and a vacuum stellarator.
[1]:
import sys
import os
sys.path.insert(0, os.path.abspath("."))
sys.path.append(os.path.abspath("../../../"))
[2]:
import numpy as np
import matplotlib.pyplot as plt
import desc
from desc.magnetic_fields import (
FourierCurrentPotentialField,
SplineMagneticField,
field_line_integrate,
)
from desc.grid import LinearGrid
from desc.geometry import FourierRZToroidalSurface
from desc.equilibrium import Equilibrium
from desc.objectives import (
BoundaryError,
VacuumBoundaryError,
FixBoundaryR,
FixBoundaryZ,
FixIota,
FixCurrent,
FixPressure,
FixPsi,
ForceBalance,
ObjectiveFunction,
)
from desc.profiles import PowerSeriesProfile
from desc.vmec import VMECIO
DESC version 0.11.0+15.gb84ffffd.dirty,using JAX backend, jax version=0.4.25, jaxlib version=0.4.25, dtype=float64
Using device: CPU, with 3.50 GB available memory
Solovev Tokamak
In the first example, we’ll solve for a free boundary tokamak with Solovev profiles.
We’ll start by loading in an external field, in this case an mgrid file used by VMEC. The external field can also be given directly by a coilset, or from a current potential on a winding surface, or several other representations. See desc.magnetic_fields and desc.coils for more.
[3]:
# need to specify currents in the coil circuits when using mgrid, just like in VMEC
extcur = [
3.884526409876309e06,
-2.935577123737952e05,
-1.734851853677043e04,
6.002137016973160e04,
6.002540940490887e04,
-1.734993103183817e04,
-2.935531536308510e05,
-3.560639108717275e05,
-6.588434719283084e04,
-1.154387774712987e04,
-1.153546510755219e04,
-6.588300858364606e04,
-3.560589388468855e05,
]
ext_field = SplineMagneticField.from_mgrid(
r"../../../tests/inputs/mgrid_solovev.nc", extcur=extcur
)
For our initial guess, we’ll use a circular torus of approximately the right major and minor radius.
[4]:
pres = PowerSeriesProfile([1.25e-1, 0, -1.25e-1])
iota = PowerSeriesProfile([-4.9e-1, 0, 3.0e-1])
surf = FourierRZToroidalSurface(
R_lmn=[4.0, 1.0],
modes_R=[[0, 0], [1, 0]],
Z_lmn=[-1.0],
modes_Z=[[-1, 0]],
NFP=1,
)
eq_init = Equilibrium(M=10, N=0, Psi=1.0, surface=surf, pressure=pres, iota=iota)
eq_init.solve();
Building objective: force
Precomputing transforms
Building objective: lcfs R
Building objective: lcfs Z
Building objective: fixed-Psi
Building objective: fixed-pressure
Building objective: fixed-iota
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
Number of parameters: 75
Number of objectives: 242
Starting optimization
Using method: lsq-exact
Optimization terminated successfully.
`ftol` condition satisfied.
Current function value: 5.892e-12
Total delta_x: 3.553e-01
Iterations: 13
Function evaluations: 20
Jacobian evaluations: 14
Start of solver
Total (sum of squares): 1.110e-01,
Maximum absolute Force error: 8.028e+04 (N)
Minimum absolute Force error: 6.175e+00 (N)
Average absolute Force error: 2.561e+04 (N)
Maximum absolute Force error: 2.193e-02 (normalized)
Minimum absolute Force error: 1.687e-06 (normalized)
Average absolute Force error: 6.995e-03 (normalized)
R boundary error: 0.000e+00 (m)
Z boundary error: 0.000e+00 (m)
Fixed-Psi error: 0.000e+00 (Wb)
Fixed-pressure profile error: 0.000e+00 (Pa)
Fixed-iota profile error: 0.000e+00 (dimensionless)
End of solver
Total (sum of squares): 5.892e-12,
Maximum absolute Force error: 5.591e+00 (N)
Minimum absolute Force error: 5.482e-04 (N)
Average absolute Force error: 1.515e-01 (N)
Maximum absolute Force error: 1.527e-06 (normalized)
Minimum absolute Force error: 1.497e-10 (normalized)
Average absolute Force error: 4.137e-08 (normalized)
R boundary error: 0.000e+00 (m)
Z boundary error: 0.000e+00 (m)
Fixed-Psi error: 0.000e+00 (Wb)
Fixed-pressure profile error: 0.000e+00 (Pa)
Fixed-iota profile error: 0.000e+00 (dimensionless)
[5]:
eq2 = eq_init.copy()
Next we’ll set up our constraints, which in this case simply fix the profiles, flux, and equilibrium constraint
[6]:
constraints = (
ForceBalance(eq=eq2),
FixIota(eq=eq2),
FixPressure(eq=eq2),
FixPsi(eq=eq2),
)
Solving a free boundary equilibrium is just like any other optimization problem. In this case our objective is to minimize boundary error, which is done by the BoundaryError objective.
Specifically, it attempts to minimize the residual in the free boundary MHD boundary conditions:
\(\mathbf{B} \cdot \mathbf{n} = 0\)
\(B^2_{in} + p - B^2_{out} = 0\)
If a sheet current is known to exist on the plasma surface (such as if the pressure at the edge is nonzero), this can be modelled by passing a FourierCurrentPotentialField to the BoundaryError objective, in which case a third equation must be satisfied:
\(\mu_0 \nabla \Phi = \mathbf{n} \times (\mathbf{B}_{out} - \mathbf{B}_{in})\)
Where \(\Phi\) is the current potential on the surface.
The current potential is represented as a FourierCurrentPotentialField which is a subclass of the standard FourierRZToroidalSurface. To include it as part of the free boundary calculation, we set the equilibrium surface to be an instance of FourierCurrentPotentialField by converting the existing surface:
[7]:
eq2.surface = FourierCurrentPotentialField.from_surface(eq2.surface, M_Phi=4)
[8]:
# For a standard free boundary solve, we set field_fixed=True. For single stage optimization, we would set to False
objective = ObjectiveFunction(BoundaryError(eq=eq2, field=ext_field, field_fixed=True))
[9]:
# we know this is a pretty simple shape so we'll only use |m| <= 2
R_modes = (
eq2.surface.R_basis.modes[np.max(np.abs(eq2.surface.R_basis.modes), 1) > 2, :],
)
Z_modes = eq2.surface.Z_basis.modes[np.max(np.abs(eq2.surface.Z_basis.modes), 1) > 2, :]
bdry_constraints = (
FixBoundaryR(eq=eq2, modes=R_modes),
FixBoundaryZ(eq=eq2, modes=Z_modes),
)
eq2, out = eq2.optimize(
objective,
constraints + bdry_constraints,
optimizer="proximal-lsq-exact",
verbose=3,
options={},
)
Building objective: Boundary error
Precomputing transforms
Timer: Precomputing transforms = 523 ms
Timer: Objective build = 921 ms
Building objective: force
Precomputing transforms
Timer: Precomputing transforms = 36.5 ms
Timer: Objective build = 129 ms
Building objective: lcfs R
Building objective: lcfs Z
Building objective: fixed-Psi
Building objective: fixed-pressure
Building objective: fixed-iota
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: Proximal projection build = 3.33 sec
Building objective: fixed-iota
Building objective: fixed-pressure
Building objective: fixed-Psi
Building objective: lcfs R
Building objective: lcfs Z
Timer: Objective build = 442 ms
Timer: Linear constraint projection build = 935 ms
Number of parameters: 11
Number of objectives: 123
Starting optimization
Using method: proximal-lsq-exact
Iteration Total nfev Cost Cost reduction Step norm Optimality
0 1 8.196e+00 4.298e+01
1 2 4.468e-01 7.749e+00 4.189e+05 7.023e+00
2 3 4.007e-02 4.068e-01 2.163e+05 1.015e+00
3 4 2.704e-02 1.303e-02 1.194e+05 7.623e-01
4 5 2.825e-03 2.421e-02 2.450e+05 2.390e-02
5 6 1.416e-03 1.409e-03 5.143e+04 3.464e-01
6 7 3.975e-05 1.376e-03 8.907e+03 6.368e-02
7 9 1.180e-05 2.795e-05 9.331e+02 4.660e-03
8 10 5.046e-06 6.753e-06 1.882e+03 7.956e-03
9 12 4.372e-06 6.741e-07 2.097e+03 6.915e-04
10 14 2.833e-06 1.539e-06 4.014e+02 2.733e-03
11 15 2.832e-06 1.077e-09 2.952e+02 1.520e-04
12 16 2.746e-06 8.516e-08 4.066e+01 1.476e-04
13 17 2.722e-06 2.397e-08 1.094e+02 2.932e-05
Optimization terminated successfully.
`ftol` condition satisfied.
Current function value: 2.722e-06
Total delta_x: 5.171e+02
Iterations: 13
Function evaluations: 17
Jacobian evaluations: 14
Timer: Solution time = 1.31 min
Timer: Avg time per step = 5.63 sec
Start of solver
Total (sum of squares): 8.196e+00,
Maximum absolute Boundary normal field error: 1.420e-02 (T*m^2)
Minimum absolute Boundary normal field error: 3.630e-08 (T*m^2)
Average absolute Boundary normal field error: 8.348e-03 (T*m^2)
Maximum absolute Boundary normal field error: 1.040e-02 (normalized)
Minimum absolute Boundary normal field error: 2.659e-08 (normalized)
Average absolute Boundary normal field error: 6.114e-03 (normalized)
Maximum absolute Boundary magnetic pressure error: 3.201e-01 (T^2*m^2)
Minimum absolute Boundary magnetic pressure error: 1.924e-01 (T^2*m^2)
Average absolute Boundary magnetic pressure error: 2.487e-01 (T^2*m^2)
Maximum absolute Boundary magnetic pressure error: 6.868e-01 (normalized)
Minimum absolute Boundary magnetic pressure error: 4.128e-01 (normalized)
Average absolute Boundary magnetic pressure error: 5.336e-01 (normalized)
Maximum absolute Boundary field jump error: 4.757e-01 (T*m^2)
Minimum absolute Boundary field jump error: 4.749e-01 (T*m^2)
Average absolute Boundary field jump error: 4.753e-01 (T*m^2)
Maximum absolute Boundary field jump error: 3.484e-01 (normalized)
Minimum absolute Boundary field jump error: 3.478e-01 (normalized)
Average absolute Boundary field jump error: 3.481e-01 (normalized)
Maximum absolute Force error: 5.591e+00 (N)
Minimum absolute Force error: 5.482e-04 (N)
Average absolute Force error: 1.515e-01 (N)
Maximum absolute Force error: 1.527e-06 (normalized)
Minimum absolute Force error: 1.497e-10 (normalized)
Average absolute Force error: 4.137e-08 (normalized)
Fixed-iota profile error: 0.000e+00 (dimensionless)
Fixed-pressure profile error: 0.000e+00 (Pa)
Fixed-Psi error: 0.000e+00 (Wb)
R boundary error: 0.000e+00 (m)
Z boundary error: 0.000e+00 (m)
End of solver
Total (sum of squares): 2.722e-06,
Maximum absolute Boundary normal field error: 2.246e-04 (T*m^2)
Minimum absolute Boundary normal field error: 5.507e-08 (T*m^2)
Average absolute Boundary normal field error: 1.036e-04 (T*m^2)
Maximum absolute Boundary normal field error: 1.645e-04 (normalized)
Minimum absolute Boundary normal field error: 4.033e-08 (normalized)
Average absolute Boundary normal field error: 7.585e-05 (normalized)
Maximum absolute Boundary magnetic pressure error: 2.053e-04 (T^2*m^2)
Minimum absolute Boundary magnetic pressure error: 8.168e-06 (T^2*m^2)
Average absolute Boundary magnetic pressure error: 1.018e-04 (T^2*m^2)
Maximum absolute Boundary magnetic pressure error: 4.404e-04 (normalized)
Minimum absolute Boundary magnetic pressure error: 1.752e-05 (normalized)
Average absolute Boundary magnetic pressure error: 2.183e-04 (normalized)
Maximum absolute Boundary field jump error: 5.668e-04 (T*m^2)
Minimum absolute Boundary field jump error: 5.253e-05 (T*m^2)
Average absolute Boundary field jump error: 3.155e-04 (T*m^2)
Maximum absolute Boundary field jump error: 4.151e-04 (normalized)
Minimum absolute Boundary field jump error: 3.847e-05 (normalized)
Average absolute Boundary field jump error: 2.310e-04 (normalized)
Maximum absolute Force error: 2.052e+01 (N)
Minimum absolute Force error: 8.049e-04 (N)
Average absolute Force error: 6.218e-01 (N)
Maximum absolute Force error: 5.607e-06 (normalized)
Minimum absolute Force error: 2.199e-10 (normalized)
Average absolute Force error: 1.698e-07 (normalized)
Fixed-iota profile error: 0.000e+00 (dimensionless)
Fixed-pressure profile error: 0.000e+00 (Pa)
Fixed-Psi error: 0.000e+00 (Wb)
R boundary error: 0.000e+00 (m)
Z boundary error: 0.000e+00 (m)
To check our solution, we can compare to a high resolution free boundary VMEC run, and we see we get extremely good agreement:
[10]:
VMECIO.plot_vmec_comparison(eq2, "../../../tests/inputs/wout_solovev_freeb.nc");
We can plot the normal magnetic field error (the plotting function automatically will add the plasma current’s contribution), and we can see that the normal field is very small for the final solution.
[11]:
desc.plotting.plot_2d(eq2, "B*n", field=ext_field)
[11]:
(<Figure size 432.324x432.324 with 2 Axes>,
<Axes: title={'center': '$\\mathbf{B} \\cdot \\hat{n} ~(\\mathrm{T})$'}, xlabel='$\\zeta$', ylabel='$\\theta$'>)
We can examine what the surface current looks like by plotting it below, and we see it is indeed quite small, only a few tens of Amps, compared to the plasma current which is ~1000x larger. In this case we could probably get equivalent results without including the sheet current term.
[12]:
desc.plotting.plot_2d(eq2.surface, "K");
[13]:
desc.plotting.plot_1d(eq2, "current");
Vacuum Stellarator
We’ll again use an mgrid file for our background field:
[10]:
extcur = [4700.0, 1000.0]
ext_field = SplineMagneticField.from_mgrid(
"../../../tests/inputs/mgrid_test.nc", extcur=extcur
)
For our initial guess, we’ll again use a circular torus of approximately the right major and minor radius.
[11]:
surf = FourierRZToroidalSurface(
R_lmn=[0.70, 0.10],
modes_R=[[0, 0], [1, 0]],
Z_lmn=[-0.10],
modes_Z=[[-1, 0]],
NFP=5,
)
eq_init = Equilibrium(M=6, N=6, Psi=-0.035, surface=surf)
eq_init.solve();
Building objective: force
Precomputing transforms
Compiling objective function and derivatives: ['force']
Number of parameters: 375
Number of objectives: 2450
Starting optimization
Using method: lsq-exact
`gtol` condition satisfied.
Current function value: 1.304e-18
Total delta_x: 2.127e-01
Iterations: 65
Function evaluations: 86
Jacobian evaluations: 66
Start of solver
Total (sum of squares): 1.158e-02,
Maximum absolute Force error: 9.655e+03 (N)
Minimum absolute Force error: 2.406e-13 (N)
Average absolute Force error: 3.130e+03 (N)
Maximum absolute Force error: 7.075e-03 (normalized)
Minimum absolute Force error: 1.763e-19 (normalized)
Average absolute Force error: 2.293e-03 (normalized)
R boundary error: 0.000e+00 (m)
Z boundary error: 0.000e+00 (m)
Fixed-Psi error: 0.000e+00 (Wb)
Fixed-pressure profile error: 0.000e+00 (Pa)
Fixed-current profile error: 0.000e+00 (A)
End of solver
Total (sum of squares): 1.304e-18,
Maximum absolute Force error: 4.408e-04 (N)
Minimum absolute Force error: 1.875e-08 (N)
Average absolute Force error: 2.920e-05 (N)
Maximum absolute Force error: 3.230e-10 (normalized)
Minimum absolute Force error: 1.374e-14 (normalized)
Average absolute Force error: 2.140e-11 (normalized)
R boundary error: 0.000e+00 (m)
Z boundary error: 0.000e+00 (m)
Fixed-Psi error: 0.000e+00 (Wb)
Fixed-pressure profile error: 0.000e+00 (Pa)
Fixed-current profile error: 0.000e+00 (A)
[12]:
eq2 = eq_init.copy()
And again we’ll set up our constraints.
[13]:
constraints = (
ForceBalance(eq=eq2),
FixCurrent(eq=eq2),
FixPressure(eq=eq2),
FixPsi(eq=eq2),
)
The BoundaryError objective we just used uses the virtual casing principle to compute the plasma component of the magnetic field, required to compute the boundary error. If we know we’re solving a vacuum equilibrium, we can skip this calculating since we know the plasma component of the field is zero. This is done in the VacuumBoundaryError objective, which is much more efficient for vacuum equilibria.
[14]:
objective = ObjectiveFunction(
VacuumBoundaryError(eq=eq2, field=ext_field, field_fixed=True)
)
For the optimization, we’ll do a “multigrid” approach where we first optimize the low order modes, and then the higher ones. For a simple problem like this it probably isn’t necessary, but for higher resolution and more complicated shaping this is much more robust.
[15]:
for k in [2, 4]:
# get modes where |m|, |n| > k
R_modes = (
eq2.surface.R_basis.modes[np.max(np.abs(eq2.surface.R_basis.modes), 1) > k, :],
)
Z_modes = eq2.surface.Z_basis.modes[
np.max(np.abs(eq2.surface.Z_basis.modes), 1) > k, :
]
# fix those modes
bdry_constraints = (
FixBoundaryR(eq=eq2, modes=R_modes),
FixBoundaryZ(eq=eq2, modes=Z_modes),
)
# optimize
eq2, out = eq2.optimize(
objective,
constraints + bdry_constraints,
optimizer="proximal-lsq-exact",
verbose=3,
options={},
)
Building objective: Vacuum boundary error
Precomputing transforms
Timer: Precomputing transforms = 237 ms
Timer: Objective build = 1.00 sec
Building objective: force
Precomputing transforms
Timer: Precomputing transforms = 226 ms
Timer: Objective build = 535 ms
Timer: Proximal projection build = 3.86 sec
Timer: Linear constraint projection build = 2.78 sec
Compiling objective function and derivatives: ['Vacuum boundary error']
Timer: Objective compilation time = 5.46 sec
Timer: Jacobian compilation time = 7.49 sec
Timer: Total compilation time = 12.9 sec
Compiling objective function and derivatives: ['force']
Timer: Objective compilation time = 2.65 sec
Timer: Jacobian compilation time = 6.61 sec
Timer: Total compilation time = 9.27 sec
Number of parameters: 25
Number of objectives: 1250
Starting optimization
Using method: proximal-lsq-exact
Iteration Total nfev Cost Cost reduction Step norm Optimality
0 1 2.181e+00 1.085e+02
1 2 3.555e-01 1.825e+00 7.803e-02 2.755e+01
2 3 3.467e-02 3.209e-01 8.711e-02 4.919e+00
3 4 1.655e-03 3.302e-02 3.865e-02 7.491e-01
4 5 9.987e-04 6.558e-04 3.826e-02 9.365e-01
5 6 3.014e-05 9.686e-04 1.011e-02 7.756e-02
6 7 5.440e-06 2.470e-05 2.929e-03 1.115e-02
7 8 5.231e-06 2.090e-07 1.264e-03 2.378e-03
8 9 5.225e-06 6.178e-09 1.345e-04 6.218e-05
Optimization terminated successfully.
`ftol` condition satisfied.
Current function value: 5.225e-06
Total delta_x: 9.407e-02
Iterations: 8
Function evaluations: 9
Jacobian evaluations: 9
Timer: Solution time = 3.57 min
Timer: Avg time per step = 23.8 sec
Start of solver
Total (sum of squares): 2.181e+00,
Maximum absolute Boundary normal field error: 1.121e-02 (T*m^2)
Minimum absolute Boundary normal field error: 2.618e-18 (T*m^2)
Average absolute Boundary normal field error: 4.520e-03 (T*m^2)
Maximum absolute Boundary normal field error: 1.016e-01 (normalized)
Minimum absolute Boundary normal field error: 2.374e-17 (normalized)
Average absolute Boundary normal field error: 4.099e-02 (normalized)
Maximum absolute Boundary magnetic pressure error: 6.633e-02 (T^2*m^2)
Minimum absolute Boundary magnetic pressure error: 3.862e-02 (T^2*m^2)
Average absolute Boundary magnetic pressure error: 5.659e-02 (T^2*m^2)
Maximum absolute Boundary magnetic pressure error: 3.817e-01 (normalized)
Minimum absolute Boundary magnetic pressure error: 2.223e-01 (normalized)
Average absolute Boundary magnetic pressure error: 3.256e-01 (normalized)
Maximum absolute Force error: 4.408e-04 (N)
Minimum absolute Force error: 1.875e-08 (N)
Average absolute Force error: 2.920e-05 (N)
Maximum absolute Force error: 3.230e-10 (normalized)
Minimum absolute Force error: 1.374e-14 (normalized)
Average absolute Force error: 2.140e-11 (normalized)
Fixed-current profile error: 0.000e+00 (A)
Fixed-pressure profile error: 0.000e+00 (Pa)
Fixed-Psi error: 0.000e+00 (Wb)
R boundary error: 0.000e+00 (m)
Z boundary error: 0.000e+00 (m)
End of solver
Total (sum of squares): 5.225e-06,
Maximum absolute Boundary normal field error: 1.401e-04 (T*m^2)
Minimum absolute Boundary normal field error: 4.884e-18 (T*m^2)
Average absolute Boundary normal field error: 4.658e-05 (T*m^2)
Maximum absolute Boundary normal field error: 1.270e-03 (normalized)
Minimum absolute Boundary normal field error: 4.428e-17 (normalized)
Average absolute Boundary normal field error: 4.224e-04 (normalized)
Maximum absolute Boundary magnetic pressure error: 5.364e-05 (T^2*m^2)
Minimum absolute Boundary magnetic pressure error: 4.749e-08 (T^2*m^2)
Average absolute Boundary magnetic pressure error: 1.130e-05 (T^2*m^2)
Maximum absolute Boundary magnetic pressure error: 3.087e-04 (normalized)
Minimum absolute Boundary magnetic pressure error: 2.733e-07 (normalized)
Average absolute Boundary magnetic pressure error: 6.502e-05 (normalized)
Maximum absolute Force error: 1.271e+01 (N)
Minimum absolute Force error: 9.713e-05 (N)
Average absolute Force error: 6.251e-01 (N)
Maximum absolute Force error: 9.316e-06 (normalized)
Minimum absolute Force error: 7.117e-11 (normalized)
Average absolute Force error: 4.580e-07 (normalized)
Fixed-current profile error: 0.000e+00 (A)
Fixed-pressure profile error: 0.000e+00 (Pa)
Fixed-Psi error: 0.000e+00 (Wb)
R boundary error: 0.000e+00 (m)
Z boundary error: 0.000e+00 (m)
Building objective: force
Precomputing transforms
Timer: Precomputing transforms = 227 ms
Timer: Objective build = 544 ms
Timer: Proximal projection build = 2.80 sec
Timer: Linear constraint projection build = 1.80 sec
Compiling objective function and derivatives: ['Vacuum boundary error']
Timer: Objective compilation time = 32.8 ms
Timer: Jacobian compilation time = 242 ms
Timer: Total compilation time = 277 ms
Compiling objective function and derivatives: ['force']
Timer: Objective compilation time = 1.95 sec
Timer: Jacobian compilation time = 5.75 sec
Timer: Total compilation time = 7.70 sec
Number of parameters: 81
Number of objectives: 1250
Starting optimization
Using method: proximal-lsq-exact
Iteration Total nfev Cost Cost reduction Step norm Optimality
0 1 5.225e-06 1.268e-02
1 4 2.795e-07 4.945e-06 2.374e-03 7.767e-03
2 6 1.979e-07 8.159e-08 5.633e-04 9.221e-04
3 7 1.674e-07 3.048e-08 1.183e-03 2.527e-03
4 8 1.375e-07 2.990e-08 1.005e-03 2.034e-03
5 9 1.182e-07 1.932e-08 9.773e-04 1.963e-03
6 11 1.044e-07 1.384e-08 4.943e-04 4.162e-04
7 12 9.596e-08 8.409e-09 9.903e-04 1.652e-03
8 13 8.477e-08 1.119e-08 9.955e-04 2.082e-03
9 14 7.476e-08 1.001e-08 1.022e-03 1.965e-03
10 15 6.883e-08 5.929e-09 1.076e-03 2.280e-03
11 16 6.335e-08 5.472e-09 1.105e-03 2.296e-03
12 17 5.837e-08 4.981e-09 1.093e-03 1.924e-03
13 18 5.560e-08 2.776e-09 1.063e-03 1.652e-03
14 19 5.440e-08 1.199e-09 1.044e-03 1.276e-03
15 20 5.070e-08 3.702e-09 2.622e-04 3.939e-04
16 21 4.985e-08 8.483e-10 5.142e-04 3.355e-04
17 23 4.946e-08 3.884e-10 2.651e-04 4.193e-04
Optimization terminated successfully.
`ftol` condition satisfied.
Current function value: 4.946e-08
Total delta_x: 1.214e-02
Iterations: 17
Function evaluations: 23
Jacobian evaluations: 18
Timer: Solution time = 3.60 min
Timer: Avg time per step = 12.0 sec
Start of solver
Total (sum of squares): 5.225e-06,
Maximum absolute Boundary normal field error: 1.401e-04 (T*m^2)
Minimum absolute Boundary normal field error: 4.884e-18 (T*m^2)
Average absolute Boundary normal field error: 4.658e-05 (T*m^2)
Maximum absolute Boundary normal field error: 1.270e-03 (normalized)
Minimum absolute Boundary normal field error: 4.428e-17 (normalized)
Average absolute Boundary normal field error: 4.224e-04 (normalized)
Maximum absolute Boundary magnetic pressure error: 5.364e-05 (T^2*m^2)
Minimum absolute Boundary magnetic pressure error: 4.749e-08 (T^2*m^2)
Average absolute Boundary magnetic pressure error: 1.130e-05 (T^2*m^2)
Maximum absolute Boundary magnetic pressure error: 3.087e-04 (normalized)
Minimum absolute Boundary magnetic pressure error: 2.733e-07 (normalized)
Average absolute Boundary magnetic pressure error: 6.502e-05 (normalized)
Maximum absolute Force error: 1.271e+01 (N)
Minimum absolute Force error: 9.713e-05 (N)
Average absolute Force error: 6.251e-01 (N)
Maximum absolute Force error: 4.529e-05 (normalized)
Minimum absolute Force error: 3.460e-10 (normalized)
Average absolute Force error: 2.227e-06 (normalized)
Fixed-current profile error: 0.000e+00 (A)
Fixed-pressure profile error: 0.000e+00 (Pa)
Fixed-Psi error: 0.000e+00 (Wb)
R boundary error: 0.000e+00 (m)
Z boundary error: 0.000e+00 (m)
End of solver
Total (sum of squares): 4.946e-08,
Maximum absolute Boundary normal field error: 2.936e-05 (T*m^2)
Minimum absolute Boundary normal field error: 4.537e-18 (T*m^2)
Average absolute Boundary normal field error: 2.703e-06 (T*m^2)
Maximum absolute Boundary normal field error: 2.662e-04 (normalized)
Minimum absolute Boundary normal field error: 4.114e-17 (normalized)
Average absolute Boundary normal field error: 2.451e-05 (normalized)
Maximum absolute Boundary magnetic pressure error: 3.467e-05 (T^2*m^2)
Minimum absolute Boundary magnetic pressure error: 5.448e-09 (T^2*m^2)
Average absolute Boundary magnetic pressure error: 4.034e-06 (T^2*m^2)
Maximum absolute Boundary magnetic pressure error: 1.995e-04 (normalized)
Minimum absolute Boundary magnetic pressure error: 3.135e-08 (normalized)
Average absolute Boundary magnetic pressure error: 2.321e-05 (normalized)
Maximum absolute Force error: 2.695e+01 (N)
Minimum absolute Force error: 1.346e-05 (N)
Average absolute Force error: 1.437e+00 (N)
Maximum absolute Force error: 9.601e-05 (normalized)
Minimum absolute Force error: 4.794e-11 (normalized)
Average absolute Force error: 5.118e-06 (normalized)
Fixed-current profile error: 0.000e+00 (A)
Fixed-pressure profile error: 0.000e+00 (Pa)
Fixed-Psi error: 0.000e+00 (Wb)
R boundary error: 0.000e+00 (m)
Z boundary error: 0.000e+00 (m)
And we see that the boundary has changed quite a lot:
[16]:
desc.plotting.plot_surfaces(eq2);
Because this is a vacuum equilibrium, we can verify the free boundary solution by tracing field lines directly from the external field:
[17]:
nturns = 100 # how many turns
nplanes = 6 # how many places within 1 field period do we want output
# for starting locations we'll pick positions on flux surfaces on the outboard midplane
grid = LinearGrid(rho=np.linspace(0, 1, 9))
r0 = eq2.compute("R", grid=grid)["R"]
z0 = eq2.compute("Z", grid=grid)["Z"]
phis = (
np.linspace(0, 2 * np.pi / eq2.NFP, nplanes, endpoint=False)
+ np.arange(0, nturns)[:, None] * 2 * np.pi / eq2.NFP
).flatten()
[18]:
rs, zs = field_line_integrate(r0, z0, phis, ext_field)
[19]:
zs = zs.reshape((nturns, nplanes, -1))
rs = rs.reshape((nturns, nplanes, -1))
if eq2.Psi < 0: # field lines are traced backwards when toroidal field < 0
rs, zs = rs[:, :: int(np.sign(eq2.Psi))], zs[:, :: int(np.sign(eq2.Psi))]
rs, zs = np.roll(rs, 1, 1), np.roll(zs, 1, 1)
[20]:
fig, ax = desc.plotting.plot_surfaces(eq2)
ax = ax.flatten()
for i in range(nplanes):
ax[i].scatter(rs[:, i, :], zs[:, i, :], s=5)
