Physics Controls

Rayleigh solves the MHD equations in spherical geometry under the Boussinesq and anelastic approximations. 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 non-dimensional anelastic formulations (based on a polytropic reference state) are provided as part of Rayleigh.

In this section, we present the equation sets solved when running in each of these three modes, and discuss the relevant control parameters for each mode. We also discuss the boundary conditions available in Rayleigh and those namelist variables that can be used to modify the code’s behavior in any of these three modes.

Anelastic Mode (dimensional)

Example Input: Rayleigh/input_examples/main_input_sun

When run in dimensional, anelastic mode, reference_type=2 must be specified in the Reference_Namelist. In that case, Rayleigh solves the following form of the MHD equations:

\[\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}} \;\;\; &\mathrm{Momentum}\\ % % \hat{\rho}(r)\,\hat{T}(r)\left[\frac{\partial \Theta}{\partial t} +\boldsymbol{v}\cdot\boldsymbol{\nabla}\Theta \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 &\mathrm{Thermal\; Energy}\\ % 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) &\mathrm{Induction} \\ % % \mathcal{D}_{ij} &= 2\hat{\rho}(r)\,\nu(r)\left[e_{ij}-\frac{1}{3}\boldsymbol{\nabla}\cdot\boldsymbol{v}\right] &\mathrm{Viscous\; Stress\; Tensor}\\ % % \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] &\mathrm{Viscous\; Heating} \\ % % \boldsymbol{\nabla}\cdot\left(\hat{\rho}(r)\,\boldsymbol{v}\right)&=0 &\mathrm{Solenoidal\; Mass\; Flux}\\ \boldsymbol{\nabla}\cdot\boldsymbol{B}&=0 &\mathrm{Solenoidal\; Magnetic\; Field}\end{aligned}\end{split}\]

Here, \(\hat{\rho}\) and \(\hat{T}\) are the reference-state density and temperature respectively. \(g\) is the gravitational acceleration, \(c_P\) is the specific heat at constant pressure, and \(\Omega_0\) is the frame rotation rate. The velocity field vector is denoted by \(\boldsymbol{v}\), the magnetic field vector by \(\boldsymbol{B}\), and the pressure by \(P\). The thermal anomoly is denoted by \(\Theta\) and should be interpreted is as entropy \(s\) in this formulation. The thermal variables satisfy the linearized equation of state

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

The kinematic viscosity, thermal diffusivity, and magnetic diffusivity are given by \(\nu\), \(\kappa\), and \(\eta\) respectively. Finally, \(Q(r)\) is an internal heating function; it might represent radiative heating or heating due to nuclear fusion, for instance.

When running in anelastic mode, the reference_type variable in the Reference_Namelist must be set to 2.

Moreover, certain variables in the Reference_Namelist and the Transport_Namelist must be specified. The Reference_Namelist variables are described in Table table_anelastic and the Transport_Namelist variables are described in Table table_anelastic_trans. Default values indicated in brackets.

Table. Anelastic.

Variables in the Reference_Namelist that must be specified when running in dimensional anelastic mode. In addition, reference_type=2 must also be specified.

Variable [Default value]	Description
poly_n [0]	polytropic index (\(P\propto\rho^n\))
poly_Nrho [0]	number of density scaleheights spanning the domain
poly_mass [0]	mass interior to \(rmin\)
poly_rho_i [0]	density at rmin, \(\rho(r=rmin)\)
pressure_specific_heat [0]	specific heat at constant pressure
angular_velocity [1.0]	cyclic frequency of the rotating frame

Table. Anelastic Transport.

Variables in the Transport_Namelist that must be specified when running in dimensional anelastic mode. In addition, reference_type=2 must also be specified in the Reference_Namelist.

Variable [Default value]	Description
nu_top [1.0]	kinematic viscosity at rmax, \(\nu(rmax)\)
nu_type [1]	determines whether \(\nu\) is constant with radius (1) or varies with density (2)
nu_power [0.0]	exponent in : \(\nu(r) = \left( \frac{\rho (r)}{\rho(r=rmax)} \right)^ {nu\_power}\); use with nu_type=2
kappa_top [1.0]	thermal diffusivity at rmax, \(\kappa(rmax)\)
kappa_type [1]	determines whether \(\kappa\) is constant with radius (1) or varies with density (2)
kappa_power [0.0]	exponent in : \(\kappa(r) = \left( \frac{\ rho(r)}{\rho(r=rmax)} \right)^ {kappa\_power}\); use with kappa_type=2
eta_top [1.0]	magnetic diffusivity at rmax, \(\eta(rmax)\)
eta_type [1]	determines whether \(\eta\) is constant with radius (1) or varies with density (2)
eta_power [0.0]	exponent in : \(\eta(r) = \left( \frac{ \rho(r)}{\rho(r=rmax)} \right)^ {eta\_power}\); use with eta_type=2

The polytropic reference state is the same as that used in the benchmarks and is described in detail in Jones et al. (2011).

See the example input file main_input_sun for a an example of how to run a solar-like model using Rayleigh’s dimensional, anelastic formulation.

Boussinesq Mode (nondimensional)

Example Input: Rayleigh/input_examples/c2001_case1_input

When run in nondimensional Boussinesq mode, reference_type=1 must be specified in the Reference_Namelist. In that case, Rayleigh employs the nondimensionalization

\[\begin{split}\begin{aligned} \mathrm{Length} &\rightarrow L &\;\;\;\; \mathrm{(Shell\; Depth)} \\ \mathrm{Time} &\rightarrow \frac{L^2}{\nu} &\;\;\;\; \mathrm{(Viscous\; Timescale)}\\ \mathrm{Temperature} &\rightarrow \Delta T&\;\;\;\; \mathrm{(Temperature\; Contrast\; Across\; Shell)} \\ \mathrm{Magnetic~Field} &\rightarrow \sqrt{\rho\mu\eta\Omega_0},\end{aligned}\end{split}\]

where \(\Omega_0\) is the rotation rate of the frame, \(\rho\) is the (constant) density of the fluid, \(\mu\) is the magnetic permeability, \(\eta\) is the magnetic diffusivity, and \(\nu\) is the kinematic viscosity. After nondimensionalizing, the following nondimensional numbers appear in our equations

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

where \(\alpha\) is coefficient of thermal expansion, \(g_0\) is the gravitational acceleration at the top of the domain, and \(\kappa\) is the thermal diffusivity.

In addition, ohmic and viscous heating, which do not appear in the Boussinesq formulation, are turned off when this nondimensionalization is specified at runtime. Rayleigh solves the following equations when running in nondimensional Boussinesq mode:

\[\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}^2\boldsymbol{v} \;\;\; &\mathrm{Momentum}\\ % % \left[\frac{\partial \Theta}{\partial t} +\boldsymbol{v}\cdot\boldsymbol{\nabla}\Theta \right] &= \frac{1}{Pr}\boldsymbol{\nabla}^2\Theta &\mathrm{Thermal\; Energy}\\ % Diffusion % % \frac{\partial \boldsymbol{B}}{\partial t} &= \boldsymbol{\nabla}\times\left(\,\boldsymbol{v}\times\boldsymbol{B}\right)+\frac{1}{Pm}\boldsymbol{\nabla}^2\boldsymbol{B} &\mathrm{Induction} \\ % % % % % % \boldsymbol{\nabla}\cdot\boldsymbol{v}&=0 &\mathrm{Solenoidal\; Velocity\; Field}\\ \boldsymbol{\nabla}\cdot\boldsymbol{B}&=0 &\mathrm{Solenoidal\; Magnetic\; Field},\end{aligned}\end{split}\]

where \(r_0 \equiv rmax\). In this formulation, \(\Theta\) should be interpreted as the temperature perturbation \(T\). Those Reference_Namelist variables that must be set for this model are indicated in Table table_boussinesq.

Note that our choice for the temperature scale assumes fixed-temperature boundary conditions. We might choose to specify fixed-flux boundary conditions and/or an internal heating, in which case the meaning of \(\Delta T\) in our equation set changes, with \(\Delta T \equiv L\frac{\partial T}{\partial r}\) instead, for some fiducial value of \(\frac{\partial T}{\partial r}\). Which regard to the temperature scaling, it is up to the user to select boundary conditions appropriate for their desired values of \(\Delta T\). If \(\Delta T\) denotes the temperature contrast across the domain, then their boundary condition variables should look like:

&boundary\_conditions\_namelist
T_Top    = 0.0d0
T_Bottom = 1.0d0
fix_tvar_top = .true.
fix_tvar_bottom = .true.
/

Alternatively, if the temperature scale is determined by a gradient at one boundary, the user should ensure that the amplitude of the temperature gradient at that boundary is 1. For example:

&boundary\_conditions\_namelist
dTdr_bottom = -1.0d0
fix_dtdr_bottom = .true.
/

Boundary conditions and internal heating are discussed in §Boundary Conditions & Internal Heating. Finally, in Boussinesq mode, the namelist variables nu_type, kappa_type, and eta_type should be set to 1. Their values will be determined by Pr and Pm, instead of nu_top, kappa_top, or eta_top.

Table. Boussinesq.

Variables in the Reference_Namelist that must be specified when running in nondimensional Boussinesq mode. In addition, reference_type=1 must also be specified.

Variable	Description
Ekman_Number	The Ekman Number \(E\)
Rayleigh_Number	The Rayleigh Number \(Ra\)
Prandtl_Number	The Prandtl Number \(Pr\)
Magnetic_Prandtl_Number	The Magnetic Prandtl Number \(Pm\)
Gravity_Power	Buoyancy coefficient = \(\frac{\mathrm{Ra}}{\mathrm {Pr}}\left(\frac{r}{rmax} \right) ^\mathrm{gravity\_power}\)

Anelastic Mode (nondimensional)

Example Input: Rayleigh/input_examples/main_input_jupiter

When running in nondimensional anelastic mode, you must set reference_type=3 in the Reference_Namelist. When this parameter is set, the following nondimensionalization is used (following [HGW16]):

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

We assume a polytropic background state (similar to dimensional anelastic mode), with gravity varying as \(\frac{1}{r^2}\). We further assume that the transport coefficients \(\nu\), \(\kappa\), and \(\eta\) do not vary with radius. The results in the nondimensionalized equations (tildes used to indicated nondimensional reference-state values):

\[\begin{split}\begin{aligned} \frac{\partial \boldsymbol{v}}{\partial t} +\boldsymbol{v}\cdot\boldsymbol{\nabla}\boldsymbol{v} %advection +2\boldsymbol{\hat{z}}\times\boldsymbol{v} &= % Coriolis \mathrm{Ra}^*\frac{r_\mathrm{max}^2}{r^2}\Theta\,\boldsymbol{\hat{r}} % buoyancy +\boldsymbol{\nabla}\left(\frac{P}{\tilde{\rho}(r)}\right) % pressure +\frac{\mathrm{E}}{\mathrm{Pm}\,\tilde{\rho}}\left(\boldsymbol{\nabla}\times\boldsymbol{B}\right)\times\boldsymbol{B} % Lorentz Force +\frac{\mathrm{E}}{\tilde{\rho(r)}}\boldsymbol{\nabla}\cdot\boldsymbol{\mathcal{D}} \;\;\; &\mathrm{Momentum}\\ % % \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{\rho}(r)\,\tilde{T}(r)\,\boldsymbol{\nabla}\Theta \right] % diffusion +Q(r) % Internal heating +\frac{\mathrm{E}\,\mathrm{Di}}{\mathrm{Ra}^*}\Phi(r,\theta,\phi) +\frac{\mathrm{Di\,E^2}}{\mathrm{Pm}^2\mathrm{R}^*}\left[\boldsymbol{\nabla}\times\boldsymbol{B}\right]^2 &\mathrm{Thermal\; Energy}\\ % Ohmic Heating % % \frac{\partial \boldsymbol{B}}{\partial t} &= \boldsymbol{\nabla}\times\left(\,\boldsymbol{v}\times\boldsymbol{B}-\frac{\mathrm{E}}{\mathrm{Pm}}\boldsymbol{\nabla}\times\boldsymbol{B}\,\right) &\mathrm{Induction} \\ % % \mathcal{D}_{ij} &= 2\tilde{\rho}(r)\left[e_{ij}-\frac{1}{3}\boldsymbol{\nabla}\cdot\boldsymbol{v}\right] &\mathrm{Viscous\; Stress\; Tensor}\\ % % \Phi(r,\theta,\phi) &= 2\,\tilde{\rho}(r)\left[e_{ij}e_{ij}-\frac{1}{3}\left(\boldsymbol{\nabla}\cdot\boldsymbol{v}\right)^2\right] &\mathrm{Viscous\; Heating} \\ % % \boldsymbol{\nabla}\cdot\left(\tilde{\rho}(r)\,\boldsymbol{v}\right)&=0 &\mathrm{Solenoidal\; Mass\; Flux}\\ \boldsymbol{\nabla}\cdot\boldsymbol{B}&=0. &\mathrm{Solenoidal\; Magnetic\; Field}\end{aligned}\end{split}\]

In the equations above, Di is the dissipation number, defined by

\[\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 anomoly \(\Theta\) should be interpreted as entropy \(s\). The symbol Ra\(^*\) is the modified Rayleigh number, given by

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

Those Reference_Namelist variables that must be set for this model are indicated in Table table_anelastic_nd. As with \(\Delta T\) in the nondimensional Boussinesq mode, the user must choose boundary conditions suitable for their definition of \(\Delta s\). As with the dimensional anelastic formulation, the background state is polytropic and is described through a polytropic index and number of density scale heights.

Note: As with the Boussinesq mode, please set the variables nu_type, kappa_type, eta_type in the Transport_Namelist.

Table. Anelastic_nd.

Variables in the Reference_Namelist that must be specified when running in nondimensional anelastic mode. In addition, reference_type=3 must also be specified.

Variable	Description
Ekman_Number	The Ekman Number E
Modified_Rayleigh_Number	The Modified Rayleigh Number Ra\(^*\)
Prandtl_Number	The Prandtl Number Pr
Magnetic_Prandtl_Number	The Magnetic Prandtl Number Pm
poly_n [0]	polytropic index (\(P\propto\rho^n\))
poly_Nrho [0]	number of density scaleheights spanning the domain

Boundary Conditions & Internal Heating

Boundary conditions are controlled through the Boundary_Conditions_Namelist. All Rayleigh simulations are run with impenetrable boundaries. These boundaries may be either no-slip or stress-free (default). If you want to employ no-slip conditions at both boundaries, set no_slip_boundaries = .true.. If you wish to set no-slip conditions at only one boundary, set no_slip_top=.true. or no_slip_bottom=.true. in the Boundary_Conditions_Namelist.

By default, magnetic boundary conditions are set to match to a potential field at each boundary.

By default, the thermal anomoly \(\Theta\) is set to a fixed value at each boundary. The upper and lower boundary-values are specified by setting T_top and T_bottom respectively in the Boundary_Conditions_Namelist. Their defaults values are 1 and 0 respectively.

Alternatively, you may specify a constant value of \(\partial\Theta/\partial r\) at each boundary. This is accomplished by setting the variables fix_dTdr_top and fix_dTdr_bottom. Values of the gradient may be enforced by setting the namelist variables dTdr_top and dTdr_bottom. Both default to a value of zero.

Generic Boundary Conditions

Boundary conditions for temperature, \(T\), and the magnetic poloidal potential, \(C\), may also be set using generic spectral input. As with initial conditions (see Generic Initial Conditions) this is done by generating a generic input file using the rayleigh_spectral_input.py script. The file output from this script can then be applied using:

• fix_Tvar_top and T_top_file (overrides T_top for a constant boundary condition) to set a fixed upper \(T\) boundary condition

• fix_dTdr_top and dTdr_top_file (overrides dTdr_top) to set a fixed upper \(T\) gradient boundary condition

• fix_Tvar_bottom and T_bottom_file (overrides T_bottom) to set a fixed lower \(T\) boundary condition

• fix_dTdr_bottom and dTdr_bottom_file (overrides dTdr_bottom) to set a fixed lower \(T\) gradient boundary condition

• fix_poloidal_top and C_top_file (overrides impose_dipole_field) to set a fixed upper \(C\) boundary condition

• fix_poloidal_bottom and C_bottom_file (overrides impose_dipole_field) to set a fixed lower \(C\) boundary condition

For example, to set a \(C\) boundary condition on both boundaries with modes (l,m) = (1,0) and (1,1) set to pre-calculated values run:

rayleigh_spectral_input.py -m 1 0 2.973662220170157 -m 1 1 0.5243368809294343+0.j -o ctop_init_bc
rayleigh_spectral_input.py -m 1 0 8.496177771914736 -m 1 1 1.4981053740840984+0.j -o cbottom_init_bc

which will generate generic spectral input files ctop_init_bc and cbottom_init_bc. Set these to be used as the boundary conditions in main_input using:

&Boundary_Conditions_Namelist
fix_poloidalfield_top = .true.
fix_poloidalfield_bottom = .true.
C_top_file = 'ctop_init_bc'
C_bottom_file = 'cbottom_init_bc'
/

This can be seen being applied in tests/generic_input.

Internal Heating

The internal heating function \(Q(r)\) is activated and described by two variables in the Reference_Namelist. These are Luminosity and heating_type. Note that these values are part of the Reference_Namelist and not the Boundary_Conditions namelist. Three heating types (0,1, and 4) are fully supported at this time. Heating type zero corresponds to no heating. This is the default.

Heating_type=1: This heating type is given by :

\[\label{eq:heating} %\frac{\partial \Theta}{\partial t}=\gamma\left( 1 -\frac{\hat{\rho}(r_\mathrm{max})\,\hat{T}(r_\mathrm{max}) }{\hat{\rho}(r)\, \hat{T}(r)} \right), Q(r)= \gamma\,\hat{\rho}(r)\, \hat{T}(r)\]

where \(\gamma\) is a normalization constant defined such that

\[%4\pi r_o^2 \hat{\rho}\hat{T}\kappa(r)\frac{\partial \Theta}{\partial r}=\mathrm{Luminosity} 4\pi\int_{r=r_\mathrm{min}}^{r=r_\mathrm{max}} Q(r)\, r^2 dr = \mathrm{Luminosity}.\]

This heating profile is particularly useful for emulating radiative heating in a stellar convection zone.

Heating_type=4: This heating type corresponds a heating that is variable in radius, but constant in energy density. Namely

\[\hat{\rho}\hat{T}\frac{\partial \Theta}{\partial t}=\gamma.\]

The constant \(\gamma\) in this case is also set by enforcing Equation eq_lum.

Note: If internal heating is used in combination with fix_dTdr_top, then the value of \(\partial\Theta/\partial r\) at the upper boundary is set by Rayleigh. Any value for dTdr_top specified in main_input is ignored. This is done to ensure consistency with the internal heating and any flux passing through the lower boundary due to the use of a fixed-flux condition. To override this behavior, set adjust_dTdr_top to .false. in the Boundary_Conditions namelist.

General Physics Controls

A number of logical variables can be used to turn certain physics on (value = .true.) or off ( value = .false.). These variables are described in Table table_logicals, with default values indicated in brackets.

Table. Logicals.

Variables in the Physical_Controls_Namelist that may be specified to control run behavior (defaults indicated in brackets)

Variable [Default value]	Description
magnetism [.false.]	Turn magnetism on or off
rotation [.false.]	Turn rotation on or off (pressure is not scaled by E when off)
lorentz_forces [.true.]	Turn Lorentz forces on or off (magnetism must be .true.)
viscous_heating [.true.]	Turn viscous heating on or off (inactive in Boussinesq mode)
ohmic_heating [.true.]	Turn ohmic heating off or on (inactive in Boussinesq mode)

Initializing a Model

A Rayleigh simulation may be initialized with a random thermal and/or magnetic field, or it may be restarted from an existing checkpoint file (see §Checkpointing for a detailed discussion of checkpointing). This behavior is controlled through the initial_conditions_namelist and the init_type and magnetic_init_type variables. The init_type variable controls the behavior of the velocity and thermal fields at initialization time. Available options are:

• init_type=-1 ; read velocity and thermal fields from a checkpoint file

• init_type=1 ; Christensen et al. (2001) case 0 benchmark init ( {\(\ell=4,m=4\)} temperature mode)

• init_type=6 ; Jones et al. (2011) steady anelastic benchmark ( {\(\ell=19,m=19\)} entropy mode)

• init_type=7 ; random temperature or entropy perturbation

• init_type=8 ; user generated temperature or entropy perturbation (see Generic Initial Conditions below)

When initializing a random thermal field, all spherical harmonic modes are independently initialized with a random amplitude whose maximum possible value is determined by the namelist variable temp_amp. The mathematical form of of this random initialization is given by

\[T(r,\theta,\phi) = \sum_\ell \sum_m c_\ell^m f(r)g(\ell)\mathrm{Y}_\ell^m(\theta,\phi),\]

where the \(c_\ell^m\)’s are (complex) random amplitudes, distributed normally within the range [-temp_amp, temp_amp]. The radial amplitude \(f(r)\) is designed to taper off to zero at the boundaries and is given by

\[f(r) = \frac{1}{2}\left[1-\mathrm{cos}\left( 2\pi\frac{r-rmin}{rmax-rmin} \right) \right].\]

The amplitude function \(g(\ell)\) concentrates power in the central band of spherical harmonic modes used in the simulation. It is given by

\[g(\ell) = \mathrm{exp}\left[ - 9\left( \frac{ 2\,\ell-\ell_\mathrm{max} }{ \ell_\mathrm{max} } \right)^2 \right],\]

which is itself, admittedly, a bit random.

When initializing using a random thermal perturbation, it is important to consider whether it makes sense to separately initialize the spherically-symmetric component of the thermal field with a profile that is in conductive balance. This is almost certainly the case when running with fixed temperature conditions. The logical namelist variable conductive_profile can be used for this purpose. It’s default value is .false. (off), and its value is ignored completely when restarting from a checkpoint. To initialize a simulation with a random temperature field superimposed on a spherically-symmetric, conductive background state, something similar to the following should appear in your main_input file:

&initial_conditions_namelist
init_type=7
temp_amp = 1.0d-4
conductive_profile=.true.
/

Alternatively, you may wish to specify an ell=0 initial thermal profile that is neither random nor conductive. To create your own profile, follow the example found in Rayleigh/examples/custom_thermal_profile/custom_thermal_profile.ipynb. Then, use the following combination of input parameters in main_input:

&initial_conditions_namelist
init_type=7
temp_amp = 1.0d-4
custom_thermal_file = 'my_custom_profile.dat'
/

This will use the radial profile stored in my_custom_profile.dat for the ell=0 component of entropy/temperature Random values will be used to initialize all other modes.

Magnetic-field initialization follows a similar pattern. Available values for magnetic_input type are:

• magnetic_init_type = -1 ; read magnetic field from a checkpoint file

• magnetic_init_type = 1 ; Christensen et al. (2001) case 0 benchmark init

• magnetic_init_type = 7 ; randomized vector potential

• magnetic_init_type=8 ; user generated magnetic potential fields (see Generic Initial Conditions below)

For the randomized magnetic field, both the poloidal and toroidal vector-potential functions are given a random power distribution described by Equation eq_init. Each mode’s random amplitude is then determined by namelist variable mag_amp. This variable should be interpreted as an approximate magnetic field strength (it’s value is rescaled appropriately for the poloidal and toroidal vector potentials, which are differentiated to yield the magnetic field).

When initializing all fields from scratch, a main_input file should contain something similar to:

&initial_conditions_namelist
init_type=7
temp_amp = 1.0d-4
conductive_profile=.true.  ! Not always necessary (problem dependent) ...
magnetic_init_type=7
mag_amp = 1.0d-1
/

Generic Initial Conditions

The user can input any initial conditions from data files generated by a python routine “rayleigh_spectral_input.py”, which can be called as a script or imported as a python class.

The available generic initial conditions options are

&initial_conditions_namelist
init_type=8
T_init_file = '<filename>'  !! Temperature
W_init_file = '<filename>'  !! Poloidal velocity potential
Z_init_file = '<filename>'  !! Toroidal velocity potential
P_init_file = '<filename>'  !! `Pressure` potential
magneic_init_type=8
C_init_file = '<filename>'  !! Poloidal magnetic potential
A_init_file = '<filename>'  !! Toroidal magnetic potential

/

where T_init_file is a user generated initial temperature field and <filename> is the name of the file generated by the python script. If T_init_file is not specified the initial field will be zero by default. The same for the other fields. Fields T, W, Z, and P are only initialized from the file if init_type=8. Fields C and A are only initialized from file if magnetic_init_type=8.

To generate a generic initial condition input file, for example, if a user wanted to specify a single mode in that input file then they could just run the script:

rayleigh_spectral_input.py -m 0 0 0 1.+0.j -o example

to specify (n,l,m) = (0,0,0) to have a coefficient 1.+0.j and output it to the file example.

This could also be done using the python as a module. In a python shell this would look like:

from rayleigh_spectral_input import *
si = SpectralInput()
si.add_mode(1., n=0, l=0, m=0)
si.write('example')

For a more complicated example, e.g. the hydrodynamic benchmark from Christensen et al. 2001, the user can specify functions of theta, phi and radius that the python will convert to spectral:

rayleigh_spectral_input.py -ar 0.35 -sd 1.0 -nt 96 -nr 64 -o example \
 -e 'import numpy as np; x = 2*radius - rmin - rmax;
 rmax*rmin/radius - rmin + 210*0.1*(1 - 3*x*x + 3*(x**4) -
 x**6)*(np.sin(theta)**4)*np.cos(4*phi)/np.sqrt(17920*np.pi)'

in “script” mode.

Alternatively, in “module” mode in a python shell:

from rayleigh_spectral_input import *
si = SpectralInput(n_theta=96, n_r=64)
rmin, rmax = radial_extents(aspect_ratio=0.35, shell_depth=1.0)
def func(theta, phi, radius):
   x = 2*radius - rmin - rmax
   return rmax*rmin/radius - rmin + 210*0.1*(1 - 3*x*x + 3*(x**4) - x**6)*(np.sin(theta)**4)*np.cos(4*phi)/np.sqrt(17920*np.pi)
si.transform_from_rtp_function(func, aspect_ratio=0.35, shell_depth=1.0)
si.write('example')

The above commands will generate a file called example which can be called by

&initial_conditions_namelist
init_type=8
T_init_file = 'example'

Note that these two examples will have produced different data formats - the first one sparse (listing only the mode specified) and the second one dense (listing all modes).

For more examples including magnetic potentials see tests/generic_input.