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]:
# from desc import set_device
# set_device("gpu")
[3]:
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
An NVIDIA GPU may be present on this machine, but a CUDA-enabled jaxlib is not installed. Falling back to cpu.
DESC version 0.12.1+206.g91d82b558.dirty,using JAX backend, jax version=0.4.31, jaxlib version=0.4.31, dtype=float64
Using device: CPU, with 4.19 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.
[4]:
# 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.
[5]:
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: 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
Number of parameters: 75
Number of objectives: 242
Starting optimization
Using method: lsq-exact
Optimization terminated successfully.
`ftol` condition satisfied.
Current function value: 9.452e-13
Total delta_x: 3.499e-01
Iterations: 39
Function evaluations: 47
Jacobian evaluations: 40
==============================================================================================================
Start --> End
Total (sum of squares): 1.110e-01 --> 9.452e-13,
Maximum absolute Force error: 8.028e+04 --> 1.357e+00 (N)
Minimum absolute Force error: 6.175e+00 --> 6.633e-05 (N)
Average absolute Force error: 2.561e+04 --> 6.681e-02 (N)
Maximum absolute Force error: 2.193e-02 --> 3.708e-07 (normalized)
Minimum absolute Force error: 1.687e-06 --> 1.812e-11 (normalized)
Average absolute Force error: 6.995e-03 --> 1.825e-08 (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 (~)
==============================================================================================================
[6]:
eq2 = eq_init.copy()
Next we’ll set up our constraints, which in this case simply fix the profiles, flux, and equilibrium constraint
[7]:
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\)
[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 = 248 ms
Timer: Objective build = 539 ms
Building objective: force
Precomputing transforms
Timer: Precomputing transforms = 41.6 ms
Timer: Objective build = 137 ms
Timer: Proximal projection build = 3.58 sec
Building objective: fixed iota
Building objective: fixed pressure
Building objective: fixed Psi
Building objective: lcfs R
Building objective: lcfs Z
Timer: Objective build = 467 ms
Timer: Linear constraint projection build = 857 ms
Number of parameters: 5
Number of objectives: 82
Timer: Initializing the optimization = 4.94 sec
Starting optimization
Using method: proximal-lsq-exact
Iteration Total nfev Cost Cost reduction Step norm Optimality
0 1 5.805e+00 3.235e+01
1 2 5.547e-01 5.250e+00 4.782e-01 7.006e+00
2 3 1.808e-02 5.366e-01 4.231e-01 9.222e-01
3 4 8.857e-03 9.228e-03 2.740e-01 7.800e-01
4 5 4.005e-05 8.817e-03 1.083e-01 4.927e-02
5 6 1.432e-06 3.862e-05 3.014e-02 1.157e-03
6 8 1.379e-06 5.230e-08 2.673e-03 1.295e-04
7 10 1.376e-06 3.046e-09 1.343e-03 4.880e-05
8 11 1.376e-06 3.546e-10 3.614e-04 1.101e-05
9 12 1.376e-06 1.617e-10 9.251e-05 9.530e-06
10 32 1.376e-06 0.000e+00 0.000e+00 9.530e-06
Warning: A bad approximation caused failure to predict improvement.
Current function value: 1.376e-06
Total delta_x: 5.907e-01
Iterations: 10
Function evaluations: 32
Jacobian evaluations: 10
Timer: Solution time = 1.64 min
Timer: Avg time per step = 8.94 sec
==============================================================================================================
Start --> End
Total (sum of squares): 5.805e+00 --> 1.376e-06,
Maximum absolute Boundary normal field error: 1.426e-02 --> 1.170e-04 (T*m^2)
Minimum absolute Boundary normal field error: 3.630e-08 --> 5.431e-08 (T*m^2)
Average absolute Boundary normal field error: 8.363e-03 --> 6.473e-05 (T*m^2)
Maximum absolute Boundary normal field error: 1.044e-02 --> 8.572e-05 (normalized)
Minimum absolute Boundary normal field error: 2.659e-08 --> 3.978e-08 (normalized)
Average absolute Boundary normal field error: 6.125e-03 --> 4.740e-05 (normalized)
Maximum absolute Boundary magnetic pressure error: 3.201e-01 --> 2.053e-04 (T^2*m^2)
Minimum absolute Boundary magnetic pressure error: 1.924e-01 --> 5.639e-06 (T^2*m^2)
Average absolute Boundary magnetic pressure error: 2.487e-01 --> 1.043e-04 (T^2*m^2)
Maximum absolute Boundary magnetic pressure error: 6.868e-01 --> 4.404e-04 (normalized)
Minimum absolute Boundary magnetic pressure error: 4.128e-01 --> 1.210e-05 (normalized)
Average absolute Boundary magnetic pressure error: 5.336e-01 --> 2.237e-04 (normalized)
Maximum absolute Force error: 1.357e+00 --> 1.999e+01 (N)
Minimum absolute Force error: 6.633e-05 --> 7.638e-03 (N)
Average absolute Force error: 6.681e-02 --> 6.325e-01 (N)
Maximum absolute Force error: 3.708e-07 --> 5.461e-06 (normalized)
Minimum absolute Force error: 1.812e-11 --> 2.087e-09 (normalized)
Average absolute Force error: 1.825e-08 --> 1.728e-07 (normalized)
Fixed iota profile error: 0.000e+00 --> 0.000e+00 (dimensionless)
Fixed pressure profile error: 0.000e+00 --> 0.000e+00 (Pa)
Fixed Psi error: 0.000e+00 --> 0.000e+00 (Wb)
R boundary error: 0.000e+00 --> 0.000e+00 (m)
Z boundary error: 0.000e+00 --> 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 576.162x576.162 with 2 Axes>,
<Axes: title={'center': '$\\mathbf{B} \\cdot \\hat{n} ~(\\mathrm{T})$'}, xlabel='$\\zeta$', ylabel='$\\theta$'>)
If a sheet current is known (or suspected) to exist on the plasma surface (such as if the pressure at the edge is nonzero), this can be modelled by making the equilibrium surface into a FourierCurrentPotentialField.
\(\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:
[12]:
eq3 = eq2.copy()
eq3.surface = FourierCurrentPotentialField.from_surface(eq3.surface, M_Phi=4)
[13]:
constraints = (
ForceBalance(eq=eq3),
FixIota(eq=eq3),
FixPressure(eq=eq3),
FixPsi(eq=eq3),
)
objective = ObjectiveFunction(BoundaryError(eq=eq3, field=ext_field, field_fixed=True))
[14]:
R_modes = eq3.surface.R_basis.modes[np.max(np.abs(eq3.surface.R_basis.modes), 1) > 2, :]
Z_modes = eq3.surface.Z_basis.modes[np.max(np.abs(eq3.surface.Z_basis.modes), 1) > 2, :]
bdry_constraints = (
FixBoundaryR(eq=eq3, modes=R_modes),
FixBoundaryZ(eq=eq3, modes=Z_modes),
)
eq3, out = eq3.optimize(
objective,
constraints + bdry_constraints,
optimizer="proximal-lsq-exact",
verbose=3,
options={},
)
Building objective: Boundary error
Precomputing transforms
Timer: Precomputing transforms = 315 ms
Timer: Objective build = 413 ms
Building objective: force
Precomputing transforms
Timer: Precomputing transforms = 38.5 ms
Timer: Objective build = 134 ms
Timer: Proximal projection build = 2.39 sec
Building objective: fixed iota
Building objective: fixed pressure
Building objective: fixed Psi
Building objective: lcfs R
Building objective: lcfs Z
Timer: Objective build = 419 ms
Timer: Linear constraint projection build = 749 ms
Number of parameters: 11
Number of objectives: 123
Timer: Initializing the optimization = 3.61 sec
Starting optimization
Using method: proximal-lsq-exact
Iteration Total nfev Cost Cost reduction Step norm Optimality
0 1 8.049e-06 4.399e-03
1 7 7.959e-06 8.966e-08 1.069e+02 1.255e-03
2 8 7.796e-06 1.632e-07 5.337e+01 1.049e-03
3 9 7.765e-06 3.109e-08 1.151e+02 3.622e-04
4 10 7.722e-06 4.248e-08 2.062e+01 1.530e-04
Optimization terminated successfully.
`ftol` condition satisfied.
Current function value: 7.722e-06
Total delta_x: 1.187e+02
Iterations: 4
Function evaluations: 10
Jacobian evaluations: 5
Timer: Solution time = 40.4 sec
Timer: Avg time per step = 8.08 sec
==============================================================================================================
Start --> End
Total (sum of squares): 8.049e-06 --> 7.722e-06,
Maximum absolute Boundary normal field error: 1.170e-04 --> 8.649e-05 (T*m^2)
Minimum absolute Boundary normal field error: 5.431e-08 --> 5.433e-08 (T*m^2)
Average absolute Boundary normal field error: 6.473e-05 --> 4.596e-05 (T*m^2)
Maximum absolute Boundary normal field error: 1.090e-04 --> 8.057e-05 (normalized)
Minimum absolute Boundary normal field error: 5.060e-08 --> 5.061e-08 (normalized)
Average absolute Boundary normal field error: 6.030e-05 --> 4.282e-05 (normalized)
Maximum absolute Boundary magnetic pressure error: 2.053e-04 --> 1.989e-04 (T^2*m^2)
Minimum absolute Boundary magnetic pressure error: 5.639e-06 --> 7.496e-06 (T^2*m^2)
Average absolute Boundary magnetic pressure error: 1.043e-04 --> 1.051e-04 (T^2*m^2)
Maximum absolute Boundary magnetic pressure error: 9.116e-04 --> 8.833e-04 (normalized)
Minimum absolute Boundary magnetic pressure error: 2.504e-05 --> 3.329e-05 (normalized)
Average absolute Boundary magnetic pressure error: 4.631e-04 --> 4.666e-04 (normalized)
Maximum absolute Boundary field jump error: 7.271e-04 --> 6.503e-04 (T*m^2)
Minimum absolute Boundary field jump error: 1.556e-05 --> 3.899e-05 (T*m^2)
Average absolute Boundary field jump error: 3.056e-04 --> 2.913e-04 (T*m^2)
Maximum absolute Boundary field jump error: 6.774e-04 --> 6.058e-04 (normalized)
Minimum absolute Boundary field jump error: 1.449e-05 --> 3.632e-05 (normalized)
Average absolute Boundary field jump error: 2.847e-04 --> 2.713e-04 (normalized)
Maximum absolute Force error: 1.999e+01 --> 1.998e+01 (N)
Minimum absolute Force error: 7.638e-03 --> 4.495e-03 (N)
Average absolute Force error: 6.325e-01 --> 6.315e-01 (N)
Maximum absolute Force error: 1.130e-05 --> 1.129e-05 (normalized)
Minimum absolute Force error: 4.319e-09 --> 2.542e-09 (normalized)
Average absolute Force error: 3.576e-07 --> 3.571e-07 (normalized)
Fixed iota profile error: 0.000e+00 --> 0.000e+00 (dimensionless)
Fixed pressure profile error: 0.000e+00 --> 0.000e+00 (Pa)
Fixed Psi error: 0.000e+00 --> 0.000e+00 (Wb)
R boundary error: 0.000e+00 --> 0.000e+00 (m)
Z boundary error: 0.000e+00 --> 0.000e+00 (m)
==============================================================================================================
We can see that including the sheet current makes very little difference in the final result:
[15]:
VMECIO.plot_vmec_comparison(eq3, "../../../tests/inputs/wout_solovev_freeb.nc");
We can see that the normal field error decreased slighlty, though not by much:
[16]:
desc.plotting.plot_2d(eq3, "B*n", field=ext_field)
[16]:
(<Figure size 576.162x576.162 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.
[17]:
desc.plotting.plot_2d(eq3.surface, "K");
[18]:
desc.plotting.plot_1d(eq3, "current");
Vacuum Stellarator
We’ll again use an mgrid file for our background field:
[19]:
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.
[20]:
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
Building objective: lcfs R
Building objective: lcfs Z
Building objective: fixed Psi
Building objective: fixed pressure
Building objective: fixed current
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
Number of parameters: 375
Number of objectives: 2450
Starting optimization
Using method: lsq-exact
`gtol` condition satisfied.
Current function value: 1.514e-18
Total delta_x: 2.119e-01
Iterations: 66
Function evaluations: 72
Jacobian evaluations: 67
==============================================================================================================
Start --> End
Total (sum of squares): 1.158e-02 --> 1.514e-18,
Maximum absolute Force error: 9.655e+03 --> 4.722e-04 (N)
Minimum absolute Force error: 2.406e-13 --> 2.098e-07 (N)
Average absolute Force error: 3.130e+03 --> 3.149e-05 (N)
Maximum absolute Force error: 7.075e-03 --> 3.460e-10 (normalized)
Minimum absolute Force error: 1.763e-19 --> 1.537e-13 (normalized)
Average absolute Force error: 2.293e-03 --> 2.307e-11 (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 current profile error: 0.000e+00 --> 0.000e+00 (A)
Fixed sheet current error: 0.000e+00 --> 0.000e+00 (~)
==============================================================================================================
[21]:
eq2 = eq_init.copy()
And again we’ll set up our constraints.
[22]:
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.
[23]:
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.
[24]:
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 = 179 ms
Timer: Objective build = 549 ms
Building objective: force
Precomputing transforms
Timer: Precomputing transforms = 41.1 ms
Timer: Objective build = 147 ms
Timer: Proximal projection build = 4.49 sec
Building objective: fixed current
Building objective: fixed pressure
Building objective: fixed Psi
Building objective: lcfs R
Building objective: lcfs Z
Timer: Objective build = 645 ms
Timer: Linear constraint projection build = 1.61 sec
Number of parameters: 25
Number of objectives: 1250
Timer: Initializing the optimization = 6.80 sec
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.556e-01 1.825e+00 7.801e-02 2.755e+01
2 3 3.469e-02 3.209e-01 8.714e-02 4.919e+00
3 4 1.656e-03 3.303e-02 3.863e-02 7.457e-01
4 5 9.891e-04 6.670e-04 3.820e-02 9.394e-01
5 6 2.787e-05 9.612e-04 9.914e-03 7.601e-02
6 7 5.401e-06 2.247e-05 2.768e-03 9.852e-03
7 8 5.231e-06 1.694e-07 1.276e-03 2.415e-03
8 9 5.225e-06 6.348e-09 1.319e-04 5.564e-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 = 2.22 min
Timer: Avg time per step = 14.8 sec
==============================================================================================================
Start --> End
Total (sum of squares): 2.181e+00 --> 5.225e-06,
Maximum absolute Boundary normal field error: 1.121e-02 --> 1.400e-04 (T*m^2)
Minimum absolute Boundary normal field error: 2.618e-18 --> 4.884e-18 (T*m^2)
Average absolute Boundary normal field error: 4.520e-03 --> 4.658e-05 (T*m^2)
Maximum absolute Boundary normal field error: 1.016e-01 --> 1.270e-03 (normalized)
Minimum absolute Boundary normal field error: 2.374e-17 --> 4.428e-17 (normalized)
Average absolute Boundary normal field error: 4.099e-02 --> 4.223e-04 (normalized)
Maximum absolute Boundary magnetic pressure error: 6.633e-02 --> 5.365e-05 (T^2*m^2)
Minimum absolute Boundary magnetic pressure error: 3.862e-02 --> 5.023e-08 (T^2*m^2)
Average absolute Boundary magnetic pressure error: 5.659e-02 --> 1.130e-05 (T^2*m^2)
Maximum absolute Boundary magnetic pressure error: 3.817e-01 --> 3.088e-04 (normalized)
Minimum absolute Boundary magnetic pressure error: 2.223e-01 --> 2.890e-07 (normalized)
Average absolute Boundary magnetic pressure error: 3.256e-01 --> 6.502e-05 (normalized)
Maximum absolute Force error: 4.722e-04 --> 1.270e+01 (N)
Minimum absolute Force error: 2.098e-07 --> 2.579e-04 (N)
Average absolute Force error: 3.149e-05 --> 6.249e-01 (N)
Maximum absolute Force error: 3.460e-10 --> 9.309e-06 (normalized)
Minimum absolute Force error: 1.537e-13 --> 1.890e-10 (normalized)
Average absolute Force error: 2.307e-11 --> 4.579e-07 (normalized)
Fixed current profile error: 0.000e+00 --> 0.000e+00 (A)
Fixed pressure profile error: 0.000e+00 --> 0.000e+00 (Pa)
Fixed Psi error: 0.000e+00 --> 0.000e+00 (Wb)
R boundary error: 0.000e+00 --> 0.000e+00 (m)
Z boundary error: 0.000e+00 --> 0.000e+00 (m)
==============================================================================================================
Timer: Objective build = 11.5 ms
Timer: Proximal projection build = 2.42 sec
Building objective: lcfs R
Building objective: lcfs Z
Timer: Objective build = 355 ms
Timer: Linear constraint projection build = 826 ms
Number of parameters: 81
Number of objectives: 1250
Timer: Initializing the optimization = 3.65 sec
Starting optimization
Using method: proximal-lsq-exact
Iteration Total nfev Cost Cost reduction Step norm Optimality
0 1 5.225e-06 1.269e-02
1 4 2.895e-07 4.935e-06 2.339e-03 8.111e-03
2 6 2.144e-07 7.510e-08 5.449e-04 5.586e-04
3 7 1.965e-07 1.798e-08 1.200e-03 2.779e-03
4 8 1.684e-07 2.812e-08 1.103e-03 2.737e-03
5 9 1.464e-07 2.191e-08 1.049e-03 2.575e-03
6 10 1.240e-07 2.243e-08 1.046e-03 2.456e-03
7 12 1.069e-07 1.707e-08 5.082e-04 5.178e-04
8 13 1.014e-07 5.531e-09 9.959e-04 1.956e-03
9 14 9.185e-08 9.564e-09 9.848e-04 1.787e-03
10 15 8.490e-08 6.949e-09 9.919e-04 1.724e-03
11 16 7.945e-08 5.452e-09 1.003e-03 1.676e-03
12 17 7.517e-08 4.279e-09 1.018e-03 1.642e-03
13 18 7.218e-08 2.991e-09 1.037e-03 1.625e-03
14 19 6.997e-08 2.210e-09 1.055e-03 1.609e-03
15 20 6.899e-08 9.848e-10 1.059e-03 1.581e-03
16 21 6.392e-08 5.064e-09 2.679e-04 1.106e-04
17 22 6.369e-08 2.372e-10 5.300e-04 4.139e-04
18 23 6.309e-08 5.974e-10 1.310e-04 2.840e-05
Optimization terminated successfully.
`ftol` condition satisfied.
Current function value: 6.309e-08
Total delta_x: 1.353e-02
Iterations: 18
Function evaluations: 23
Jacobian evaluations: 19
Timer: Solution time = 2.24 min
Timer: Avg time per step = 7.09 sec
==============================================================================================================
Start --> End
Total (sum of squares): 5.225e-06 --> 6.309e-08,
Maximum absolute Boundary normal field error: 1.400e-04 --> 3.109e-05 (T*m^2)
Minimum absolute Boundary normal field error: 4.884e-18 --> 4.705e-18 (T*m^2)
Average absolute Boundary normal field error: 4.658e-05 --> 3.024e-06 (T*m^2)
Maximum absolute Boundary normal field error: 1.270e-03 --> 2.819e-04 (normalized)
Minimum absolute Boundary normal field error: 4.428e-17 --> 4.266e-17 (normalized)
Average absolute Boundary normal field error: 4.223e-04 --> 2.742e-05 (normalized)
Maximum absolute Boundary magnetic pressure error: 5.365e-05 --> 3.712e-05 (T^2*m^2)
Minimum absolute Boundary magnetic pressure error: 5.023e-08 --> 3.760e-08 (T^2*m^2)
Average absolute Boundary magnetic pressure error: 1.130e-05 --> 5.133e-06 (T^2*m^2)
Maximum absolute Boundary magnetic pressure error: 3.088e-04 --> 2.136e-04 (normalized)
Minimum absolute Boundary magnetic pressure error: 2.890e-07 --> 2.164e-07 (normalized)
Average absolute Boundary magnetic pressure error: 6.502e-05 --> 2.954e-05 (normalized)
Maximum absolute Force error: 1.270e+01 --> 1.445e+01 (N)
Minimum absolute Force error: 2.579e-04 --> 2.773e-04 (N)
Average absolute Force error: 6.249e-01 --> 9.437e-01 (N)
Maximum absolute Force error: 9.309e-06 --> 1.059e-05 (normalized)
Minimum absolute Force error: 1.890e-10 --> 2.032e-10 (normalized)
Average absolute Force error: 4.579e-07 --> 6.915e-07 (normalized)
Fixed current profile error: 0.000e+00 --> 0.000e+00 (A)
Fixed pressure profile error: 0.000e+00 --> 0.000e+00 (Pa)
Fixed Psi error: 0.000e+00 --> 0.000e+00 (Wb)
R boundary error: 0.000e+00 --> 0.000e+00 (m)
Z boundary error: 0.000e+00 --> 0.000e+00 (m)
==============================================================================================================
And we see that the boundary has changed quite a lot:
[25]:
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:
[26]:
fig, ax = desc.plotting.plot_surfaces(eq2)
# for starting locations we'll pick positions on flux surfaces on the outboard midplane
grid_trace = LinearGrid(rho=np.linspace(0, 1, 9))
r0 = eq2.compute("R", grid=grid_trace)["R"]
z0 = eq2.compute("Z", grid=grid_trace)["Z"]
fig, ax = desc.plotting.poincare_plot(
ext_field, r0, z0, ntransit=200, NFP=eq2.NFP, ax=ax
);
[ ]: