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("../../../"))
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 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
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.
[5]:
# 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
)
/CODES/DESC/desc/utils.py:572: UserWarning: mgrid does not appear to contain vector potential information. Vector potential will not be computable.
warnings.warn(msg, err)
For our initial guess, we’ll use a circular torus of approximately the right major and minor radius.
[6]:
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
`gtol` condition satisfied. (gtol=1.00e-08)
Current function value: 5.135e-13
Total delta_x: 3.499e-01
Iterations: 38
Function evaluations: 45
Jacobian evaluations: 39
==============================================================================================================
Start --> End
Total (sum of squares): 6.012e-02 --> 5.135e-13,
Maximum absolute Force error: 8.028e+04 --> 1.357e+00 (N)
Minimum absolute Force error: 6.175e+00 --> 7.141e-05 (N)
Average absolute Force error: 2.561e+04 --> 6.696e-02 (N)
Maximum absolute Force error: 1.614e-02 --> 2.728e-07 (normalized)
Minimum absolute Force error: 1.242e-06 --> 1.436e-11 (normalized)
Average absolute Force error: 5.148e-03 --> 1.346e-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 (~)
==============================================================================================================
[7]:
eq2 = eq_init.copy()
Next we’ll set up our constraints, which in this case simply fix the profiles, flux, and equilibrium constraint
[8]:
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\)
[9]:
# 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))
[10]:
# 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 = 19.1 ms
Timer: Objective build = 922 ms
Building objective: force
Precomputing transforms
Timer: Precomputing transforms = 40.2 ms
Timer: Objective build = 55.7 ms
Timer: Objective build = 1.01 ms
Timer: Eq Update LinearConstraintProjection build = 87.6 ms
Timer: Proximal projection build = 3.28 sec
Building objective: fixed iota
Building objective: fixed pressure
Building objective: fixed Psi
Building objective: lcfs R
Building objective: lcfs Z
Timer: Objective build = 231 ms
Timer: LinearConstraintProjection build = 1.51 sec
Number of parameters: 5
Number of objectives: 82
Timer: Initializing the optimization = 5.11 sec
Starting optimization
Using method: proximal-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.140e+00
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.145e+00 2.498e+00
1 2 3.010e-01 2.844e+00 4.780e-01 7.608e-01
2 3 1.039e-02 2.906e-01 4.236e-01 1.078e-01
3 4 4.647e-03 5.743e-03 2.706e-01 7.546e-02
4 5 1.870e-05 4.629e-03 1.076e-01 4.508e-03
5 6 8.599e-07 1.784e-05 1.242e-02 4.459e-05
6 7 8.496e-07 1.028e-08 6.442e-03 6.004e-05
7 9 8.489e-07 7.712e-10 6.524e-04 1.284e-06
8 10 8.475e-07 1.340e-09 1.709e-04 7.322e-07
Optimization terminated successfully.
`ftol` condition satisfied. (ftol=1.00e-02)
Current function value: 8.475e-07
Total delta_x: 5.909e-01
Iterations: 8
Function evaluations: 10
Jacobian evaluations: 9
Timer: Solution time = 35.7 sec
Timer: Avg time per step = 3.97 sec
==============================================================================================================
Start --> End
Total (sum of squares): 3.145e+00 --> 8.475e-07,
Maximum absolute Boundary normal field error: 1.426e-02 --> 1.698e-04 (T*m^2)
Minimum absolute Boundary normal field error: 3.630e-08 --> 5.452e-08 (T*m^2)
Average absolute Boundary normal field error: 8.363e-03 --> 7.633e-05 (T*m^2)
Maximum absolute Boundary normal field error: 8.957e-03 --> 1.067e-04 (normalized)
Minimum absolute Boundary normal field error: 2.281e-08 --> 3.426e-08 (normalized)
Average absolute Boundary normal field error: 5.255e-03 --> 4.796e-05 (normalized)
Maximum absolute Boundary magnetic pressure error: 3.201e-01 --> 2.174e-04 (T^2*m^2)
Minimum absolute Boundary magnetic pressure error: 1.924e-01 --> 8.599e-06 (T^2*m^2)
Average absolute Boundary magnetic pressure error: 2.487e-01 --> 1.093e-04 (T^2*m^2)
Maximum absolute Boundary magnetic pressure error: 5.055e-01 --> 3.433e-04 (normalized)
Minimum absolute Boundary magnetic pressure error: 3.038e-01 --> 1.358e-05 (normalized)
Average absolute Boundary magnetic pressure error: 3.928e-01 --> 1.725e-04 (normalized)
Maximum absolute Force error: 1.357e+00 --> 2.900e+01 (N)
Minimum absolute Force error: 7.141e-05 --> 3.927e-04 (N)
Average absolute Force error: 6.696e-02 --> 8.898e-01 (N)
Maximum absolute Force error: 2.728e-07 --> 5.831e-06 (normalized)
Minimum absolute Force error: 1.436e-11 --> 7.896e-11 (normalized)
Average absolute Force error: 1.346e-08 --> 1.789e-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:
[11]:
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.
[12]:
desc.plotting.plot_2d(eq2, "B*n", field=ext_field);
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:
[13]:
eq3 = eq2.copy()
eq3.surface = FourierCurrentPotentialField.from_surface(eq3.surface, M_Phi=4)
[14]:
constraints = (
ForceBalance(eq=eq3),
FixIota(eq=eq3),
FixPressure(eq=eq3),
FixPsi(eq=eq3),
)
objective = ObjectiveFunction(BoundaryError(eq=eq3, field=ext_field, field_fixed=True))
[15]:
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,
ftol=1e-4,
# make the equilibrium solve subproblem ftol a bit lower, as we don't expect big changes in
# the equilibrium, so we need to resolve the equilibrium more accurately to capture
# the small changes we will see during this optimization
options={"solve_options": {"ftol": 1e-4}},
)
Building objective: Boundary error
Precomputing transforms
Timer: Precomputing transforms = 420 ms
Timer: Objective build = 453 ms
Building objective: force
Precomputing transforms
Timer: Precomputing transforms = 38.9 ms
Timer: Objective build = 50.7 ms
Timer: Objective build = 809 us
Timer: Eq Update LinearConstraintProjection build = 2.89 sec
Timer: Proximal projection build = 17.3 sec
Building objective: fixed iota
Building objective: fixed pressure
Building objective: fixed Psi
Building objective: lcfs R
Building objective: lcfs Z
Timer: Objective build = 193 ms
Timer: LinearConstraintProjection build = 1.58 sec
Number of parameters: 11
Number of objectives: 123
Timer: Initializing the optimization = 19.1 sec
Starting optimization
Using method: proximal-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.770e+00
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 5.212e-06 5.853e-04
1 7 5.036e-06 1.762e-07 6.905e+01 4.411e-04
2 8 5.004e-06 3.188e-08 1.502e+02 6.830e-04
3 9 4.955e-06 4.845e-08 3.654e+01 1.791e-04
4 10 4.939e-06 1.601e-08 2.705e+01 1.490e-04
5 11 4.920e-06 1.965e-08 6.971e+01 2.561e-04
6 13 4.910e-06 9.589e-09 1.764e+01 1.225e-04
7 14 4.904e-06 5.917e-09 1.524e+01 1.123e-04
8 15 4.893e-06 1.085e-08 2.918e+01 9.735e-05
9 16 4.876e-06 1.683e-08 5.601e+01 1.141e-04
10 18 4.876e-06 1.913e-10 2.852e+01 2.519e-04
11 19 4.871e-06 5.291e-09 8.660e+00 1.155e-04
12 20 4.869e-06 2.085e-09 1.655e+01 1.210e-04
13 21 4.867e-06 1.787e-09 1.556e+01 1.107e-04
14 22 4.866e-06 1.516e-09 1.637e+01 1.292e-04
15 23 4.864e-06 1.430e-09 4.122e+00 6.452e-05
Optimization terminated successfully.
`xtol` condition satisfied. (xtol=1.00e-06)
Current function value: 4.864e-06
Total delta_x: 3.784e+02
Iterations: 15
Function evaluations: 23
Jacobian evaluations: 16
Timer: Solution time = 1.53 min
Timer: Avg time per step = 5.75 sec
==============================================================================================================
Start --> End
Total (sum of squares): 5.216e-06 --> 4.864e-06,
Maximum absolute Boundary normal field error: 1.699e-04 --> 1.136e-04 (T*m^2)
Minimum absolute Boundary normal field error: 5.452e-08 --> 5.464e-08 (T*m^2)
Average absolute Boundary normal field error: 7.630e-05 --> 4.477e-05 (T*m^2)
Maximum absolute Boundary normal field error: 1.359e-04 --> 9.093e-05 (normalized)
Minimum absolute Boundary normal field error: 4.363e-08 --> 4.372e-08 (normalized)
Average absolute Boundary normal field error: 6.106e-05 --> 3.582e-05 (normalized)
Maximum absolute Boundary magnetic pressure error: 2.165e-04 --> 2.184e-04 (T^2*m^2)
Minimum absolute Boundary magnetic pressure error: 8.356e-06 --> 9.080e-06 (T^2*m^2)
Average absolute Boundary magnetic pressure error: 1.083e-04 --> 1.088e-04 (T^2*m^2)
Maximum absolute Boundary magnetic pressure error: 7.066e-04 --> 7.128e-04 (normalized)
Minimum absolute Boundary magnetic pressure error: 2.728e-05 --> 2.964e-05 (normalized)
Average absolute Boundary magnetic pressure error: 3.535e-04 --> 3.552e-04 (normalized)
Maximum absolute Boundary field jump error: 7.593e-04 --> 6.001e-04 (T*m^2)
Minimum absolute Boundary field jump error: 2.238e-05 --> 1.107e-05 (T*m^2)
Average absolute Boundary field jump error: 3.134e-04 --> 2.967e-04 (T*m^2)
Maximum absolute Boundary field jump error: 6.076e-04 --> 4.802e-04 (normalized)
Minimum absolute Boundary field jump error: 1.791e-05 --> 8.858e-06 (normalized)
Average absolute Boundary field jump error: 2.507e-04 --> 2.374e-04 (normalized)
Maximum absolute Force error: 2.900e+01 --> 1.733e+01 (N)
Minimum absolute Force error: 3.927e-04 --> 2.058e-03 (N)
Average absolute Force error: 8.898e-01 --> 5.552e-01 (N)
Maximum absolute Force error: 1.205e-05 --> 7.202e-06 (normalized)
Minimum absolute Force error: 1.632e-10 --> 8.552e-10 (normalized)
Average absolute Force error: 3.698e-07 --> 2.307e-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)
==============================================================================================================
If you would like to store the objective values before and after equilibrium solve, you can access them through out.
[16]:
# info["Objective values"] is a dictionary of the objective values
for key in out["Objective values"].keys():
if isinstance(out["Objective values"][key], dict):
for subkey in out["Objective values"][key].keys():
print(key, subkey, out["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(out["Objective values"][key], list):
for i in range(len(out["Objective values"][key])):
for subkey in out["Objective values"][key][i].keys():
print(key, subkey, out["Objective values"][key][i][subkey])
Total (sum of squares) f 4.864107276434656e-06
Total (sum of squares) f0 5.216459016961856e-06
Boundary Error: Boundary normal field error: {'f_max': Array(0.00011363, dtype=float64), 'f_min': Array(5.46424341e-08, dtype=float64), 'f_mean': Array(4.47664356e-05, dtype=float64), 'f_max_norm': Array(9.09288654e-05, dtype=float64), 'f_min_norm': Array(4.37243783e-08, dtype=float64), 'f_mean_norm': Array(3.58216943e-05, dtype=float64), 'f0_max': Array(0.00016989, dtype=float64), 'f0_min': Array(5.45234388e-08, dtype=float64), 'f0_mean': Array(7.63045645e-05, dtype=float64), 'f0_max_norm': Array(0.00013594, dtype=float64), 'f0_min_norm': Array(4.36291594e-08, dtype=float64), 'f0_mean_norm': Array(6.10582179e-05, dtype=float64)}
Boundary Error: Boundary magnetic pressure error: {'f_max': Array(0.00021837, dtype=float64), 'f_min': Array(9.08049005e-06, dtype=float64), 'f_mean': Array(0.00010881, dtype=float64), 'f_max_norm': Array(0.00071282, dtype=float64), 'f_min_norm': Array(2.964156e-05, dtype=float64), 'f_mean_norm': Array(0.00035519, dtype=float64), 'f0_max': Array(0.00021648, dtype=float64), 'f0_min': Array(8.35570923e-06, dtype=float64), 'f0_mean': Array(0.00010829, dtype=float64), 'f0_max_norm': Array(0.00070664, dtype=float64), 'f0_min_norm': Array(2.72756487e-05, dtype=float64), 'f0_mean_norm': Array(0.0003535, dtype=float64)}
Boundary Error: Boundary field jump error: {'f_max': Array(0.00060006, dtype=float64), 'f_min': Array(1.10699802e-05, dtype=float64), 'f_mean': Array(0.00029669, dtype=float64), 'f_max_norm': Array(0.00048016, dtype=float64), 'f_min_norm': Array(8.85809735e-06, dtype=float64), 'f_mean_norm': Array(0.00023741, dtype=float64), 'f0_max': Array(0.00075926, dtype=float64), 'f0_min': Array(2.23795296e-05, dtype=float64), 'f0_mean': Array(0.00031336, dtype=float64), 'f0_max_norm': Array(0.00060756, dtype=float64), 'f0_min_norm': Array(1.79078959e-05, dtype=float64), 'f0_mean_norm': Array(0.00025075, dtype=float64)}
Force error: f_max 17.327988137271074
Force error: f0_max 29.001406963641454
Force error: f_min 0.0020575743372581346
Force error: f0_min 0.00039269495826235875
Force error: f_mean 0.5551535430390399
Force error: f0_mean 0.8898339041013144
Force error: f_max_norm 7.201948803400666e-06
Force error: f0_max_norm 1.2053716018507558e-05
Force error: f_min_norm 8.551797772905136e-10
Force error: f0_min_norm 1.6321392664598554e-10
Force error: f_mean_norm 2.307358109504866e-07
Force error: f0_mean_norm 3.698374081341599e-07
Fixed iota profile error: f 0.0
Fixed iota profile error: f0 0.0
Fixed pressure profile error: f 0.0
Fixed pressure profile error: f0 0.0
Fixed Psi error: f 0.0
Fixed Psi error: f0 0.0
R boundary error: f 0.0
R boundary error: f0 0.0
Z boundary error: f 0.0
Z boundary error: f0 0.0
We can see that including the sheet current makes very little difference in the final result:
[17]:
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:
[18]:
desc.plotting.plot_2d(eq3, "B*n", field=ext_field);
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.
[19]:
desc.plotting.plot_2d(eq3.surface, "K");
[20]:
desc.plotting.plot_1d(eq3, "current");
Vacuum Stellarator
We’ll again use an mgrid file for our background field:
[21]:
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.
[22]:
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=8, N=4, 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: 452
Number of objectives: 2754
Starting optimization
Using method: lsq-exact
`gtol` condition satisfied. (gtol=1.00e-08)
Current function value: 1.048e-13
Total delta_x: 1.743e-01
Iterations: 13
Function evaluations: 20
Jacobian evaluations: 14
==============================================================================================================
Start --> End
Total (sum of squares): 1.940e-02 --> 1.048e-13,
Maximum absolute Force error: 9.735e+03 --> 8.287e-02 (N)
Minimum absolute Force error: 8.771e-14 --> 8.860e-07 (N)
Average absolute Force error: 3.135e+03 --> 6.443e-03 (N)
Maximum absolute Force error: 9.130e-03 --> 7.772e-08 (normalized)
Minimum absolute Force error: 8.227e-20 --> 8.310e-13 (normalized)
Average absolute Force error: 2.940e-03 --> 6.043e-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 current profile error: 0.000e+00 --> 0.000e+00 (A)
Fixed sheet current error: 0.000e+00 --> 0.000e+00 (~)
==============================================================================================================
[23]:
eq2 = eq_init.copy()
And again we’ll set up our constraints.
[24]:
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.
[25]:
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.
[26]:
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 = 742 ms
Timer: Objective build = 1.68 sec
Building objective: force
Precomputing transforms
Timer: Precomputing transforms = 64.7 ms
Timer: Objective build = 82.2 ms
Timer: Objective build = 1.00 ms
Timer: Eq Update LinearConstraintProjection build = 137 ms
Timer: Proximal projection build = 4.21 sec
Building objective: fixed current
Building objective: fixed pressure
Building objective: fixed Psi
Building objective: lcfs R
Building objective: lcfs Z
Timer: Objective build = 357 ms
Timer: LinearConstraintProjection build = 1.95 sec
Number of parameters: 25
Number of objectives: 1122
Timer: Initializing the optimization = 6.63 sec
Starting optimization
Using method: proximal-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.125e+00
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.555e+00 2.586e+00
1 2 5.639e-01 2.991e+00 7.921e-02 7.699e-01
2 3 5.480e-02 5.091e-01 9.051e-02 1.804e-01
3 4 2.191e-03 5.261e-02 4.101e-02 2.399e-02
4 5 1.409e-03 7.817e-04 3.945e-02 1.895e-02
5 6 3.663e-05 1.373e-03 1.054e-02 2.662e-03
6 7 8.220e-06 2.841e-05 4.111e-03 9.039e-04
7 8 6.639e-06 1.582e-06 6.899e-04 1.565e-05
8 9 6.637e-06 1.714e-09 2.233e-04 3.278e-06
Optimization terminated successfully.
`ftol` condition satisfied. (ftol=1.00e-02)
Current function value: 6.637e-06
Total delta_x: 9.407e-02
Iterations: 8
Function evaluations: 9
Jacobian evaluations: 9
Timer: Solution time = 1.41 min
Timer: Avg time per step = 9.40 sec
==============================================================================================================
Start --> End
Total (sum of squares): 3.555e+00 --> 6.637e-06,
Maximum absolute Boundary normal field error: 1.125e-02 --> 1.312e-04 (T*m^2)
Minimum absolute Boundary normal field error: 2.618e-18 --> 4.862e-18 (T*m^2)
Average absolute Boundary normal field error: 4.524e-03 --> 4.620e-05 (T*m^2)
Maximum absolute Boundary normal field error: 1.154e-01 --> 1.346e-03 (normalized)
Minimum absolute Boundary normal field error: 2.685e-17 --> 4.987e-17 (normalized)
Average absolute Boundary normal field error: 4.641e-02 --> 4.739e-04 (normalized)
Maximum absolute Boundary magnetic pressure error: 6.635e-02 --> 3.724e-05 (T^2*m^2)
Minimum absolute Boundary magnetic pressure error: 3.862e-02 --> 1.966e-08 (T^2*m^2)
Average absolute Boundary magnetic pressure error: 5.659e-02 --> 1.086e-05 (T^2*m^2)
Maximum absolute Boundary magnetic pressure error: 4.888e-01 --> 2.743e-04 (normalized)
Minimum absolute Boundary magnetic pressure error: 2.845e-01 --> 1.448e-07 (normalized)
Average absolute Boundary magnetic pressure error: 4.168e-01 --> 7.999e-05 (normalized)
Maximum absolute Force error: 8.287e-02 --> 2.182e+00 (N)
Minimum absolute Force error: 8.860e-07 --> 1.898e-06 (N)
Average absolute Force error: 6.443e-03 --> 8.664e-02 (N)
Maximum absolute Force error: 7.772e-08 --> 2.047e-06 (normalized)
Minimum absolute Force error: 8.310e-13 --> 1.780e-12 (normalized)
Average absolute Force error: 6.043e-09 --> 8.126e-08 (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 = 1.14 ms
Timer: Objective build = 1.03 ms
Timer: Eq Update LinearConstraintProjection build = 130 ms
Timer: Proximal projection build = 1.62 sec
Building objective: lcfs R
Building objective: lcfs Z
Timer: Objective build = 139 ms
Timer: LinearConstraintProjection build = 1.67 sec
Number of parameters: 81
Number of objectives: 1122
Timer: Initializing the optimization = 3.55 sec
Starting optimization
Using method: proximal-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.189e-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.637e-06 9.404e-04
1 4 3.016e-07 6.335e-06 2.330e-03 3.010e-04
2 6 1.862e-07 1.154e-07 5.654e-04 3.287e-05
3 7 1.504e-07 3.583e-08 1.221e-03 1.328e-04
4 8 1.079e-07 4.249e-08 1.053e-03 9.543e-05
5 10 9.110e-08 1.682e-08 4.751e-04 1.908e-05
6 11 8.096e-08 1.014e-08 9.489e-04 7.552e-05
7 12 6.796e-08 1.300e-08 9.628e-04 8.097e-05
8 13 5.647e-08 1.149e-08 9.788e-04 8.564e-05
9 14 4.619e-08 1.028e-08 9.929e-04 8.750e-05
10 15 3.720e-08 8.994e-09 1.003e-03 8.655e-05
11 16 2.966e-08 7.532e-09 1.010e-03 8.369e-05
12 17 2.368e-08 5.979e-09 1.014e-03 7.970e-05
13 18 1.918e-08 4.505e-09 1.015e-03 7.511e-05
14 19 1.594e-08 3.239e-09 1.011e-03 7.055e-05
15 20 1.395e-08 1.993e-09 1.003e-03 6.735e-05
16 21 1.361e-08 3.347e-10 9.961e-04 6.671e-05
17 22 8.356e-09 5.257e-09 2.571e-04 5.401e-06
18 24 8.233e-09 1.233e-10 1.298e-04 1.333e-06
19 25 8.079e-09 1.535e-10 2.585e-04 5.022e-06
20 27 7.970e-09 1.093e-10 1.314e-04 1.373e-06
21 28 7.837e-09 1.324e-10 2.601e-04 5.047e-06
22 30 7.736e-09 1.019e-10 1.320e-04 1.350e-06
23 31 7.620e-09 1.159e-10 2.614e-04 4.954e-06
24 33 7.523e-09 9.677e-11 1.327e-04 1.297e-06
25 34 7.429e-09 9.388e-11 2.636e-04 4.848e-06
26 36 7.343e-09 8.571e-11 1.341e-04 1.269e-06
27 37 7.269e-09 7.401e-11 2.679e-04 4.802e-06
28 38 7.181e-09 8.824e-11 2.724e-04 4.822e-06
29 40 7.107e-09 7.409e-11 1.391e-04 1.235e-06
30 41 7.058e-09 4.938e-11 2.833e-04 4.823e-06
Optimization terminated successfully.
`ftol` condition satisfied. (ftol=1.00e-02)
Current function value: 7.058e-09
Total delta_x: 1.448e-02
Iterations: 30
Function evaluations: 41
Jacobian evaluations: 31
Timer: Solution time = 1.43 min
Timer: Avg time per step = 2.76 sec
==============================================================================================================
Start --> End
Total (sum of squares): 6.637e-06 --> 7.058e-09,
Maximum absolute Boundary normal field error: 1.312e-04 --> 4.753e-06 (T*m^2)
Minimum absolute Boundary normal field error: 4.862e-18 --> 4.436e-18 (T*m^2)
Average absolute Boundary normal field error: 4.620e-05 --> 1.367e-06 (T*m^2)
Maximum absolute Boundary normal field error: 1.346e-03 --> 4.876e-05 (normalized)
Minimum absolute Boundary normal field error: 4.987e-17 --> 4.551e-17 (normalized)
Average absolute Boundary normal field error: 4.739e-04 --> 1.403e-05 (normalized)
Maximum absolute Boundary magnetic pressure error: 3.724e-05 --> 3.241e-06 (T^2*m^2)
Minimum absolute Boundary magnetic pressure error: 1.966e-08 --> 1.617e-08 (T^2*m^2)
Average absolute Boundary magnetic pressure error: 1.086e-05 --> 7.005e-07 (T^2*m^2)
Maximum absolute Boundary magnetic pressure error: 2.743e-04 --> 2.387e-05 (normalized)
Minimum absolute Boundary magnetic pressure error: 1.448e-07 --> 1.191e-07 (normalized)
Average absolute Boundary magnetic pressure error: 7.999e-05 --> 5.160e-06 (normalized)
Maximum absolute Force error: 2.182e+00 --> 5.495e+00 (N)
Minimum absolute Force error: 1.898e-06 --> 7.700e-05 (N)
Average absolute Force error: 8.664e-02 --> 1.766e-01 (N)
Maximum absolute Force error: 2.047e-06 --> 5.153e-06 (normalized)
Minimum absolute Force error: 1.780e-12 --> 7.222e-11 (normalized)
Average absolute Force error: 8.126e-08 --> 1.656e-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:
[27]:
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.
[28]:
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
);
Note that while one can see a continuum of flux surfaces in the vacuum field, the one that the equilibrium free boundary solution will converge is determined by the value of eq2.Psi, the net enclosed toroidal magnetic flux. For example, if we were to change eq2.Psi to a smaller value and run free-boundary, the new equilibrium’s boundary would converge to an interior flux surface of the original, one corresponding to the smaller net enclosed toroidal flux.