Overview of Diagnostic Outputs in Rayleigh

The purpose of this document is to describe Rayleigh’s internal menu system used for specifying diagnostic outputs. Rayleigh’s design includes an onboard diagnostics package that allows a user to output a variety of system quantities as the run evolves. These include system state variables, such as velocity and entropy, as well as derived quantities, such as the vector components of the Lorentz force and the kinetic energy density. Each diagnostic quantity is requested by adding its associated menu number to the main_input file. Radial velocity, for instance, has menu code 1, \(\theta\)-component of velocity has menu code 2, etc.

A few points to keep in mind are

• This document is intended to describe the diagnostics output menu only. A complete description of Rayleigh’s diagnostic package is provided in Rayleigh/doc/Diagnostic_Plotting.pdf. A more in-depth description of the anelastic and Boussinesq modes available in Rayleigh is provided in Rayleigh/doc/user_guide.pdf.

• A number of output methods may be used to output any system diagnostic. No diagnostic is linked to a particular output method. The same diagnostic might be output in volume-averaged, azimuthally-averaged, and fully 3-D form, for instance.

• You may notice a good deal of redundancy in the available outputs. For instance, the azimuthal velocity, \(v_\phi\), and its zonal average, \(\overline{v_\phi}\), are both available as outputs. Were the user to output both of these in an azimuthally-averaged format, the result would be the same. 3-D output, however, would not yield the same result. This redundancy has been added to help with post-processing calculations in which it can be useful to have all data products in a similar format.

• Given the degree of redundancy found in the list below, you may be surprised to notice that several values are not available for output at all. Some of these are best added as custom-user diagnostics and may be included in a future release. Many, however, may be obtained by considering either the sum, or difference, of those outputs already available.

Definitions 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 Equation Sets Solved by 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 + \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) = \mathcal{D}_{ij}e_{ij} &= 2\,\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\,\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 \\ \mathrm{f}_{14}(r)&\rightarrow 0\; &c_{10}&\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} \\ \mathrm{f}_{14}(r)&\rightarrow \frac{d\hat{S}}{dr }&c_{10}&\rightarrow L_*.\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. \(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}^*}\\ \mathrm{f}_{14}(r)&\rightarrow \frac{d\tilde{S}}{dr }&c_{10}&\rightarrow L_*.\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 + v_r\frac{d\hat{S}}{dr}\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}\]