Physics models
TORAX provides a modular framework for incorporating various physics models. This section summarizes the presently implemented models for geometry, transport, sources, and neoclassical physics. Extended physics features are a focus of the short-term development roadmap. Specific plans are summarized following each subsection.
Magnetic geometry
TORAX presently supports four geometry models for determining the metric coefficients and flux surface averaged geometry variables required for the transport equations.
CHEASE: This model utilizes equilibrium data from the |chease| fixed boundary Grad-Shafranov equilibrium code. Users provide a CHEASE output file, and TORAX extracts the relevant geometric quantities.
FBT: This model utilizes equilibrium data from the |fbt| free boundary Grad-Shafranov equilibrium code. Users provide FBT output files, and TORAX extracts the relevant geometric quantities.
EQDSK: This model utilizes equilibrium data provided in the EQDSK format. Users provide EQDSK formatted files. TORAX extracts the relevant geometric quantities and carries out the necessary flux surface averaging based on traced contours of the 2D poloidal flux data.
Circular: For testing and demonstration purposes, TORAX includes a simplified circular geometry model. This model assumes a circular plasma cross-section and corrects for elongation to approximate the metric coefficients.
Using \(\psi\) and the magnetic geometry terms, the toroidal current density is calculated as follows:
where \(V\) is the volume enclosed by a flux surface.
The safety factor \(q\) is calculated as follows:
With \(q(\rho=0)=\frac{2B_0}{\mu_0 j_{tot}(\rho=0) R_0}\).
For simulations of tokamak scenarios with dynamic equilibra, sequences of geometry data can be provided. TORAX interpolates geometry data at each end of each timestep interval, using both sets of values in the theta solver method. Time-derivative terms are calculated based on finite differences. Generalization to enable self-consistent geometry calculations by models within the solver loop is planned.
Plasma composition, initial and prescribed conditions
Presently, TORAX accommodates a single main ion species and single impurity species, which can be comprised of time-dependent mixtures of ions with fractional abundances summing to 1. This is useful for example for simulating isotope mixes. Based on the ion symbols and fractional abundances, the average mass of each species is determined. The average charge state of each ion in each mixture is determined by Mavrin polynomial fits to ADAS atomic data |mavrin|. For the temperature ranges of interest in the tokamak core, these are well approximated as 1D functions of electron temperature. All ions with atomic numbers below Carbon are assumed to be fully ionized.
The impurity and main ion densities are constrained by the plasma effective charge, \(Z_\mathrm{eff}\), which is a user-provided 2D array in both time and space, as well as quasineutrality.
\(n_i\), and \(n_{imp}\), are solved from the following system of equations, where \(Z_\mathrm{eff}\) and the electron density are known, and \(Z_\mathrm{imp}\) is the average impurity charge of the impurity mixture, with the average charge state for each ion determined from the Mavrin polynomials.
Initial conditions for the evolving profiles \(T_i\), \(T_e\), \(n_e\), and \(\psi\) are user-configurable.
Additionally, if any of the \(T_{i,e}\), \(n_e\) equations are set to be non-evolving (i.e., not evolved by the PDE solver), then time-dependent prescribed profiles are user-configurable.
For the poloidal flux \(\psi(\hat{\rho})\), the user specifies if the initial condition is based on a prescribed current profile, \(j_\mathrm{tor}=j_0(1-\hat{\rho}^2)^\nu\) (with \(j_0\) scaled to match \(I_p\), and \(\nu\) is user-configurable), or from the \(\psi\) provided in a geometry file. The prescribed current profile option is always used for the circular geometry. The total current \(I_p\) can be user-configured or determined by the geometry data. In the latter case, the \(\psi\) provided by the geometry data can still be used, but is scaled by the ratio of the desired \(I_p\) and the geometry-provided \(I_p\).
Transport models
Turbulent transport determines the values of the transport coefficients (\(\chi_i\), \(\chi_e\), \(D_e\), \(V_e\)) in the TORAX equation summary, and is a key ingredient in core transport simulations. Theory-based turbulent transport models provide the largest source of nonlinearity in the PDE system. TORAX currently offers five transport models:
Constant: This simple model sets all transport coefficients to constant, user-configurable values. While not physically realistic, it can be useful for testing purposes.
CGM: The critical gradient model (CGM) is a simple theory-based model, capturing the basic feature of tokamak turbulent transport, critical temperature gradient transport thresholds. The model is a simple way to introduce transport coefficient nonlinearity, and is mainly used for testing purposes.
\[\begin{split}\chi_i = \begin{cases} \chi_{GB} \text{C} (\frac{R}{L_{Ti}} - \frac{R}{L_{Ti,\textit{crit}}})^{\\alpha} & \text{if } \frac{R}{L_{Ti}} \ge \frac{R}{L_{Ti,\textit{crit}}} \\ \chi_{min} & \text{if } \frac{R}{L_{Ti}} < \frac{R}{L_{Ti,\textit{crit}}} \end{cases}\end{split}\]with the GyroBohm scaling factor
\[\chi_{GB} = \frac{(A_i m_p)^{1/2}}{(eB_0)^2} \frac{(T_i k_B)^{3/2}}{R_\textit{maj}}\]and the |guo-romanelli| ion-temperature-gradient (ITG) mode critical gradient formula.
\[R/L_{Ti,crit} = \frac{4}{3}(1 + T_i/T_e)(1 + 2|\hat{s}|/q)\]where \(\chi_\textit{min}\) is a user-configurable minimum allowed \(\chi\), \(L_{Ti}\equiv-\frac{T_i}{\nabla T_i}\) is the ion temperature gradient length, \(A_i\) is the main ion atomic mass number, \(m_p\) the proton mass, \(e\) the electron charge, \(B_0\) the magnetic field on axis, and \(R_\mathrm{maj}\) the major radius. The stiffness factor \(C\) and the exponent \(\alpha\) are user-configurable model parameters.
Regarding additional transport coefficient outputs, the electron heat conductivity, \(\chi_e\), and particle diffusivity, \(D_e\), are scaled to \(\chi_i\) using user-configurable model parameters. The particle convection velocity, \(V_e\), is user-defined.
Bohm-GyroBohm: A widely used semi-empirical model summing terms proportional to Bohm and gyro-Bohm scaling factors |bohm-gyrobohm|.
The heat diffusivities for electrons and ions are given by:
\[\chi_e = \alpha_{e, \text{B}} \chi_{e, \text{B}} + \alpha_{e, \text{gB}} \chi_{e, \text{gB}}\]\[\chi_i = \alpha_{i, \text{B}} \chi_{i, \text{B}} + \alpha_{i, \text{gB}} \chi_{i, \text{gB}}\]where \(\alpha_{s, \text{B}}\) and \(\alpha_{s, \text{gB}}\) are the coefficients for the Bohm and gyro-Bohm contribution for species \(s\) respectively. These are given by:
\[\chi_{e, \text{B}} = 0.5 \chi_{i, \text{B}} = \frac{R_\text{min} q^2}{e B_\text{0} n_e} \left| \frac{\partial p_e}{\partial \rho_{\text{tor}}} \right|\]\[\chi_{e, \text{gB}} = 2 \chi_{i, \text{gB}} = \frac{\sqrt{T_e}}{B_\text{0}^2} \left| \frac{\partial T_e}{\partial \rho_{\text{tor}}} \right|\]where \(R_\text{min}\) is the minor radius, \(q\) is the safety factor, \(e\) is the elementary charge, \(B_\text{0}\) is the toroidal magnetic field at the magnetic axis, \(n_e\) is the electron density, \(\rho_{\text{tor}}\) is the (unnormalized) toroidal flux coordinate, \(p_e\) is the electron pressure, and \(T_e\) is the electron temperature.
The electron diffusivity is given by |garzotti2003|:
\[D_e = \eta \frac{\chi_e \chi_i}{\chi_e + \chi_i}\]where \(\eta\) is a weighting factor given by:
\[\eta = c_1 + (c_2 - c_1) \rho_{\text{tor}}\]where \(c_1\) and \(c_2\) are user-defined parameters.
There is little discussion in the literature about setting the electron convectivity from the Bohm/gyro-Bohm model. Following RAPTOR’s
vpdn_chiescalmethod, in TORAX, we set the electron convectivity proportional to the diffusivity,\[v_e = c_v D_e\]where \(c_v\) is a user-defined parameter.
The default values for the model parameters are as follows:
\(\alpha_{e,i,\text{B}} = 8e^{-5}\)
\(\alpha_{e,i,\text{gB}} = 5e^{-6}\)
\(c_1 = 1.0\)
\(c_2 = 0.3\)
\(c_v = -0.1\)
Please note that the Bohm-GyroBohm model TORAX implementation is presently experimental and subject to ongoing verification against established simulations.
QLKNN: This is a ML-surrogate model trained on a large dataset of the |qualikiz| quasilinear gyrokinetic code. TORAX presently employs the recently released open source |qlknn_7_11| model, based on QuaLiKiz v2.8.1 data. Its training data contains a combination of a 11D input hypercube for the plasma bulk region, and a 7D input hypercube focusing on the L-mode edge region near the LCFS. The model provides separate ouputs for ion-temperature-gradient (ITG), trapped-electron-mode (TEM), and electron-temperature-gradient (ETG) mode turbulent fluxes. The NNs take as input local plasma parameters, such as normalized gradients of temperature and density (both electron and impurity density gradients), temperature ratios, safety factor, magnetic shear, main ion dilution, and normalized collisionality, and outputs turbulent fluxes for ion and electron heat and particle transport. The QLKNN model is significantly faster than direct gyrokinetic simulations, enabling fast and accurate simulation within its range of validity. The ability to seamlessly couple ML-surrogate models is a key TORAX feature. TORAX is also coupled to the QLKNN-hyper-10D model |qlknn10d|, including dedicated JAX inference code written in |flax|.
QuaLiKiz: TORAX can be configured to use the |qualikiz| quasilinear gyrokinetic transport model itself. Since QuaLiKiz is an external code (written in Fortran), both |qualikiz| and its associated |qualikiz-pythontools| must be installed separately. The path to the QuaLiKiz executable must be set in a
TORAX_QLK_EXEC_PATHenvironment variable. If this environment variable is not set, then the default is~/qualikiz/QuaLiKiz. See above links for installation instructions. QuaLiKiz and TORAX exchange data via file I/O on a temporary directory. Since transport model calls are ostensibly carried out within JAX-compiled functions, running with QuaLiKiz requires disabling JAX compilation by settingTORAX_COMPILATION_ENABLED=False. While other solutions exist, this is the simplest and most straightforward approach. In any case QuaLiKiz becomes the simulation bottleneck and limits the overall simulation speed regardless of JAX compilation in the rest of the stack. Furthermore, QuaLiKiz must be run with thelinearsolver, ideally with theuse_predictor_correctorsolver option, sincenewton_raphsonrequires autodiff which is not supported for QuaLiKiz. Running with QuaLiKiz is not a typical workflow due to its computational expense (2-3 orders of magnitude slower than with QLKNN). Its use-cases are:
Evaluating ML-surrogates against their ground truth, i.e. for QLKNN.
Carrying out higher-fidelity simulations to verify faster simulations carried out with ML-surrogates. It is a major advantage to be able to use TORAX as the same framework for both ML-surrogate and high-fidelity simulations.
For all transport models, optional spatial smoothing of the transport coefficients using a Gaussian convolution kernel is implemented, to improve solver convergence rates, an issue which can arise with stiff transport coefficients such as from QLKNN. Furthermore, for all transport models, the user can set inner (towards the center) and/or outer (towards the edge) radial zones where the transport coefficients are prescribed to fixed values.
An edge-transport-barrier, or pedestal, is set up in TORAX through an adaptive source term which sets a desired value (pedestal height) of \(T_e\), \(T_i\) and \(n_e\), at a user-configurable location (pedestal width). Two different variants are available, one setting the pedestal pressure and temperature ratios, and the other setting the pedestal temperatures directly.
In the TORAX roadmap, coupling to additional transport models is envisaged, including to additional semi-empirical models and H-mode confinement scaling law adaptive models, as well as more ML-surrogates of theory-based models, both for core turbulence and pedestal predictions. A more physically consistent approach for setting up the pedestal will be implemented by incorporating adaptive transport coefficients in the pedestal region, as opposed to an adaptive local source/sink term.
Neoclassical physics
TORAX employs the Sauter model |sauter99| to calculate the bootstrap current density, \(j_{bs}\), and the neoclassical conductivity, \(\sigma_{||}\), used in the current diffusion equation. The Sauter model is a widely-used analytical formulation that provides a relatively fast and differentiable approximation for these neoclassical quantities.
Future work can incorporate more recent neoclassical physics parameterizations, and also set neoclassical transport coefficients themselves. This can be of importance for ion heat transport in the inner core. When extending TORAX to include impurity transport, incorporating fast analytical neoclassical models for heavy impurity transport will be of great importance.
Rotation Physics
The radial electric field (\(E_r\)), which drives \(E \times B\) plasma rotation, is a crucial factor in turbulent transport, as \(E \times B\) shear can suppress turbulence.
The radial electric field \(E_r\) is determined from the radial force balance equation for the main ions:
where \(P_i\) is the main ion pressure, \(n_i\) is the main ion density, \(Z_i\) is the main ion charge number, \(e\) is the elementary charge, \(v_{\phi}\) is the toroidal rotation velocity, \(v_{\theta}\) is the poloidal rotation velocity, \(B_{\theta}\) is the poloidal magnetic field, and \(B_{\phi}\) is the toroidal magnetic field. The derivatives are with respect to a midplane-averaged radial coordinate.
The poloidal velocity \(v_{\theta}\) is calculated using neoclassical formulas. Specifically, it implements Equation 33 from |kim1991|. The formula used is:
where \(k_{neo}\) is a neoclassical coefficient, \(dT_i/dr\) is the radial gradient of the ion temperature, and \(B_{total}\) is the total magnetic field.
The neoclassical coefficient \(k_{neo}\) is based on Equation (6.135) from |hinton1976|. This coefficient depends on the normalized ion collisionality (\(\nu_{i}^{*}\)) and the inverse aspect ratio (\(\epsilon\)). The limits of this formula are approximately 1.17 in the banana regime (\(\nu_{i}^{*} \rightarrow 0\)) and -2.1 in the Pfirsch-Schluter regime (\(\nu_{i}^{*} \rightarrow \infty\)).
The normalized ion collisionality \(\nu_{i}^{*}\) is calculated based on |sauter99|, and depends on the safety factor (\(q\)), geometry, ion density (\(n_i\)), ion temperature (\(T_i\)), effective charge number (\(Z_{eff}\)), and the ion-ion Coulomb logarithm (\(\log \Lambda_{ii}\)).
The \(E \times B\) velocity (\(v_{E \times B}\)) is derived from the radial electric field \(E_r\) and the total magnetic field as follows:
Rotation effects are currently disabled by default in transport models. They can be enabled through the transport model configuration.
Rotation in Transport Models:
TGLFNN-ukaea: The TGLFNN-ukaea model includes the \(E \times B\) shearing rate as an input feature, directly informing the model’s turbulent transport predictions. To enable this, set the use_rotation parameter to True in the transport model configuration.
QLKNN: The QLKNN transport model incorporates a “rotation rule” that modifies turbulent fluxes (specifically for ITG and TEM modes) based on E×B shear. The modification combines two physical effects:
Waltz rule applied to poloidal/pressure contributions (from \(\nabla P_i\) and \(v_\theta B_\phi\)): Simple suppression \(f_{waltz} = -\alpha\), where \(\alpha\) is set by the
shear_suppression_alphaparameter (default 1.0). See |waltz1998|.Victor rule applied to toroidal contributions (from \(-v_\phi B_\theta\)): A fitted model that can produce both suppression and enhancement depending on local plasma parameters, due to competing stabilizing (E×B shear) and destabilizing (parallel velocity shear) effects. See |qlknn10d|.
The combined scaling factor is:
\[f_{rot} = \text{clip}\left(1 + f_{waltz}\frac{|\gamma_{pp}|}{\gamma_{max}} + f_{victor}\frac{|\gamma_{tor}|}{\gamma_{max}}, 0\right)\]where \(\gamma_{pp}\) and \(\gamma_{tor}\) are the E×B shearing rates from poloidal/pressure and toroidal contributions respectively.
The application of the rotation rule is controlled by the
rotation_modeconfiguration parameter. Options forrotation_modeare:off: No rotation correction is applied.half_radius: The rotation correction is only applied to the outer half of the radius (\(\hat{\rho} > 0.5\)).full_radius: The rotation correction is applied everywhere.
Tuning Rotation Effects:
Two parameters are available to fine-tune the impact of rotation, both of them default to 1.0:
rotation_multiplier: Located in the transport model configs, this parameter scales the \(E \times B\) shear term.poloidal_velocity_multiplier: Found under theneoclassicalconfiguration, this parameter directly scales the poloidal velocity term.
Edge models
TORAX supports coupling to reduced edge/divertor models to provide physically consistent boundary conditions for the core transport solver. Currently, only the Extended Lengyel model is supported.
Extended Lengyel Model
The Extended Lengyel model describes the 1D parallel heat and particle transport in the Scrape-Off Layer (SOL) and divertor. It extends the classic Lengyel model by including cross-field transport in the divertor, power and momentum loss due to neutral ionization close to the divertor target and turbulent broadening of the upstream heat flux channel. The implementation follows |body2025|.
The model relates upstream quantities (e.g. power crossing separatrix \(P_{SOL}\), separatrix density \(n_{e,sep}\)) and divertor impurity content, to downstream quantities such as target temperature \(T_{e,t}\). The separatrix temperature is calculated from a 2-point model.
Modes of Operation:
Forward Mode: Given the impurity mix and upstream conditions, calculate the target temperature \(T_{e,t}\).
Inverse Mode: Given a desired target temperature \(T_{e,t}\) (e.g., 5 eV to represent detachment onset), calculate the required concentration of a given mix of seeded impurity species.
Physics Features:
Impurity Radiation: Radiative cooling rates \(L_z(T_e)\) are calculated using polynomial fits from |mavrin2017|, which are valid at low edge temperatures (down to ~1 eV).
Enrichment: Impurity enrichment (\(c_{div}/c_{core}\)) can be specified manually or calculated using the empirical scaling from |kallenbach2024|: \(E \propto Z^{-0.5} p_0^{-0.4} (E_{ion,z}/E_{ion,D})^{-5.8}\).
Coupling to Core:
When enabled, the edge model updates the following TORAX runtime parameters for the next time step. This is done via explicit coupling: the parameters at time \(t\) are passed to the edge model, which calculates the updated parameters at time \(t+\Delta t\).
Boundary Conditions: The separatrix electron temperature calculated by the model sets the core \(T_e\) boundary condition. \(T_i\) boundary condition is set via a configured ratio.
Impurity Composition:
In Inverse Mode, the required seeded impurity concentration updates the core impurity profile (scaled by enrichment factor).
The model also enforces consistency for fixed background impurities between core and edge.
Sources
The source terms in the TORAX equation summary are comprised of a summation of individual source/sink terms. Each of these terms can be configured to be either:
Implicit: Where needed in the theta-method, the source term is calculated based on the current guess for the state at \(t+\Delta t\).
Explicit: The source term is always calculated based on the state of the system at the beginning of the timestep, even if the solver \(\theta>0\). This makes the PDE system less nonlinear, avoids the incorporation of the source in the residual Jacobian if solving with Newton-Raphson, and leads to a single source calculation per timestep.
Explicit treatment is less accurate, but can be justified and computationally beneficial for sources with complex but slow-evolving physics. Furthermore, explicit source calculations do not need to be JAX-compatible, since explicit sources are an input into the PDE solver, and do not require JIT compilation. Conversely, implicit treatment can be important for accurately resolving the impact of fast-evolving source terms.
All sources can optionally be set to zero, prescribed with explicit values or calculated with a dedicated physics-based model. However, the code modular structure facilitates easy coupling of additional source models in future work. Specifics of source models currently implemented in TORAX follow:
Ion-electron heat exchange
The collisional heat exchange power density is calculated as
where \(A_i\) is the atomic mass number of the main ion species, \(m_p\) is the proton mass, and \(\tau_e\) is the electron collision time, given by:
where \(\epsilon_0\) is the permittivity of free space, \(m_e\) is the electron mass, \(e\) is the elementary charge, and \(\ln \Lambda_{ei}\) is the Coulomb logarithm for electron-ion collisions given by:
\(Q_{ei}\) is added to the electron heat sources, meaning that positive \(Q_{ei}\) with \(T_i>T_e\) heats the electrons. Conversely, \(-Q_{ei}\) is added to the ion heat sources.
Fusion power
TORAX optionally calculates the fusion power density assuming a 50-50 deuterium-tritium (D-T) fuel mixture using the |bosch-hale| parameterization for the D-T fusion reactivity \(\langle \sigma v \rangle\):
where \(E_{fus} = 17.6\) MeV is the energy released per fusion reaction, \(n_i\) is the ion density, and \(\langle \sigma v \rangle\) is given by:
with
and
where \(T_i\) is the ion temperature in keV, \(m_rc^2\) is the reduced mass of the D-T pair. The values of \(m_rc^2\), the Gamov constant \(B_G\), and the constants \(C_1\) through \(C_7\) are provided in the Bosch-Hale paper.
TORAX partitions the fusion power between ions and electrons using the parameterized fusion particle slowing down model of Mikkelsen, which neglects the slowing down time itself.
Ohmic power
The Ohmic power density, \(P_\mathrm{ohm}\), arising from resistive dissipation of the plasma current, is calculated as:
where \(j_\mathrm{tor}\) is the toroidal current density, and \(R_\mathrm{maj}\) is the major radius. The loop voltage \(\frac{\partial \psi}{\partial t}\) is calculated according to the \(\psi\) equation in the TORAX equation summary. \(P_\mathrm{ohm}\) is then included as a source term in the electron heat transport equation.
Auxiliary Heating and Current Drive
For prescribing generic non-physics-based auxiliary heating and current drive sources, TORAX provides built-in Gaussian formulations of a generic ion and electron heat source, and a generic current drive source, with time-dependent user configurable locations, Gaussian width, amplitude, and fractional heating of ions and electrons.
More sophisticated physics-based models and/or ML-surrogates of specific auxiliary heating and current drive sources can be coupled modularly to TORAX, enhancing the fidelity of the simulation. By setting these as explicit sources, these can also come from external codes (not necessarily JAX compatible) coupled to TORAX in larger workflows.
Available physics-based models and/or ML-surrogates are listed below.
Electron-Cyclotron Heating and Current Drive
The electron-cyclotron current drive can be calculated from the heating power density, \(Q_\mathrm{EC}(\rho) [Wm^{-3}]\), and a dimensionless EC current drive efficiency profile, \(\zeta_\mathrm{EC}(\rho)\) . The current drive parallel to the magnetic field, in \([Am^{-2}]\), is then given by:
where \(\epsilon_0\) is the vacuum permittivity, \(F = B_\phi R\) is the flux function, \(e\) is the elementary charge, \(R_\mathrm{maj}\) is the device major radius, \(T_e\) is the electron temperature in joules, and \(n_e\) is the electron density per cubic meter. This relationship is based on the Lin-Liu model |lin-liu|. The derivation can be found here.
Particle Sources
Similar to auxiliary heating and current drive, particle sources can also be configured using either prescribed formulas. Presently, TORAX provides three built-in formula-based particle sources for the \(n_e\) equation:
Gas Puff: An exponential function with configurable parameters models the ionization of neutral gas injected from the plasma edge.
Pellet Injection: A Gaussian function approximates the deposition of particles from pellets injected into the plasma core. The time-dependent configuration parameter feature allows either a continuous approximation or discrete pellets to be modelled.
Generic particle source: An additional Gaussian function which can be configured to model arbitrary particle sources, e.g. to mock-up an NBI source.
Radiation
Bremsstrahlung
Uses the model from Wesson, John, and David J. Campbell. Tokamaks. Vol. 149.
An optional correction for relativistic effects from Stott PPCF 2005 can be
enabled with the flag use_relativistic_correction.
When the Mavrin impurity radiation model is also active, the bremsstrahlung source automatically excludes the impurity bremsstrahlung component (using only the main-ion contribution to \(Z_\text{eff}\)) to avoid double-counting, since the Mavrin model already accounts for impurity bremsstrahlung via higher-fidelity fits to ADAS data.
Cyclotron Radiation
Uses the total radiation power from |albajar2001|, with a deposition profile from |artaud2018|. The Albajar model includes a parameterization of the temperature profile which in TORAX is fit by simple grid search for computational efficiency.
Impurity Radiation
The following models are available:
Set the impurity radiation to be a constant fraction of the total external input power.
Polynomial fits to ADAS data from Mavrin |mavrin|. Provides radiative cooling rates for the following impurity species:
Helium
Lithium
Beryllium
Carbon
Nitrogen
Oxygen
Neon
Argon
Krypton
Xenon
Tungsten
These cooling curves are multiplied by the electron density and impurity densities to obtain the impurity radiation power density.
The valid temperature range of the fit is [0.1-100] keV.
Ion Cyclotron Resonance Heating
Presently this source is implemented for a SPARC specific ICRH scenario.
A core Ion Cyclotron Range of Frequencies (ICRF) heating surrogate model trained on TORIC ICRH spectrum solver simulations is used to provide power profiles for Helium-3, Tritium (via its second harmonic) and electrons |wallace2024|.
A “Stix distribution” [Stix, Nuc. Fus. 1975] is used to model the non-thermal Helium-3 distribution based on an analytic solution to the Fokker-Planck equation to estimate the birth energy of Helium-3.
TORAX partitions the Helium-3 power between ions and electrons using the parameterized model of Mikkelsen, as for Fusion Power.
It is assumed that all tritium heating goes to ions and all electron heating goes to electrons.
MHD models
Currently only a sawtooth model is implemented, although the TORAX config and internal APIs are designed to accommodate other models in the future.
Sawtooth Model
Sawteeth are periodic oscillations in the core plasma temperature, density, and current caused by the growth and subsequent rapid reconnection of an m=1, n=1 kink instability within the plasma volume where the safety factor, \(q\), drops below unity. The sawtooth crash is triggered by a state-dependent critical magnetic shear at the \(q=1\) surface.
The TORAX sawtooth model comprises two components:
Trigger Model: determines the conditions under which a sawtooth crash is initiated.
Redistribution Model: Modifies the plasma profiles (temperature, density, poloidal flux) following a crash to simulate the rapid transport event.
Currently only simple Trigger and Redistribution models are implemented.
The simple trigger model checks for the following conditions at each time
step:
Existence of a q=1 surface: The safety factor profile q must drop below 1.
Minimum radius: The normalized radius of the q=1 surface (\(\hat{\rho}_{q=1}\)) must be greater than a specified minimum value (
minimum_radius). This prevents spurious triggers very close to the magnetic axis.Critical magnetic shear: The magnetic shear (\(\hat{s}\)) at the \(q=1\) surface must exceed a critical value (
s_critical).
The simple redistribution model mplements a simplified redistribution
process based on conserving particle number, energy, and current within a
mixing radius.
Mixing Radius Calculation: The mixing radius (\(\hat{\rho}_{mix}\)) is calculated as
mixing_radius_multiplier* \(\hat{\rho}_{q=1}\).Profile Flattening: Profiles (\(T_i\), \(T_e\), \(n_e\), \(j_{tot}\)) within the mixing radius are partially flattened.
Inside the q=1 surface (\(\hat{\rho} < \hat{\rho}_{q=1}\)), a linear profile is created, where the value at \(\hat{\rho}=0\) is set to a multiple
flattening_factorof the value at \(\hat{\rho}_{q=1}\).Between the \(q=1\) surface and the mixing radius (\(\hat{\rho}_{q=1} \le \hat{\rho} < \hat{\rho}_{mix}\)), the profile is linearly interpolated between the flattened value at \(\hat{\rho}_{q=1}\) and the original profile value at \(\hat{\rho}_{mix}\).
The values at \(\hat{\rho}_{q=1}\) are set by conservation laws.
Conservation: The flattening process conserves the total number of particles (\(\int n_e dV\)), total electron and ion thermal energy (\(\int \frac{3}{2} n_{e,i} T_{e,i} dV\)), and total current (\(\int j_{tot} dS\)) within the mixing radius. The values at \(\hat{\rho}_{q=1}\) after flattening are adjusted to ensure these conservation laws are met. The poloidal flux profile is then recalculated to be consistent with the flattened current profile, and the pre-crash poloidal flux boundary condition.
Derived Quantities: After redistribution, derived quantities like ion densities (\(n_i\), \(n_{imp}\)), impurity charge states (\(Z_{imp}\)), safety factor (\(q\)), magnetic shear (\(\hat{s}\)), and sources, are recalculated based on the modified profiles.
Simulation Step During Crash
When a sawtooth crash is triggered:
The standard PDE time step is skipped.
The redistribution model is applied to modify the
core_profiles.A short, fixed time step (user-configurable
crash_step_duration) is taken.Derived quantities are recalculated based on the redistributed profiles.
The simulation then always proceeds to the next regular PDE time step. Subsequent sawtooth crashes are not allowed, to prevent continuous consecutive crashes if the trigger condition is still met immediately after redistribution.
See the sawtooth configuration details for a summary of the input parameters.
The current sawtooth model is a simplified representation. Future development plans include: Implementing more sophisticated trigger models (e.g., based on the Porcelli model). Developing more physically detailed redistribution models (e.g., Kadomtsev reconnection models).