.. raw:: latex \clearpage .. _physics: 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: .. math:: \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} Here, :math:`\hat{\rho}` and :math:`\hat{T}` are the reference-state density and temperature respectively. :math:`g` is the gravitational acceleration, :math:`c_P` is the specific heat at constant pressure, and :math:`\Omega_0` is the frame rotation rate. The velocity field vector is denoted by :math:`\boldsymbol{v}`, the magnetic field vector by :math:`\boldsymbol{B}`, and the pressure by :math:`P`. The thermal anomoly is denoted by :math:`\Theta` and should be interpreted is as entropy :math:`s` in this formulation. The thermal variables satisfy the linearized equation of state .. math:: \frac{P}{\hat{P}}= \frac{T}{\hat{T}} + \frac{\rho}{\hat{\rho}} The kinematic viscosity, thermal diffusivity, and magnetic diffusivity are given by :math:`\nu`, :math:`\kappa`, and :math:`\eta` respectively. Finally, :math:`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: .. centered:: **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 | | | (:math:`P\propto\rho^n`) | +-----------------------------------+-----------------------------------+ | poly_Nrho [0] | number of density scaleheights | | | spanning the domain | +-----------------------------------+-----------------------------------+ | poly_mass [0] | mass interior to :math:`rmin` | +-----------------------------------+-----------------------------------+ | poly_rho_i [0] | density at rmin, | | | :math:`\rho(r=rmin)` | +-----------------------------------+-----------------------------------+ | pressure_specific_heat [0] | specific heat at constant | | | pressure | +-----------------------------------+-----------------------------------+ | angular_velocity [1.0] | cyclic frequency of the rotating | | | frame | +-----------------------------------+-----------------------------------+ .. _table_anelastic_trans: .. centered:: **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, | | | :math:`\nu(rmax)` | +-----------------------------------+-----------------------------------+ | nu_type [1] | determines whether :math:`\nu` is | | | constant with radius (1) or | | | varies with density (2) | +-----------------------------------+-----------------------------------+ | nu_power [0.0] | exponent in : | | | :math:`\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, | | | :math:`\kappa(rmax)` | +-----------------------------------+-----------------------------------+ | kappa_type [1] | determines whether :math:`\kappa` | | | is constant with radius (1) or | | | varies with density (2) | +-----------------------------------+-----------------------------------+ | kappa_power [0.0] | exponent in : | | | :math:`\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, | | | :math:`\eta(rmax)` | +-----------------------------------+-----------------------------------+ | eta_type [1] | determines whether :math:`\eta` | | | is constant with radius (1) or | | | varies with density (2) | +-----------------------------------+-----------------------------------+ | eta_power [0.0] | exponent in : | | | :math:`\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. .. raw:: latex \clearpage .. _physics_boussinesq_nondimensional: 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 .. math:: \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} where :math:`\Omega_0` is the rotation rate of the frame, :math:`\rho` is the (constant) density of the fluid, :math:`\mu` is the magnetic permeability, :math:`\eta` is the magnetic diffusivity, and :math:`\nu` is the kinematic viscosity. After nondimensionalizing, the following nondimensional numbers appear in our equations .. math:: \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} where :math:`\alpha` is coefficient of thermal expansion, :math:`g_0` is the gravitational acceleration at the top of the domain, and :math:`\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: .. math:: \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} where :math:`r_0 \equiv rmax`. In this formulation, :math:`\Theta` should be interpreted as the temperature perturbation :math:`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 :math:`\Delta T` in our equation set changes, with :math:`\Delta T \equiv L\frac{\partial T}{\partial r}` instead, for some fiducial value of :math:`\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 :math:`\Delta T`. If :math:`\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 §\ :ref:`boundary_conditions`. 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: .. centered:: **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 :math:`E` | +-----------------------------------+-----------------------------------+ | Rayleigh_Number | The Rayleigh Number :math:`Ra` | +-----------------------------------+-----------------------------------+ | Prandtl_Number | The Prandtl Number :math:`Pr` | +-----------------------------------+-----------------------------------+ | Magnetic_Prandtl_Number | The Magnetic Prandtl Number | | | :math:`Pm` | +-----------------------------------+-----------------------------------+ | Gravity_Power | Buoyancy coefficient = | | | :math:`\frac{\mathrm{Ra}}{\mathrm | | | {Pr}}\left(\frac{r}{rmax} \right) | | | ^\mathrm{gravity\_power}` | +-----------------------------------+-----------------------------------+ .. raw:: latex \clearpage 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 :cite:`Heimpel_etal2016`): .. math:: \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} We assume a polytropic background state (similar to dimensional anelastic mode), with gravity varying as :math:`\frac{1}{r^2}`. We further assume that the transport coefficients :math:`\nu`, :math:`\kappa`, and :math:`\eta` do not vary with radius. The results in the nondimensionalized equations (tildes used to indicated nondimensional reference-state values): .. math:: \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} In the equations above, Di is the dissipation number, defined by .. math:: \mathrm{Di}= \frac{g_o\,\mathrm{L}}{c_\mathrm{P}\,T_o}, where :math:`g_o` and :math:`T_o` are the gravitational acceleration and temperature at the outer boundary respectively. Once more, the thermal anomoly :math:`\Theta` should be interpreted as entropy :math:`s`. The symbol Ra\ :math:`^*` is the modified Rayleigh number, given by .. math:: \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 :math:`\Delta T` in the nondimensional Boussinesq mode, the user must choose boundary conditions suitable for their definition of :math:`\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: .. centered:: **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\ :math:`^*` | +-----------------------------------+-----------------------------------+ | Prandtl_Number | The Prandtl Number Pr | +-----------------------------------+-----------------------------------+ | Magnetic_Prandtl_Number | The Magnetic Prandtl Number Pm | +-----------------------------------+-----------------------------------+ | poly_n [0] | polytropic index | | | (:math:`P\propto\rho^n`) | +-----------------------------------+-----------------------------------+ | poly_Nrho [0] | number of density scaleheights | | | spanning the domain | +-----------------------------------+-----------------------------------+ .. _boundary_conditions: 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 :math:`\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 :math:`\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, :math:`T`, and the magnetic poloidal potential, :math:`C`, may also be set using generic spectral input. As with initial conditions (see :ref:`sec:generic_ic`) 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 :math:`T` boundary condition - **fix_dTdr_top** and **dTdr_top_file** (overrides **dTdr_top**) to set a fixed upper :math:`T` gradient boundary condition - **fix_Tvar_bottom** and **T_bottom_file** (overrides **T_bottom**) to set a fixed lower :math:`T` boundary condition - **fix_dTdr_bottom** and **dTdr_bottom_file** (overrides **dTdr_bottom**) to set a fixed lower :math:`T` gradient boundary condition - **fix_poloidal_top** and **C_top_file** (overrides **impose_dipole_field**) to set a fixed upper :math:`C` boundary condition - **fix_poloidal_bottom** and **C_bottom_file** (overrides **impose_dipole_field**) to set a fixed lower :math:`C` boundary condition For example, to set a :math:`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 :math:`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 : .. math:: \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 :math:`\gamma` is a normalization constant defined such that .. _eq_lum: .. math:: %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 .. math:: \hat{\rho}\hat{T}\frac{\partial \Theta}{\partial t}=\gamma. The constant :math:`\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 :math:`\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: .. centered:: **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 §\ :ref:`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 ( {:math:`\ell=4,m=4`} temperature mode) - init_type=6 ; Jones et al. (2011) steady anelastic benchmark ( {:math:`\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 .. _eq_init: .. math:: T(r,\theta,\phi) = \sum_\ell \sum_m c_\ell^m f(r)g(\ell)\mathrm{Y}_\ell^m(\theta,\phi), where the :math:`c_\ell^m`\ ’s are (complex) random amplitudes, distributed normally within the range [-temp_amp, temp_amp]. The radial amplitude :math:`f(r)` is designed to taper off to zero at the boundaries and is given by .. math:: f(r) = \frac{1}{2}\left[1-\mathrm{cos}\left( 2\pi\frac{r-rmin}{rmax-rmin} \right) \right]. The amplitude function :math:`g(\ell)` concentrates power in the central band of spherical harmonic modes used in the simulation. It is given by .. math:: 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 / .. _sec:generic_ic: 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 = '' !! Temperature W_init_file = '' !! Poloidal velocity potential Z_init_file = '' !! Toroidal velocity potential P_init_file = '' !! `Pressure` potential magneic_init_type=8 C_init_file = '' !! Poloidal magnetic potential A_init_file = '' !! Toroidal magnetic potential / where `T_init_file` is a user generated initial temperature field and 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`.