Rayleigh User Guide

Rayleigh_User_Guide

User Manual:

Open the PDF directly: View PDF .
Page Count: 22

Rayleigh User Guide

Version 0.9.1

Nicholas Featherstone

April 30, 2018

Overview

Rayleigh solves the MHD equations, in a rotating frame, within spherical shells, using the anelastic or Boussi-

nesq approximations. Derivatives in Rayleigh are calculated using a spectral transform scheme. Spherical har-

monics are used as basis functions in the horizontal direction. Chebyshev polynomials are employed in radius.

Time-stepping is accomplished used the semi-implicit Crank-Nicolson method for the linear terms, and the Adams-

Bashforth method for the nonlinear terms. Both methods are second-order in time.

This document serves as a guide to installation and running Rayleigh. Rayleigh’s diagnostics package is dis-

cussed in the companion document Rayleigh/doc/Diagnostics Plotting.html,pdf

Rayleigh is written by Nicholas Featherstone, with National-Science-Foundation support through the Geody-

namo Working Group of the Computational Infrastructure for Geodynamics (CIG; PI: Louise Kellogg).

The CIG Geodynamo Working Group Members are: Jonathon Aurnou, Benjamin Brown, Bruce Buﬀet, Nicholas

Featherstone, Gary Glatzmaier, Eric Heien, Moritz Heimpel ,Lorraine Hwang, Louise Kellogg, Hiroaki Matsui,

Peter Olson, Krista Soderlund, Sabine Stanley, Rakesh Yadav

Rayleigh’s implementation of the pseudo-spectral algorithm and its parallel design would not have been pos-

sible without earlier work by Gary Glatzmaier and Thomas Clune, described in:

1. Glatzmaier, G.A., 1984, J. Comp. Phys., 55, 461

2. Clune, T.C., Elliott, J.R., Miesch, M.S. & Toomre, J., 1999, Parallel Comp.,25, 361

Contents

1 Compiling and Installing Rayleigh 3

1.1 Third-PartyDependencies........................................... 3

1.2 Compilation................................................... 3

1.3 VerifyingYourInstallation........................................... 3

2 Running the Code 4

2.1 Preparation................................................... 4

2.2 Code Execution and Load-Balancing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4

2.2.1 LoadBalancing............................................. 5

2.2.2 Specifying Resolution & Domain Bounds . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

2.3 Controlling Run Length & Time Stepping . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

3 Running a Benchmark 7

4 Physics Controls 9

4.1 AnelasticMode(dimensional)......................................... 9

4.2 Boussinesq Mode (nondimensional) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

4.3 AnelasticMode(nondimensional)....................................... 13

4.4 Boundary Conditions & Internal Heating . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

4.4.1 InternalHeating ............................................ 14

4.5 GeneralPhysicsControls ........................................... 14

4.6 InitializingaModel............................................... 14

5 Checkpointing 17

5.1 StandardCheckpoints ............................................. 17

5.1.1 Generating Standard Checkpoints . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

5.1.2 Restarts From Standard Checkpoints . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

5.2 Quicksaves ................................................... 18

5.2.1 1.2.1 Generating Quicksaves . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

5.2.2 1.2.2 Restarting from Quicksaves . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

5.3 1.3CheckpointLogs .............................................. 19

6 Diagnostic Outputs 19

7 I/O Redirection 20

8 Ensemble Mode 21

9 Acknowledging Rayleigh 22

1 Compiling and Installing Rayleigh

A detailed explanation of the installation process may be found in the root directory of the code repository at:

Rayleigh/INSTALL.

We provide an abbreviated version of those instructions here.

1.1 Third-Party Dependencies

In order to compile Rayleigh, you will need to have MPI (Message Passing Interface) installed along with a Fortran

2003-compliant compiler. Rayleigh has been successfully compiled with GNU, Intel, and IBM compilers (PGI has

not been tested yet). Rayleigh’s conﬁgure script provides native support for the Intel and GNU compilers. See

Rayleigh\INSTALL for an example of conﬁguration using the IBM compiler.

Rayleigh depends on the following third party libraries:

1. BLAS (Basic Linear Algebra Subprograms)

2. LAPACK (Linear Algebra PACKage)

3. FFTW 3.x (Fastest Fourier Transform in the West)

You will need to install these libraries before compiling Rayleigh. If you plan to run Rayleigh on Intel processors,

we suggest installing Intel’s Math Kernel Library (MKL) in lieu of installing these libraries individually. The Math

Kernel Library provides optimized versions of BLAS, LAPACK, and FFTW. It has been tuned, by Intel, for optimal

performance on Intel processors. At the time of this writing, MKL is provided free of charge. You may ﬁnd it here.

1.2 Compilation

Rayleigh is compiled using the standard Linux installation scheme of conﬁgure/make/make-install. From within the

Rayleigh directory, run these commands:

1. ./conﬁgure – See Rayleigh/INSTALL or run ./conﬁgure --help to see relevant options.

2. make – This compiles the code

3. make install – This places the Rayleigh executables in Rayleigh/bin. If you would like to place them in (say)

/home/my rayleigh, run conﬁgure as: ./conﬁgure –preﬁx=/home/my rayleigh.

For most builds, two executables will be created: rayleigh.opt and rayleigh.dbg. Use them as follows:

1. When running production jobs, use rayleigh.opt.

2. If you encounter an unexpected crash and would like to report the error, rerun the job with rayleigh.dbg. This

version of the code is compiled with debugging symbols. It will (usually) produces meaningful error messages

in place of the gibberish that is output when rayleigh.opt crashes.

If conﬁgure detects the Intel compiler, you will be presented with a number of choices for the vectorization option. If

you select all, rayleigh.opt will not be created. Instead, rayleigh.sse, rayleigh.avx, etc. will be placed in Rayleigh/bin.

This is useful if running on a machine with heterogeneous node architectures (e.g., Pleiades). If you are not running

on such a machine, pick the appropriate vectorization level, and rayleigh.opt will be compiled using that vectorization

level.

1.3 Verifying Your Installation

Rayleigh comes with a benchmarking mode that help you verify that the installation is performing correctly. If you

are running Rayleigh for the ﬁrst time, or running on a new machine, follow along with the example in §3, and verify

that you receive an accurate benchmark report before running a custom model.

2 Running the Code

Whenever you run a new simulation, a similar series of steps must be performed. A summary of the typical Rayleigh

work ﬂow is:

1. Create a unique directory for storing simulation output

2. Create a main input ﬁle

3. Copy or soft link the Rayleigh executable into the simulation directory

4. Modify main input as desired

5. Run the code

6. Examine output and restart simulation as necessary

2.1 Preparation

Each simulation run using Rayleigh should have its own directory. The code is run from within that directory, and

any output is stored in various subdirectories created by Rayleigh at run time. Wherever you create your simulation

directory, ensure that you have suﬃcient space to store the output.

Do not run Rayleigh from within the source code directory.

Do not cross the beams: no running two models from within the same directory.

After you create your run directory, you will want to copy (cp) or soft link (ln -s ) the executable from Rayleigh/bin

to your run directory. Soft-linking is recommended; if you recompile the code, the executable remains up-to-date. If

running on an IBM machine, copy the script named Rayleigh/etc/make dirs to your run directory and execute the

script. This will create the directory structure expected by Rayleigh for its outputs. This step is unnecessary when

compiling with the Intel or GNU compilers.

Next, you must create a main input ﬁle. This ﬁle contains the information that describes how your simulation is

run. Rayleigh always looks for a ﬁle named main input in the directory that it is launched from. Copy one of the

sample input ﬁles from the Rayleigh/etc/input examples/ into your run directory, and rename it to main input. The

ﬁle named benchmark diagnostics input can be used to generate output for the diagnostics plotting tutorial (see §6).

Finally, Rayleigh has some OpenMP-related logic that is still in development. We do not support Rayleigh’s OpenMP

mode at this time, but on some systems, it can be important to explicitly disable OpenMP in order to avoid tripping

any OpenMP ﬂags used by external libraries, such as Intel’s MKL. Please be sure and run the following command

before executing Rayleigh. This command should be precede each call to Rayleigh.

export OMP_NUM_THREADS=1 (bash)

setenv OMP_NUM_THREADS 1 (c-shell)

2.2 Code Execution and Load-Balancing

Rayleigh is parallelized using MPI and a 2-D domain decomposition. The 2-D domain decomposition means that we

envision the MPI Ranks as being distributed in rows and columns. The number of MPI ranks within a row is nprow

and the number of MPI ranks within a column is npcol. When Rayleigh is run with N MPI ranks, the following

constraint must be satisﬁed:

N = npcol ×nprow.(1)

If this constraint is not satisﬁed , the code will print an error message and exit. The values of nprow and npcol can

be speciﬁed in main input or on the command line via the syntax:

mpiexec -np 8 ./rayleigh.opt -nprow 4 -npcol 2

2.2.1 Load Balancing

Rayleigh’s performance is sensitive to the values of nprow and npcol, as well as the number of radial grid points

Nrand latitudinal grid points Nθ. If you examine the main input ﬁle, you will see that it is divided into Fortran

namelists. The ﬁrst namelist is the problemsize namelist. Within this namelist, you will see a place to specify nprow

and npcol. Edit main input so that nprow and npcol agree with the N you intend to use (or use the command-line

syntax mentioned above). The dominate eﬀect on parallel scalability is the number of messages sent per iteration.

For optimal message counts, nprow and npcol should be as close to one another in value as possible.

1. N = nprow ×npcol.

2. nprow and npcol should be equal or within a factor of two of one another.

The value of nprow determines how spherical harmonics are distributed across processors. Spherical harmonics are

distributed in high-m/low-mpairs, where mis the azimuthal wavenumber. Each process is responsible for all `-values

associated with those m’s contained in memory.

The value of npcol determines how radial levels are distributed across processors. Radii are distributed uniformly

across processes in contiguous chunks. Each process is responsible for a range of radii ∆r.

The number of spherical harmonic degrees N`is deﬁned by

N`=2

3Nθ(2)

For optimal load-balancing, nprow should divide evenly into Nrand npcol should divide evenly into the number of

high-m/low-mpairs (i.e., N`/2). Both nprow and npcol must be at least 2.

In summary,

1. nprow ≥2.

2. npcol ≥2.

3. n×npcol =Nr(for integer n).

4. k×npcol =1

3Nθ(for integer k).

2.2.2 Specifying Resolution & Domain Bounds

As discussed, the number of radial grid points is denoted by Nr, and the number of θgrid points by Nθ. The number

of grid points in the φdirection is always Nφ= 2 ×Nθ.Nrand Nθmay each be deﬁned in the problemsize namelist

of main input:

&problemsize_namelist

n_r = 48

n_theta = 96

Nrand Nθmay also be speciﬁed at the command line (overriding the values in main input) via:

mpiexec -np 8 ./rayleigh.opt -nr 48 -ntheta 96

If desired, the maximal spherical harmonic degree `max ≡N`−1 can be speciﬁed in lieu of Nθ. The example above

may equivalently be written:

&problemsize_namelist

n_r = 48

l_max = 63

The radial domain bounds are determined by the namelist variables rmin (the lower radial boundary) and rmax

(the upper radial boundary):

&problemsize_namelist

rmin = 1.0

rmax = 2.0

Alternatively, the user may specify the shell depth (rmax −rmin) and aspect ratio (rmin/rmax) in lieu of rmin

and rmax. The preceding example may then be written as:

&problemsize_namelist

aspect_ratio = 0.5

shell_depth = 1.0

Note that the interpretation of rmin and rmax depends on whether your simulation is dimensional or nondimensional.

We discuss these alternative formulations in §4

2.3 Controlling Run Length & Time Stepping

A simulation’s runtime and time-step size can be controlled using the temporal controls namelist. The length

of time for which a simulation runs before completing is controlled by the namelist variable max time minutes.

The maximum number of time steps that a simulation will run for is determined by the value of the namelist

max iterations. The simulation will complete when it has run for max time minutes minutes or when it has run

for max iterations time steps – whichever occurs ﬁrst.

Time-step size in Rayleigh is controlled by the Courant-Friedrichs-Lewy condition (CFL; as determined by the ﬂuid

velocity and Alfv´en speed). A safety factor of cﬂmax is applied to the maximum time step determined by the CFL.

Time-stepping is adaptive. An additional variable cﬂmin is used to determine if the time step should be increased.

The user may also specify the maximum allowed time-step size through the namelist variable max time step. The

minimum allowable time-step size is controlled through the variable min time step. If the CFL condition is less

than this value, the simulation will exit.

Let ∆tbe the current time-step size, and let tCFL be the maximum time-step size as determined by the CFL limit.

The following logic is employed by Rayleigh when calculating the time-step size:

•IF {∆t≥cﬂmax ×tCFL }THEN {∆tis set to cﬂmax ×tCFL }.

•IF {∆t≤cﬂmin ×tCFL }THEN {∆tis set to cﬂmax ×tCFL }.

•IF{tCFL ≥max time step }THEN {∆tis set to max time step }

•IF{tCFL ≤min time step }THEN {Rayleigh Exits }

The default values for these variables are:

&temporal_controls_namelist

max_iterations = 1000000

max_time_minutes = 1d8

cflmax = 0.6d0

cflmin = 0.4d0

max_time_step = 1.0d0

min_time_step = 1.0d-13

3 Running a Benchmark

Rayleigh has been programmed with internal testing suite so that its results may be compared against benchmarks

described in the following two papers:

1. Christensen, U.R., et al. 2001, A Numerical Dynamo Benchmark,PEPI, 128, 25

2. Jones, C.A., et al., 2011, Anelastic-Convective-Driven Dynamo Benchmarks,Icarus, 216, 120

We recommend running a benchmark whenever running Rayleigh on a new machine for the ﬁrst time, or after

recompiling the code. The Christensen et al. (2001) reference describes two Boussinesq tests that Rayleigh’s results

may be compared against. The Jones et al. (2011) reference describes anelastic tests. Rayleigh has been tested

successfully against two benchmarks from each of these papers. Input ﬁles for these diﬀerent tests are enumerated

in Table 1 below. In addition to the input ﬁles listed in Table 1, input examples appropriate for use as a template

for new runs are provided with the input suﬃx (as opposed to the minimal suﬃx. These input ﬁles still have

benchmark mode active. Be sure to turn this ﬂag oﬀ if not running a benchmark.

Important: If you are not running a benchmark, but only wish to modify an existing benchmark-input ﬁle, delete

the line containing the text “benchmark mode=X.” When benchmark mode is active, custom inputs, such as Rayleigh

number, are overridden and reset to their benchmark-appropriate values.

We suggest using the c2001 case0 minimal input ﬁle for installation veriﬁcation. Algorithmically, there

is little diﬀerence between the MHD, non-MHD, Boussinesq, and anelastic modes of Rayleigh. As a result, when

installing the code on a new machine, it is normally suﬃcient to run the cheapest benchmark, case 0 from Christensen

2001.

To run this benchmark, create a directory from within which to run your benchmark, and follow along with the

commands below. Modify the directory structure a each step as appropriate:

1. mkdir path to my sim

2. cd path to my sim

3. cp path to rayleigh/Rayleigh/etc/input examples/c2001 case0 minimal main input

4. cp path to rayleigh/Rayleigh/bin/rayleigh.opt rayleigh.opt (or use ln -s in lieu of cp)

5. mpiexec -np N./rayleigh.opt -nprow X-npcol Y-nr R-ntheta T

For the value N, select the number of cores you wish to run with. For this short test, 32 cores is more than suﬃcient.

Even with only four cores, the lower-resolution test suggested below will only take around half an hour. The values

Xand Yare integers that describe the process grid. They should both be at least 2, and must satisfy the expression

N=X×Y. (3)

Some suggested combinations are {N,X,Y}={32,4,8},{16,4,4},{8,2,4},{4,2,2}. The values Rand Tdenote the

number of radial and latitudinal collocation points respectively. Select either {R,T}={48,64}or {R,T}={64,96}.

The lower-resolution case takes about 3 minutes to run on 32 Intel Haswell cores. The higher-resolution case takes

about 12 minutes to run on 32 Intel Haswell cores.

Once your simulation has run, examine the ﬁle path to my sim/Benchmark Reports/00025000. You should see

output similar to that presented in Tables 2 or 3. Your numbers may diﬀer slightly, but all values should have a %

diﬀerence of less than 1. If this condition is satisﬁed, your installation is working correctly.

Paper Benchmark Input File

Christensen et al. 2001 Case 0 Rayleigh/etc/input examples/c2001 case0 minimal

Christensen et al. 2001 Case 1(MHD) Rayleigh/etc/input examples/c2001 case1 minimal

Jones et al. 2011 Steady Hydro Rayleigh/etc/input examples/j2011 steady hydro minimal

Jones et al. 2011 Steady MHD Rayleigh/etc/input examples/j2011 steady MHD minimal

Table 1: Benchmark-input examples useful for verifying Rayleigh’s installation. Those from Christensen et al.

(2001) are Boussinesq. Those from Jones et al. (2011) are anelastic.

Observable Measured Suggested % Diﬀerence Std. Dev.

Kinetic Energy 58.347827 58.348000 -0.000297 0.000000

Temperature 0.427416 0.428120 -0.164525 0.000090

Vphi -10.118053 -10.157100 -0.384434 0.012386

Drift Frequency 0.183272 0.182400 0.477962 0.007073

Table 2: Rayleigh benchmark report for Christensen et al. (2001) case 0 when run with nr=48 and ntheta=64. Run

time was approximately 3 minutes when run on 32 Intel Haswell cores.

Run command: mpiexec -np 32 ./rayleigh.opt -nprow 4 -npcol 8 -nr 48 -ntheta 64

Observable Measured Suggested % Diﬀerence Std. Dev.

Kinetic Energy 58.347829 58.348000 -0.000294 0.000000

Temperature 0.427786 0.428120 -0.077927 0.000043

Vphi -10.140183 -10.157100 -0.166551 0.005891

Drift Frequency 0.182276 0.182400 -0.067994 0.004877

Table 3: Rayleigh benchmark report for Christensen et al. (2001) case 0 when run with nr=64 and ntheta=96. Run

time was approximately 12 minutes when run on 32 Intel Haswell cores.

Run command: mpiexec -np 32 ./rayleigh.opt -nprow 4 -npcol 8 -nr 64 -ntheta 96

4 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.

4.1 Anelastic Mode (dimensional)

Example Input: Rayleigh/etc/input examples/main input sun

When run in dimensional, anelastic mode, reference type=2 must be speciﬁed in the Reference Namelist. In that

case, Rayleigh solves the following form of the MHD equations:

ˆρ(r)∂v

∂t +v·∇v+ 2Ω0ˆz×v=ˆρ(r)

g(r)Θ ˆr+ ˆρ(r)∇P

ˆρ(r)+1

4π(∇×B)×B+∇·DMomentum

ˆρ(r)ˆ

T(r)∂Θ

∂t +v·∇Θ=∇·hˆρ(r)ˆ

T(r)κ(r)∇Θi+Q(r) + Φ(r, θ, φ) + η(r)

4π[∇×B]2Thermal Energy

∂B

∂t =∇×(v×B−η(r)∇×B) Induction

Dij = 2ˆρ(r)ν(r)eij −1

3∇·vViscous Stress Tensor

Φ(r, θ, φ) = 2 ˆρ(r)ν(r)eij eij −1

3(∇·v)2Viscous Heating

∇·(ˆρ(r)v) = 0 Solenoidal Mass Flux

∇·B= 0 Solenoidal Magnetic Field

Here, ˆρand ˆ

Tare the reference-state density and temperature respectively. gis the gravitational acceleration, cP

is the speciﬁc heat at constant pressure, and Ω0is the frame rotation rate. The velocity ﬁeld vector is denoted by

v, the magnetic ﬁeld vector by B, and the pressure by P. The thermal anomoly is denoted by Θ and should be

interpreted is as entropy sin this formulation. The thermal variables satisfy the linearized equation of state

P=T

T+ρ

ˆρ(4)

The kinematic viscosity, thermal diﬀusivity, and magnetic diﬀusivity are given by ν,κ, and η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 speciﬁed. The

Reference Namelist variables are described in Table 4 and the Transport Namelist variables are described in Table

5. Default values indicated in brackets.

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 ﬁle main input sun for a an example of how to run a solar-like model using Rayleigh’s

dimensional, anelastic formulation.

Variable Description

poly n [0] polytropic index (P∝ρ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, ρ(r=rmin)

pressure speciﬁc heat [0] speciﬁc heat at constant pressure

angular velocity [1.0] cyclic frequency of the rotating frame

Table 4: Variables in the Reference Namelist that must be speciﬁed when running in dimensional anelastic mode.

In addition, reference type=2 must also be speciﬁed.

Variable Description

nu top [1.0] kinematic viscosity at rmax, ν(rmax)

nu type [1] determines whether νis constant with radius (1) or varies with density (2)

nu power [0.0] exponent in : ν(r) = ρ(r)

ρ(r=rmax)nu power; use with nu type=2

kappa top [1.0] thermal diﬀusivity at rmax, κ(rmax)

kappa type [1] determines whether κis constant with radius (1) or varies with density (2)

kappa power [0.0] exponent in : κ(r) = ρ(r)

ρ(r=rmax)kappa power; use with kappa type=2

eta top [1.0] magnetic diﬀusivity at rmax, η(rmax)

eta type [1] determines whether ηis constant with radius (1) or varies with density (2)

eta power [0.0] exponent in : η(r) = ρ(r)

ρ(r=rmax)eta power; use with eta type=2

Table 5: Variables in the Transport Namelist that must be speciﬁed when running in dimensional anelastic mode.

In addition, reference type=2 must also be speciﬁed in the Reference Namelist.

4.2 Boussinesq Mode (nondimensional)

Example Input: Rayleigh/etc/input examples/c2001 case1 input

When run in nondimensional Boussinesq mode, reference type=1 must be speciﬁed in the Reference Namelist. In

that case, Rayleigh employs the nondimensionalization

Length →L(Shell Depth)

Time →L2

ν(Viscous Timescale)

Temperature →∆T(Temperature Contrast Across Shell)

Magnetic Field →pρµηΩ0,

where Ω0is the rotation rate of the grame, ρis the (constant) density of the ﬂuid, µis the magnetic permeability, η

is the magnetic diﬀusivity, and νis the kinematic viscosity. After nondimensionalizing, the following nondimensional

numbers appear in our equations

P r =ν

κPrandtl Number

P m =ν

ηMagnetic Prandtl Number

E=ν

Ω0L2Ekman Number

Ra =αg0∆T L3

νκ Rayleigh Number,

where αis coeﬃcient of thermal expansion, g0is the gravitational acceleration at the top of the domain, and κis

the thermal diﬀusivity.

In addition, ohmic and viscous heating, which do not appear in the Boussinesq formulation, are turned oﬀ when this

nondimensionalization is speciﬁed at runtime. Rayleigh solves the following equations when running in nondimen-

sional Boussinesq mode:

∂v

∂t +v·∇v+2

Eˆz×v=Ra

P r r

ron

Θˆr−1

E∇P+1

E P m (∇×B)×B+∇2vMomentum

∂Θ

∂t +v·∇Θ=1

P r ∇2Θ Thermal Energy

∂B

∂t =∇×(v×B) + 1

P m ∇2BInduction

∇·v= 0 Solenoidal Velocity Field

∇·B= 0 Solenoidal Magnetic Field,

where r0≡rmax. In this formulation, Θ should be interpreted as the temperature perturbation T. Those Refer-

ence Namelist variables that must be set for this model are indicated in Table 6.

Note that our choice for the temperature scale assumes ﬁxed-temperature boundary conditions. We might choose to

specify ﬁxed-ﬂux boundary conditions and/or an internal heating, in which case the meaning of ∆Tin our equation

set changes, with ∆T≡L∂T

∂r instead, for some ﬁducial value of ∂T

∂r . Which regard to the temperature scaling, it is up

to the user to select boundary conditions appropriate for their desired values of ∆T. If ∆Tdenotes 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.

Variable Description

Ekman Number The Ekman Number E

Rayleigh Number The Rayleigh Number Ra

Prandtl Number The Prandtl Number P r

Magnetic Prandtl Number The Magnetic Prandtl Number P m

Gravity Power Buoyancy coeﬃcient = Ra

Pr r

rmax gravity power

Table 6: Variables in the Reference Namelist that must be speciﬁed when running in nondimensional Boussinesq

mode. In addition, reference type=1 must also be speciﬁed.

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 §4.4. 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.

4.3 Anelastic Mode (nondimensional)

Example Input: Rayleigh/etc/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 Heimpel et al., 2016, Nat. Geo, 9, 19):

Length →L(Shell Depth)

Time →1

Ω0

(Rotational Timescale)

Temperature →To≡ˆ

T(rmax) (Reference −State Temperature at Upper Boundary)

Density →ρo≡ˆρ(rmax) (Reference −State Density at Upper Boundary)

Entropy →∆s(Entropy Constrast Across Shell)

Magnetic Field →pρoµηΩ0.

We assume a polytropic background state (similar to dimensional anelastic mode), with gravity varying as 1

r2. We fur-

ther assume that the transport coeﬃcients ν,κ, and ηdo not vary with radius. The results in the nondimensionalized

equations (tildes used to indicated nondimensional reference-state values):

∂v

∂t +v·∇v+ 2ˆz×v= Ra∗r2

max

r2Θˆr+∇P

˜ρ(r)+E

Pm ˜ρ(∇×B)×B+E

ρ(r)

∇·DMomentum

˜ρ(r)˜

T(r)∂Θ

∂t +v·∇Θ=E

Pr∇·h˜ρ(r)˜

T(r)∇Θi+Q(r) + E Di

Ra∗Φ(r, θ, φ) + Di E2

Pm2R∗[∇×B]2Thermal Energy

∂B

∂t =∇×v×B−E

Pm∇×BInduction

Dij = 2˜ρ(r)eij −1

3∇·vViscous Stress Tensor

Φ(r, θ, φ) = 2 ˜ρ(r)eij eij −1

3(∇·v)2Viscous Heating

∇·(˜ρ(r)v) = 0 Solenoidal Mass Flux

∇·B= 0.Solenoidal Magnetic Field

In the equations above, Di is the dissipation number, deﬁned by

Di = goL

cPTo

,(5)

where goand Toare the gravitational acceleration and temperature at the outer boundary respectively. Once more,

the thermal anomoly Θ should be interpreted as entropy s. The symbol Ra∗is the modiﬁed Rayleigh number, given

Ra∗=go

cPΩ2

∆s

L(6)

Those Reference Namelist variables that must be set for this model are indicated in Table 7. As with ∆Tin the

nondimensional Boussinesq mode, the user must choose boundary conditions suitable for their deﬁnition of ∆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 Trans-

port Namelist.

4.4 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

Variable Description

Ekman Number The Ekman Number E

Modiﬁed Rayleigh Number The Modiﬁed Rayleigh Number Ra∗

Prandtl Number The Prandtl Number Pr

Magnetic Prandtl Number The Magnetic Prandtl Number Pm

poly n [0] polytropic index (P∝ρn)

poly Nrho [0] number of density scaleheights spanning the domain

Table 7: Variables in the Reference Namelist that must be speciﬁed when running in nondimensional anelastic

mode. In addition, reference type=3 must also be speciﬁed.

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.

Magnetic boundary conditions are set to match to a potential ﬁeld at each boundary. There are no supported

alternatives at this time.

By default, the thermal anomoly Θ is set to a ﬁxed value at each boundary. The upper and lower boundary-values

are speciﬁed 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 ∂Θ/∂r at each boundary. This is accomplished by setting the

variables ﬁx dTdr top and ﬁx 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.

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

Q(r) = γˆρ(r)ˆ

T(r) (7)

where γis a normalization constant deﬁned such that

4πZr=rmax

r=rmin

Q(r)r2dr = Luminosity.(8)

This heating proﬁle 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

ˆρˆ

T∂Θ

∂t =γ. (9)

The constant γin this case is also set by enforcing Equation 8.

4.5 General Physics Controls

A number of logical variables can be used to turn certain physics on (value = .true.) or oﬀ ( value = .false.). These

variables are described in Table 8, with default values indicated in brackets.

4.6 Initializing a Model

A Rayleigh simulation may be initialized with a random thermal and/or magnetic ﬁeld, or it may be restarted from

an existing checkpoint ﬁle (see §5 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 ﬁelds at initialization time. Available options are:

Variable Description

magnetism [.false.] Turn magnetism on or oﬀ

rotation [.false.] Turn rotation on or oﬀ (pressure is not scaled by E when oﬀ)

lorentz forces [.true.] Turn Lorentz forces on or oﬀ (magnetism must be .true.)

viscous heating [.true.] Turn viscous heating on or oﬀ (inactive in Boussinesq mode)

ohmic heating [.true.] Turn ohmic heating oﬀ or on (inactive in Boussinesq mode)

Table 8: Variables in the Physical Controls Namelist that may be speciﬁed to control run behavior (defaults

indicated in brackets)

•init type=-1 ; read velocity and thermal ﬁelds from a checkpoint ﬁle

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

•init type=6 ; Jones et al. (2011) steady anelastic benchmark ( {`= 19, m = 19}entropy mode)

•init type=7 ; random temperature or entropy perturbation

When initializing a random thermal ﬁeld, 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, θ, φ) = X

`f(r)g(`)Ym

`(θ, φ),(10)

where the cm

`’s are (complex) random amplitudes, distributed normally within the range [-temp amp, temp amp].

The radial amplitude f(r) is designed to taper oﬀ to zero at the boundaries and is given by

f(r) = 1

21−cos 2πr−rmin

rmax −rmin.(11)

The amplitude function g(`) concentrates power in the central band of spherical harmonic modes used in the simu-

lation. It is given by

g(`) = exp "−92`−`max

`max 2#,(12)

which is itself, admittedly, a bit random.

When initializing using a random thermal perturbation, it is important to consider whether it makes sense to sepa-

rately initialize the spherically-symmetric component of the thermal ﬁeld with a proﬁle that is in conductive balance.

This is almost certainly the case when running with ﬁxed temperature conditions. The logical namelist variable

conductive proﬁle can be used for this purpose. It’s default value is .false. (oﬀ), and its value is ignored com-

pletely when restarting from a checkpoint. To initialize a simulation with a random temperature ﬁeld superimposed

on a spherically-symmetric, conductive background state, something similar to the following should appear in your

main input ﬁle:

&initial_conditions_namelist

init_type=7

temp_amp = 1.0d-4

conductive_profile=.true.

Magnetic-ﬁeld initialization follows a similar pattern. Available values for magnetic input type are:

•magnetic init type = -1 ; read magnetic ﬁeld from a checkpoint ﬁle

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

•magnetic init type = 7 ; randomized vector potential

For the randomized magnetic ﬁeld, both the poloidal and toroidal vector-potential functions are given a random

power distribution described by Equation 10. Each mode’s random amplitude is then determined by namelist

variable mag amp. This variable should be interpreted as an approximate magnetic ﬁeld strength (it’s value is

rescaled appropriately for the poloidal and toroidal vector potentials, which are diﬀerentiated to yield the magnetic

ﬁeld).

When initializing all ﬁelds from scratch, a main input ﬁle 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

Filename Contents

00010000 W Poloidal Stream function (at time step 10,000)

00010000 Z Toroidal Stream function

00010000 P Pressure

00010000 S Entropy

00010000 C Poloidal Vector Potential

00010000 A Toroidal Vector Potential

00010000 WAB Adams-Bashforth (A-B) terms for radial momentum (W) equation

00010000 ZAB A-B terms for radial vorticity (Z) equation

00010000 PAB A-B terms for horizontal divergence of momentum (dWdr) equation

00010000 SAB A-B terms for Entropy equation

00010000 CAB A-B terms for C-equation

00010000 AAB A-B terms for A-equation

00010000 grid etc grid and time-stepping info

Table 9: Example checkpoint ﬁle collection for a time step 10,000. File contents are indicated.

5 Checkpointing

We refer to saved states in Rayleigh as checkpoints. A single checkpoint consists of 13 ﬁles when magnetism is

activated, and 9 ﬁles when magnetism is turned oﬀ. A checkpoint written at time step Xcontains all information

needed to advance the system to time step X+1. Checkpoint ﬁlenames end with a suﬃx indicating the contents of

the ﬁle (see Table 9). Each checkpoint ﬁlename possess a preﬁx as well. Files belonging to the same checkpoint share

the same preﬁx. A checkpoint ﬁle collection, written at time step 10,000 would look like that shown in Table 9.

These ﬁles contain all information needed to advance the system state from time step 10,000 to time step 10,001.

Checkpoints come in two ﬂavors, denoted by two diﬀerent preﬁx conventions: standard checkpoints and quick-

saves. This section discusses how to generate and restart from both types of checkpoints.

5.1 Standard Checkpoints

Standard checkpoints are intended to serve as regularly spaced restart points for a given run. These ﬁles begin with

an 8-digit preﬁx indicating the time step at which the checkpoint was created.

5.1.1 Generating Standard Checkpoints

The frequency with which standard checkpoints are generated can be controlled by modifying the checkpoint interval

variable in the temporal controls namelist. For example, if you want to generate a checkpoint once every 50,000

time steps, you would modify your main input ﬁle to read:

&temporal_controls_namelist

checkpoint_interval = 50000 ! Checkpoint every 50,000 time steps

The default value of checkpoint interval is 1,000,000, which is typically much larger than what you will use in practice.

5.1.2 Restarts From Standard Checkpoints

Restarting from a checkpoint is accomplished by ﬁrst assigning a value of -1 to the variables init type and/or

magnetic init type in the initial conditions namelist. In addition, the time step from which you wish to

restart from should be speciﬁed using the restart iter variable in the initial conditions namelist. The example

below will restart both the magnetic and hydrodynamic variables using the set of checkpoint ﬁles beginning with the

preﬁx 00005000.

&initial_conditions_namelist

init_type = -1 !Restart magnetic and hydro variables from time step 5,000

magnetic_init_type = -1

restart_iter = 5000

When both values are set to -1, hydrodynamic and magnetic variables are read from the relevant checkpoint ﬁles.

Alternatively, magnetic and hydrodynamic variables may be initialized separately. This allows you to add magnetism

to an already equilibrated hydrodynamic case, for instance. The example below will initialize the system with a

random magnetic ﬁeld, but it will read the hydrodynamic information (W,Z,S,P) from a checkpoint created at time

step 50,000:

&initial_conditions_namelist

init_type = -1 ! Restart hydro from time step 5,000

magnetic_init_type = 7 ! Add a random magnetic field

restart_iter = 5000

In addition to specifying the checkpoint time step manually, you can tell Rayleigh to simply restart using the last

checkpoint written by assigning a value of zero to init type:

&initial_conditions_namelist

init_type = -1

magnetic_init_type = 7

restart_iter = 0 ! Restart using the most recent checkpoint

In this case, Rayleigh reads the last checkpoint ﬁle (found within the Checkpoints directory) to determine which

checkpoint it reads.

5.2 Quicksaves

In practice, Rayleigh checkpoints are used for two purposes: (1) guarding against unexpected crashes and (2)

supplementing a run’s record with a series of restart points. While standard checkpoints may serve both purposes,

the frequency with which checkpoints are written for purpose (1) is often much higher than that needed for purpose

(2). This means that a lot of data culling is performed at the end of a run or, if disk space is a particularly limiting

factor, during the run itself. For this reason, Rayleigh has a quicksave checkpointing scheme in addition to the

standard scheme. Quicksaves can be written with high-cadence, but require low storage due to the rotating reuse of

quicksave ﬁles.

5.2.1 1.2.1 Generating Quicksaves

The cadence with which quicksaves are written can be speciﬁed by setting the quicksave interval variable in the

temporal controls namelist. Alternatively, the elapsed wall time (in minutes) that passes between quicksaves

may be controlled by specifying the quicksave minutes variable. If both quicksave interval and quicksave minutes

are speciﬁed, quicksave minutes takes precedence.

What distinguishes quicksaves from standard checkpoints is that only a speciﬁed number of quicksaves exist on the

disk at any given time. That number is determined by the value of num quicksaves. Quicksave ﬁles begin with the

preﬁx quicksave XX, where XX is a 2-digit code, ranging from 1 through num quicksaves and indicates the quicksave

number. Consider the following example:

&temporal_controls_namelist

checkpoint_interval = 35000 ! Generate a standard checkpoint once every 35,000 time steps

quicksave_interval = 10000 ! Generate a quicksave once every 10,000 time steps

num_quicksaves = 2 ! Keep only two quicksaves on disk at a time

At time step 10,000, a set of checkpoint ﬁles beginning with preﬁx quicksave 01 will be generated. At time step

20,000, a set of checkpoint ﬁles beginning with preﬁx quicksave 02 will be generated. Following that, at time step

30,000, another checkpoint will be generated, but it will overwrite the existing quicksave 01 ﬁles. At time step 40,000,

the quicksave 02 ﬁles will be overwritten, and so forth. Because the num quicsaves was set to 2, ﬁlenames with

preﬁx quicksave 03 will never be generated.

Note that checkpoints beginning with an 8-digit preﬁx (e.g., 00035000) are still written to disk regularly and are not

aﬀected by the quicksave checkpointing. On time steps where a quicksave and a standard checkpoint would both

be written, only the standard checkpoint is written. Thus, at time step 70,000 in the example above, a standard

checkpoint would be written, and the ﬁles beginning with quicksave 01 would remain unaltered.

5.2.2 1.2.2 Restarting from Quicksaves

Restarting from quicksave XX may be accomplished by specifying the value of restart iter to be -XX (i.e., the negative

of the quicksave you wish to restart from). The following example shows how to restart the hydrodynamic variables

from quicksave 02, while also initializing a random magnetic ﬁeld. “‘

&initial_conditions_namelist

init_type = -1 ! Restart hydro variables from a checkpoint

magnetic_init_type = 7 ! Initialize a random magnetic field

restart_iter = -2 ! Restart from quicksave number 2

Note that the ﬁle last checkpoint contains the number of last checkpoint written. This might be a quicksave or a

standard checkpoint. Specifying a value of zero for restart iter thus works with quicksaves and standard checkpoints

alike.

5.3 1.3 Checkpoint Logs

When checkpoints are written, the number of the most recent checkpoint is appended to a ﬁle named checkpoint log,

found in the Checkpoints directory. The checkpoint log can be used to identify the time step number of a quicksave

ﬁle that otherwise has no identifying information. While this information is also contained in the grid etc ﬁle, those

are written in unformatted binary and cumbersome to access from the terminal command line.

An entry in the log of ”00050000 02” means that a checkpoint was written at time step 50,000 to quicksave 02. An

entry lacking a two-digit number indicates that a standard checkpoint was written at that time step. The most

recent entry in the checkpoint log always comes at the end of the ﬁle.

6 Diagnostic Outputs

Rayleigh comes bundled with an in-situ diagnostics package that allows the user to sample a simulation in a variety

of ways, and at user-speciﬁed intervals throughout a run. This package is comprised of roughly 17,000 lines of code

(about half of the Rayleigh code base), and it is complex enough that we describe it in two other documents. We

refer the user to :

1. The diagnostics plotting manual, provided in three formats:

•Rayleigh/etc/analysis/Diagnostics Plotting.ipynb (Jupyter Python notebook format; recommended for

interactive use)

•Rayleigh/doc/Diagnostics Plotting.html (recommended for optimal viewing; generated from the .ipynb

ﬁle)

•Rayleigh/doc/Diagnostics Plotting.pdf (same content as .html and .ipynb, but formatting quality is infe-

rior)

2. Rayleigh/doc/rayleigh output variables.pdf – This companion document provides the output menu system

referred to in the main diagnostics documentation.

A number of stand-alone Python plotting examples may also be found in the Rayleigh/etc/analysis directory.

7 I/O Redirection

Rayleigh writes all text output (e.g., error messages, iteration counter, etc.) to stdout by default. Diﬀerent computing

centers handle stdout in diﬀerent ways, but typically one of two path is taken. On some machines, a log ﬁle is created

immediately and updated continuously as the simulation runs. On other machines, stdout is buﬀered on-node and

written to disk only when the run has terminated.

There are situations where it can be advantageous to have a regularly updated log ﬁle whose update frequency may

be controlled. This feature exists in Rayleigh and may be accessed by assigning values to stdout ﬂush interval

and stdout ﬁle in the io controls namelist.

&io_controls_namelist

stdout_flush_interval = 1000

stdout_file = ‘routput’

Set stdout ﬁle to the name of a ﬁle that will contain Rayleigh’s text output. In the example above, a ﬁle named

routput will be appear in the simulation directory and will be updated periodically throughout the run. The variable

stdout ﬂush interval determines how many lines of text are buﬀered before they are ﬂushed to routput. Rayleigh

prints time-step information during each time step, and so setting this variable to a relatively large number (e.g.,

100+) prevents excessive disk access from occurring throughout the run. In the example above, a text buﬀer ﬂush

will occur once 1000 lines of text have been accumulated.

Changes in the time-step size and self-termination of the run will also force a text-buﬀer ﬂush. Unexpected crashes

and sudden termination by the system job scheduler do not force a buﬀer ﬂush. Note that the default value of

stdout ﬁle is ‘noﬁle’. If this value is speciﬁed, output will directed to normal stdout.

8 Ensemble Mode

Rayleigh can also be used to run multiple simulations under the umbrella of a single executable. This functionality is

particularly useful for running parameter space studies, which often consist of mulitple, similarly-sized simulations,

in one shot. Moreover, as some queuing systems favor large jobs over small jobs, an ensemble mode is useful for

advancing multiple small simulations through the queue in a reasonable timeframe.

Running Rayleigh in ensemble mode is relatively straightforward. To begin with, create a directory for each simulation

as you normally would, and place an appropriately modiﬁed main input into each directory. These directories should

all reside within the same parent directory. Within that parent directory, you should place a copy of the Rayleigh

executable (or a softlink). In addition, you should create a text ﬁle named run list that contains the name of each

simulation directory, one name per line. An ensemble job may then be executed by calling Rayleigh with nruns

command line ﬂag as:

user@machinename ~/runs/ $ mpiexec -np Y ./rayleigh.opt -nruns X

Here, Y is the total number of cores needed by all X simulations listed in run list.

Example: Suppose you wish to run three simulations at once from within a parent directory named ensemble and

that the simulation directories are named run1, run2, and run3. When performing an ls from within ensemble, you

should see 5 items.

user@machinename ~/runs/ $ cd ensemble

user@machinename ~/runs/ensemble $ ls

rayleigh.opt run1 run2 run3 run_list

In this example, the contents of run list should be the local names of your ensemble run-directories, namely run1,

run2, and run3.

user@machinename ~runs/ensemble $ more run_list

run1

run2

run3

<-- place an empty line here

Note that some Fortran implementations will not read the last line in run list unless it ends in a newline character.

Avoid unexpected crashes by hitting ”enter” following your ﬁnal entry in run list.

Before running Rayleigh, make sure you know how many cores each simulation needs by examining the main input

ﬁles:

user@machinename ~runs/ensemble $ head run1/main_input

&problemsize_namelist

n_r = 128

n_theta = 192

nprow = 16

npcol = 16

user@machinename ~runs/ensemble $ head run2/main_input

&problemsize_namelist

n_r = 128

n_theta = 384

nprow = 32

npcol = 16

user@machinename ~runs/ensemble $ head run3/main_input

&problemsize_namelist

n_r = 64

n_theta = 192

nprow = 16

npcol = 16

In this example, we need a total of 1024 cores (256+512+256) to execute three simulations, and so the relevant call

to Rayleigh would be:

user@machinename ~/runs/ $ mpiexec -np 1024 ./rayleigh.opt -nruns 3

Closing Notes: When running in ensemble mode, it is strongly recommended that you redirect standard output for

each simulation to a text ﬁle (see §7). Otherwise, all simulations write to the same default (machine-dependent) log

ﬁle, making it diﬃcult to read. Moreover, some machines such as NASA Pleiades will terminate a run if the log ﬁle

becomes too long. This is easy to do when multiple simulations are writing to the same ﬁle.

Finally, The ﬂags -nprow and -npcol are ignored when -nruns is speciﬁed. The row and column conﬁguration for

all simulations needs to be speciﬁed in their respective main input ﬁles instead.

9 Acknowledging Rayleigh

We politely ask that you acknowledge the author, the National Science Foundation, and CIG in any work enabled by

Rayleigh. This includes work deriving from results generated by Rayleigh. It also includes work that derives from

the use of software that incorporates original or modiﬁed Rayleigh source code.

Sample acknowledgement text can always be found in Rayleigh/ACKNOWLEDGE.

You are reading the documentation for Rayleigh version 0.9.0. A formal paper describing the numerics and paral-

lelization methods is forthcoming. In the interim, please cite the following two references in any academic publications

that present results enabled by Rayleigh:

1. Featherstone, N.A. & Hindman, B.W., 2016, Astrophys. J.,818, 32, DOI: 10.3847/0004-637X/818/1/32

2. Matsui et al., 2016, Geochemistry, Geophysics, Geosystems,17,1586, DOI: 10.1002/2015GC006159

Rayleigh User Guide

Rayleigh_User_Guide

Navigation menu

Versions of this User Manual:

Views

Navigation