"""Contains functionality to export equilibria to eqdsk format and files."""
import numpy as np
from pathlib import Path
from scipy.interpolate import interp1d
import scipy.integrate as sint
from torx.units import Quantity
from torx.decorators import autodoc_function
from .numerical_m import NumericalEquilibrium
from .io.eqdsk_io_m import write_eqdsk_file
from torx.geometry import Polygon2D
[docs]
@autodoc_function
def export_to_eqdsk(equi: NumericalEquilibrium, filepath: Path,
dims: tuple=(129, 129), name: str="TorX",
shot: str="00000", time: str="0.0",
comment: str="TorX", dry: bool=False):
"""
Export the equilibrium to eqdsk format.
Calculates estimations for poloidal current function and pressure based on
simple formulas without sophisticated constraints. Maps the poloidal flux
on an RZ matrix with bounds given by the limiter polygon and number of
points provided by dims. Creates a safety factor profile on a uniform
psi grid via field line tracing (may take some time). Returns the
dictionary that was written to file.
Can specify dry mode True to only get the content dictionary without
writing the file.
Note
----
For future improvements of the physics output, one might consider
using the actual poloidal current for F, Eq. (6.13) in [1]_, and
provide a boundary condition for the pressure at the separatrix.
References
----------
.. [1] J. P. Freidberg, "Ideal MHD", Cambridge University Press
(2014).
"""
assert isinstance(equi, NumericalEquilibrium), \
f"Equi must be of type numerical (was {type(equi)})!"
# NOTE: This routine uses mixed normalization, i.e. lengths and fields
# are normalized while poloidal fluxes are in Wb. At the end
# conversion to SI is performed.
meter = equi.R0.to("m").m
tesla = equi.B0.to("T").m
mu0 = Quantity(1, "mu0").to_base_units().m
nw = dims[0]
nh = dims[1]
rmaxis = equi.axis_r_norm.values
zmaxis = equi.axis_z_norm.values
rcentr = rmaxis
bcentr = equi.magfield_absolute(equi.axis_r_norm,
equi.axis_z_norm).values
simag = equi.psi_axis
sibry = equi.psi_separatrix
boundary = equi.get_boundary_polygon()
limiter = equi.get_limiter_polygon(npoints=boundary.x_points.size)
rmax = np.max(limiter.x_points)
# Use limiter polygon as "computational box"
rdim = np.abs(np.max(limiter.x_points) - np.min(limiter.x_points))
zdim = np.abs(np.max(limiter.y_points) - np.min(limiter.y_points))
rleft = np.min(limiter.x_points)
zmid = 0.5 * zdim
# Create RZ grid and calculate flux
R = np.linspace(rleft, rleft + rdim, nw)
Z = np.linspace(zmid - 0.5 * zdim, zmid + 0.5 * zdim, nh)
psirz = equi.psi_interpolator(x=R, y=Z, dx=0, dy=0, grid=True)
# Create flux grid and calculate quantities
R_interp = np.linspace(rcentr, rmax, nw)
psi_interp = equi.psi_interpolator(x=R_interp, y=equi.axis_z_norm, \
dx=0, dy=0, grid=False)
psi_u = np.linspace(simag, sibry, nw, endpoint=False)
R_u = interp1d(psi_interp, R_interp, kind="cubic",
fill_value="extrapolate")(psi_u)
# NOTE: The innermost point sits on the magnetic axis; safety_factor
# offsets it internally to avoid the singular near-axis trace.
r_u = R_u - equi.axis_r_norm.values
Z_u = np.ones_like(R_u) * equi.axis_z_norm.values
# NOTE: We have to convert to SI here since we add terms of different
# normalization.
lhs = equi.gradshaf_star(R_u, zmaxis) / meter**2
fpol = equi.magfield_component_toroidal(R_u, Z_u).values * R_u \
* tesla * meter
ffprim = fpol * np.gradient(fpol, psi_u)
# Eq. (6.15) inverted, see the docstring reference above
pprime = -(lhs + ffprim) / R_u**2 / meter**2 / mu0
# NOTE: We have no information about the pressure in the equi and the
# result from integrating pprime is determined up to a constant.
# We shift the profile such that there is no negative pressure at
# the outermost radial location for now.
pres = sint.cumulative_trapezoid(pprime, psi_u, initial=0)
pres = pres + np.abs(np.min(pres))
# Trace field lines for each psi to get the safety factor. Also get
# the trace result for further processing
qpsi, ode_sol = equi.safety_factor(R_u, Z_u, return_sol=True)
# Calculate flux surface area element and total current
theta = np.linspace(0.0, 2.0 * np.pi, nw, endpoint=False)
fs_area = np.array([Polygon2D(sol(theta)[0,:],
sol(theta)[1,:]).signed_area() \
for sol in ode_sol])
# Eq. (6.16), see the docstring reference above
current = -sint.trapezoid(lhs / R_u / meter, fs_area * meter**2) / mu0
# Boundary and limiter
nbbbs = boundary.x_points.size
rbbbs = boundary.x_points * meter
zbbbs = boundary.y_points * meter
limitr = limiter.x_points.size
rlim = limiter.x_points * meter
zlim = limiter.y_points * meter
# Converted lengths and magfields to SI
rmaxis = rmaxis * meter
zmaxis = zmaxis * meter
rcentr = rcentr * meter
rdim = rdim * meter
zdim = zdim * meter
rleft = rleft * meter
zmid = zmid * meter
bcentr = bcentr * tesla
content = dict(name=name, shot=shot, time=time, comment=comment,
nw=nw, nh=nh, rdim=rdim, zdim=zdim, rcentr=rcentr,
rleft=rleft, zmid=zmid, rmaxis=rmaxis, zmaxis=zmaxis,
simag=simag, sibry=sibry, bcentr=bcentr,
current=current, fpol=fpol, pres=pres, ffprim=ffprim,
pprime=pprime, psirz=psirz, qpsi=qpsi,
nbbbs=nbbbs, limitr=limitr, rbbbs=rbbbs,
zbbbs=zbbbs, rlim=rlim, zlim=zlim)
if not dry:
write_eqdsk_file(filepath, content)
return content
