The Stellarator Reduced MHD Model

The JOREK stellarator model (reduced MHD) features two distinct models:

  1. Model 180 does not evolve the plasma in time but only calculates the initial conditions from VMEC/GVEC.
  2. Model 183 actually evolves the plasma in time (keeping the same hard-coded compiling parameters as model 180).

The models include some hard-coded switches, such as the inclusion of the parallel velocity $v_{\parallel}$.

A detailed description of the original JOREK stellarator model implementations can be found in the following literature:

  1. N. Nikulsin. Models and methods for nonlinear magnetohydrodynamic simulations of stellarators (2021), PhD thesis, TUM.
  2. R. Ramasamy Equilibrium and initial value problem simulation studies of nonlinear magnetohydrodynamics in stellarators (2022), PhD thesis, TUM. Here, particulary Ch. 7 is important.
  3. N. Nikulsin et al. JOREK3D: An extension of the JOREK nonlinear MHD code to stellarators, Physics of Plasmas 29.6 (2022).
  4. R. Ramasamy et al. Nonlinear MHD modeling of soft beta limits in W7-AS, Nucl. Fusion 64 086030 (2024).

The stellarator model equations

Let us split the magnetic field into the vacuum part $\nabla\chi$ and the part that is generated by currents in the plasma $\vec B_p$: \(\vec B = \nabla\chi + \vec B_p.\) Then, Clebsch coordinates $(\psi_v,\beta_v,\chi)$ can be constructed for the vacuum field, and by definition must satisfy \(\nabla\chi = \nabla\psi_v\times\nabla\beta_v.\) Now, if we consider just the plasma-induced field $\vec B_p$, then its vector potential can be written as: \(\vec A_p = A_\chi\nabla\chi + A_{\psi_v}\nabla\psi_v,\) where gauge freedom was used to eliminate the $\nabla\beta_v$ component. Thus, letting $A_\chi = \Psi$ and $A_{\psi_v} = \Omega$, the total magnetic field is \(\vec B = \nabla\chi + \nabla\Psi\times\nabla\chi + \nabla\Omega\times\nabla\psi_v.\)

As shown in [1], any velocity field can be written as: \(\newcommand{\vpar}{v_\parallel} \newcommand{\pgrad}{\nabla^\perp} \vec v = \frac{\nabla\Phi\times\nabla\chi}{B_v^2} + \vpar\vec B + \pgrad\zeta,\) where $B_v = \Vert\nabla\chi\Vert$ and $\nabla^\perp = \nabla - B_v^{-2}\nabla\chi\nabla\chi\cdot\nabla$.

After inserting these ansatzes into the resistive MHD equations, as shown on the Base Fluid Models page, we project Faraday's law on $\nabla\psi_v$ and $\nabla\chi$, and apply the following projection operators to the momentum equation: \(\begin{gathered} \nabla\chi\cdot\nabla\times(B_v^{-2} \\ \vec B\cdot \\ \nabla\chi\cdot\nabla\times(B_v^{-2}\nabla\times \end{gathered}\) This results in the following full MHD equations: \(\newcommand{\tderiv}[1]{\frac{\partial#1}{\partial t}} \newcommand{\llderiv}{\partial^\parallel} \newcommand{\psderiv}{\partial_{\psi_v}} \newcommand{\pLap}{\Delta^\perp} \begin{aligned} &\left[\psi_v,\tderiv{\Psi}\right] = \left[\frac{[\Psi,\Phi]-\llderiv\Phi}{B_v},\psi_v\right] - \frac{F_v}{B_v}\left[\Omega,\frac{[\psi_v,\Phi]}{B_v}\right]_{\psi_v} + \llderiv(F_v\psderiv\zeta) + \left[(\zeta,\Psi),\psi_v\right] - \frac{F_v}{B_v}\left[\Omega,(\zeta,\psi_v)\right]_{\psi_v} + \frac{1}{B_v}\nabla\cdot(\eta\nabla\psi_v\times\vec j), \\ \\ &\left[\tderiv{\Omega},\psi_v\right] = -\left[\Omega,\frac{[\psi_v,\Phi]}{B_v}\right] + \left[\psi_v,\frac{[\Omega,\Phi]}{B_v}\right] - 2(B_v,\zeta) - B_v\pLap\zeta - \left[\Omega,(\zeta,\psi_v)\right] + \left[\psi_v,(\zeta,\Omega)\right] + \frac{1}{B_v}\nabla\cdot(\eta\nabla\chi\times\vec j), \\ \\ &\nabla\cdot\left[\frac{\rho}{B_v^2}\pgrad\tderiv{\Phi} + \rho\tderiv{}(\vpar\pgrad\Psi) - \frac{\rho}{B_v}\tderiv{}(\vpar\llderiv\Omega)\nabla\psi_v\right] - B_v\left[\frac{\rho}{B_v^2},\tderiv{\zeta}\right] = \frac{B_v}{2}\left[\frac{\rho}{B_v^2},v^2\right] + B_v\left[\frac{\rho\omega^\chi}{B_v^4},\Phi\right] \\ &- \nabla\cdot(\rho \vpar\vec\omega^\perp) + B_v\left[\frac{\rho \vpar\omega^\chi}{B_v^2},\Psi\right] + B_v\left[\frac{\rho \vpar\omega^{\psi_v}}{B_v^2},\Omega\right] + \nabla\cdot\left(\frac{\rho\omega^\chi}{B_v^2}\pgrad\zeta - \frac{P}{B_v^2}\pgrad\Phi - P\vpar\pgrad\Psi + \frac{P\vpar\llderiv\Omega}{B_v}\nabla\psi_v\right) \\ &- B_v\left[\frac{\rho\vpar\vec\omega\cdot\nabla\Omega}{B_v^2},\psi_v\right] + B_v\left[\frac{P}{B_v^2},\zeta\right] + \nabla\cdot\left(\frac{B^\chi}{B_v^2}\vec j - \frac{j^\chi}{B_v^2}\vec B\right) + B_v\left[\frac{1}{B_v^2},p\right] + \nabla\cdot(\mu_\perp\pgrad\pLap\Phi), \\ \\ &\rho\left(\tderiv{\Phi},\Psi\right) - \rho\frac{F_v}{B_v}\llderiv\Omega\psderiv\tderiv{\Phi} + \rho B^2\tderiv{\vpar} + \frac{\rho\vpar}{2}\tderiv{B^2} + \rho B_v\left[\tderiv{\zeta},\Psi\right] + \rho F_v\left[\tderiv{\zeta},\Omega\right]_{\psi_v} = -\frac{\rho B_v}{2}\llderiv v^2 \\ &- \frac{\rho B_v}{2}\left[v^2,\Psi\right] - \frac{\rho F_v}{2}\left[v^2,\Omega\right]_{\psi_v} - \frac{\rho\omega^\chi}{B_v}\llderiv\Phi - \frac{\rho\omega^\chi}{B_v}\left[\Phi,\Psi\right] - \frac{\rho\omega^\chi F_v}{B_v^2}\left[\Phi,\Omega\right]_{\psi_v} + \frac{\rho B^\chi}{B_v^2}\vec\omega\cdot\nabla\Phi - \rho\vec\omega\cdot(\nabla\zeta\times\nabla\chi) \\ &+ \rho\omega^\chi(\Psi,\zeta) - \rho\vec\omega\cdot\nabla\Omega F_v\psderiv\zeta + \rho\omega^{\psi_v}(\Omega,\zeta) - \vec v\cdot\vec B~P - B_v\llderiv p - B_v\left[p,\Psi\right] - F_v\left[p,\Omega\right]_{\psi_v} + \nabla\cdot(\mu_\parallel\nabla_\perp\vpar + \mu_{\parallel,\parallel}\nabla_\parallel\vpar), \\ \\ &B_v\left[\frac{\rho}{B_v^2},\tderiv{\Phi}\right] + B_v\left[\rho\tderiv{\vpar},\Psi\right] + B_v\left[\rho \vpar,\tderiv{\Psi}\right] - B_v\left[\frac{\rho}{B_v}\tderiv{}(\vpar\llderiv\Omega),\psi_v\right] + \nabla\cdot\left(\rho\pgrad\tderiv{\zeta}\right) \\ &= -\nabla\cdot\left(\frac{\rho}{2}\pgrad v^2 + \frac{\rho\omega^\chi}{B_v^2}\pgrad\Phi - \frac{\rho \vpar B^\chi}{B_v^2}\nabla\chi\times\vec\omega + \frac{\rho \vpar\omega^\chi}{B_v^2}\nabla\chi\times\vec B\right) + B_v\left[\frac{\rho\omega^\chi}{B_v^2},\zeta\right] - B_v\left[\frac{P}{B_v^2},\Phi\right] \\ &- B_v\left[P\vpar,\Psi\right] + B_v\left[\frac{P\vpar\llderiv\Omega}{B_v},\psi_v\right] - \nabla\cdot\left(P\pgrad\zeta + \frac{B^\chi}{B_v^2}\nabla\chi\times\vec j - \frac{j^\chi}{B_v^2}\nabla\chi\times\vec B\right) - \pLap p + \nabla\cdot(\mu_\zeta\pgrad\pLap\zeta), \\ \\ &\tderiv{\rho} = - B_v\left[\frac{\rho}{B_v^2},\Phi\right] - B_v\llderiv(\rho \vpar) - B_v[\rho \vpar,\Psi] - F_v[\rho \vpar,\Omega]_{\psi_v} - \nabla\cdot(\rho\pgrad\zeta) + P, \\ \\ &\tderiv{p} = - \frac{1}{B_v}\left[p,\Phi\right] - \vpar B_v\llderiv p - \vpar B_v\left[p,\Psi\right] - \vpar F_v \left[p,\Omega\right]_{\psi_v} - (\zeta,p) - \gamma p B_v\left[\frac{1}{B_v^2},\Phi\right] - \gamma p B_v\llderiv \vpar - \gamma p B_v\left[\vpar,\Psi\right] \\ &- \gamma p F_v\left[\vpar,\Omega\right]_{\psi_v} - \gamma p\pLap\zeta + \nabla\cdot\vec q + (\gamma-1)(S_e + \eta j^2), \end{aligned}\) where $[f,g] = \nabla\chi\cdot(\nabla f\times\nabla g)/B_v$ and $[f,g]{\psi_v} = \nabla\psi_v\cdot(\nabla f\times\nabla g)/\Vert\nabla\psi_v\Vert$ are Poisson brackets, $(f,g) = \pgrad f\cdot\pgrad g$ is an inner product, $\nabla\parallel = B^{-2}\vec B\vec B\cdot\nabla$ and $\nabla_\perp = \nabla - \nabla_\parallel$. No further assumptions have been made in this derivation beyond resistive MHD being valid, but we have added generic viscosity terms to the projected momentum equations. These generic viscosity terms do not accurately model viscous effects in a plasma, but they are often useful in avoiding numerical instabilities.

Reduced MHD

The equations that are actually implemented in JOREK are the reduced MHD equations, and one further step is required to obtain them from the full MHD equations above. Namely, we set $\Omega = 0$ and $\zeta = 0$ and drop the evolution equations for $\Omega$ and $\zeta$ from the full set of equations above. To conserve energy, we must also drop the components of $\vec j$ perpendicular to $\nabla\chi$ from the $\Psi$ evolution equation. The resulting reduced MHD equations are as follows: \(\newcommand{\red}[1]{{\color{red}#1}} \newcommand{\blue}[1]{{\color{blue}#1}} \begin{aligned} B_v\tderiv{\Psi} &= \llderiv\Phi - \left[\Psi,\Phi\right] - \eta\frac{j^\chi}{B_v} \blue{+ \eta^{(h)} B_v\pLap\left(\frac{j^\chi}{B_v^2}\right)}, \\ \nabla\cdot\left[\frac{\rho}{B_v^2}\pgrad\tderiv{\Phi} \red{+ \rho\tderiv{}(\vpar\pgrad\Psi)}\right] &= \frac{B_v}{2}\left[\frac{\rho}{B_v^2},v^2\right] + B_v\left[\frac{\rho\omega^\chi}{B_v^4},\Phi\right] - \nabla\cdot(\rho \vpar\vec\omega^\perp) + B_v\left[\frac{\rho \vpar\omega^\chi}{B_v^2},\Psi\right] \\ &- \nabla\cdot\left(\frac{P}{B_v^2}\pgrad\Phi \red{+ P\vpar\pgrad\Psi}\right) - \nabla\cdot\left(\frac{j^\chi}{B_v^2}\vec B\right) + B_v\left[\frac{1}{B_v^2},p\right] + \nabla\cdot(\mu_\perp\pgrad\pLap\Phi) \blue{- \mu_\perp^{(h)}\Delta^2\pLap\Phi}, \\ \red{\rho\left(\tderiv{\Phi},\Psi\right)} + \rho B^2\tderiv{\vpar} + \frac{\rho\vpar}{2}\tderiv{B^2} &= -\frac{\rho B_v}{2}\llderiv v^2 - \frac{\rho B_v}{2}\left[v^2,\Psi\right] \red{- \frac{\rho\omega^\chi}{B_v}\llderiv\Phi - \frac{\rho\omega^\chi}{B_v}\left[\Phi,\Psi\right] + \rho\vec\omega\cdot\nabla\Phi} - (\red{(\Phi,\Psi)} + \vpar B^2)P \\ &- B_v\llderiv p - B_v\left[p,\Psi\right] + \nabla\cdot(\mu_\parallel\nabla_\perp\vpar + \mu_{\parallel,\parallel}\nabla_\parallel\vpar) \blue{- \mu_\parallel^{(h)}\Delta^2\vpar}, \\ \tderiv{\rho} &= - B_v\left[\frac{\rho}{B_v^2},\Phi\right] - B_v\llderiv(\rho \vpar) - B_v[\rho \vpar,\Psi] + P, \\ \tderiv{p} &= - \frac{1}{B_v}\left[p,\Phi\right] - \vpar B_v\llderiv p - \vpar B_v\left[p,\Psi\right] - \gamma p B_v\left[\frac{1}{B_v^2},\Phi\right] - \gamma p B_v\llderiv \vpar - \gamma p B_v\left[\vpar,\Psi\right] \\ &+ \nabla\cdot\vec q + (\gamma-1)(S_e + \eta j^2) \blue{- \kappa_\perp^{(h)}\Delta^2 T}, \end{aligned}\) where $j^\chi = \nabla\chi\cdot\vec j = -\nabla\cdot(B_v^2\pgrad\Psi)/\mu_0$ and $\omega^\chi = \nabla\chi\cdot\vec\omega = -\pLap\Phi \red{- \nabla\cdot(B_v^2\vpar\pgrad\Psi)}$.

These equations, without the terms highlighted in red, are implemented in model 183. A further reduction is possible by setting $\vpar = 0$ and dropping the $\vpar$ evolution equation from the set of equations above. This is achieved by setting with_vpar = .false. in the hard-coded parameters before compiling JOREK. Conversely, when with_vpar = .true., all of the equations above are solved and no assumptions on $\vpar$ are made. Note that none of the red terms would have contributed anyway if with_vpar = .false..

We have also added hyperresistivity, hyperviscosity and hyperconductivity terms (highlighted in blue) to the original equations. These terms do not model any physical effects, but are sometimes necessary to suppress numerical instabilities. If these terms are used, one should make sure that they are small enough not to affect the physics.

The vacuum field

Dommaschk potentials

TO DO

Finite element Laplace solve

TO DO

The GVEC import

TO DO