Bootstrap Current Self-Consistency
This tutorial demonstrates how to optimize a quasi-symmetric equilibrium to have a self-consistent bootstrap current profile. This is performed by minimizing the difference between the toroidal currents \(\langle J \cdot B \rangle\) computed from the MHD equilibrium and from the Redl formula. The Redl formula is only valid in the limit of perfect quasi-symmetry, so this procedure will not work for other configurations that are not quasi-symmetric.
There are two methods that can be used, and both will be shown:
Optimize the current profile for self-consistency
Iteratively solve the equilibrium with new current profiles
These methods should be equivalent, although one might be faster than the other depending on the particular problem.
[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!
[2]:
# 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.
[ ]:
# 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)
[3]:
import numpy as np
np.set_printoptions(linewidth=np.inf)
import matplotlib.pyplot as plt
plt.rcParams["font.size"] = 14
from desc.compat import rescale
from desc.equilibrium import EquilibriaFamily
from desc.examples import get
from desc.grid import LinearGrid
from desc.objectives import (
BootstrapRedlConsistency,
FixAtomicNumber,
FixBoundaryR,
FixBoundaryZ,
FixCurrent,
FixElectronDensity,
FixElectronTemperature,
FixIonTemperature,
FixPsi,
ForceBalance,
ObjectiveFunction,
)
from desc.plotting import plot_1d
from desc.profiles import PowerSeriesProfile
An NVIDIA GPU may be present on this machine, but a CUDA-enabled jaxlib is not installed. Falling back to cpu.
DESC version 0.13.0+702.ge6b8a02dc.dirty,using JAX backend, jax version=0.4.33, jaxlib version=0.4.33, dtype=float64
Using device: CPU, with 48.53 GB available memory
As an example, we will reproduce the QA results from Landreman et al. (2022).
We will start with the “precise QA” example equilibrium, scaled to the ARIES-CS reactor size.
[4]:
eq0 = get("precise_QA")
eq0 = rescale(eq0, L=("R0", 10), B=("B0", 5.86))
This equilibrium has the vacuum profiles \(p = 0 ~\text{Pa}\) and \(\frac{2\pi}{\mu_0} I = 0 ~\text{A}\). Calculating the bootstrap current requires knowledge of the temperature and density profiles for each species in the plasma. We replace the pressure profile with the following kinetic profiles corresponding to \(\langle\beta\rangle=2.5\%\):
\(n_e = n_i = 2.38\times10^{20} (1 - \rho^{10}) ~\text{m}^{-3}\)
\(T_e = T_i = 9.45\times10^{3} (1 - \rho^{2}) ~\text{eV}\)
The temperature profiles must be given for both ions and electrons, but only the electron density profile is specified. The ion density profile is given by the effective atomic number \(Z_{eff}\) as \(n_i = n_e / Z_{eff}\). The plasma pressure will then be computed as
\(p = e (n_e T_e + n_i T_i)\).
[5]:
eq0.pressure = None # must remove the pressure profile before setting kinetic profiles
eq0.atomic_number = PowerSeriesProfile([1])
eq0.electron_density = PowerSeriesProfile(params=[1, -1], modes=[0, 10]) * 2.38e20
eq0.electron_temperature = PowerSeriesProfile(params=[1, -1], modes=[0, 2]) * 9.45e3
eq0.ion_temperature = PowerSeriesProfile(params=[1, -1], modes=[0, 2]) * 9.45e3
# the existing current profile is the vacuum case eq0.current = PowerSeriesProfile([0])
We need to re-solve the equilibrium force balance with the new profiles.
[6]:
eq0, _ = eq0.solve(objective="force", optimizer="lsq-exact", verbose=3)
Building objective: force
Precomputing transforms
Timer: Precomputing transforms = 552 ms
Timer: Objective build = 1.39 sec
Building objective: lcfs R
Building objective: lcfs Z
Building objective: fixed Psi
Building objective: fixed current
Building objective: fixed electron density
Building objective: fixed electron temperature
Building objective: fixed ion temperature
Building objective: fixed atomic number
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 = 872 ms
Timer: Linear constraint projection build = 5.06 sec
Number of parameters: 856
Number of objectives: 5346
Timer: Initializing the optimization = 7.37 sec
Starting optimization
Using method: lsq-exact
Iteration Total nfev Cost Cost reduction Step norm Optimality
0 1 8.133e-03 4.417e-02
1 4 5.389e-04 7.594e-03 9.344e-02 1.577e-02
2 5 1.230e-04 4.159e-04 1.292e-01 7.423e-03
3 7 5.177e-05 7.123e-05 9.991e-02 5.108e-03
4 8 1.029e-05 4.148e-05 8.118e-02 1.355e-03
5 10 1.955e-06 8.335e-06 3.413e-02 3.363e-04
6 12 1.040e-06 9.146e-07 1.924e-02 2.155e-04
7 14 8.142e-07 2.259e-07 8.804e-03 7.844e-05
8 15 7.655e-07 4.865e-08 1.597e-02 9.577e-05
9 16 6.887e-07 7.684e-08 1.585e-02 1.683e-04
10 17 6.137e-07 7.495e-08 1.630e-02 1.110e-04
11 18 5.330e-07 8.071e-08 1.683e-02 1.205e-04
12 19 4.294e-07 1.036e-07 1.729e-02 9.057e-05
13 21 3.795e-07 4.991e-08 8.852e-03 2.436e-05
14 22 3.550e-07 2.456e-08 1.724e-02 8.022e-05
15 23 3.160e-07 3.896e-08 1.711e-02 8.099e-05
16 24 2.829e-07 3.305e-08 1.708e-02 8.322e-05
17 25 2.548e-07 2.815e-08 1.704e-02 8.526e-05
18 26 2.313e-07 2.349e-08 1.700e-02 8.822e-05
19 27 2.121e-07 1.921e-08 1.696e-02 9.137e-05
20 28 1.967e-07 1.544e-08 1.692e-02 9.341e-05
21 29 1.841e-07 1.254e-08 1.683e-02 9.284e-05
22 30 1.736e-07 1.047e-08 1.665e-02 8.847e-05
23 31 1.648e-07 8.820e-09 1.633e-02 7.896e-05
24 32 1.579e-07 6.903e-09 1.573e-02 6.173e-05
25 33 1.546e-07 3.286e-09 1.478e-02 4.372e-05
26 34 1.434e-07 1.127e-08 3.668e-03 2.439e-06
27 35 1.433e-07 1.202e-10 6.968e-03 1.317e-05
28 36 1.420e-07 1.261e-09 1.766e-03 7.229e-07
Optimization terminated successfully.
`ftol` condition satisfied.
Current function value: 1.420e-07
Total delta_x: 5.309e-01
Iterations: 28
Function evaluations: 36
Jacobian evaluations: 29
Timer: Solution time = 1.33 min
Timer: Avg time per step = 2.76 sec
==============================================================================================================
Start --> End
Total (sum of squares): 8.133e-03 --> 1.420e-07,
Maximum absolute Force error: 6.965e+07 --> 5.438e+05 (N)
Minimum absolute Force error: 2.392e-01 --> 5.723e+00 (N)
Average absolute Force error: 5.326e+06 --> 3.297e+04 (N)
Maximum absolute Force error: 1.576e-02 --> 1.230e-04 (normalized)
Minimum absolute Force error: 5.413e-11 --> 1.295e-09 (normalized)
Average absolute Force error: 1.205e-03 --> 7.460e-06 (normalized)
R boundary error: 0.000e+00 --> 1.807e-15 (m)
Z boundary error: 0.000e+00 --> 3.294e-16 (m)
Fixed Psi error: 0.000e+00 --> 0.000e+00 (Wb)
Fixed current profile error: 0.000e+00 --> 0.000e+00 (A)
Fixed electron density profile error: 0.000e+00 --> 0.000e+00 (m^-3)
Fixed electron temperature profile error: 0.000e+00 --> 0.000e+00 (eV)
Fixed ion temperature profile error: 0.000e+00 --> 0.000e+00 (eV)
Fixed atomic number profile error: 0.000e+00 --> 0.000e+00 (dimensionless)
Fixed sheet current error: 0.000e+00 --> 0.000e+00 (~)
==============================================================================================================
Now we have our initial equilibrium, which does not have a self-consistent bootstrap current:
[7]:
fig, ax = plot_1d(eq0, "<J*B> Redl", linecolor="b", lw=2, label="Redl")
fig, ax = plot_1d(eq0, "<J*B>", linecolor="r", lw=2, label="MHD", ax=ax)
ax.legend(loc="best")
ax.set_title("Initial Equilibrium");
We need to create a grid on which to evaluate the boootstrap current self-consistency. The bootstrap current is a radial profile, but the grid must have finite poloidal and toroidal resolution to accurately compute flux surface quantities. The Redl formula is undefined where the kinetic profiles vanish, so in our example we do not include points at \(\rho=0\) or \(\rho=1\).
[8]:
grid = LinearGrid(
M=eq0.M_grid,
N=eq0.N_grid,
NFP=eq0.NFP,
sym=eq0.sym,
rho=np.linspace(1 / eq0.L_grid, 1, eq0.L_grid) - 1 / (2 * eq0.L_grid),
)
Our current profile will be represented as a power series of the form:
\(I = c_0 + c_1 \rho + c_2 \rho^2 + \mathcal{O}(\rho^3)\)
Physically, the current should vanish on the magnetic axis so \(c_0 = 0\). And in order for the MHD equilibrium to be analytic, it should scale as \(\mathcal{O}(\rho^2)\) near the magnetic axis so \(c_1 = 0\) also. However, the Redl bootstrap current formula scales as \(\mathcal{O}(\sqrt{\rho})\) near the magnetic axis. This is incorrect, because the drift-kinetic equation from the Redl formula does not account for finite orbit width effects that become important near the axis.
Typically, we use even power series with sym=True for all equilibrium profiles to give the desired analycity conditions. For bootstrap current optimizations, it is recommended to use the full power series with sym=False while also enforcing \(c_0 = c_1 = 0\). This prevents getting good self-consistency near the magnetic axis, but allows for good agreement throughout the rest of the plasma volume and results in high quality equilibria overall.
[9]:
eq0.current = PowerSeriesProfile(np.zeros((eq0.L + 1,)), sym=False)
1. Optimization
In this method, we will optimize the current profile to minimize the self-consistency errors evaluated by the BootstrapRedlConsistency objective. This objective requires the helicity, which for QA is \((M, N) = (1, 0)\).
In this example we will only optimize the current profile, so all other profiles and the plasma boundary are constrained to be fixed. It is recommended to use a very small value for gtol when optimizing the bootstrap current.
[10]:
eq1 = eq0.copy()
[11]:
objective = ObjectiveFunction(
BootstrapRedlConsistency(eq=eq1, grid=grid, helicity=(1, 0)),
)
constraints = (
FixAtomicNumber(eq=eq1),
FixBoundaryR(eq=eq1),
FixBoundaryZ(eq=eq1),
FixCurrent(eq=eq1, indices=[0, 1]), # fix c_0=c_1=0 current profile coefficients
FixElectronDensity(eq=eq1),
FixElectronTemperature(eq=eq1),
FixIonTemperature(eq=eq1),
FixPsi(eq=eq1),
ForceBalance(eq=eq1),
)
eq1, _ = eq1.optimize(
objective=objective,
constraints=constraints,
optimizer="proximal-lsq-exact",
maxiter=10,
verbose=3,
)
Building objective: Bootstrap current self-consistency (Redl)
Precomputing transforms
Timer: Precomputing transforms = 842 ms
Timer: Objective build = 1.53 sec
Building objective: force
Precomputing transforms
Timer: Precomputing transforms = 59.3 ms
Timer: Objective build = 83.9 ms
Timer: Proximal projection build = 6.72 sec
Building objective: fixed atomic number
Building objective: lcfs R
Building objective: lcfs Z
Building objective: fixed current
Building objective: fixed electron density
Building objective: fixed electron temperature
Building objective: fixed ion temperature
Building objective: fixed Psi
Timer: Objective build = 262 ms
Timer: Linear constraint projection build = 2.79 sec
Number of parameters: 7
Number of objectives: 16
Timer: Initializing the optimization = 9.83 sec
Starting optimization
Using method: proximal-lsq-exact
Iteration Total nfev Cost Cost reduction Step norm Optimality
0 1 1.135e-02 1.416e-01
1 2 3.300e-05 1.132e-02 1.753e+07 6.683e-03
2 3 1.268e-06 3.173e-05 3.493e+07 6.175e-04
3 4 8.387e-07 4.289e-07 6.939e+07 4.352e-04
4 5 7.104e-07 1.283e-07 2.068e+07 1.602e-04
5 6 6.900e-07 2.039e-08 1.932e+06 8.287e-05
6 7 6.831e-07 6.911e-09 6.306e+05 3.310e-05
7 8 6.803e-07 2.763e-09 9.909e+05 2.269e-05
Optimization terminated successfully.
`ftol` condition satisfied.
Current function value: 6.803e-07
Total delta_x: 1.305e+08
Iterations: 7
Function evaluations: 8
Jacobian evaluations: 8
Timer: Solution time = 2.49 min
Timer: Avg time per step = 18.6 sec
==============================================================================================================
Start --> End
Total (sum of squares): 1.135e-02 --> 6.803e-07,
Maximum absolute Bootstrap current self-consistency error: 2.861e+06 --> 3.816e+04 (T A m^-2)
Minimum absolute Bootstrap current self-consistency error: 1.687e+05 --> 3.592e+03 (T A m^-2)
Average absolute Bootstrap current self-consistency error: 1.804e+06 --> 1.237e+04 (T A m^-2)
Maximum absolute Bootstrap current self-consistency error: 2.143e-01 --> 2.858e-03 (normalized)
Minimum absolute Bootstrap current self-consistency error: 1.263e-02 --> 2.690e-04 (normalized)
Average absolute Bootstrap current self-consistency error: 1.351e-01 --> 9.267e-04 (normalized)
Fixed atomic number profile error: 0.000e+00 --> 0.000e+00 (dimensionless)
R boundary error: 0.000e+00 --> 1.752e-18 (m)
Z boundary error: 0.000e+00 --> 2.173e-19 (m)
Fixed current profile error: 0.000e+00 --> 0.000e+00 (A)
Fixed electron density profile error: 0.000e+00 --> 0.000e+00 (m^-3)
Fixed electron temperature profile error: 0.000e+00 --> 0.000e+00 (eV)
Fixed ion temperature profile error: 0.000e+00 --> 0.000e+00 (eV)
Fixed Psi error: 0.000e+00 --> 0.000e+00 (Wb)
Maximum absolute Force error: 5.438e+05 --> 1.563e+06 (N)
Minimum absolute Force error: 5.723e+00 --> 4.336e+00 (N)
Average absolute Force error: 3.297e+04 --> 7.910e+04 (N)
Maximum absolute Force error: 1.230e-04 --> 3.537e-04 (normalized)
Minimum absolute Force error: 1.295e-09 --> 9.810e-10 (normalized)
Average absolute Force error: 7.460e-06 --> 1.790e-05 (normalized)
==============================================================================================================
When plotting the bootstrap current profiles, we see the MHD equilibrium now has very good agreement with the Redl formula.
[12]:
fig, ax = plot_1d(eq1, "<J*B> Redl", linecolor="b", lw=2, label="Redl")
fig, ax = plot_1d(eq1, "<J*B>", linecolor="r", lw=2, label="MHD", ax=ax)
ax.legend(loc="best")
ax.set_title("Method 1: Optimization");
2. Iterative Solves
In this method, we iteratively solve the equilibrium with updated guesses for the current profile. The current profile is computed such that the parallel current is consistent with the Redl formula, according to Equation C3 in Landreman & Catto (2012). This is the same approach as STELLOPT VBOOT with SFINCS, and it usually converges in only a few iterations.
[13]:
eq2 = eq0.copy()
fam2 = EquilibriaFamily(eq2)
[14]:
niters = 3
for k in range(niters):
eq2 = eq2.copy()
# compute new guess for the current profile, consistent with Redl formula
data = eq2.compute("current Redl", grid)
current = grid.compress(data["current Redl"])
rho = grid.compress(data["rho"])
# fit the current profile to a power series, with c_0=c_1=0
XX = np.fliplr(np.vander(rho, eq2.L + 1)[:, :-2])
eq2.c_l = np.pad(np.linalg.lstsq(XX, current, rcond=None)[0], (2, 0))
# re-solve the equilibrium
eq2, _ = eq2.solve(objective="force", optimizer="lsq-exact", verbose=3)
fam2.append(eq2)
Building objective: force
Precomputing transforms
Timer: Precomputing transforms = 60.0 ms
Timer: Objective build = 81.2 ms
Building objective: lcfs R
Building objective: lcfs Z
Building objective: fixed Psi
Building objective: fixed current
Building objective: fixed electron density
Building objective: fixed electron temperature
Building objective: fixed ion temperature
Building objective: fixed atomic number
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 = 326 ms
Timer: Linear constraint projection build = 376 ms
Number of parameters: 856
Number of objectives: 5346
Timer: Initializing the optimization = 807 ms
Starting optimization
Using method: lsq-exact
Iteration Total nfev Cost Cost reduction Step norm Optimality
0 1 4.316e-03 7.380e-02
1 2 6.062e-04 3.709e-03 1.787e-01 1.457e-02
2 3 8.345e-05 5.228e-04 1.160e-01 5.285e-03
3 5 4.531e-06 7.892e-05 2.265e-02 6.459e-04
4 7 2.080e-06 2.451e-06 1.281e-02 2.150e-04
5 9 1.825e-06 2.557e-07 6.214e-03 2.159e-05
6 10 1.824e-06 1.059e-09 1.104e-02 8.529e-05
7 11 1.714e-06 1.098e-07 3.190e-03 8.634e-06
8 12 1.690e-06 2.351e-08 6.080e-03 2.522e-05
9 13 1.686e-06 4.460e-09 1.223e-02 7.575e-05
10 14 1.642e-06 4.433e-08 3.476e-03 6.591e-06
11 15 1.632e-06 9.743e-09 6.202e-03 1.806e-05
Optimization terminated successfully.
`ftol` condition satisfied.
Current function value: 1.632e-06
Total delta_x: 2.174e-01
Iterations: 11
Function evaluations: 15
Jacobian evaluations: 12
Timer: Solution time = 32.4 sec
Timer: Avg time per step = 2.70 sec
==============================================================================================================
Start --> End
Total (sum of squares): 4.316e-03 --> 1.632e-06,
Maximum absolute Force error: 6.338e+07 --> 1.411e+06 (N)
Minimum absolute Force error: 1.530e+02 --> 4.526e+01 (N)
Average absolute Force error: 4.041e+06 --> 9.982e+04 (N)
Maximum absolute Force error: 1.434e-02 --> 3.192e-04 (normalized)
Minimum absolute Force error: 3.463e-08 --> 1.024e-08 (normalized)
Average absolute Force error: 9.143e-04 --> 2.259e-05 (normalized)
R boundary error: 0.000e+00 --> 1.752e-18 (m)
Z boundary error: 0.000e+00 --> 2.170e-19 (m)
Fixed Psi error: 0.000e+00 --> 0.000e+00 (Wb)
Fixed current profile error: 0.000e+00 --> 2.982e-08 (A)
Fixed electron density profile error: 0.000e+00 --> 0.000e+00 (m^-3)
Fixed electron temperature profile error: 0.000e+00 --> 0.000e+00 (eV)
Fixed ion temperature profile error: 0.000e+00 --> 0.000e+00 (eV)
Fixed atomic number profile error: 0.000e+00 --> 0.000e+00 (dimensionless)
Fixed sheet current error: 0.000e+00 --> 0.000e+00 (~)
==============================================================================================================
Building objective: force
Precomputing transforms
Timer: Precomputing transforms = 61.9 ms
Timer: Objective build = 83.2 ms
Building objective: lcfs R
Building objective: lcfs Z
Building objective: fixed Psi
Building objective: fixed current
Building objective: fixed electron density
Building objective: fixed electron temperature
Building objective: fixed ion temperature
Building objective: fixed atomic number
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 = 236 ms
Timer: Linear constraint projection build = 346 ms
Number of parameters: 856
Number of objectives: 5346
Timer: Initializing the optimization = 673 ms
Starting optimization
Using method: lsq-exact
Iteration Total nfev Cost Cost reduction Step norm Optimality
0 1 2.297e-04 1.434e-02
1 2 3.563e-06 2.262e-04 5.152e-02 7.821e-04
2 3 1.217e-06 2.345e-06 1.572e-02 9.091e-05
3 5 1.150e-06 6.753e-08 7.105e-03 2.866e-05
4 7 1.142e-06 7.311e-09 3.896e-03 5.647e-06
Optimization terminated successfully.
`ftol` condition satisfied.
Current function value: 1.142e-06
Total delta_x: 5.435e-02
Iterations: 4
Function evaluations: 7
Jacobian evaluations: 5
Timer: Solution time = 10.0 sec
Timer: Avg time per step = 2.01 sec
==============================================================================================================
Start --> End
Total (sum of squares): 2.297e-04 --> 1.142e-06,
Maximum absolute Force error: 1.791e+07 --> 1.301e+06 (N)
Minimum absolute Force error: 6.347e+00 --> 2.539e+00 (N)
Average absolute Force error: 9.465e+05 --> 8.112e+04 (N)
Maximum absolute Force error: 4.052e-03 --> 2.944e-04 (normalized)
Minimum absolute Force error: 1.436e-09 --> 5.745e-10 (normalized)
Average absolute Force error: 2.142e-04 --> 1.835e-05 (normalized)
R boundary error: 0.000e+00 --> 2.711e-20 (m)
Z boundary error: 0.000e+00 --> 6.986e-21 (m)
Fixed Psi error: 0.000e+00 --> 0.000e+00 (Wb)
Fixed current profile error: 0.000e+00 --> 2.107e-08 (A)
Fixed electron density profile error: 0.000e+00 --> 0.000e+00 (m^-3)
Fixed electron temperature profile error: 0.000e+00 --> 0.000e+00 (eV)
Fixed ion temperature profile error: 0.000e+00 --> 0.000e+00 (eV)
Fixed atomic number profile error: 0.000e+00 --> 0.000e+00 (dimensionless)
Fixed sheet current error: 0.000e+00 --> 0.000e+00 (~)
==============================================================================================================
Building objective: force
Precomputing transforms
Timer: Precomputing transforms = 61.2 ms
Timer: Objective build = 81.8 ms
Building objective: lcfs R
Building objective: lcfs Z
Building objective: fixed Psi
Building objective: fixed current
Building objective: fixed electron density
Building objective: fixed electron temperature
Building objective: fixed ion temperature
Building objective: fixed atomic number
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 = 225 ms
Timer: Linear constraint projection build = 388 ms
Number of parameters: 856
Number of objectives: 5346
Timer: Initializing the optimization = 712 ms
Starting optimization
Using method: lsq-exact
Iteration Total nfev Cost Cost reduction Step norm Optimality
0 1 7.983e-06 2.380e-03
1 2 1.727e-06 6.255e-06 2.300e-02 5.136e-04
2 3 1.441e-06 2.859e-07 1.791e-02 3.753e-04
3 5 1.183e-06 2.584e-07 1.051e-02 5.586e-05
4 7 1.139e-06 4.430e-08 4.935e-03 9.875e-06
5 9 1.132e-06 6.331e-09 2.533e-03 2.903e-06
Optimization terminated successfully.
`ftol` condition satisfied.
Current function value: 1.132e-06
Total delta_x: 2.154e-02
Iterations: 5
Function evaluations: 9
Jacobian evaluations: 6
Timer: Solution time = 15.1 sec
Timer: Avg time per step = 2.52 sec
==============================================================================================================
Start --> End
Total (sum of squares): 7.983e-06 --> 1.132e-06,
Maximum absolute Force error: 3.109e+06 --> 1.440e+06 (N)
Minimum absolute Force error: 1.218e+02 --> 7.819e+00 (N)
Average absolute Force error: 1.859e+05 --> 8.172e+04 (N)
Maximum absolute Force error: 7.035e-04 --> 3.258e-04 (normalized)
Minimum absolute Force error: 2.757e-08 --> 1.769e-09 (normalized)
Average absolute Force error: 4.207e-05 --> 1.849e-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 current profile error: 0.000e+00 --> 0.000e+00 (A)
Fixed electron density profile error: 0.000e+00 --> 0.000e+00 (m^-3)
Fixed electron temperature profile error: 0.000e+00 --> 0.000e+00 (eV)
Fixed ion temperature profile error: 0.000e+00 --> 0.000e+00 (eV)
Fixed atomic number profile error: 0.000e+00 --> 0.000e+00 (dimensionless)
Fixed sheet current error: 0.000e+00 --> 0.000e+00 (~)
==============================================================================================================
We can plot the current profile at each iteration to visualize how it changed:
[15]:
fig, ax = plot_1d(fam2[0], "current", linecolor="k", lw=2, label="0")
fig, ax = plot_1d(fam2[1], "current", linecolor="g", lw=2, label="1", ax=ax)
fig, ax = plot_1d(fam2[2], "current", linecolor="b", lw=2, label="2", ax=ax)
fig, ax = plot_1d(fam2[3], "current", linecolor="r", lw=2, label="3", ax=ax)
ax.legend(loc="best");
With this method the MHD equilibrium also has very good agreement with the Redl formula.
[16]:
fig, ax = plot_1d(eq2, "<J*B> Redl", linecolor="b", lw=2, label="Redl")
fig, ax = plot_1d(eq2, "<J*B>", linecolor="r", lw=2, label="MHD", ax=ax)
ax.legend(loc="best")
ax.set_title("Method 2: Iterative Solves");
Comparison
Even though both methods give good self-consistency for the bootstrap current, they do result in slightly different coefficients for the current profile:
[17]:
print(eq1.c_l)
print(eq2.c_l)
[ 0.00000000e+00 0.00000000e+00 -6.30160671e+03 4.08854743e+06 1.56695263e+07 -6.53752503e+07 9.52881073e+07 -5.75991584e+07 1.02645085e+07]
[ 0.00000000e+00 0.00000000e+00 4.45156650e+05 -1.55399717e+06 4.20378770e+07 -1.24735194e+08 1.64818692e+08 -9.79245665e+07 1.92579542e+07]
In this example, the first method of optimization gave better self-consistency but was noticeably slower than the second method of iterative solves.