Underlying Physics

Rayleigh solves the MHD equations in spherical geometry under the Boussinesq and anelastic approximations. This section will provide a basic overview of those equations as well as the mathematical approach Rayleigh uses to solve them.

Notation and Conventions

Vector and Tensor Notation

All vector quantities are represented in bold italics. Components of a vector are indicated in non-bold italics, along with a subscript indicating the direction associated with that component. Unit vectors are written in lower-case, bold math font and are indicated by the use of a hat character. For example, a vector quantity \(\boldsymbol{a}\) would represented as

(1)\[ \boldsymbol{a} = a_r\boldsymbol{\hat{a}}+a_\theta\boldsymbol{\hat{\theta}}+a_\phi\boldsymbol{\hat{\phi}}.\]

The symbols (\(\boldsymbol{\hat{r}}\), \(\boldsymbol{\hat{\theta}}\), \(\boldsymbol{\hat{\phi}}\)) indicate the unit vectors in the (\(r\),\(\theta\),\(\phi\)) directions, and (\(a_r\), \(a_\theta\), \(a_\phi\)) indicate the components of \(\boldsymbol{a}\) along those directions respectively.

Vectors may be written in lower case, as with the velocity field \(\boldsymbol{v}\), or in upper case as with the magnetic field \(\boldsymbol{B}\). Tensors are indicated by bold, upper-case, script font, as with the viscous stress tensor \(\boldsymbol{\mathcal{D}}\). Tensor components are indicated in non-bold, and with directional subscripts (i.e., \(\mathcal{D}_{r\theta}\)).

Reference-State Values

The hat notation is also used to indicate reference-state quantities. These quantities are scalar, and they are not written in bold font. They vary only in radius and have no \(\theta\)-dependence or \(\phi\)-dependence. The reference-state density is indicated by \(\hat{\rho}\) and the reference-state temperature by \(\hat{T}\), for instance.

Averaged and Fluctuating Values

Most of the output variables have been decomposed into a zonally-averaged value, and a fluctuation about that average. The average is indicated by an overbar, such that

(2)\[ \overline{a}\equiv \frac{1}{2\pi}\int_{0}^{2\pi} a(r,\theta,\phi)\, \mathrm{d}\phi.\]

Fluctations about that average are indicated by a prime superscript, such that

(3)\[ a'(r,\theta,\phi)\equiv a(r,\theta,\phi)-\overline{a}(r,\theta)\]

Finally, some quantities are averaged over the full sphere. These are indicated by a double-zero subscript (i.e. \(\ell=0,\,m=0\)), such that

(4)\[a_{00}\equiv \frac{1}{4\pi}\int_{0}^{2\pi}\int_{0}^{\pi} a(r,\theta,\phi)\, r\mathrm{sin}\,\theta\mathrm{d}\theta\mathrm{d}\phi.\]

The System of Equations Solved in Rayleigh

Rayleigh solves the Boussinesq or anelastic MHD equations in spherical geometry. Both the equations that Rayleigh solves and its diagnostics can be formulated either dimensionally or nondimensionally. A nondimensional Boussinesq formulation, as well as dimensional and nondimensional anelastic formulations (based on a polytropic reference state) are provided as part of Rayleigh. The user may employ alternative formulations via the custom Reference-state interface. To do so, they must specify the functions \(\mathrm{f}_i\) and the constants \(c_i\) in Equations (5)-(11) at input time (in development).

The general form of the momentum equation solved by Rayleigh is given by

(5)\[\begin{split} \mathrm{f}_1(r)\left[\frac{\partial \boldsymbol{v}}{\partial t} + \boldsymbol{v}\cdot\boldsymbol{\nabla}\boldsymbol{v} %advection + c_1\boldsymbol{\hat{z}}\times\boldsymbol{v} \right] =\ % Coriolis &c_2\,\mathrm{f}_2(r)\Theta\,\boldsymbol{\hat{r}} % buoyancy - c_3\,\mathrm{f}_1(r)\boldsymbol{\nabla}\left(\frac{P}{\mathrm{f}_1(r)}\right) % pressure \\ &+ c_4\left(\boldsymbol{\nabla}\times\boldsymbol{B}\right)\times\boldsymbol{B} % Lorentz Force + c_5\boldsymbol{\nabla}\cdot\boldsymbol{\mathcal{D}},\end{split}\]

where the stress tensor \(\mathcal{D}\) is given by

(6)\[ \mathcal{D}_{ij} = 2\mathrm{f}_1(r)\,\mathrm{f}_3(r)\left[e_{ij} - \frac{1}{3}\left(\boldsymbol{\nabla}\cdot\boldsymbol{v}\right)\delta_{ij}\right].\]

Here \(e_{ij}\) and \(\delta_{ij}\) refer to the standard rate-of-strain tensor and Kronecker delta, respectively.

The velocity field is denoted by \(\boldsymbol{v}\), the thermal anomoly by \(\Theta\), the pressure by \(P\), and the magnetic field by \(\boldsymbol{B}\). All four of these quantities (eight, if you count the three components each for \(\boldsymbol{v}\) and \(\boldsymbol{B}\)) are 3-dimensional functions of position, in contrast to the 1-dimensional functions of radius \(\mathrm{f}_i(r)\). The velocity and magnetic fields are subject to the constraints

(7)\[ \boldsymbol{\nabla}\cdot\left[\mathrm{f}_1(r)\,\boldsymbol{v}\right] = 0\]

and

(8)\[ \boldsymbol{\nabla}\cdot\boldsymbol{B}=0,\]

respectively. The evolution of \(\Theta\) is described by

(9)\[\begin{split}\mathrm{f}_1(r)\,\mathrm{f}_4(r)\left[\frac{\partial \Theta}{\partial t} + \boldsymbol{v}\cdot\boldsymbol{\nabla}\Theta + c_{11}\,\mathrm{f}_{14}(r)v_r\right] =\ c_6\,\boldsymbol{\nabla}\cdot\left[\mathrm{f}_1(r)\,\mathrm{f}_4(r)\,\mathrm{f}_5(r)\,\boldsymbol{\nabla}\Theta \right] \\ +\ c_{10}\,\mathrm{f}_6(r) + c_8\,\Phi(r,\theta,\phi) + c_9\,\mathrm{f}_7(r)|\boldsymbol{\nabla}\times\boldsymbol{B}|^2,\end{split}\]

where the viscous heating \(\Phi\) is given by

(10)\[\begin{split} \Phi(r,\theta,\phi) = c_5\,\mathcal{D}_{ij}e_{ij} &= 2\,c_5\,\mathrm{f}_1(r)\mathrm{f}_3(r)\left[e_{ij}e_{ij} - \frac{1}{3}\left(\boldsymbol{\nabla}\cdot\boldsymbol{v}\right)^2\right] \\ &= 2\,c_5\,\mathrm{f}_1(r)\mathrm{f}_3(r)\left[e_{ij} - \frac{1}{3}\left(\boldsymbol{\nabla}\cdot\boldsymbol{v}\right)\delta_{ij}\right]^2.\end{split}\]

Finally, the evolution of \(\boldsymbol{B}\) is described by the induction equation

(11)\[ \frac{\partial \boldsymbol{B}}{\partial t} = \boldsymbol{\nabla}\times\left[\boldsymbol{v}\times\boldsymbol{B} - c_7\,\mathrm{f}_7(r)\boldsymbol{\nabla}\times\boldsymbol{B}\,\right].\]

Note that when Rayleigh actually solves the equations, the following additional derivative functions are used:

\[\begin{split}\mathrm{f}_8(r) &= \frac{d\ln{\mathrm{f}_1}}{dr}\\ \mathrm{f}_9(r) &= \frac{d^2\ln{\mathrm{f}_1}}{dr^2}\\ \mathrm{f}_{10}(r) &= \frac{d\ln{\mathrm{f}_4}}{dr}\\ \mathrm{f}_{11}(r) &= \frac{d\ln{\mathrm{f}_3}}{dr}\\ \mathrm{f}_{12}(r) &= \frac{d\ln{\mathrm{f}_5}}{dr}\\ \mathrm{f}_{13}(r) &= \frac{d\ln{\mathrm{f}_7}}{dr}.\end{split}\]

When supplying a custom reference state, the user may specify the six derivative functions “by hand.” If the user fails to do so, Rayleigh will compute the required derivatives (only if the user supplies the function whose derivative is to be taken) from the function’s Chebyshev coefficients.

Note that equations (5)-(11) could have been formulated in other ways. For instance, we could combine \(\mathrm{f}_1\) and \(\mathrm{f}_3\) into a single function in Equation (10). The form of the equations presented here has been chosen to reflect that actually used in the code, which was originally written dimensionally.

We now describe the nondimensional Boussinesq, and dimensional/nondimensional anelastic formulations used in the code.

Nondimensional Boussinesq Formulation of the MHD Equations

Rayleigh can be run using a nondimensional, Boussinesq formulation of the MHD equations (reference_type=1). The nondimensionalization employed is as follows:

\[\begin{split}\begin{aligned} \mathrm{Length} &\rightarrow L &\;\;\;\; \mathrm{(Shell\; Depth)} \\ \mathrm{Time} &\rightarrow \frac{L^2}{\nu_o} &\;\;\;\; \mathrm{(Viscous\; Timescale)}\\ \mathrm{Temperature} &\rightarrow \Delta T&\;\;\;\; \mathrm{(Temperature\; Contrast\; Across\; Shell)} \\ \mathrm{Magnetic\; Field} &\rightarrow \sqrt{\hat{\rho}\mu\eta_o\Omega_0} \\ \mathrm{Reduced\; Pressure} &\rightarrow \nu_o\Omega_0&\;\;\; (\mathrm{[Thermodynamic\; Pressure]}/\hat{\rho}),\end{aligned}\end{split}\]

where \(\Omega_0\) is the rotation rate of the frame, \(\hat{\rho}\) is the (constant) density of the fluid, \(\eta_o\) is the magnetic diffusivity at the top of the domain (i.e., at \(r=r_o\)), \(\nu_o\) is the kinematic viscosity at the top of the domain, and \(\mu\) is the magnetic permeability. Note that in Gaussian units for vacuum, \(\mu=4\pi\). After nondimensionalizing, the following nondimensional numbers appear in our equations:

\[\begin{split}\begin{aligned} Pr &=\frac{\nu_o}{\kappa_o} &\;\;\;\;\;\; \mathrm{Prandtl\; Number} \\ Pm &=\frac{\nu_o}{\eta_o} &\;\;\;\;\;\; \mathrm{Magnetic\; Prandtl\; Number} \\ E &=\frac{\nu_o}{\Omega_0\,L^2} &\;\;\;\;\;\; \mathrm{Ekman\; Number} \\ Ra &=\frac{\alpha g_o \Delta T\,L^3}{\nu_o\kappa_o} &\;\;\;\;\;\; \mathrm{Rayleigh\; Number},\end{aligned}\end{split}\]

where \(\alpha\) is coefficient of thermal expansion, \(g_o\) is the gravitational acceleration at the top of the domain, and \(\kappa\) is the thermal diffusivity. Adopting this nondimensionalization is equivalent to assigning the following to the functions \(\mathrm{f}_i(r)\) and the constants \(c_i\):

\[\begin{split}\begin{aligned} \mathrm{f}_1(r) &\rightarrow 1\; &c_1 &\rightarrow \frac{2}{E} \\ \mathrm{f}_2(r) &\rightarrow \left(\frac{r}{r_o}\right)^n \; &c_2 &\rightarrow \frac{Ra}{Pr} \\ \mathrm{f}_3(r) &\rightarrow \tilde{\nu}(r)\; &c_3 &\rightarrow \frac{1}{E}\\ \mathrm{f}_4(r) &\rightarrow 1\; &c_4 &\rightarrow \frac{1}{E\,Pm} \\ \mathrm{f}_5(r) &\rightarrow \tilde{\kappa}(r)\; &c_5 &\rightarrow 1 \\ \mathrm{f}_6(r) &\rightarrow 0\; &c_6 &\rightarrow \frac{1}{Pr} \\ \mathrm{f}_7(r) &\rightarrow \tilde{\eta}(r)\; &c_7 &\rightarrow \frac{1}{Pm} \\ &\vdots &c_8&\rightarrow 0\\ &\vdots &c_9&\rightarrow 0 \\ &\vdots &c_{10}&\rightarrow 0 \\ \mathrm{f}_{14}(r)&\rightarrow 0\; &c_{11}&\rightarrow 0.\end{aligned}\end{split}\]

Here the tildes denote nondimensional radial profiles, e.g., \(\tilde{\nu}(r) = \nu(r)/\nu_o\).

Our choice of \(\mathrm{f}_{14}(r)\rightarrow 0\) sets the default atmosphere in non-dimensional Boussinesq to be neutrally stable. For other choices (i.e., convectively stable or unstable), one must use the custom-reference-state framework.

Our choice of \(\mathrm{f}_2(r)\) allows gravity to vary with radius based on the value of the exponent \(n\), which has a default value of \(0\) in the code. Note also that our definition of \(Ra\) assumes fixed-temperature boundary conditions. We might specify fixed-flux boundary conditions and/or an internal heating through a suitable choice \(c_{10}\mathrm{f}_6(r)\), in which case the meaning of \(Ra\) in our equation set changes, with \(Ra\) denoting a flux Rayleigh number instead. In addition, ohmic and viscous heating, which do not appear in the Boussinesq formulation, are turned off when this nondimensionalization is specified at runtime. When these substitutions are made, Equations (5)-(11) transform as follows.

\[\begin{split}\begin{aligned} \left[\frac{\partial \boldsymbol{v}}{\partial t} + \boldsymbol{v}\cdot\boldsymbol{\nabla}\boldsymbol{v} %advection + \frac{2}{E}\boldsymbol{\hat{z}}\times\boldsymbol{v} \right] &= % Coriolis \frac{Ra}{Pr}\left(\frac{r}{r_o}\right)^n\Theta\,\boldsymbol{\hat{r}} % buoyancy - \frac{1}{E}\boldsymbol{\nabla}P % pressure + \frac{1}{E\,Pm}\left(\boldsymbol{\nabla}\times\boldsymbol{B}\right)\times\boldsymbol{B} % Lorentz Force + \boldsymbol{\nabla}\cdot\boldsymbol{\mathcal{D}}& \\ % % \left[\frac{\partial \Theta}{\partial t} + \boldsymbol{v}\cdot\boldsymbol{\nabla}\Theta \right] &= \frac{1}{Pr}\boldsymbol{\nabla}\cdot\left[\tilde{\kappa}(r)\boldsymbol{\nabla}\Theta\right] \\ % Diffusion % % \frac{\partial \boldsymbol{B}}{\partial t} &= \boldsymbol{\nabla}\times\left[\boldsymbol{v}\times\boldsymbol{B} - \frac{1}{Pm}\tilde{\eta}(r)\boldsymbol{\nabla}\times\boldsymbol{B}\right]\\ \mathcal{D}_{ij} &= 2\tilde{\nu}(r)e_{ij} \\ % % \boldsymbol{\nabla}\cdot\boldsymbol{v}&=0\\ \boldsymbol{\nabla}\cdot\boldsymbol{B}&=0 \end{aligned}\end{split}\]

Here \(\Theta\) refers to the temperature (perturbation from the background) and \(P\) to the reduced pressure (ratio of the thermodynamic pressure to the constant density).

Dimensional Anelastic Formulation of the MHD Equations

When run in dimensional, anelastic mode (cgs units; reference_type=2 ), the following values are assigned to the functions \(\mathrm{f}_i\) and the constants \(c_i\):

\[\begin{split}\begin{aligned} \mathrm{f}_1(r) &\rightarrow \hat{\rho}(r)\; &c_1 &\rightarrow 2\Omega_0 \\ \mathrm{f}_2(r) &\rightarrow \frac{\hat{\rho}(r)}{c_P}g(r)\; &c_2 &\rightarrow 1 \\ \mathrm{f}_3(r) &\rightarrow \nu(r)\; &c_3 &\rightarrow 1\\ \mathrm{f}_4(r) &\rightarrow \hat{T}(r)\; &c_4 &\rightarrow \frac{1}{4\pi} \\ \mathrm{f}_5(r) &\rightarrow \kappa(r)\; &c_5 &\rightarrow 1 \\ \mathrm{f}_6(r) &\rightarrow \frac{Q(r)}{L_*}\; &c_6 &\rightarrow 1 \\ \mathrm{f}_7(r) &\rightarrow \eta(r)\; &c_7 &\rightarrow 1 \\ &\vdots &c_8&\rightarrow 1\\ &\vdots &c_9&\rightarrow \frac{1}{4\pi} \\ &\vdots &c_{10}&\rightarrow L_* \\ \mathrm{f}_{14}(r)&\rightarrow \frac{d\hat{S}}{dr }&c_{11}&\rightarrow 1.\end{aligned}\end{split}\]

Here \(\hat{\rho}(r)\), \(\hat{T}(r)\), and \(d\hat{S}/dr\) are the spherically symmetric, time-independent reference-state density, temperature, and entropy gradient, respectively. The thermal variables satisfy the linearized equation of state

\[\frac{P}{\hat{P}}= \frac{T}{\hat{T}} + \frac{\rho}{\hat{\rho}}\]

\(g(r)\) is the gravitational acceleration, \(c_P\) is the specific heat at constant pressure, and \(\Omega_0\) is the frame rotation rate. The viscous, thermal, and magnetic diffusivities (also assumed to be spherically symmetric and time-independent) are given by \(\nu(r)\), \(\kappa(r)\), and \(\eta(r)\), respectively. Note that the entropy gradient term \(f_{14}(r)v_r\) is only used in Equation (9) if advect_reference_state=.true.. Finally, \(Q(r)\) is an internal heating function; it might represent radiative heating or heating due to nuclear fusion, for instance. In our convention, the volume integral of \(\mathrm{f}_6(r)\) equals unity, and \(c_{10}\) equals the luminosity or heating_integral \(L_*\) specified in the main_input file. When using a custom reference state, this allows easy adjustment of the luminosity using the override_constants formalism, e.g.,

override_constants(10) = T

ra_constants(10) = 3.846d33

specified in the in the reference_namelist.

Note that in the anelastic formulation, the thermal variable \(\Theta\) is interpreted as the entropy perturbation, rather than the temperature perturbation. When these substitutions are made, Equations (5)-(11) transform as follows.

\[\begin{split}\begin{aligned} \hat{\rho}(r)\left[\frac{\partial \boldsymbol{v}}{\partial t} +\boldsymbol{v}\cdot\boldsymbol{\nabla}\boldsymbol{v} %advection +2\Omega_0\boldsymbol{\hat{z}}\times\boldsymbol{v} \right] =\; % Coriolis &\frac{\hat{\rho}(r)}{c_P}g(r)\Theta\,\boldsymbol{\hat{r}} % buoyancy +\hat{\rho}(r)\boldsymbol{\nabla}\left(\frac{P}{\hat{\rho}(r)}\right) % pressure \\ &+\frac{1}{4\pi}\left(\boldsymbol{\nabla}\times\boldsymbol{B}\right)\times\boldsymbol{B} % Lorentz Force +\boldsymbol{\nabla}\cdot\boldsymbol{\mathcal{D}}\\ % % \hat{\rho}(r)\,\hat{T}(r)\left[\frac{\partial \Theta}{\partial t} +\boldsymbol{v}\cdot\boldsymbol{\nabla}\Theta + v_r\frac{d\hat{S}}{dr}\right] =\; &\boldsymbol{\nabla}\cdot\left[\hat{\rho}(r)\,\hat{T}(r)\,\kappa(r)\,\boldsymbol{\nabla}\Theta \right] % diffusion +Q(r) % Internal heating \\ &+\Phi(r,\theta,\phi) +\frac{\eta(r)}{4\pi}\left[\boldsymbol{\nabla}\times\boldsymbol{B}\right]^2\\ % Ohmic Heating % % \frac{\partial \boldsymbol{B}}{\partial t} =\; &\boldsymbol{\nabla}\times\left[\,\boldsymbol{v}\times\boldsymbol{B}-\eta(r)\boldsymbol{\nabla}\times\boldsymbol{B}\,\right] \\ % % \mathcal{D}_{ij} =\; &2\hat{\rho}(r)\,\nu(r)\left[e_{ij}-\frac{1}{3}\left(\boldsymbol{\nabla}\cdot\boldsymbol{v}\right)\delta_{ij}\right] \\ % % \Phi(r,\theta,\phi) =\; &2\,\hat{\rho}(r)\nu(r)\left[e_{ij}e_{ij}-\frac{1}{3}\left(\boldsymbol{\nabla}\cdot\boldsymbol{v}\right)^2\right] \\ % % \boldsymbol{\nabla}\cdot\left[\hat{\rho}(r)\,\boldsymbol{v}\right] =\; &0 \\ \boldsymbol{\nabla}\cdot\boldsymbol{B} =\; &0. \end{aligned}\end{split}\]

Nondimensional Anelastic MHD Equations

To run in nondimensional anelastic mode, you must set reference_type=3 in the Reference_Namelist. The reference state is assumed to be polytropic with a \(\frac{1}{r^2}\) profile for gravity. When this mode is active, the following nondimensionalization is used (following Heimpel et al., 2016, Nat. Geo., 9, 19 ):

\[\begin{split}\begin{aligned} \mathrm{Length} &\rightarrow L \equiv r_o - r_i &\;\;\;\; \mathrm{(Shell\; Depth)} \\ \mathrm{Time} &\rightarrow \frac{1}{\Omega_0} &\;\;\;\; \mathrm{(Rotational\; Timescale)}\\ \mathrm{Temperature} &\rightarrow T_o\equiv\hat{T}(r_o)&\;\;\;\; \mathrm{(Reference\; Temperature\; at\; Upper\; Boundary)} \\ \mathrm{Density} &\rightarrow \rho_o\equiv\hat{\rho}(r_o)&\;\;\;\; \mathrm{(Reference\; Density\; at\; Upper\; Boundary)} \\ \mathrm{Entropy} &\rightarrow \Delta{s}&\;\;\;\; \mathrm{(Entropy\; Constrast\; Across\; Shell)} \\ \mathrm{Magnetic~Field} &\rightarrow \sqrt{\hat{\rho}_o\mu\eta_o\Omega_0} \\ \mathrm{Pressure} &\rightarrow \rho_oL^2\Omega_0^2.\end{aligned}\end{split}\]

When run in this mode, Rayleigh employs a polytropic background state, with an assumed \(\frac{1}{r^2}\) variation in gravity. These choices result in the functions \(\mathrm{f}_i\) and the constants \(c_i\) (tildes indicate nondimensional reference-state variables):

\[\begin{split}\begin{aligned} \mathrm{f}_1(r) &\rightarrow \tilde{\rho}(r)\; &c_1 &\rightarrow 2 \\ \mathrm{f}_2(r) &\rightarrow \tilde{\rho}(r)\frac{r_\mathrm{max}^2}{r^2}\; &c_2 &\rightarrow \mathrm{Ra}^* \\ \mathrm{f}_3(r) &\rightarrow \tilde{\nu}(r)\; &c_3 &\rightarrow 1\\ \mathrm{f}_4(r) &\rightarrow \tilde{T}(r)\; &c_4 &\rightarrow \frac{\mathrm{E}}{\mathrm{Pm}} \\ \mathrm{f}_5(r) &\rightarrow \tilde{\kappa}(r)\; &c_5 &\rightarrow \mathrm{E} \\ \mathrm{f}_6(r) &\rightarrow \frac{\tilde{Q}(r)}{L_*}; &c_6 &\rightarrow \frac{\mathrm{E}}{\mathrm{Pr}} \\ \mathrm{f}_7(r) &\rightarrow \tilde{\eta}(r) \; &c_7 &\rightarrow \frac{\mathrm{E}}{\mathrm{Pm}} \\ &\vdots &c_8&\rightarrow \frac{\mathrm{E}\,\mathrm{Di}}{\mathrm{Ra}^*}\\ &\vdots &c_9&\rightarrow \frac{\mathrm{E}^2\,\mathrm{Di}}{\mathrm{Pm}^2\mathrm{Ra}^*}\\ &\vdots &c_{10}&\rightarrow L_* \\ \mathrm{f}_{14}(r)&\rightarrow 0&c_{11}&\rightarrow 0.\end{aligned}\end{split}\]

As in the Boussinesq case, the nondimensional diffusivities are defined according to, e.g., \(\tilde{\nu}(r) \equiv \nu(r)/\nu_o\). The nondimensional heating \(\tilde{Q}(r)\) is defined such that its volume integral equals the nondimensional luminosity or heating_integral set in the main_input file. As in the dimensional anelastic case, the volume integral of \(\mathrm{f}_6(r)\) equals unity, and \(\mathrm{c}_{10} = L_*\). The unit for luminosity in this nondimensionalization (to get a dimensional luminosity from the nondimensional \(L_*\)) is \(\rho_oL^3T_o\Delta s\Omega_0\).

Two new nondimensional numbers appear in our equations, in addition to those defined for the Boussinesq case. \(\mathrm{Di}\), the dissipation number, is defined by

(12)\[ \mathrm{Di}= \frac{g_o\,\mathrm{L}}{c_\mathrm{P}\,T_o},\]

where \(g_o\) and \(T_o\) are the gravitational acceleration and temperature at the outer boundary respectively. Once more, the thermal anomaly \(\Theta\) should be interpreted as (nondimensional) entropy. The symbol \(\mathrm{Ra}^*\) is the modified Rayleigh number, given by

(13)\[\mathrm{Ra}^*=\frac{g_o}{c_\mathrm{P}\Omega_0^2}\frac{\Delta s}{L}\]

We thus arrive at the following nondimensionalized equations:

\[\begin{split}\begin{aligned} \tilde{\rho}(r)\left[\frac{\partial \boldsymbol{v}}{\partial t} + \boldsymbol{v}\cdot\boldsymbol{\nabla}\boldsymbol{v} %advection + 2\boldsymbol{\hat{z}}\times\boldsymbol{v}\right] =\; % Coriolis &\mathrm{Ra}^*\tilde{\rho}(r)\left(\frac{r_\mathrm{max}^2}{r^2}\right)\Theta\,\boldsymbol{\hat{r}} % buoyancy + \tilde{\rho}(r)\boldsymbol{\nabla}\left(\frac{P}{\tilde{\rho}(r)}\right) % pressure \\ &+ \frac{\mathrm{E}}{\mathrm{Pm}}\left(\boldsymbol{\nabla}\times\boldsymbol{B}\right)\times\boldsymbol{B} % Lorentz Force + \mathrm{E}\boldsymbol{\nabla}\cdot\boldsymbol{\mathcal{D}}\\ % % \tilde{\rho}(r)\,\tilde{T}(r)\left[\frac{\partial \Theta}{\partial t} + \boldsymbol{v}\cdot\boldsymbol{\nabla}\Theta\right] =\; &\frac{\mathrm{E}}{\mathrm{Pr}}\boldsymbol{\nabla}\cdot\left[\tilde{\kappa}(r)\tilde{\rho}(r)\,\tilde{T}(r)\,\boldsymbol{\nabla}\Theta \right] % diffusion + \tilde{Q}(r) % Internal heating \\ &+ \frac{\mathrm{E}\,\mathrm{Di}}{\mathrm{Ra}^*}\Phi(r,\theta,\phi) + \frac{\mathrm{Di\,E^2}}{\mathrm{Pm}^2\mathrm{Ra}^*}\tilde{\eta}(r)|\boldsymbol{\nabla}\times\boldsymbol{B}|^2 \\ % Ohmic Heating % % \frac{\partial \boldsymbol{B}}{\partial t} =\; &\boldsymbol{\nabla}\times\left[\,\boldsymbol{v}\times\boldsymbol{B}-\frac{\mathrm{E}}{\mathrm{Pm}}\tilde{\eta}(r)\boldsymbol{\nabla}\times\boldsymbol{B}\,\right] \\ % % \mathcal{D}_{ij} =\; &2\tilde{\rho}(r)\tilde{\nu}(r)\left[e_{ij} - \frac{1}{3}\boldsymbol{\nabla}\cdot\boldsymbol{v}\right] \\ % % \Phi(r,\theta,\phi) =\; &2\tilde{\rho}(r)\tilde{\nu}(r)\left[e_{ij}e_{ij} - \frac{1}{3}\left(\boldsymbol{\nabla}\cdot\boldsymbol{v}\right)^2\right] \\ % % \boldsymbol{\nabla}\cdot\left[\tilde{\rho}(r)\boldsymbol{v}\right]=\; &0 \\ \boldsymbol{\nabla}\cdot\boldsymbol{B} =\; &0. \end{aligned}\end{split}\]

The Streamfunction Formulation

The velocity field in Rayleigh is evolved subject to the solenoidal constraint

(14)\[ \boldsymbol{\nabla}\cdot\left[\mathrm{f}_1(r)\,\boldsymbol{v}\right] = 0.\]

This is accomplished by casting \(\mathrm{f}_1\boldsymbol{v}\) in terms of streamfunctions such that

(15)\[ \mathrm{f_1}\,\boldsymbol{v} = \boldsymbol{\nabla}\times\boldsymbol{\nabla}\times\left( W\,\boldsymbol{\hat{r}}\right)+\boldsymbol{\nabla}\times\left( Z\,\boldsymbol{\hat{r}}\right),\]

where W and Z are referred to as the poloidal and toroidal stream functions respectively. Rather than evolving the three components of \(\boldsymbol{v}\) directly, the momentum equations are cast in terms of these variables before advancing the timestep. The velocity components are related to the streamfunctions via the relations:

(16)\[ \mathrm{f_1}v_r = - \frac{1}{r^2\mathrm{sin}\theta}\frac{\partial}{\partial\theta}\left(\mathrm{sin}\theta\frac{\partial W}{\partial\theta} \right)-\frac{1}{r^2\mathrm{sin}^2\theta}\frac{\partial^2 W}{\partial\phi^2},\]

(17)\[ \mathrm{f_1}\,v_{\theta} = \frac{1}{r}\frac{\partial^2 W}{\partial r\partial\theta}+ \frac{1}{r\mathrm{sin}\theta}\frac{\partial Z}{\partial\phi},\]

and

(18)\[ \mathrm{f_1}v_{\phi} = \frac{1}{r\mathrm{sin}\theta}\frac{\partial^2 W}{\partial r\partial\phi} - \frac{1}{r}\frac{\partial Z}{\partial\theta}.\]

When the velocity field and streamfunctions are projected onto a spherical harmonic basis, two additional useful relations are given by

(19)\[ \left[\mathrm{f_1}\,v_r\right]_\ell^m = \frac{\ell(\ell+1)}{r^2}W_\ell^m\]

and

(20)\[ \left[ \left\{\boldsymbol{\nabla}\times\left(\mathrm{f_1}\,\boldsymbol{v}\right)\right\}_r \right]_\ell^m = \frac{\ell(\ell+1)}{r^2}Z_\ell^m.\]

A similar decomposition is performed on the magnetic field to ensure it remains divergence free. In that case, the magnetic field is projected onto flux functions such that

(21)\[ \boldsymbol{B} = \boldsymbol{\nabla}\times\boldsymbol{\nabla}\times\left( C\,\boldsymbol{\hat{r}}\right)+\boldsymbol{\nabla}\times\left( A\,\boldsymbol{\hat{r}}\right),\]

where C and A are the poloidal and toroidal flux functions respectively. Similar to the velocity field, the components of \(\boldsymbol{B}\) satisfy

(22)\[ B_r = - \frac{1}{r^2\mathrm{sin}\theta}\frac{\partial}{\partial\theta}\left(\mathrm{sin}\theta\frac{\partial C}{\partial\theta} \right)-\frac{1}{r^2\mathrm{sin}^2\theta}\frac{\partial^2 C}{\partial\phi^2},\]

(23)\[ B_{\theta} = \frac{1}{r}\frac{\partial^2 C}{\partial r\partial\theta}+ \frac{1}{r\mathrm{sin}\theta}\frac{\partial A}{\partial\phi},\]

(24)\[ B_{\phi} = \frac{1}{r\mathrm{sin}\theta}\frac{\partial^2 C}{\partial r\partial\phi} - \frac{1}{r}\frac{\partial A}{\partial\theta},\]

(25)\[ \left[B_r\right]_\ell^m = \frac{\ell(\ell+1)}{r^2}C_\ell^m,\]

and

(26)\[ \left[ \left\{\boldsymbol{\nabla}\times\boldsymbol{B}\right\}_r \right]_\ell^m = \frac{\ell(\ell+1)}{r^2}A_\ell^m.\]

The Pseudospectral Approach

Section needed.

Parallelization and Performance

Section needed.