"""Contains Dommaschk equilibrium."""
from pathlib import Path
import numpy as np
import xarray as xr
from torx import Quantity, make_xarray, genex_domm_resources_dir
from torx.units_handling_m import check_units
from torx.equilibrium.equilibrium_m import EquilibriumBaseClass
from torx.geometry.polygon_3d_m import Polygon3D
from torx.autogrid_decorators_m import autogrid_method
from torx.autodoc_decorators_m import autodoc_class
from ._dommaschk_initialization_m import _init_CD_CN, _init_Imn
from ._dommaschk_fitting_coefficients_m import DommaschkFittingCoefficients
[docs]
@autodoc_class
class DommaschkEquilibrium(EquilibriumBaseClass):
"""
Dommaschk type stellarator equilibrium.
Analytic representation of the vacuum magnetic field of stellarator
configurations, as described in [1].
[1] W. Dommaschk, "Representations for vacuum potentials in
stellarators", Computer Physics Communications 40, pg. 203 (1986).
"""
[docs]
def __init__(
self,
filepath: str="",
x0: float=1.0,
y0: float=0.0,
phi0: float=0.0,
num_field_periods: int=5,
l_pol: int=4,
m_tor_consecutive: int=1,
fitting_coef= None,
bndry_from_file: bool=False,
outer_bndry_file: str="",
outer_bndry_var: str="",
inner_bndry_file: str="",
inner_bndry_var: str="",
xmin: float= 0.80,
xmax: float= 1.20,
ymin: float=-0.20,
ymax: float= 0.20,
B0: Quantity=Quantity("1T"),
L_ref: float=1.0,
**kwargs
):
"""
Initialize the DommaschkEquilibrium.
Arguments to the initializer match the parameters supplied in PARALLAX.
You can either manually specify the arguments,
or pass a **kwargs_dictionary generated from a parameter file.
The B0 parameter sets the value of the toroidal magnetic field at the
normalized axis location (x0, y0, phi0).
The L_ref parameter sets the length normalization in meters. Note that
in 3D, the magnetic axis might not always be located at
R = L_ref meters.
"""
self.filepath = filepath
self.x0 = x0
self.y0 = y0
self.phi0 = phi0
self.num_field_periods = num_field_periods
self.l_pol = l_pol
self.m_tor_consecutive = m_tor_consecutive
self.bndry_from_file = bndry_from_file
self.outer_bndry_file = outer_bndry_file
self.outer_bndry_var = outer_bndry_var
self.inner_bndry_file = inner_bndry_file
self.inner_bndry_var = inner_bndry_var
self.xmin = xmin
self.xmax = xmax
self.ymin = ymin
self.ymax = ymax
self.B0 = B0
self._norm_length = Quantity(L_ref, "m")
self.m = xr.DataArray(np.arange(self.m_tor_consecutive + 1) \
* self.num_field_periods, dims=("m"))
if type(fitting_coef) == type(None):
self.fitting_coef = DommaschkFittingCoefficients(l_pol,
m_tor_consecutive)
elif type(fitting_coef) == DommaschkFittingCoefficients:
assert fitting_coef.l_pol == l_pol
assert fitting_coef.m_tor_consecutive == m_tor_consecutive
self.fitting_coef = fitting_coef
else:
self.fitting_coef = DommaschkFittingCoefficients(l_pol,
m_tor_consecutive,
fitting_coef)
_init_CD_CN(self)
_init_Imn(self)
self._init_B_norm()
self._init_boundary()
def _init_B_norm(self):
"""
Determine the magnetic field normalization.
Uses the non-normalized value of btor evaluated at the location of the
magnetic axis. This is used to normalize all subsequent magnetic
field calculations.
Set the normalization factor to 1 to calculate the non-normalized value
of btor.
"""
self._B_norm_internal = 1.0
btor_unnormalized = self.magfield_component_toroidal(self.x0, self.y0,
phi=self.phi0)
self._B_norm_internal = btor_unnormalized
def _init_boundary(self):
"""
Load the inner and/or outer boundary from file, if specified.
If no boundary is given, construct one from the box limits.
"""
default_polygon_file = genex_domm_resources_dir \
/ "flux_polygons_dommaschk_default.nc"
if self.outer_bndry_file == "default":
self.outer_bndry_file = default_polygon_file
if self.inner_bndry_file == "default":
self.inner_bndry_file = default_polygon_file
if self.bndry_from_file:
outer_boundary_filepath = Path(self.filepath) / Path(self.outer_bndry_file)
self.outer_bndry = Polygon3D.initialize_from_file(
outer_boundary_filepath,
self.outer_bndry_var,
num_field_periods=self.num_field_periods)
inner_boundary_filepath = Path(self.filepath) / Path(self.inner_bndry_file)
self.inner_bndry = Polygon3D.initialize_from_file(
inner_boundary_filepath,
self.inner_bndry_var,
num_field_periods=self.num_field_periods) \
if self.inner_bndry_file != "" else None
else:
# Construct the outer boundary polygon using the box limits in a
# counterclockwise direction
vertices = np.array([[[self.xmin, self.xmax, self.xmax, self.xmin],
[self.ymin, self.ymin, self.ymax, self.ymax]]])
self.outer_bndry = Polygon3D(
polygon_vertices=vertices,
phi_array=np.array(0.0),
num_field_periods=self.num_field_periods)
self.inner_bndry = None
@property
def axis_r_norm(self):
"""Magnetic axis radial position in normalized units."""
return make_xarray(self.x0)
@property
def axis_z_norm(self):
"""Magnetic axis vertical position in normalized units."""
return make_xarray(self.y0)
@property
def axis_phi(self):
"""Magnetic axis toroidal position."""
return make_xarray(self.phi0)
@property
def B0(self) -> Quantity:
"""Magnetic field normalization."""
return self._norm_magfield
@B0.setter
def B0(self, value: Quantity):
"""Set the magnetic field normalization."""
check_units(value, {"[mass]":1, "[time]":-2, "[current]":-1}, "B0")
self._norm_magfield = value
@property
def R0(self) -> Quantity:
"""Length normalization."""
return self._norm_length
@R0.setter
def R0(self, value: Quantity):
"""Set the length normalization."""
check_units(value, {"[length]":1}, "R0")
self._norm_length = value
[docs]
@autogrid_method
def magfield_component_r(self, r_norm, z_norm, **kwargs) \
-> xr.DataArray:
"""
Return the radial component of the vacuum magnetic field B.
Calculated according to B_R = dV/dR.
"""
if "phi" not in kwargs:
raise TypeError("Missing required keyword argument: 'phi'")
phi = kwargs["phi"]
R = self.convert_length_to_normalized(r_norm)
Z = self.convert_length_to_normalized(z_norm)
logR = np.log(R)
cos_m_phi = np.cos(self.m * phi)
sin_m_phi = np.sin(self.m * phi)
I_mn = self.I_mn_coef * Z**self.I_mn_power
dCD_mk_dR = (self.dCD_dR_r_coef \
* R**self.dC_dR_r_power).sum(dim="terms") \
+ (self.dCD_dR_lnr_coef \
* R**self.dC_dR_lnr_power * logR).sum(dim="terms")
dD_ml_dR = (I_mn.isel(l=slice(1, None)) * dCD_mk_dR).sum(dim="k")
dCN_mk_dR = (self.dCN_dR_r_coef \
* R**self.dC_dR_r_power).sum(dim="terms") \
+ (self.dCN_dR_lnr_coef \
* R**self.dC_dR_lnr_power * logR).sum(dim="terms")
dN_ml_1_dR = (I_mn.isel(l=slice(0, -1)) * dCN_mk_dR).sum(dim="k")
# Eq. 12 (differentiated with respect to R)
dV_ml_dR = ( self.fitting_coef[0, :, :] * cos_m_phi \
+ self.fitting_coef[1, :, :] * sin_m_phi \
) * dD_ml_dR \
+ ( self.fitting_coef[2, :, :] * cos_m_phi \
+ self.fitting_coef[3, :, :] * sin_m_phi \
) * dN_ml_1_dR
Br = dV_ml_dR.sum(("m", "l")) / self._B_norm_internal
return Br.assign_attrs(norm=self.B0)
[docs]
@autogrid_method
def magfield_component_z(self, r_norm, z_norm, **kwargs) \
-> xr.DataArray:
"""
Return the vertical component of the vacuum magnetic field B.
Calculated according to B_Z = dV/dZ.
"""
if "phi" not in kwargs:
raise TypeError("Missing required keyword argument: 'phi'")
phi = kwargs["phi"]
R = self.convert_length_to_normalized(r_norm)
Z = self.convert_length_to_normalized(z_norm)
logR = np.log(R)
cos_m_phi = np.cos(self.m * phi)
sin_m_phi = np.sin(self.m * phi)
dI_mn_dZ = self.dI_dZ_coef * Z**self.dI_dZ_power
CD_mk = (self.CD_r_coef * R**self.C_r_power).sum(dim="terms") \
+ (self.CD_lnr_coef * R**self.C_lnr_power * logR).sum(dim="terms")
dD_ml_dZ = (dI_mn_dZ.isel(l=slice(1, None)) * CD_mk).sum(dim="k")
CN_mk = (self.CN_r_coef * R**self.C_r_power).sum(dim="terms") \
+ (self.CN_lnr_coef * R**self.C_lnr_power * logR).sum(dim="terms")
dN_ml_1_dZ = (dI_mn_dZ.isel(l=slice(0, -1)) * CN_mk).sum(dim="k")
# Eq. 12 (differentiated with respect to Z)
dV_ml_dZ = ( self.fitting_coef[0, :, :] * cos_m_phi \
+ self.fitting_coef[1, :, :] * sin_m_phi \
) * dD_ml_dZ \
+ ( self.fitting_coef[2, :, :] * cos_m_phi \
+ self.fitting_coef[3, :, :] * sin_m_phi \
) * dN_ml_1_dZ
Bz = dV_ml_dZ.sum(("m", "l")) / self._B_norm_internal
return Bz.assign_attrs(norm=self.B0)
[docs]
@autogrid_method
def magfield_component_toroidal(self, r_norm, z_norm, **kwargs) \
-> xr.DataArray:
"""
Return the poloidal component of the vacuum magnetic field B.
Calculated according to B_phi = 1/R dV/dphi.
"""
if "phi" not in kwargs:
raise TypeError("Missing required keyword argument: 'phi'")
phi = kwargs["phi"]
R = self.convert_length_to_normalized(r_norm)
Z = self.convert_length_to_normalized(z_norm)
logR = np.log(R)
dsin_m_phi_dphi = self.m * np.cos(self.m * phi)
dcos_m_phi_dphi = -self.m * np.sin(self.m * phi)
I_mn = self.I_mn_coef * Z**self.I_mn_power
CD_mk = (self.CD_r_coef * R**self.C_r_power).sum(dim="terms") \
+ (self.CD_lnr_coef * R**self.C_lnr_power * logR).sum(dim="terms")
D_ml = (I_mn.isel(l=slice(1, None)) * CD_mk).sum(dim="k")
CN_mk = (self.CN_r_coef * R**self.C_r_power).sum(dim="terms") \
+ (self.CN_lnr_coef * R**self.C_lnr_power * logR).sum(dim="terms")
N_ml_1 = (I_mn.isel(l=slice(0, -1)) * CN_mk).sum(dim="k")
# Eq. 12 (differentiated with respect to phi)
dV_ml_dphi = ( self.fitting_coef[0, :, :] * dcos_m_phi_dphi \
+ self.fitting_coef[1, :, :] * dsin_m_phi_dphi \
) * D_ml \
+ ( self.fitting_coef[2, :, :] * dcos_m_phi_dphi \
+ self.fitting_coef[3, :, :] * dsin_m_phi_dphi \
) * N_ml_1
# The leading 1.0 is the derivative of the phi term in Eq. 1
Btor = (1.0 + dV_ml_dphi.sum(("m", "l"))) / (R * self._B_norm_internal)
return Btor.assign_attrs(norm=self.B0)
[docs]
def get_boundary_polygon(self, **kwargs):
"""Return the boundary polygon."""
return {"outer":self.outer_bndry, "inner":self.inner_bndry}
[docs]
def normalized_flux_surface_label(self, r_norm, z_norm, **kwargs):
"""Return the normalized flux surface label."""
raise NotImplementedError
