Systèmes de coordonnées magnétiques dans les équilibres de tokamak#
Ce tutoriel compare quatre systèmes de coordonnées de flux utilisés en fusion magnétique : PEST, Boozer, Hamada et les coordonnées à arc égal (Equal-arc).
Que sont les coordonnées de flux ?#
Dans un tokamak, les lignes de champ magnétique reposent sur des surfaces toroïdales emboîtées appelées surfaces de flux. Un système de coordonnées de flux \((\psi, \theta, \varphi)\) vérifie :
\(\psi\) étiquette les surfaces de flux, par exemple \(\psi_{\rm norm} = 0\) sur l’axe et \(1\) à la LCFS
\(\varphi\) est l’angle toroïdal standard
\(\theta\) est un angle poloïdal défini pour avoir des propriétés utiles
Les quatre systèmes ne diffèrent que par le choix de \(\theta\) :
Système |
Définition de l’angle poloïdal \(\theta\) |
Propriété clé |
|---|---|---|
PEST |
\(\partial(\mathbf{B}\cdot\nabla\theta^*) / \partial\theta^* = 0\) : lignes de champ droites |
Coordonnées à lignes de champ droites les plus simples |
Boozer |
PEST + jacobien fonction du flux : \(\sqrt{g} = \sqrt{g}(\psi)\) |
\(\mathbf{J}\cdot\nabla\varphi = {\rm const}(\psi)\) ; utilisé pour les sorties de codes stellarator/tokamak |
Hamada |
Aire égale depuis l’axe magnétique |
\(\mathbf{J}\cdot\nabla\theta = {\rm const}(\psi)\) ; simplifie les matrices de stabilité MHD |
Arc égal (Equal-arc) |
La longueur d’arc le long de la surface de flux est uniforme en \(\theta_{\rm ea}\) |
Robuste numériquement ; construction simple |
Pourquoi ces choix comptent-ils ?#
Lignes de champ droites (PEST, Boozer, Hamada) : un mode de nombre toroïdal \(n\) et de nombre poloïdal \(m\) satisfait \(q = m/n\) à la résonance. Cela simplifie l’analyse de Fourier des modes MHD.
Boozer : la dérive \(1/B^2\) est purement radiale, ce qui simplifie la théorie du transport néoclassique et l’équation cinétique.
Hamada : les lignes de courant sont également droites, ce qui est nécessaire pour certains codes de stabilité MHD.
Arc égal (Equal-arc) : pratique pour les grilles numériques ; décrit bien la frontière du plasma.
Mise en place : équilibre analytique de Solov’ev#
Nous utilisons l’équilibre de Solov’ev de Cerfon & Freidberg (2010), mis à l’échelle d’une taille de tokamak de référence (\(R_0 \approx 1{,}86\) m).
[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 → de taille tokamak de référence, R0≈1.86 m)
eq = solovev_iter_like(scale=0.3)
print(f"Équilibre : 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"Axe magnétique : R={Rmaxis:.3f} m, Z={Zmaxis:.3f} m")
Équilibre : R0=1.86 m, a=0.60 m, B0=5.3 T
κ=1.70, δ=0.33, q0=1.50
Axe magnétique : R=1.946 m, Z=0.000 m
[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 surfaces, ntheta=181 points poloïdaux
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"Maillage PEST construit : {ns} surfaces, {ntheta} points poloïdaux")
print(f"plage S : [{S[1]:.3f}, {S[-1]:.3f}]")
print(f"plage q : [{q_iS[1]:.2f}, {q_iS[-1]:.2f}]")
Maillage PEST construit : 40 surfaces, 181 points poloïdaux
plage S : [0.028, 0.972]
plage q : [1.49, 2.39]
Section 3 : coordonnées à arc égal (Equal-arc)#
L’angle à arc égal \(\theta_{\rm ea}\) paramètre chaque surface de flux de sorte que la longueur d’arc le long de la surface soit proportionnelle à \(\theta_{\rm ea}\). C’est la transformation de coordonnées non triviale la plus simple.
[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"Maillage Equal-arc construit : shape {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"Segments de longueur d’arc sur la surface {i_mid}: moyenne={ds.mean()*100:.2f} cm, écart-type={ds.std()*100:.3f} cm")
print(f" Uniformité (écart-type/moyenne) : {ds.std()/ds.mean():.4f}")
Maillage Equal-arc construit : shape (40, 181)
Segments de longueur d’arc sur la surface 20: moyenne=1.27 cm, écart-type=0.000 cm
Uniformité (écart-type/moyenne) : 0.0001
[4]:
# Plot equal-arc coordonnées
fig, axes = plt.subplots(1, 2, figsize=(12, 5))
# Left: R-Z view of flux surfaces 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='Axe')
ax.set_xlabel('R (m)')
ax.set_ylabel('Z (m)')
ax.set_title('Maillage Equal-arc (vue 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='Idéal (uniforme)')
ax.set_xlabel(r'$\theta_{\rm ea}$ (rad)')
ax.set_ylabel('Incrément de longueur d’arc / moyenne')
ax.set_title('Uniformité de longueur d’arc en coordonnées 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 : les segments de longueur d’arc sont uniformes à ~1 % près (discrétisation d’interpolation)")
Equal-arc : les segments de longueur d’arc sont uniformes à ~1 % près (discrétisation d’interpolation)
Section 4 : coordonnées PEST à lignes de champ droites#
En coordonnées PEST, le facteur de sécurité \(q(\psi)\) vérifie :
\[q(\psi) = \frac{\mathbf{B}\cdot\nabla\varphi}{\mathbf{B}\cdot\nabla\theta^*} = {\rm const. \; sur \; chaque \; surface}\]
Cela signifie que les lignes de champ apparaissent comme des droites dans le plan \((\theta^*, \varphi)\), de pente \(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='q(S) analytique')
ax.set_xlabel('S = √ψ_norm')
ax.set_ylabel('Facteur de sécurité q')
ax.set_title('Profil du facteur de sécurité : PEST vs analytique')
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 (m)')
ax.set_ylabel('Z (m)')
ax.set_title('Maillage PEST (vue R-Z)')
ax.set_aspect('equal')
plt.tight_layout()
plt.savefig('pest_coords.png', dpi=100, bbox_inches='tight')
plt.show()
[6]:
# Demonstrate straight lignes de champ 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('Lignes de champ dans PEST $(\\theta^*, \\varphi)$ space — droites de pente $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()
Section 5 : coordonnées de Boozer#
Les coordonnées de Boozer étendent PEST en imposant en plus que le jacobien \(\sqrt{g}\) soit une grandeur de surface de flux, dépendant seulement de \(\psi\) et non de \(\theta\). Cela garantit \(\mathbf{J}\cdot\nabla\varphi = {\rm const}(\psi)\).
Construction : \(\theta_B = \theta^* + \lambda(\psi, \theta^*)\) où
\[\frac{\partial\lambda}{\partial\theta^*} = \frac{\langle\sqrt{g}\rangle}{\sqrt{g}} - 1\]
avec \(\langle\cdot\rangle\) la moyenne sur surface de flux \(\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"Maillage Boozer construit : shape {R_B.shape}")
print(f"Correction max |λ| : {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 surfaces
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('Correction d’angle 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 (m)')
ax.set_ylabel('Z (m)')
ax.set_title('Maillage Boozer (vue R-Z)')
ax.set_aspect('equal')
plt.tight_layout()
plt.savefig('boozer_coords.png', dpi=100, bbox_inches='tight')
plt.show()
Maillage Boozer construit : shape (40, 181)
Correction max |λ| : 0.5429 rad
Section 6 : coordonnées de Hamada#
Les coordonnées de Hamada exigent que les lignes de champ et les lignes de densité de courant soient droites :
\[\mathbf{J}\cdot\nabla\psi = 0, \quad \mathbf{J}\cdot\nabla\theta_H = 0, \quad \mathbf{J}\cdot\nabla\varphi_H = {\rm const}(\psi)\]
Pour les équilibres axisymétriques, cela revient à rendre l’angle poloïdal à aire égale : l’aire balayée depuis l’axe magnétique jusqu’à l’angle \(\theta_H\) est proportionnelle à \(\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"Maillage Hamada construit : shape {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 surfaces
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='Idéal (linéaire)')
ax.set_xlabel(r'$\theta_H / \pi$ (Hamada)')
ax.set_ylabel('Aire cumulée / total')
ax.set_title('Hamada : l’aire cumulée est linéaire in $\\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 (m)')
ax.set_ylabel('Z (m)')
ax.set_title('Maillage Hamada (vue R-Z)')
ax.set_aspect('equal')
plt.tight_layout()
plt.savefig('hamada_coords.png', dpi=100, bbox_inches='tight')
plt.show()
Maillage Hamada construit : shape (40, 181)
Section 7 : panneau de comparaison - les quatre systèmes de coordonnées#
Ici, nous superposons les quatre grilles dans le plan R-Z pour le même ensemble de surfaces de flux.
[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 surfaces (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 (m)', fontsize=11)
ax.set_ylabel('Z (m)', fontsize=11)
ax.set_title(f'{name} coordonnées', fontsize=13, fontweight='bold')
ax.set_aspect('equal')
# Add annotation
ann = {'PEST': 'Lignes de champ B droites',
'Equal-arc': 'Longueur d’arc uniforme',
'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('Comparaison des systèmes de coordonnées magnétiques\n(lignes colorées = surfaces de flux, lignes fines = maillage poloïdal)',
fontsize=14, y=1.01)
plt.savefig('coords_comparison.png', dpi=100, bbox_inches='tight')
plt.show()
Section 8 : interprétation physique#
8.1 Couplage des modes de Fourier#
Dans les coordonnées à lignes de champ droites (PEST, Boozer, Hamada), une perturbation avec un seul nombre de mode toroïdal \(n\) résonne sur la surface où \(q = m/n\). La décomposition de Fourier en \(\theta\) possède une structure de modes propre.
En revanche, avec l’angle géométrique ou l’angle à arc égal (Equal-arc), un mode unique en \((m,n)\) apparaît comme plusieurs composantes de Fourier et couple différentes valeurs de \(m\) : c’est le problème de couplage de Fourier.
[10]:
# Demonstrate: Spectre de Fourier 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, # shape (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'Étiquette de flux normalisée $S$', fontsize=10)
ax.set_ylabel('Mode poloïdal m', fontsize=10)
ax.set_yticks(np.arange(0, m_show, 2))
ax.set_title(f'R(θ) Spectre de Fourier: {name}', fontsize=11)
# Annotation: PEST devrait être 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('Contenu de Fourier : PEST vs coordonnées Equal-arc', fontsize=12)
plt.tight_layout()
plt.show()
8.2 Tableau de synthèse#
[11]:
summary = """
╔══════════════╦═══════════════════════╦══════════════════════════════╦══════════════════════════╗
║ Coordonnée ║ Définition de θ ║ Avantage principal ║ Utilisé dans ║
╠══════════════╬═══════════════════════╬══════════════════════════════╬══════════════════════════╣
║ PEST ║ lignes B droites ║ Couplage Fourier minimal ║ codes MHD génériques ║
║ ║ B·∇θ*/B·∇φ = q(ψ) ║ q = m/n à la résonance ║ ║
╠══════════════╬═══════════════════════╬══════════════════════════════╬══════════════════════════╣
║ Boozer ║ PEST + √g = √g(ψ) ║ dérive 1/B² purement radiale ║ sorties de codes ║
║ ║ (jacobien uniforme) ║ théorie néoclassique épurée ║ ║
╠══════════════╬═══════════════════════╬══════════════════════════════╬══════════════════════════╣
║ Hamada ║ aire égale depuis axe ║ J·∇θ = const, stabilité MHD ║ codes de stabilité ║
║ ║ ∝ aire poloïdale ║ matrices simplifiées ║ ║
╠══════════════╬═══════════════════════╬══════════════════════════════╬══════════════════════════╣
║ Equal-arc ║ longueur arc uniforme ║ construction simple ║ grilles numériques, FEM ║
║ ║ ds/dθ_ea = const ║ bonne résolution du bord ║ ║
╚══════════════╩═══════════════════════╩══════════════════════════════╩══════════════════════════╝
"""
print(summary)
╔══════════════╦═══════════════════════╦══════════════════════════════╦══════════════════════════╗
║ Coordonnée ║ Définition de θ ║ Avantage principal ║ Utilisé dans ║
╠══════════════╬═══════════════════════╬══════════════════════════════╬══════════════════════════╣
║ PEST ║ lignes B droites ║ Couplage Fourier minimal ║ codes MHD génériques ║
║ ║ B·∇θ*/B·∇φ = q(ψ) ║ q = m/n à la résonance ║ ║
╠══════════════╬═══════════════════════╬══════════════════════════════╬══════════════════════════╣
║ Boozer ║ PEST + √g = √g(ψ) ║ dérive 1/B² purement radiale ║ sorties de codes ║
║ ║ (jacobien uniforme) ║ théorie néoclassique épurée ║ ║
╠══════════════╬═══════════════════════╬══════════════════════════════╬══════════════════════════╣
║ Hamada ║ aire égale depuis axe ║ J·∇θ = const, stabilité MHD ║ codes de stabilité ║
║ ║ ∝ aire poloïdale ║ matrices simplifiées ║ ║
╠══════════════╬═══════════════════════╬══════════════════════════════╬══════════════════════════╣
║ Equal-arc ║ longueur arc uniforme ║ construction simple ║ grilles numériques, FEM ║
║ ║ ds/dθ_ea = const ║ bonne résolution du bord ║ ║
╚══════════════╩═══════════════════════╩══════════════════════════════╩══════════════════════════╝
8.3 Relation entre les systèmes#
Les quatre systèmes sont reliés par des transformations d’angle de la forme \(\theta_{\rm new} = \theta^* + f(\psi, \theta^*)\) :
Arc égal (Equal-arc) -> paramétrisation à longueur d’arc égale (purement géométrique)
PEST -> enroulement uniforme des lignes de champ (nécessite l’intégration de \(B_R, B_Z\))
Boozer -> PEST + correction périodique pour aplatir le jacobien
Hamada -> transformation à aire égale (utilise l’aire enfermée et revient à rendre \(\sqrt{g}\) proportionnel à l’aire totale de la surface)
Pour un tokamak, Boozer et Hamada sont souvent très similaires, car les deux grandeurs moyennées sur surface de flux, le jacobien et l’aire, sont gouvernées par le même équilibre de pression.
[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"Comparaison des représentations d’angle poloïdal sur la surface S={S[i_surf]:.3f}")
print(f" Facteur de sécurité 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='Identité (géom = coord)')
ax.set_xlabel(r'Angle coordonné $\theta / \pi$')
ax.set_ylabel(r'Angle géométrique $\theta_{\rm geom} / \pi$')
ax.set_title('Angle géométrique vs angle coordonné pour chaque système\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("Les quatre systèmes ne diffèrent que par la distribution de θ autour de la surface de flux.")
Comparaison des représentations d’angle poloïdal sur la surface S=0.474
Facteur de sécurité q = 1.638
Les quatre systèmes ne diffèrent que par la distribution de θ autour de la surface de flux.