# 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:

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

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

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

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:

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]):

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):

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

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

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 :

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

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

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

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

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

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.
/
```

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.