Магнитные координатные системы в равновесиях токамака#
Этот учебник сравнивает четыре потоковые координатные системы, используемые в магнитном синтезе: PEST, Boozer, Hamada и Equal-arc.
Что такое координаты потока?#
В токамаке магнитные силовые линии лежат на вложенных тороидальных поверхностях, называемых магнитными поверхностями. Потоковая координатная система \((\psi, \theta, \varphi)\) устроена так, что:
\(\psi\) помечает магнитные поверхности (например, \(\psi_{\rm norm} = 0\) на оси, \(1\) на LCFS)
\(\varphi\) - стандартный тороидальный угол
\(\theta\) - полоидальный угол, выбранный для полезных свойств
Четыре системы отличаются только выбором \(\theta\):
Система |
Определение полоидального угла \(\theta\) |
Ключевое свойство |
|---|---|---|
PEST |
\(\partial(\mathbf{B}\cdot\nabla\theta^*) / \partial\theta^* = 0\): прямые силовые линии |
Самые простые координаты с прямыми силовыми линиями |
Boozer |
PEST + якобиан зависит только от потоковой координаты: \(\sqrt{g} = \sqrt{g}(\psi)\) |
\(\mathbf{J}\cdot\nabla\varphi = {\rm const}(\psi)\); используется в кодах стеллараторов и выходных данных VMEC |
Hamada |
Равноплощадный полоидальный угол |
\(\mathbf{J}\cdot\nabla\theta = {\rm const}(\psi)\); упрощает матрицы устойчивости MHD |
Equal-arc |
Длина дуги вдоль магнитной поверхности равномерна по \(\theta_{\rm ea}\) |
Численно устойчива; проста в построении |
Почему этот выбор важен?#
Прямые силовые линии (PEST, Boozer, Hamada): мода с тороидальным числом \(n\) и полоидальным числом \(m\) удовлетворяет \(q = m/n\) в резонансе. Это упрощает Фурье-анализ MHD-мод.
Boozer: дрейф \(1/B^2\) чисто радиален, что упрощает неоклассическую теорию транспорта и кинетическое уравнение.
Hamada: линии тока также прямые, что требуется некоторым MHD-кодам устойчивости.
Equal-arc: практична для численных сеток; хорошо разрешает границу плазмы.
Настройка: аналитическое равновесие Соловьёва#
Используется равновесие Соловьёва Cerfon & Freidberg (2010), масштабированное до типичного размера токамака (\(R0 \approx 1.86\) м).
[1]:
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.cm as cm
from matplotlib.gridspec import GridSpec
from pyna.toroidal.equilibrium.Solovev import solovev_iter_like
from pyna.toroidal.coords.PEST import build_PEST_mesh
from pyna.toroidal.coords.EqualArc import build_equal_arc_mesh
from pyna.toroidal.coords.Hamada import build_Hamada_mesh
from pyna.toroidal.coords.Boozer import build_Boozer_mesh
# Create equilibrium (scale=0.3 → tokamak-sized, R0≈1.86 m)
eq = solovev_iter_like(scale=0.3)
print(f"Равновесие: R0={eq.R0:.2f} m, a={eq.a:.2f} m, B0={eq.B0:.1f} T")
print(f"κ={eq.kappa:.2f}, δ={eq.delta:.2f}, q0={eq.q0:.2f}")
Rmaxis, Zmaxis = eq.magnetic_axis
print(f"Магнитная ось: R={Rmaxis:.3f} m, Z={Zmaxis:.3f} m")
Равновесие: R0=1.86 м, a=0.60 м, B0=5.3 Тл
κ=1.70, δ=0.33, q0=1.50
Магнитная ось: R=1.946 м, Z=0.000 м
[2]:
# Build the background grid for the equilibrium
nR, nZ = 150, 150
R_grid = np.linspace(0.3 * eq.R0, 1.5 * eq.R0, nR)
Z_grid = np.linspace(-eq.a * eq.kappa * 1.3, eq.a * eq.kappa * 1.3, nZ)
Rg, Zg = np.meshgrid(R_grid, Z_grid, indexing='ij')
BR_grid, BZ_grid = eq.BR_BZ(Rg, Zg)
BPhi_grid = eq.Bphi(Rg)
psi_norm_grid = eq.psi(Rg, Zg)
# Build PEST mesh
# ns=40 radial поверхностей, ntheta=181 полоидальных точек
ns = 40
ntheta = 181
S, TET, R_mesh, Z_mesh, q_iS = build_PEST_mesh(
R_grid, Z_grid, BR_grid, BZ_grid, BPhi_grid, psi_norm_grid,
Rmaxis, Zmaxis, ns=ns, ntheta=ntheta
)
print(f"Сетка PEST построена: {ns} поверхностей, {ntheta} полоидальных точек")
print(f"Диапазон S: [{S[1]:.3f}, {S[-1]:.3f}]")
print(f"Диапазон q: [{q_iS[1]:.2f}, {q_iS[-1]:.2f}]")
Сетка PEST построена: 40 поверхностей, 181 полоидальная точка
Диапазон S: [0.028, 0.972]
Диапазон q: [1.49, 2.39]
Раздел 3: координаты Equal-arc#
Угол равной дуги (Equal-arc) \(\theta_{\rm ea}\) параметризует каждую магнитную поверхность так, что длина дуги вдоль поверхности пропорциональна \(\theta_{\rm ea}\). Это простейшее нетривиальное преобразование координат.
[3]:
# Build equal-arc mesh
_, TET_ea, R_ea, Z_ea = build_equal_arc_mesh(S, TET, R_mesh, Z_mesh, n_theta=ntheta)
print(f"Сетка equal-arc построена: форма {R_ea.shape}")
# Verify arc length uniformity on a mid-radius surface
i_mid = ns // 2
R_s = R_ea[i_mid, :]
Z_s = Z_ea[i_mid, :]
R_loop = R_s[:-1]; Z_loop = Z_s[:-1] # drop endpoint duplicate
R_closed = np.append(R_loop, R_loop[0])
Z_closed = np.append(Z_loop, Z_loop[0])
ds = np.sqrt(np.diff(R_closed)**2 + np.diff(Z_closed)**2)
print(f"Сегменты длины дуги на поверхности {i_mid}: mean={ds.mean()*100:.2f} cm, std={ds.std()*100:.3f} cm")
print(f" Равномерность (СКО/среднее): {ds.std()/ds.mean():.4f}")
Сетка equal-arc построена: форма (40, 181)
Сегменты длины дуги на поверхности 20: среднее=1.27 см, std=0.000 см
Равномерность (СКО/среднее): 0.0001
[4]:
# Plot equal-arc координаты
fig, axes = plt.subplots(1, 2, figsize=(12, 5))
# Left: R-Z view of flux поверхностей with equal-arc mesh lines
ax = axes[0]
# Plot psi_norm contours as background
cs = ax.contour(Rg, Zg, psi_norm_grid, levels=np.linspace(0, 1, 11),
colors='lightgray', linewidths=0.5)
# Plot every 5th surface in equal-arc
stride = max(1, ns // 10)
colors = cm.plasma(np.linspace(0.2, 0.9, ns // stride + 1))
for k, i in enumerate(range(1, ns, stride)):
ax.plot(R_ea[i, :], Z_ea[i, :], color=colors[k], lw=1.0)
# Plot a few poloidal lines (constant θ_ea)
for j in range(0, ntheta-1, ntheta // 8):
ax.plot(R_ea[1:, j], Z_ea[1:, j], 'b-', lw=0.5, alpha=0.5)
ax.plot(Rmaxis, Zmaxis, 'r+', ms=10, mew=2, label='Axis')
ax.set_xlabel('R (м)')
ax.set_ylabel('Z (м)')
ax.set_title('Сетка equal-arc (вид R-Z)')
ax.set_aspect('equal')
ax.legend()
# Right: arc-length increments as function of theta_ea
ax = axes[1]
for k, i in enumerate(range(2, ns-1, stride)):
R_s = R_ea[i, :-1]; Z_s = Z_ea[i, :-1]
R_c = np.append(R_s, R_s[0]); Z_c = np.append(Z_s, Z_s[0])
ds_s = np.sqrt(np.diff(R_c)**2 + np.diff(Z_c)**2)
ds_s_norm = ds_s / ds_s.mean()
ax.plot(TET_ea[:-1], ds_s_norm, color=colors[k], lw=0.8,
label=f'S={S[i]:.2f}' if k % 2 == 0 else None)
ax.axhline(1.0, color='k', ls='--', lw=1, label='Идеал (равномерно)')
ax.set_xlabel(r'$\theta_{\rm ea}$ (rad)')
ax.set_ylabel('Приращение длины дуги / среднее')
ax.set_title('Равномерность длины дуги в координатах equal-arc')
ax.set_xlim(0, 2*np.pi)
ax.legend(fontsize=8, ncol=2)
plt.tight_layout()
plt.savefig('equal_arc_coords.png', dpi=100, bbox_inches='tight')
plt.show()
print("Equal-arc: сегменты длины дуги равномерны в пределах ~1% (дискретизация интерполяции)")
Equal-arc: сегменты длины дуги равномерны в пределах ~1% (дискретизация интерполяции)
Раздел 4: PEST-координаты с прямыми силовыми линиями#
В координатах PEST коэффициент безопасности \(q(\psi)\) удовлетворяет:
\[q(\psi) = \frac{\mathbf{B}\cdot\nabla\varphi}{\mathbf{B}\cdot\nabla\theta^*} = {\rm const}\]
Это означает, что силовые линии выглядят как прямые линии в плоскости \((\theta^*, \varphi)\) с наклоном \(1/q\).
[5]:
# Show q(S) profile from PEST integration
q_valid = q_iS[1:]
S_valid = S[1:]
# Compare with analytic equilibrium q
psi_vals = S_valid**2
q_analytic = eq.q_profile(psi_vals, n_theta=256)
fig, axes = plt.subplots(1, 2, figsize=(12, 5))
ax = axes[0]
ax.plot(S_valid, q_valid, 'b-o', ms=3, label='PEST q(S)')
finite = np.isfinite(q_analytic)
ax.plot(S_valid[finite], q_analytic[finite], 'r--', label='Analytic q(S)')
ax.set_xlabel('S = √ψ_norm')
ax.set_ylabel('Коэффициент безопасности q')
ax.set_title('Профиль q: PEST против аналитики')
ax.legend()
ax.grid(True, alpha=0.3)
# PEST mesh in R-Z
ax = axes[1]
cs = ax.contour(Rg, Zg, psi_norm_grid, levels=np.linspace(0, 1, 11),
colors='lightgray', linewidths=0.5)
stride_s = max(1, ns // 10)
colors_pest = cm.viridis(np.linspace(0.2, 0.9, ns // stride_s + 1))
for k, i in enumerate(range(1, ns, stride_s)):
ax.plot(R_mesh[i, :], Z_mesh[i, :], color=colors_pest[k], lw=1.0)
for j in range(0, ntheta-1, ntheta // 8):
ax.plot(R_mesh[1:, j], Z_mesh[1:, j], 'g-', lw=0.5, alpha=0.5)
ax.plot(Rmaxis, Zmaxis, 'r+', ms=10, mew=2)
ax.set_xlabel('R (м)')
ax.set_ylabel('Z (м)')
ax.set_title('Сетка PEST (вид R-Z)')
ax.set_aspect('equal')
plt.tight_layout()
plt.savefig('pest_coords.png', dpi=100, bbox_inches='tight')
plt.show()
[6]:
# Demonstrate straight field lines in PEST theta-phi space
# In PEST: a field line traces theta* = phi/q + const
fig, ax = plt.subplots(figsize=(8, 5))
phi_line = np.linspace(0, 4 * np.pi, 200)
surface_indices = [ns // 5, ns // 3, ns // 2, 2 * ns // 3, 4 * ns // 5]
colors_fl = cm.Set1(np.linspace(0, 0.8, len(surface_indices)))
for k, i_s in enumerate(surface_indices):
if i_s >= ns or not np.isfinite(q_iS[i_s]):
continue
q_s = q_iS[i_s]
# Field line: theta* = phi / q (PEST straight field line)
theta_fl = (phi_line / q_s) % (2 * np.pi)
ax.plot(phi_line / np.pi, theta_fl / np.pi, color=colors_fl[k],
label=f'S={S[i_s]:.2f}, q={q_s:.2f}', lw=1.5)
ax.set_xlabel(r'$\varphi / \pi$')
ax.set_ylabel(r'$\theta^* / \pi$ (PEST)')
ax.set_title('Силовые линии в PEST $(\\theta^*, \\varphi)$ space — прямые линии с наклоном $1/q$')
ax.legend(fontsize=9)
ax.set_xlim(0, 4)
ax.set_ylim(0, 2)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('pest_fieldlines.png', dpi=100, bbox_inches='tight')
plt.show()
Раздел 5: координаты Boozer#
Координаты Boozer расширяют PEST, дополнительно требуя, чтобы якобиан \(\sqrt{g}\) был поверхностной величиной (зависел только от \(\psi\), а не от \(\theta\)). Это обеспечивает \(\mathbf{J}\cdot\nabla\varphi = {\rm const}(\psi)\).
Построение: \(\theta_B = \theta^* + \lambda(\psi, \theta^*)\), где
\[\frac{\partial\lambda}{\partial\theta^*} = \frac{\langle\sqrt{g}\rangle}{\sqrt{g}} - 1\]
а \(\langle\cdot\rangle\) обозначает среднее по магнитной поверхности \(\frac{1}{2\pi}\int_0^{2\pi}\cdot\,d\theta^*\).
[7]:
# Build Boozer mesh
_, TET_B, R_B, Z_B, lambda_corr = build_Boozer_mesh(
S, TET, R_mesh, Z_mesh, q_iS, equilibrium=eq, n_theta=ntheta
)
print(f"Сетка Boozer построена: форма {R_B.shape}")
print(f"Макс. поправка |λ|: {np.nanmax(np.abs(lambda_corr)):.4f} rad")
# Show the angle correction lambda(psi, theta*)
fig, axes = plt.subplots(1, 2, figsize=(12, 5))
ax = axes[0]
# lambda correction as a function of theta for several поверхностей
for k, i_s in enumerate(range(2, ns-1, max(1, ns // 8))):
lam = lambda_corr[i_s, :]
if np.any(np.isfinite(lam)):
ax.plot(TET / np.pi, lam, label=f'S={S[i_s]:.2f}',
color=cm.plasma(0.1 + 0.8 * i_s / ns))
ax.set_xlabel(r'$\theta^* / \pi$ (PEST)')
ax.set_ylabel(r'$\lambda = \theta_B - \theta^*$ (rad)')
ax.set_title('Угловая поправка Boozer $\\lambda(\\psi, \\theta^*)$')
ax.legend(fontsize=8, ncol=2)
ax.grid(True, alpha=0.3)
ax = axes[1]
cs = ax.contour(Rg, Zg, psi_norm_grid, levels=np.linspace(0, 1, 11),
colors='lightgray', linewidths=0.5)
stride_s = max(1, ns // 10)
colors_b = cm.inferno(np.linspace(0.2, 0.9, ns // stride_s + 1))
for k, i in enumerate(range(1, ns, stride_s)):
ax.plot(R_B[i, :], Z_B[i, :], color=colors_b[k], lw=1.0)
for j in range(0, ntheta-1, ntheta // 8):
ax.plot(R_B[1:, j], Z_B[1:, j], 'm-', lw=0.5, alpha=0.5)
ax.plot(Rmaxis, Zmaxis, 'r+', ms=10, mew=2)
ax.set_xlabel('R (м)')
ax.set_ylabel('Z (м)')
ax.set_title('Сетка Boozer (вид R-Z)')
ax.set_aspect('equal')
plt.tight_layout()
plt.savefig('boozer_coords.png', dpi=100, bbox_inches='tight')
plt.show()
Сетка Boozer построена: форма (40, 181)
Макс. поправка |λ|: 0.5429 рад
Раздел 6: координаты Hamada#
Координаты Hamada требуют, чтобы и силовые линии, и линии плотности тока были прямыми:
\[\mathbf{J}\cdot\nabla\psi = 0, \quad \mathbf{J}\cdot\nabla\theta_H = 0, \quad \mathbf{J}\cdot\nabla\varphi_H = {\rm const}(\psi)\]
Для осесимметричных равновесий это сводится к тому, что полоидальный угол становится равноплощадным: площадь, заметенная от магнитной оси до угла \(\theta_H\), пропорциональна \(\theta_H\).
[8]:
# Build Hamada mesh
_, TET_H, R_H, Z_H = build_Hamada_mesh(S, TET, R_mesh, Z_mesh, n_theta=ntheta)
print(f"Сетка Hamada построена: форма {R_H.shape}")
# Verify equal-area property
R_ax = R_mesh[0, 0]
Z_ax = Z_mesh[0, 0]
fig, axes = plt.subplots(1, 2, figsize=(12, 5))
ax = axes[0]
# Show cumulative area vs theta_H for several поверхностей
for k, i_s in enumerate(range(2, ns-2, max(1, ns // 8))):
R_s = R_H[i_s, :-1]; Z_s = Z_H[i_s, :-1]
R_c = np.append(R_s, R_s[0]); Z_c = np.append(Z_s, Z_s[0])
dA = 0.5 * ((R_c[:-1] - R_ax) * (Z_c[1:] - Z_ax) -
(R_c[1:] - R_ax) * (Z_c[:-1] - Z_ax))
A_cum = np.cumsum(dA)
A_total = A_cum[-1]
if abs(A_total) > 1e-10:
ax.plot(TET_H[:-1] / np.pi, A_cum / A_total,
color=cm.cool(0.1 + 0.8 * i_s / ns),
label=f'S={S[i_s]:.2f}')
# Ideal line
ax.plot([0, 2], [0, 1], 'k--', lw=2, label='Ideal (linear)')
ax.set_xlabel(r'$\theta_H / \pi$ (Hamada)')
ax.set_ylabel('Накопленная площадь / всего')
ax.set_title('Hamada: накопленная площадь линейна по $\\theta_H$')
ax.legend(fontsize=8, ncol=2)
ax.grid(True, alpha=0.3)
ax = axes[1]
cs = ax.contour(Rg, Zg, psi_norm_grid, levels=np.linspace(0, 1, 11),
colors='lightgray', linewidths=0.5)
colors_h = cm.cool(np.linspace(0.2, 0.9, ns // stride_s + 1))
for k, i in enumerate(range(1, ns, stride_s)):
ax.plot(R_H[i, :], Z_H[i, :], color=colors_h[k], lw=1.0)
for j in range(0, ntheta-1, ntheta // 8):
ax.plot(R_H[1:, j], Z_H[1:, j], 'c-', lw=0.5, alpha=0.5)
ax.plot(Rmaxis, Zmaxis, 'r+', ms=10, mew=2)
ax.set_xlabel('R (м)')
ax.set_ylabel('Z (м)')
ax.set_title('Сетка Hamada (вид R-Z)')
ax.set_aspect('equal')
plt.tight_layout()
plt.savefig('hamada_coords.png', dpi=100, bbox_inches='tight')
plt.show()
Сетка Hamada построена: форма (40, 181)
Раздел 7: панель сравнения - все четыре координатные системы#
Здесь все четыре сетки накладываются в плоскости R-Z для одного и того же набора магнитных поверхностей.
[9]:
fig = plt.figure(figsize=(16, 14))
gs = GridSpec(2, 2, figure=fig, hspace=0.3, wspace=0.3)
meshes = [
('PEST', TET, R_mesh, Z_mesh, 'viridis', 'g'),
('Equal-arc', TET_ea, R_ea, Z_ea, 'plasma', 'b'),
('Boozer', TET_B, R_B, Z_B, 'inferno', 'm'),
('Hamada', TET_H, R_H, Z_H, 'cool', 'c'),
]
stride_s = max(1, ns // 8)
stride_t = ntheta // 12
for idx, (name, tet, R_m, Z_m, cmap, pline_color) in enumerate(meshes):
ax = fig.add_subplot(gs[idx // 2, idx % 2])
# Flux surface contours as background
ax.contour(Rg, Zg, psi_norm_grid, levels=np.linspace(0.05, 0.95, 10),
colors='lightgray', linewidths=0.4)
# Mesh поверхностей (constant S)
surface_colors = plt.get_cmap(cmap)(np.linspace(0.2, 0.9, ns // stride_s + 1))
for k, i in enumerate(range(1, ns, stride_s)):
ax.plot(R_m[i, :], Z_m[i, :], color=surface_colors[k], lw=1.2)
# Poloidal lines (constant theta)
for j in range(0, len(tet)-1, stride_t):
ax.plot(R_m[1:, j], Z_m[1:, j], color=pline_color,
lw=0.6, alpha=0.6)
ax.plot(Rmaxis, Zmaxis, 'r+', ms=12, mew=2)
ax.set_xlabel('R (м)', fontsize=11)
ax.set_ylabel('Z (м)', fontsize=11)
ax.set_title(f'{name} координаты', fontsize=13, fontweight='bold')
ax.set_aspect('equal')
# Add annotation
ann = {'PEST': 'Straight B field lines',
'Equal-arc': 'Uniform arc length',
'Boozer': 'Straight B + uniform Jacobian',
'Hamada': 'Straight B + equal area'}[name]
ax.text(0.05, 0.95, ann, transform=ax.transAxes, fontsize=9,
verticalalignment='top', bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.5))
fig.suptitle('Сравнение магнитных координатных систем\n(coloured lines = flux поверхностей, thin lines = poloidal mesh)',
fontsize=14, y=1.01)
plt.savefig('coords_comparison.png', dpi=100, bbox_inches='tight')
plt.show()
Раздел 8: физическая интерпретация#
8.1 Связь Фурье-мод#
В координатах с прямыми силовыми линиями (PEST, Boozer, Hamada) возмущение с одним тороидальным номером моды \(n\) резонирует на поверхности, где \(q = m/n\). Фурье-разложение по \(\theta\) имеет чистую модовую структуру.
Напротив, при геометрическом угле (или Equal-arc) одна мода в \((m, n)\) проявляется как несколько Фурье-компонент, связывая разные значения \(m\); это так называемая проблема связи Фурье.
[10]:
# Demonstrate: Фурье-спектр of R(theta) — 2D heatmap (m vs S) for PEST and equal-arc
m_show = 12 # show modes m=0..m_show-1
# Compute FFT of R(theta) - <R> at each surface for both coordinate systems
S_arr = np.linspace(0.1, 0.9, ns) # approx normalised flux label per surface
fft_PEST = np.zeros((ns, m_show))
fft_EA = np.zeros((ns, m_show))
for i in range(ns):
R_pest = R_mesh[i, :-1] - R_mesh[i, :-1].mean()
R_ea_i = R_ea[i, :-1] - R_ea[i, :-1].mean()
n_pts_pest = len(R_pest)
n_pts_ea = len(R_ea_i)
fft_p = np.abs(np.fft.rfft(R_pest))[:m_show]
fft_p = fft_p / (fft_p.max() + 1e-30)
fft_PEST[i, :len(fft_p)] = fft_p[:m_show]
fft_e = np.abs(np.fft.rfft(R_ea_i))[:m_show]
fft_e = fft_e / (fft_e.max() + 1e-30)
fft_EA[i, :len(fft_e)] = fft_e[:m_show]
fig, axes = plt.subplots(1, 2, figsize=(13, 4.5))
for ax, data, name in zip(
axes,
[fft_PEST, fft_EA],
['PEST (straight field line)', 'Equal-arc θ'],
):
# log scale heatmap: rows = radial surface index, cols = mode m
log_data = np.log10(data + 1e-6)
im = ax.imshow(
log_data.T, # форма (m_show, ns): mode vs surface
origin='lower',
aspect='auto',
cmap='hot_r',
vmin=log_data.max() - 3, # show 3 decades
vmax=log_data.max(),
extent=[S_arr[0], S_arr[-1], -0.5, m_show - 0.5],
interpolation='nearest',
)
cbar = fig.colorbar(im, ax=ax, pad=0.02, shrink=0.85)
cbar.set_label(r'$\log_{10}$ normalised FFT amplitude', fontsize=8)
ax.set_xlabel(r'Нормированная метка потока $S$', fontsize=10)
ax.set_ylabel('Полоидальная мода m', fontsize=10)
ax.set_yticks(np.arange(0, m_show, 2))
ax.set_title(f'R(θ) Фурье-спектр: {name}', fontsize=11)
# Annotation: PEST should be concentrated at m=1; equal-arc spreads
axes[0].text(0.05, 0.9, 'Spectrum concentrated\nat m=1 (ideal)',
transform=axes[0].transAxes, fontsize=8, color='white',
bbox=dict(boxstyle='round', facecolor='steelblue', alpha=0.5))
axes[1].text(0.05, 0.9, 'Spectrum spreads to\nhigher m (non-ideal)',
transform=axes[1].transAxes, fontsize=8, color='white',
bbox=dict(boxstyle='round', facecolor='tomato', alpha=0.5))
plt.suptitle('Фурье-состав: PEST против equal-arc', fontsize=12)
plt.tight_layout()
plt.show()
8.2 Сводная таблица#
[11]:
summary = """
╔══════════════╦═══════════════════════╦══════════════════════════════╦══════════════════════════╗
║ Coordinate ║ Definition of θ ║ Main advantage ║ Used in ║
╠══════════════╬═══════════════════════╬══════════════════════════════╬══════════════════════════╣
║ PEST ║ Straight B field lines║ Minimal Fourier coupling ║ PEST, GATO, ELITE codes ║
║ ║ B·∇θ*/B·∇φ = q(ψ) ║ q = m/n at resonance ║ ║
╠══════════════╬═══════════════════════╬══════════════════════════════╬══════════════════════════╣
║ Boozer ║ PEST + √g = √g(ψ) ║ 1/B² drift is purely radial ║ коды стелларатора, вывод VMEC ║
║ ║ (uniform Jacobian) ║ Cleaner neoclassical theory ║ ║
╠══════════════╬═══════════════════════╬══════════════════════════════╬══════════════════════════╣
║ Hamada ║ Equal area from axis ║ J·∇θ = const, MHD stability ║ TERPSICHORE, CASTOR ║
║ ║ ∝ swept poloidal area ║ matrices simplified ║ ║
╠══════════════╬═══════════════════════╬══════════════════════════════╬══════════════════════════╣
║ Equal-arc ║ Uniform arc length ║ Simple construction ║ Numerical grids, FEM ║
║ ║ ds/dθ_ea = const ║ Resolves boundary well ║ ║
╚══════════════╩═══════════════════════╩══════════════════════════════╩══════════════════════════╝
"""
print(summary)
╔══════════════╦═══════════════════════╦══════════════════════════════╦══════════════════════════╗
║ Координаты ║ Определение θ ║ Главное преимущество ║ Используется в ║
╠══════════════╬═══════════════════════╬══════════════════════════════╬══════════════════════════╣
║ PEST ║ Прямые линии B ║ Минимальная связь Фурье ║ PEST, GATO, ELITE ║
║ ║ B·∇θ*/B·∇φ = q(ψ) ║ q = m/n в резонансе ║ ║
╠══════════════╬═══════════════════════╬══════════════════════════════╬══════════════════════════╣
║ Boozer ║ PEST + √g = √g(ψ) ║ дрейф 1/B² чисто радиален ║ коды stellarator, VMEC ║
║ ║ (равномерный якобиан) ║ чище неоклассическая теория ║ ║
╠══════════════╬═══════════════════════╬══════════════════════════════╬══════════════════════════╣
║ Hamada ║ Равная площадь от оси ║ J·∇θ = const, устойчивость MHD║ TERPSICHORE, CASTOR ║
║ ║ ∝ заметенной площади ║ матрицы упрощаются ║ ║
╠══════════════╬═══════════════════════╬══════════════════════════════╬══════════════════════════╣
║ Equal-arc ║ Равная длина дуги ║ Простое построение ║ численные сетки, FEM ║
║ ║ ds/dθ_ea = const ║ хорошо разрешает границу ║ ║
╚══════════════╩═══════════════════════╩══════════════════════════════╩══════════════════════════╝
8.3 Связь между системами#
Все четыре системы связаны преобразованиями угла вида \(\theta_{\rm new} = \theta^* + f(\psi, \theta^*)\):
Equal-arc -> параметризация с равной длиной дуги (чисто геометрическая)
PEST -> равномерная намотка силовой линии (требует интегрирования \(B_R, B_Z\))
Boozer -> PEST + периодическая поправка для выравнивания якобиана
Hamada -> равноплощадное преобразование (использует охваченную площадь, эквивалентно тому, что \(\sqrt{g}\) становится пропорциональным полной площади поверхности)
Для токамака Boozer и Hamada часто очень похожи, потому что обе усредненные по магнитной поверхности величины (якобиан и площадь) управляются одним и тем же балансом давления.
[12]:
# Final comparison: all four theta angles for a single field line
# Show how theta_PEST, theta_B, theta_H, theta_ea relate on the midradius surface
i_surf = ns // 2
print(f"Сравнение представлений полоидального угла на поверхности S={S[i_surf]:.3f}")
print(f" Коэффициент безопасности q = {q_iS[i_surf]:.3f}")
# For each coordinate system, compute the geometric angle (atan2(Z-Zax, R-Rax))
fig, ax = plt.subplots(figsize=(10, 5))
for name, tet_arr, R_m, Z_m, color in [
('PEST', TET, R_mesh, Z_mesh, 'blue'),
('Equal-arc', TET_ea, R_ea, Z_ea, 'green'),
('Boozer', TET_B, R_B, Z_B, 'red'),
('Hamada', TET_H, R_H, Z_H, 'purple'),
]:
R_s = R_m[i_surf, :]
Z_s = Z_m[i_surf, :]
theta_geom = np.arctan2(Z_s - Zmaxis, R_s - Rmaxis) % (2 * np.pi)
ax.plot(tet_arr / np.pi, theta_geom / np.pi, label=name, color=color, lw=1.5)
# Diagonal = geometric angle equals coordinate angle
ax.plot([0, 2], [0, 2], 'k--', lw=1, alpha=0.5, label='Identity (geom = coord)')
ax.set_xlabel(r'Координатный угол $\theta / \pi$')
ax.set_ylabel(r'Геометрический угол $\theta_{\rm geom} / \pi$')
ax.set_title('Геометрический угол vs coordinate angle for each system\n' +
f'(surface S={S[i_surf]:.2f}, q={q_iS[i_surf]:.2f})')
ax.legend()
ax.grid(True, alpha=0.3)
ax.set_xlim(0, 2)
ax.set_ylim(0, 2)
plt.tight_layout()
plt.savefig('theta_comparison.png', dpi=100, bbox_inches='tight')
plt.show()
print("Все четыре системы отличаются только распределением θ вокруг магнитной поверхности.")
Сравнение представлений полоидального угла на поверхности S=0.474
Коэффициент безопасности q = 1.638
Все четыре системы отличаются только распределением θ вокруг магнитной поверхности.