# Contents

## Contents

Nick Featherstone (January, 2018)

**NOTE:** This document can be viewed in PDF or HTML (recommended) form. It can also be run as an interactive Jupyter notebook.

# Contents#

Running a Benchmark with Sample Output

Configuring your Python environment

Overview of Rayleigh’s Diagnostic Package

Global Averages

Shell Averages

Azimuthal Averages

Simulation Slices

Spherical Harmonic Spectra

Point Probes

Modal Outputs

# I. Running a Benchmark with Sample Output#

Before you can plot data, you will need to generate data. The code samples in this notebook assume that you have run the model described by the input file found in:

rayleigh/input_examples/benchmark_diagnostics_input

This input file instructs *Rayleigh* to run the Christensen et al. (2001) hydrodynamic (case 0) benchmark. Running this model with the prescribed outputs will generate approximately 70 MB of data.

To run this model: 1. Create a directory for your simulation (e.g., **mkdir my_test_run**) 2. Copy the input file: **cp rayleigh/input_examples/benchmark_diagnostics_input my_test_run/main_input** 3. Copy or soft-link the *rayleigh* executable: **cp Rayleigh/bin/rayleigh.opt my_test_run/.** 4. Run the code: **mpiexec -np N ./rayleigh -nprow n -npcol m** (choose values of {N,n,m} such that n x m = N)

The code will run for 40,000 timesteps, or four viscous diffusion times. While it runs, *Rayleigh* will perform an in-situ analysis of the accuracy benchmark. Reports are written once every 1,000 time steps and are stored in the *Benchmark_Reports* subdirectory. Examine file 00030000 and ensure that you see similar results to those below. Your exact numbers may differ slightly, but all quantities should be under 1% difference.

Observable |
Measured |
Suggested |
% Difference |
Std. Dev. |
---|---|---|---|---|

Kinetic Energy |
58.347827 |
58.348000 |
-0.000297 |
0.000000 |

Temperature |
0.427424 |
0.428120 |
-0.162460 |
0.000101 |

Vphi |
-10.119483 |
-10.157100 |
-0.370356 |
0.013835 |

Drift Frequency |
0.183016 |
0.182400 |
0.337630 |
0.007295 |

If necessary, copy the data to the system on which you intend to conduct your analysis. Before you can plot, you will need to configure your Python environment.

# II. Configuring Your Python Environment#

Rayleigh comes packaged with a Python library (*rayleigh_diagnostics.py*) that provides data structures and methods associated with each type of diagnostic output in Rayleigh. This library relies on Numpy and is compatible with Python 3.x or 2.x (The *print* function is imported from the **future** module).

If you wish to follow along with the plotting examples described in this document, you will need to have the Numpy and Matplotlib Python packages installed. The following versions of these packages were used when creating these examples: * Matplotlib v2.0.2 * Numpy v1.13.1

Unless you are experienced at installing and managing Python packages, I recommend setting up a virtual environment for Python using Conda. You may also install the required packages manually, but the advantage of this approach is that you maintain an entirely separate version of Python and related packages for this project. Below are directions for setting up a Python/Conda environment with Intel-optimized Python packages on a Linux system (Mac and Windows work similarly).

## Conda Installation on Linux Systems#

**Step 1:** Download the appropriate Miniconda installation script from https://conda.io/miniconda.html (choose Python 3.x)

**Step 2:** Make the shell script executable via: **chmod +x Miniconda3-latest-Linux-x86_64.sh** (or similar script name)

**Step 3:** Run the installation script: **./Miniconda3-latest-Linux-x86_64.sh**

*NOTE:* The default installation directory is your home directory. This is also where Python packages for your Conda environments will be installed. Avoid installing to a disk with limited space (user home directories on HPC systems are often limited to a few GB).

*NOTE:* Unless you have a specific reason not to do so, answer “yes” to the question concerning prepending to PATH.

**Step 5:** Update your Conda: **conda update conda**

**Step 6:** Add the Intel Conda channel: **conda config –add channels intel**

**Step 7:** Create a virtual environment for Intel’s Conda distribution: **conda create -n idp intelpython3_full python=3**

*NOTE:* In this case, *idp* will be your virtual environment name. You are free to pick an alternative when running conda create.

*NOTE:* A number of Python packages will be downloaded, including Numpy and Matplotlib. The process may appear to hang at the last step. Be patient.

**Step 8:** Activate your virtual environment: **source activate idp**

**Step 9:** Verify your installation. Type **python** and then type the following commands at the prompt: 1. import numpy 2. import matplotlib

If those commands worked without error, you may close Python ( type **exit()** ). You can revert to your native environment by typing **source deactivate** (or just close the terminal). Whenever you wish to access your newly-installed Python, type **source activate idp** first, before running python.

## Preparing to Plot#

All examples in this document rely on the rayleigh_diagnostics module. This module is located in Rayleigh/post_processing, along with several standalone scripts copies from the individual sections of this document. For example, the script **plot_G_Avgs.py** contains the code from section IV below. All python files you wish to use will need to reside in either your run directory (recommended) or a directory within your PYTHONPATH.

We suggest copying all python files to your my_test_run directory: 1. cp Rayleigh/post_processing/*.py my_test_run/. 2. cp Rayleigh/post_processing/*.ipynb my_test_run/.

## The Jupyter Notebook#

This document resides in three places: 1. Rayleigh/doc/Diagnostic_Plotting.pdf 2. Rayleigh/doc/Diagnostic_Plotting.html 3. Rayleigh/post_processing/Diagnostic_Plotting.ipynb

The third file is a Jupyter notebook file. This source code was used to generate the html and pdf documents. The notebook is designed to be run from within a Rayleigh simulation directory. If you wish to follow along interactively, copy the Jupyter notebook file from Rayleigh/post_processing/ into your Rayleigh simulation directory (step 2 from *Preparing to Plot*). You can run the file in Jupyter via: 1. source activate idp 2. jupyter notebook (from within your
my_test_run directory) 3. select Diagnostic_Plotting.ipynb in the file menu that presents itself.

When finished: 1. To close the notebook, type **ctrl+c** and enter “yes” when prompted to shut down the notebook server. 2. type **source deactivate**

# III. Overview of Diagnostics in Rayleigh#

*Rayleigh’s* diagnostics package facilitates the in-situ analysis of a simulation using a variety of sampling methods. Each sampling method may be applied to a unique set of sampled quantities. Sampling methods are hereafter referred to as *output types* and sampled quantities as *output variables*.

Files of each output type are stored in a similarly-named subdirectory within the *Rayleigh* simulation directory. Output files are numbered by the time step of the final data record stored in the file. Output behavior for each simulation is controlled through the *main_input* file. For each output type, the user specifies the output variables, cadence, records-per-file, and other properties by modifying the appropriate variables in the **output_namelist** section of *main_input*.

## Basic Output Control#

Each output type in *Rayleigh* has at least three namelist variables that govern its behavior:

**{OutputType}_values**: comma-separated list of menu codes corresponding to the desired output variables

**{OutputType}_frequency**: integer value that determines how often this type of output is performed

**{OutputType}_nrec**: integer value that determines how many records are stored in each output file.

All possible output variables and their associated menu codes are described in **rayleigh/doc/rayleigh_output_variables.pdf** You may find it useful to have that document open while following along with examples in this notebook.

As an example of how these variables work, suppose that we want to occasionally output equatorial cuts (output type) of temperature, kinetic energy density, and radial velocity (output variables). At the same time, we might wish to dump full-volume averages (output type) of kinetic and magnetic energy (output variables) with a higher cadence. In that case, something similar to the following would appear in main_input:

This tells *Rayleigh* to output full-volume-averages of kinetic energy density (value code 401) and magnetic energy density (value code 1101) once every 50 time steps, with 100 records per file. Files are named based on the time step number of their final record. As a result, information from time steps 50, 100, 150, …, 4950, 5000 will be stored in the file named *G_Avgs/00005000*. Time steps 5050 through 10,000 will stored in *G_Avgs/00010000*, and so on.

For the equatorial cuts, *Rayleigh* will output radial velocity (code 1), the kinetic energy density (code 401) and temperature (code 501) in the equatorial plane once every 2,500 time steps, storing two time steps per file. Data from time steps 2,500 and 5,000 will be stored in *Equatorial_Slices/00005000*. Data from time steps 7,500 and 10,000 will be stored in *Equatorial_Slices/00010000* , and so on.

*This general organizational scheme for output was adapted from that developed by Thomas Clune for the ASH code.*

## Positional Output Control#

Many of *Rayleigh’s* output types allow the user to specify a set of gridpoints at which to sample the simulation. A user can, for example, output spherical surfaces sampled at arbitrary radii, or a meridional plane sampled at a specific longitude. This behavior is controlled through additional namelist variables; we refer to these variables as positional specifiers. In the sections that follow, positional specifiers associated with a given output type, if any, will be defined.

Positional specifiers are either *indicial* or *normalized*. In the *main_input* file, indicial specifiers can be assigned a comma-separated list of grid indices on which to perform the output. For example,

shellslice_levels = 1, 32, 64, 128

instructs *Rayleigh* to output shell slices at { radius[1], radius[32], radius[64], radius[128]}. Note that radius[1] is the outer boundary.

While useful in some situations, specifying indices can lead to confusion if a simulations resolution needs to be changed at some point during a model’s evolution. For example if the radial grid initially had 128 points, index 128 would correspond to the lower boundary. If the resolution were to double, index 128 would correspond to mid-shell.

For this reason, all positional specifiers may also be written in normalized form. Instead of integers, the normalized specifier is assigned a comma separated list of real values in the range [0,1]. The value of zero corresponds to the lowest-value grid coordinate (e.g., the inner radial boundary or theta=0 pole). The value 1 corresponds to the maximal coordinate (e.g., the outer radial boundary or theta=pi pole). A value of 0.5 corresponds to mid-domain. Normalized coordinates are indicated by adding *_nrm* to the indicial specifier’s name. For example,

shellslice_levels_nrm= 0, 0.5, 0.95

instructs *Rayleigh* to output shell slices at the lower boundary, mid-shell, and slightly below the upper boundary. *Rayleigh* does not interpolate, but instead picks the grid coordinate closest to each specified normalized coordinate.

We recommend using normalized coordinates to avoid inconsistencies between restarts. They also overcome difficulties associated with the non-uniform nature of the radial and theta grids wherein grid points cluster near the boundaries.

**Positional Ranges** Ranges of coordinates can be specified using shorthand, if desired. The inclusive coordinate range [X,Y] is indicated by a positive/negative number pair appearing in the indicial or normalized coordinate list. Multiple ranges can be specified within a list. For example,

shellslice_levels = 1,10,-15, 16, 20,-25, 128

would instruct *Rayleigh* to output shell slices at radial indices = { 1, 10, 11, 12, 13, 14, 15, 16, 20, 21, 22, 23, 24, 25, 128}

Similarly,

shellslice_levels_nrm = 0,-0.5, 1.0

instructs *Rayleigh* to output shells at all radii in the lower half of the domain, and at the outer boundary.

# IV. Global Averages#

**Summary:** Full-volume averages of requested output variables over the full, spherical shell

**Subdirectory:** G_Avgs

**main_input prefix:** globalavg

**Python Class:** G_Avgs

**Additional Namelist Variables:**

*Before proceeding, ensure that you have copied Rayleigh/post_processing/rayleigh_diagnostics.py to your simulation directory. This Python module is required for reading Rayleigh output into Python.*

Examining the *main_input* file, we see that the following output values have been denoted for the Global Averages (see *rayleigh_output_variables.pdf* for the mathematical formulae):

Menu Code |
Description |
---|---|

401 |
Full Kinetic Energy Density (KE) |

402 |
KE (radial motion) |

403 |
KE (theta motion) |

404 |
KE (phi motion) |

405 |
Mean Kinetic Energy Density (MKE) |

406 |
MKE (radial motion) |

407 |
MKE (theta motion) |

408 |
MKE (phi motion) |

409 |
Fluctuating Kinetic Energy Density (FKE) |

410 |
FKE (radial motion) |

411 |
FKE (theta motion) |

412 |
FKE (phi motion) |

In the example that follows, we will plot the time-evolution of these different contributions to the kinetic energy budget. We begin with the following preamble:

```
[ ]:
```

```
%matplotlib inline
from rayleigh_diagnostics import G_Avgs, build_file_list
import matplotlib.pyplot as plt
import numpy
```

The preamble for each plotting example will look similar to that above. We import the numpy and matplotlib.pyplot modules, aliasing the latter to *plt*. We also import two items from *rayleigh_diagnostics*: a helper function *build_file_list* and the *GlobalAverage* class.

The *G_Avgs* class is the Python class that corresponds to the full-volume averages stored in the *G_Avgs* subdirectory of each Rayleigh run.

We will use the build_file_list function in many of the examples that follow. It’s useful when processing a time series of data, as opposed to a single snapshot. This function accepts three parameters: a beginning time step, an ending time step, and a subdirectory (path). It returns a list of all files found in that directory that lie within the inclusive range [beginning time step, ending time step]. The file names are prepended with the subdirectory name, as shown below.

```
[ ]:
```

```
# Build a list of all files ranging from iteration 0 million to 1 million
files = build_file_list(0,1000000,path='G_Avgs')
print(files)
```

We can create an instance of the G_Avgs class by initializing it with a filename. The optional keyword parameter *path* is used to specify the directory. If *path* is not specified, its value will default to the subdirectory name associated with the datastructure (*G_Avgs* in this instance).

Each class was programmed with a **docstring** describing the class attributes. Once you created an instance of a rayleigh_diagnostics class, you can view its attributes using the help function as shown below.

```
[ ]:
```

```
a = G_Avgs(filename=files[0],path='') # Here, files[0]='G_Avgs/00010000'
#a= G_Avgs(filename='00010000') would yield an equivalent result
help(a)
```

Examining the docstring, we see a few important attributes that are common to the other outputs discussed in this document: 1. niter – the number of time steps in the file 2. nq – the number of output variables stored in the file 3. qv – the menu codes for those variables 4. vals – the actual data 5. time – the simulation time corresponding to each output dump

The first step in plotting a time series is to collate the data.

```
[ ]:
```

```
# Loop over all files and concatenate their data into a single array
nfiles = len(files)
for i,f in enumerate(files):
a = G_Avgs(filename=f,path='')
if (i == 0):
nq = a.nq
niter = a.niter
gavgs = numpy.zeros((niter*nfiles,nq),dtype='float64')
iters = numpy.zeros(niter*nfiles,dtype='int32')
time = numpy.zeros(niter*nfiles,dtype='float64')
i0 = i*niter
i1 = (i+1)*niter
gavgs[i0:i1,:] = a.vals
time[i0:i1] = a.time
iters[i0:i1] = a.iters
```

## The Lookup Table (LUT)#

The next step in the process is to identify where within the *gavgs* array our deisired output variables reside. Every Rayleigh file object possesses a lookup table (lut). The lookup table is a python list used to identify the index within the vals array where a particular menu code resides. For instance, the menu code for the theta component of the velocity is 2. The location of v_theta in the vals array is then stored in lut[2].

Note that you should never assume that output variables are stored in any particular order. Moreover, the lookup table is unique to each file and is likely to change during a run if you modify the output variables in between restarts. When running the benchmark, we kept a consistent set of outputs throughout the entirety of the run. This means that the lookup table did not change between outputs and that we can safely use the final file’s lookup table (or any other file’s table) to reference our data.

## Plotting Kinetic Energy#

Let’s examine the different contributions to the kinetic energy density in our models. Before we can plot, we should use the lookup table to identify the location of each quantity we are interested in plotting.

```
[ ]:
```

```
#The indices associated with our various outputs are stored in a lookup table
#as part of the GlobalAverage data structure. We define several variables to
#hold those indices here:
lut = a.lut
ke = lut[401] # Kinetic Energy (KE)
rke = lut[402] # KE associated with radial motion
tke = lut[403] # KE associated with theta motion
pke = lut[404] # KE associated with azimuthal motion
#We also grab some energies associated with the mean (m=0) motions
mke = lut[405]
mrke = lut[406] # KE associated with mean radial motion
mtke = lut[407] # KE associated with mean theta motion
mpke = lut[408] # KE associated with mean azimuthal motion
#We also output energies associated with the fluctuating/nonaxisymmetric
#motions (e.g., v- v_{m=0})
fke = lut[409]
frke = lut[410] # KE associated with mean radial motion
ftke = lut[411] #KE associated with mean theta motion
fpke = lut[412] # KE associated with mean azimuthal motion
```

To begin with, let’s plot the total, mean, and fluctuating kinetic energy density during the initial transient phase, and then during the equilibrated phase.

```
[ ]:
```

```
sizetuple=(10,3)
fig, ax = plt.subplots(ncols=2, figsize=sizetuple)
ax[0].plot(time, gavgs[:,ke], label='KE')
ax[0].plot(time, gavgs[:,mke],label='MKE')
ax[0].plot(time, gavgs[:,fke], label='FKE')
ax[0].legend(loc='center right', shadow=True)
ax[0].set_xlim([0,0.2])
ax[0].set_title('Equilibration Phase')
ax[0].set_xlabel('Time')
ax[0].set_ylabel('Energy')
ax[1].plot(time, gavgs[:,ke], label='KE')
ax[1].plot(time, gavgs[:,mke], label = 'MKE')
ax[1].plot(time,gavgs[:,fke],label='FKE')
ax[1].legend(loc='center right', shadow=True)
ax[1].set_title('Entire Time-Trace')
ax[1].set_xlabel('Time')
ax[1].set_ylabel('Energy')
saveplot = False # Plots appear in the notebook and are not written to disk (set to True to save to disk)
savefile = 'energy_trace.pdf' #Change .pdf to .png if pdf conversion gives issues
plt.tight_layout()
plt.show()
```

We can also look at the energy associated with each velocity component. Note that we log scale in the last plot. There is very little mean radial or theta kinetic energy; it is mostly phi energy.

```
[ ]:
```

```
sizetuple=(5,10)
xlims=[0,0.2]
fig, ax = plt.subplots(ncols=1, nrows=3, figsize=sizetuple)
ax[0].plot(time, gavgs[:,ke], label='KE')
ax[0].plot(time, gavgs[:,rke],label='RKE')
ax[0].plot(time, gavgs[:,tke], label='TKE')
ax[0].plot(time, gavgs[:,pke], label='PKE')
ax[0].legend(loc='center right', shadow=True)
ax[0].set_xlim(xlims)
ax[0].set_title('Total KE Breakdown')
ax[0].set_xlabel('Time')
ax[0].set_ylabel('Energy')
ax[1].plot(time, gavgs[:,fke], label='FKE')
ax[1].plot(time, gavgs[:,frke], label='FRKE')
ax[1].plot(time, gavgs[:,ftke], label='FTKE')
ax[1].plot(time, gavgs[:,fpke], label='FPKE')
ax[1].legend(loc='center right', shadow=True)
ax[1].set_xlim(xlims)
ax[1].set_title('Fluctuating KE Breakdown')
ax[1].set_xlabel('Time')
ax[1].set_ylabel('Energy')
ax[2].plot(time, gavgs[:,mke], label='MKE')
ax[2].plot(time, gavgs[:,mrke], label='MRKE')
ax[2].plot(time, gavgs[:,mtke], label='MTKE')
ax[2].plot(time, gavgs[:,mpke], label='MPKE')
ax[2].legend(loc='lower right', shadow=True)
ax[2].set_xlim(xlims)
ax[2].set_title('Mean KE Breakdown')
ax[2].set_xlabel('Time')
ax[2].set_ylabel('Energy')
ax[2].set_yscale('log')
plt.tight_layout()
plt.show()
```

# V. Shell Averages#

**Summary:** Spherical averages of requested output variables. Each output variable is stored as a 1-D function of radius.

**Subdirectory:** Shell_Avgs

**main_input prefix:** shellavg

**Python Class:** Shell_Avgs

**Additional Namelist Variables:**

The Shell-Averaged outputs are useful for examining how quantities vary as a function of radius. They are particularly useful for examining the distribution of energy as a function of radius, or the heat flux balance established by the system.

Examining the *main_input* file, we see that the following output values have been denoted for the Shell Averages (see *rayleigh_output_variables.pdf* for mathematical formulae):

Menu Code |
Description |
---|---|

1 |
Radial Velocity |

2 |
Theta Velocity |

3 |
Phi Velocity |

501 |
Temperature Perturbation |

1438 |
Radial Convective Heat Flux |

1468 |
Radial Conductive Heat Flux |

In the example that follows, we will plot the spherically-averaged velocity field as a function of radius, the mean temperature profile, and the radial heat flux. We begin with a preamble similar to that used for the Global Averages. Using the help function, we see that the Shell_Avgs data structure is similar to that of the G_Avgs. There are three important differences: * There is a radius attribute (necessary if we want to plot anything vs. radius) * The dimensionality of the values array has changed; radial index forms the first dimension. * The second dimension of the values array has a length of 4. In addition to the spherical mean, the 1st, 2nd and 3rd moments are stored in indices 0,1,2, and 3 respectively.

```
[ ]:
```

```
from rayleigh_diagnostics import Shell_Avgs, build_file_list
import matplotlib.pyplot as plt
import numpy
# Build a list of all files ranging from iteration 0 million to 1 million
files = build_file_list(0,1000000,path='Shell_Avgs')
a = Shell_Avgs(filename=files[0], path='')
help(a)
```

While it can be useful to look at instaneous snapshots of Shell Averages, it’s often useful to examine these outputs in a time-averaged sense. Let’s average of all 200 snapshots in the last file that was output. We could average over data from multiple files, but since the benchmark run achieves a nearly steady state, a single file will do in this case.

```
[ ]:
```

```
nfiles = len(files)
nr = a.nr
nq = a.nq
nmom = 4
niter = a.niter
radius = a.radius
savg=numpy.zeros((nr,nmom,nq),dtype='float64')
for i in range(niter):
savg[:,:,:] += a.vals[:,:,:,i]
savg = savg*(1.0/niter)
lut = a.lut
vr = lut[1] # Radial Velocity
vtheta = lut[2] # Theta Velocity
vphi = lut[3] # Phi Velocity
thermal = lut[501] # Temperature
eflux = lut[1440] # Convective Heat Flux (radial)
cflux = lut[1470] # Conductive Heat Flux (radial)
```

## Velocity vs. Radius#

Next, we plot the mean velocity field, and its first moment, as a function of radius. Notice that the radial and theta velocity components have a zero spherical mean. Since we are running an incompressible model, this is a good sign!

```
[ ]:
```

```
sizetuple = (7,7)
fig, ax = plt.subplots(nrows=2, ncols =1, figsize=sizetuple)
ax[0].plot(radius,savg[:,0,vr],label=r'$v_r$')
ax[0].plot(radius,savg[:,0,vtheta], label=r'$v_\theta$')
ax[0].plot(radius,savg[:,0,vphi], label=r'$v_\phi$')
ax[0].legend(shadow=True,loc='lower right')
ax[0].set_xlabel('Radius')
ax[0].set_ylabel('Velocity')
ax[0].set_title('Spherically-Averaged Velocity Components')
ax[1].plot(radius,savg[:,1,vr],label=r'$v_r$')
ax[1].plot(radius,savg[:,1,vtheta], label=r'$v_\theta$')
ax[1].plot(radius,savg[:,1,vphi], label=r'$v_\phi$')
ax[1].legend(shadow=True,loc='upper left')
ax[1].set_xlabel('Radius')
ax[1].set_ylabel('Velocity')
ax[1].set_title('Velocity Components: First Spherical Moment')
plt.tight_layout()
plt.show()
```

## Radial Temperature Profile#

We might also look at temperature …

```
[ ]:
```

```
fig, ax = plt.subplots()
ax.plot(radius,savg[:,0,thermal],label='Temperature (mean)')
ax.plot(radius,savg[:,1,thermal]*10, label='Temperature (standard dev.)')
ax.legend(shadow=True,loc='upper right')
ax.set_xlabel('Radius')
ax.set_ylabel('Temperature')
ax.set_title('Radial Temperature Profile')
plt.show()
```

## Heat Flux Contributions#

We can also examine the balance between convective and conductive heat flux. In this case, before plotting these quantities as a function of radius, we normalize them by the surface area of the sphere to form a luminosity.

```
[ ]:
```

```
fpr=4.0*numpy.pi*radius*radius
elum = savg[:,0,eflux]*fpr
clum = savg[:,0,cflux]*fpr
tlum = elum+clum
fig, ax = plt.subplots()
ax.plot(radius,elum,label='Convection')
ax.plot(radius,clum, label='Conduction')
ax.plot(radius,tlum, label='Total')
ax.set_title('Flux Balance')
ax.set_ylabel(r'Energy Flux ($\times 4\pi r^2$)')
ax.set_xlabel('Radius')
ax.legend(shadow=True)
plt.tight_layout()
plt.show()
```

# VI. Azimuthal Averages#

**Summary:** Azimuthal averages of requested output variables. Each output variable is stored as a 2-D function of radius and latitude.

**Subdirectory:** AZ_Avgs

**main_input prefix:** azavg

**Python Class:** AZ_Avgs

**Additional Namelist Variables:**

Azimuthally-Averaged outputs are particularly useful for examining a system’s mean flows (i.e., differential rotation and meridional circulation).

Examining the *main_input* file, we see that the following output values have been denoted for the Azimuthal Averages (see *rayleigh_output_variables.pdf* for mathematical formulae):

Menu Code |
Description |
---|---|

1 |
Radial Velocity |

2 |
Theta Velocity |

3 |
Phi Velocity |

201 |
Radial Mass Flux |

202 |
Theta Mass Flux |

501 |
Temperature Perturbation |

In the example that follows, we demonstrate how to plot azimuthal averages, including how to generate streamlines of mass flux. Note that since the benchmark is Boussinesq, our velocity and mass flux fields are identical. This is not the case when running an anelastic simulation.

We begin with the usual preamble and also import two helper routines used for displaying azimuthal averages.

Examining the data structure, we see that the vals array is dimensioned to account for latitudinal variation, and that we have new attributes costheta and sintheta used for referencing locations in the theta direction.

```
[ ]:
```

```
from rayleigh_diagnostics import AZ_Avgs, build_file_list, plot_azav, streamfunction
import matplotlib.pyplot as plt
import pylab
import numpy
#from azavg_util import *
files = build_file_list(30000,40000,path='AZ_Avgs')
az = AZ_Avgs(files[0],path='')
help(az)
```

Before creating our plots, let’s time-average over the last two files that were output (thus sampling the equilibrated phase).

```
[ ]:
```

nfiles = len(files)
tcount=0
for i in range(nfiles):
az=AZ_Avgs(files[i],path='')
if (i == 0):
nr = az.nr
ntheta = az.ntheta
nq = az.nq
azavg=numpy.zeros((ntheta,nr,nq),dtype='float64')
for j in range(az.niter):
azavg[:,:,:] += az.vals[:,:,:,j]
tcount+=1
azavg = azavg*(1.0/tcount) # Time steps were uniform for this run, so a simple average will suffice
lut = az.lut
vr = azavg[:,:,lut[1]]
vtheta = azavg[:,:,lut[2]]
vphi = azavg[:,:,lut[3]]
rhovr = azavg[:,:,lut[201]]
rhovtheta = azavg[:,:,lut[202]]
temperature = azavg[:,:,lut[501]]
radius = az.radius
costheta = az.costheta
sintheta = az.sintheta

Before we render, we need to do some quick post-processing: 1. Remove the spherical mean temperature from the azimuthal average. 2. Convert v_phi into omega 3. Compute the magnitude of the mass flux vector 4. Compute stream function associated with the mass flux field

```
[ ]:
```

```
#Subtrace the ell=0 component from temperature at each radius
for i in range(nr):
temperature[:,i]=temperature[:,i] - numpy.mean(temperature[:,i])
#Convert v_phi to an Angular velocity
omega=numpy.zeros((ntheta,nr))
for i in range(nr):
omega[:,i]=vphi[:,i]/(radius[i]*sintheta[:])
#Generate a streamfunction from rhov_r and rhov_theta
psi = streamfunction(rhovr,rhovtheta,radius,costheta,order=0)
#contours of mass flux are overplotted on the streamfunction PSI
rhovm = numpy.sqrt(rhovr**2+rhovtheta**2)*numpy.sign(psi)
```

**NOTE:**If you want to save any of these figures, you can mimic the saveplot logic at the bottom of this example.

```
[ ]:
```

```
# We do a single row of 3 images
# Spacing is default spacing set up by subplot
figdpi=300
sizetuple=(5.5*3,3*3)
tsize = 20 # title font size
cbfsize = 10 # colorbar font size
fig, ax = plt.subplots(ncols=3,figsize=sizetuple,dpi=figdpi)
plt.rcParams.update({'font.size': 14})
#temperature
#ax1 = f1.add_subplot(1,3,1)
units = '(nondimensional)'
plot_azav(fig,ax[0],temperature,radius,costheta,sintheta,mycmap='RdYlBu_r',boundsfactor = 2,
boundstype='rms', units=units, fontsize = cbfsize)
ax[0].set_title('Temperature',fontsize=tsize)
#Differential Rotation
#ax1 = f1.add_subplot(1,3,2)
units = '(nondimensional)'
plot_azav(fig,ax[1],omega,radius,costheta,sintheta,mycmap='RdYlBu_r',boundsfactor = 1.5,
boundstype='rms', units=units, fontsize = cbfsize)
ax[1].set_title(r'$\omega$',fontsize=tsize)
#Mass Flux
#ax1 = f1.add_subplot(1,3,3)
units = '(nondimensional)'
plot_azav(fig,ax[2],psi,radius,costheta,sintheta,mycmap='RdYlBu_r',boundsfactor = 1.5,
boundstype='rms', units=units, fontsize = cbfsize, underlay = rhovm)
ax[2].set_title('Mass Flux',fontsize = tsize)
saveplot=False
if (saveplot):
p.savefig(savefile)
else:
plt.show()
```

# VII. Simulation Slices#

## VII.1 Equatorial Slices#

**Summary:** 2-D profiles of selected output variables in the equatorial plane.

**Subdirectory:** Equatorial_Slices

**main_input prefix:** equatorial

**Python Class:** Equatorial_Slices

**Additional Namelist Variables:**

The equatorial-slice output type allows us to examine how the fluid properties vary in longitude and radius.

Examining the *main_input* file, we see that the following output values have been denoted for the Equatorial Slices (see *rayleigh_output_variables.pdf* for mathematical formulae):

Menu Code |
Description |
---|---|

1 |
Radial Velocity |

2 |
Theta Velocity |

3 |
Phi Velocity |

In the example that follows, we demonstrate how to create a 2-D plot of radial velocity in the equatorial plane (at a single time step).

We begin with the usual preamble. Examining the data structure, we see that the *vals* array is dimensioned to account for longitudinal variation, and that we have the new coordinate attribute *phi*.

```
[ ]:
```

```
from rayleigh_diagnostics import Equatorial_Slices
import numpy
import matplotlib.pyplot as plt
from matplotlib import ticker, font_manager
istring = '00040000'
es = Equatorial_Slices(istring)
tindex =1 # Grab second time index from this file
help(es)
```

```
[ ]:
```

```
################################
# Equatorial Slice
#Set up the grid
remove_mean = True # Remove the m=0 mean
nr = es.nr
nphi = es.nphi
r = es.radius/numpy.max(es.radius)
phi = numpy.zeros(nphi+1,dtype='float64')
phi[0:nphi] = es.phi
phi[nphi] = numpy.pi*2 # For display purposes, it is best to have a redunant data point at 0,2pi
#We need to generate a cartesian grid of x-y coordinates (both X & Y are 2-D)
radius_matrix, phi_matrix = numpy.meshgrid(r,phi)
X = radius_matrix * numpy.cos(phi_matrix)
Y = radius_matrix * numpy.sin(phi_matrix)
qindex = es.lut[1] # radial velocity
field = numpy.zeros((nphi+1,nr),dtype='float64')
field[0:nphi,:] =es.vals[:,:,qindex,tindex]
field[nphi,:] = field[0,:] #replicate phi=0 values at phi=2pi
#remove the mean if desired (usually a good idea, but not always)
if (remove_mean):
for i in range(nr):
the_mean = numpy.mean(field[:,i])
field[:,i] = field[:,i]-the_mean
#Plot
sizetuple=(8,5)
fig, ax = plt.subplots(figsize=(8,8))
tsize = 20 # title font size
cbfsize = 10 # colorbar font size
img = ax.pcolormesh(X,Y,field,cmap='jet')
ax.axis('equal') # Ensure that x & y axis ranges have a 1:1 aspect ratio
ax.axis('off') # Do not plot x & y axes
# Plot bounding circles
ax.plot(r[nr-1]*numpy.cos(phi), r[nr-1]*numpy.sin(phi), color='black') # Inner circle
ax.plot(r[0]*numpy.cos(phi), r[0]*numpy.sin(phi), color='black') # Outer circle
ax.set_title(r'$v_r$', fontsize=20)
#colorbar ...
cbar = plt.colorbar(img,orientation='horizontal', shrink=0.5, aspect = 15, ax=ax)
cbar.set_label('nondimensional')
tick_locator = ticker.MaxNLocator(nbins=5)
cbar.locator = tick_locator
cbar.update_ticks()
cbar.ax.tick_params(labelsize=cbfsize) #font size for the ticks
t = cbar.ax.xaxis.label
t.set_fontsize(cbfsize) # font size for the axis title
plt.tight_layout()
plt.show()
```

## VII.2 Meridional Slices#

**Summary:** 2-D profiles of selected output variables sampled in meridional planes.

**Subdirectory:** Meridional_Slices

**main_input prefix:** meridional

**Python Class:** Meridional_Slices

**Additional Namelist Variables:**

meridional_indices (indicial) : indices along longitudinal grid at which to output meridional planes.

meridional_indices_nrm (normalized) : normalized longitudinal grid coordinates at which to output

The meridional-slice output type allows us to examine how the fluid properties vary in latitude and radius.

Examining the *main_input* file, we see that the following output values have been denoted for the Meridional Slices (see *rayleigh_output_variables.pdf* for mathematical formulae):

Menu Code |
Description |
---|---|

1 |
Radial Velocity |

2 |
Theta Velocity |

3 |
Phi Velocity |

In the example that follows, we demonstrate how to create a 2-D plot of radial velocity in a meridional plane. The procedure is similar to that used to plot an azimuthal average.

We begin with the usual preamble and import the *plot_azav* helper function. Examining the data structure, we see that it is similar to the AZ_Avgs data structure. The *vals* array possesses an extra dimension relative to its AZ_Avgs counterpart to account for the multiple longitudes that may be output, we see attributes *phi* and *phi_indices* have been added to reference the longitudinal grid.

```
[ ]:
```

```
#####################################
# Meridional Slice
from rayleigh_diagnostics import Meridional_Slices, plot_azav
import numpy
import matplotlib.pyplot as plt
from matplotlib import ticker, font_manager
# Read the data
istring = '00040000'
ms = Meridional_Slices(istring)
tindex =1 # All example quantities were output with same cadence. Grab second time-index from all.
help(ms)
```

```
[ ]:
```

radius = ms.radius
costheta = ms.costheta
sintheta = ms.sintheta
phi_index = 0 # We only output one Meridional Slice
vr_ms = ms.vals[phi_index,:,:,ms.lut[1],tindex]
units = 'nondimensional'
# Plot
sizetuple=(8,5)
fig, ax = plt.subplots(figsize=(8,8))
tsize = 20 # title font size
cbfsize = 10 # colorbar font size
ax.axis('equal') # Ensure that x & y axis ranges have a 1:1 aspect ratio
ax.axis('off') # Do not plot x & y axes
plot_azav(fig,ax,vr_ms,radius,costheta,sintheta,mycmap='RdYlBu_r',boundsfactor = 4.5,
boundstype='rms', units=units, fontsize = cbfsize)
ax.set_title('Radial Velocity',fontsize=tsize)
plt.tight_layout()
plt.show()

## VII.3 Shell Slices#

**Summary:** 2-D, spherical profiles of selected output variables sampled in at discrete radii.

**Subdirectory:** Shell_Slices

**main_input prefix:** shellslice

**Python Class:** Shell_Slices

**Additional Namelist Variables:**

shellslice_levels (indicial) : indices along radial grid at which to output spherical surfaces.

shellslice_levels_nrm (normalized) : normalized radial grid coordinates at which to output spherical surfaces.

The shell-slice output type allows us to examine how the fluid properties vary on spherical surfaces.

Examining the *main_input* file, we see that the following output values have been denoted for the Shell Slices (see *rayleigh_output_variables.pdf* for mathematical formulae):

Menu Code |
Description |
---|---|

1 |
Radial Velocity |

2 |
Theta Velocity |

3 |
Phi Velocity |

In the example that follows, we demonstrate how to create a 2-D plot of the radial velocity on a Cartesian, lat-lon grid.

Plotting on a lat-lon grid is straightforward and illustrated below. The shell-slice data structure is also displayed via the help() function in the example below and contains information needed to define the spherical grid for plotting purposes.

```
[ ]:
```

```
#####################################
# Shell Slice
from rayleigh_diagnostics import Shell_Slices
import numpy
import matplotlib.pyplot as plt
from matplotlib import ticker, font_manager
# Read the data
istring = '00040000'
ss = Shell_Slices(istring)
help(ss)
ntheta = ss.ntheta
nphi = ss.nphi
costheta = ss.costheta
theta = numpy.arccos(costheta)
#help(ss)
tindex =1 # All example quantities were output with same cadence. Grab second time-index from all.
rindex = 0 # only output one radius
sizetuple=(8,8)
vr = ss.vals[:,:,rindex,ss.lut[1],tindex]
fig, ax = plt.subplots(figsize=sizetuple)
img = plt.imshow(numpy.transpose(vr), extent=[0,360,-90,90])
ax.set_xlabel( 'Longitude')
ax.set_ylabel( 'Latitude')
ax.set_title( 'Radial Velocity')
plt.tight_layout()
plt.show()
```

By running the cell below, we can plot different output quantities on a spherical surface. In the example shown here, we plot all three velocity components (\(u_r\), \(u_{\theta}\) and \(u_{\phi}\)) projected onto a spherical surface. For demonstration purposes, we illustrate each velocity component using different colormaps, at different latitudinal centers of vantage point etc. (for more details see comments within cell below as well as the Jupyter notebook titled “plot_shells.ipynb”).

Note that in order to successfully run the cell below, you also need to have the orthographic projection code “projection.py” within the same directory/folder as this notebook.

```
[ ]:
```

```
import numpy
import matplotlib.pyplot as plt
from matplotlib import gridspec
from rayleigh_diagnostics import Shell_Slices
from projection import plot_ortho
# This plots various data from a single shell_slice file.
# 3 diferent plots in 1 row and 3 columns are created.
# It's easy to hack this to work with multiple files
s1=Shell_Slices('00040000')
data = numpy.zeros((s1.nphi,s1.ntheta),dtype='float64')
costheta = s1.costheta
nrows=2
ncols=2
pltin = 9 # Size of each subimage in inches (will be square)
# number of rows and columns
nrow=1
ncol=3
#We use gridspec to set up a grid. We actually have nrow*2 rows, with every other
#row being a 'spacer' row that's 10% the height of the main rows.
#This was the simplest way I could come up with the have the color bars appear nicely.
fig = plt.figure(constrained_layout=False, figsize=(pltin*ncol,pltin*nrow*1.1))
spec = gridspec.GridSpec(ncols=ncol, nrows=nrow*2, figure=fig, height_ratios=[1,.1]*nrow, width_ratios=[1]*ncol)
plt.rcParams.update({'font.size': 16})
#quantities codes to plot -- here all three velocity components
qi = [1,2, 3]
nm = [r'u$_r$', r'u$_\theta$', r'u$_\phi$']
qinds = [qi, qi] # Quantity codes to plot
names = [nm, nm] # Names for labeling
lv = [[1]*ncol , [2]*ncol] # Shell levels to plot (top row is level 1, bottom row is level 2)
style1=['-','--',':']
style2=['-', '-', '-']
styles = [style1, style2] # Line style of grid lines
gwidth1=[0.5 , 1 , 1.5]
gwidth2=[1,1,1]
gwidths = [gwidth1, gwidth2] # width of grid lines for each image (Default: True)
hwidths1=[2.5,2.5,2.5] # Width of the horizon line or each image (Default: 2)
hwidths2=[2,2,2]
hwidths=[hwidths1,hwidths2]
cmaps1 = ["RdYlBu_r", "seismic", 'PiYG'] # A color table for each image (Default: RdYlBu_r)
cmaps2 = ["RdYlBu_r"]*4
cmaps = [cmaps1, cmaps1]
pgrids1 = [True, True, True]
pgrids2 = [True, True, True]
pgrids = [pgrids1, pgrids2] # Plot grids, or not for each image (Default: True)
latcens1 = [60, 45, 15]
latcens = [latcens1, latcens1] # Latitudinal center of vantage point (Default: 45 N)
loncens1 = [0,0,0]
loncens2 = [30,30,30] # Longitudinal center of vantage point (Default: 0)
loncens = [loncens1,loncens2]
##########################################################
# If the grid is plotted, the number of latitude lines
# for the grid can be controlled via the nlats keyword.
# Default: 9
# Note that if nlats is even, the equator will not be drawn
nlats1 = [3,5,7]
nlats2 = [4,6,8]
nlats = [nlats1, nlats1]
##############################################################################
# Similarly, the nlons keyword can be used to control longitude lines
# More precisely, it controls the number of MERIDIANS (great circles) drawn
# Default: 8
nlons1 = [4,8,12]
nlons2 = [4,12,16]
nlons = [nlons1,nlons1]
#Longitude grid-lines can be drawn in one of two ways:
# 1) Completely to the pole (polar_style = 'polar')
# 2) Truncated at the last line of latitue drawn (polar_style = 'truncated')
# Default: "truncated"
pstyle1 = ['truncated', 'polar', 'truncated']
pstyle = [pstyle1, pstyle1]
##############################################################
# We can also control the way in which the image is saturated
# via the scale_type keyword. There are three possibilities:
# 1) scale_type=['rms', a], where a is of type 'float'
# In this instance, the image bounds are -a*rms(data), +a*rms(data)
# 2) scale_type = ['abs', a]
# In this instance, the image bounds are -a*abs(data), +a*abs(data)
# 3) scale_type= ['force', [a,b]]
# In this instance, the image bounds are a,b
# 4) scale_type = [None,None]
# In this instance, the image bounds are min(image), max(image)
# Default: [None,None]
# Note that rms and abs are taken over projected image values, not input data
# (you only see half the data in the image)
scale_type1 = [['rms',2.0 ], [None,None], ['abs', 0.5]]
scale_type2 = [['force', [-1500,1500]], ['force',[-10000,10000]], ['rms',2.5]]
scale_types = [scale_type1, scale_type1]
# Number of pixels across each projected, interpolated image
# 768 is the default and seems to do a reasonable job
nyzi = 768
for j in range(ncol):
for i in range(nrow):
data[:,:] = s1.vals[:,:,lv[i][j],s1.lut[qinds[i][j]],0]
row_ind = 2*i # skip over space allowed for color bars
col_ind = j
print("ROW/COL: ", row_ind, col_ind)
ax = fig.add_subplot(spec[row_ind,col_ind])
cspec = spec[row_ind+1,col_ind]
caxis=None
plot_ortho(data,s1.costheta,fig,ax,caxis, hwidth=hwidths[i][j], gridstyle=styles[i][j],
gridwidth=gwidths[i][j], nyz=nyzi, colormap=cmaps[i][j],
plot_grid=pgrids[i][j], latcen=latcens[i][j], loncen= loncens[i][j],
pole_style=pstyle[i][j], nlats = nlats[i][j],scale_type=scale_types[i][j])
ptitle=names[i][j]+" (r_index = "+str(lv[i][j])+")"
ax.set_title(ptitle)
# You can save the plot as a figure
#plt.savefig('flows.pdf')
```

# VIII. Spherical Harmonic Spectra#

**Summary:** Spherical Harmonic Spectra sampled at discrete radii.

**Subdirectory:** Shell_Spectra

**main_input prefix:** shellspectra

**Python Classes:**

Shell_Spectra : Complete data structure associated with Shell_Spectra outputs.

PowerSpectrum : Reduced data structure – contains power spectrum of velocity and/or magnetic fields only.

**Additional Namelist Variables:**

shellspectra_levels (indicial) : indices along radial grid at which to output spectra.

shellspectra_levels_nrm (normalized) : normalized radial grid coordinates at which to output spectra.

The shell-spectra output type allows us to examine the spherical harmonic decomposition of output variables at discrete radii.

Examining the *main_input* file, we see that the following output values have been denoted for the Shell Spectra (see *rayleigh_output_variables.pdf* for mathematical formulae):

Menu Code |
Description |
---|---|

1 |
Radial Velocity |

2 |
Theta Velocity |

3 |
Phi Velocity |

Spherical harmonic spectra can be read into Python using either the **Shell_Spectra** or **PowerSpectrum** classes.

The **Shell_Spectra** class provides the full complex spectra, as a function of degree ell and azimuthal order m, for each specified output variable. It possesses an attribute named *lpower* that contains the associated power for each variable, along with its m=0 contributions separated and removed.

The **Power_Spectrum** class can be used to read a Shell_Spectra file and quickly generate a velocity or magnetic-field power spectrum. For this class to work correctly, your file must contain all three components of either the velocity or magnetic field. Other variables are ignored (use Shell_Spectrum’s lpower for those).

We illustrate how to use these two classes below. As usual, we call the help() function to display the docstrings that describe the different data structures embodied by each class.

```
[ ]:
```

```
import matplotlib.pyplot as plt
from matplotlib import ticker
import numpy
from rayleigh_diagnostics import Shell_Spectra, Power_Spectrum
istring = '00040000'
tind = 0
rind = 0
#help(ss)
vpower = Power_Spectrum(istring)
help(vpower)
power = vpower.power
fig, ax = plt.subplots(nrows=3, figsize=(6,6))
ax[0].plot(power[:,rind,tind,0])
ax[0].set_xlabel(r'Degree $\ell$')
ax[0].set_title('Velocity Power (total)')
ax[1].plot(power[:,rind,tind,1])
ax[1].set_xlabel(r'Degree $\ell$')
ax[1].set_title('Velocity Power (m=0)')
ax[2].plot(power[:,rind,tind,2])
ax[2].set_xlabel(r'Degree $\ell$')
ax[2].set_title('Velocity Power ( total - {m=0} )')
plt.tight_layout()
plt.show()
fig, ax = plt.subplots()
ss = Shell_Spectra(istring)
help(ss)
mmax = ss.mmax
lmax = ss.lmax
power_spectrum = numpy.zeros((lmax+1,mmax+1),dtype='float64')
for i in range(1,4): # i takes on values 1,2,3
qind=ss.lut[i]
complex_spectrum = ss.vals[:,:,rind,qind,tind]
power_spectrum = power_spectrum+numpy.real(complex_spectrum)**2 + numpy.imag(complex_spectrum)**2
power_spectrum = numpy.transpose(power_spectrum)
tiny = 1e-6
img=ax.imshow(numpy.log10(power_spectrum+tiny), origin='lower')
ax.set_ylabel('Azimuthal Wavenumber m')
ax.set_xlabel(r'Degree $\ell$')
ax.set_title('Velocity Power Spectrum')
#colorbar ...
cbar = plt.colorbar(img) # ,shrink=0.5, aspect = 15)
cbar.set_label('Log Power')
tick_locator = ticker.MaxNLocator(nbins=5)
cbar.locator = tick_locator
cbar.update_ticks()
cbar.ax.tick_params() #font size for the ticks
plt.show()
```

# IX. Point Probes#

**Summary:** Point-wise sampling of desired output variables.

**Subdirectory:** Point_Probes

**main_input prefix:** point_probe

**Python Class:** Point_Probes

**Additional Namelist Variables:**

point_probe_r : radial indices for point-probe output

point_probe_theta : theta indices for point-probe output

point_probe_phi : phi indices for point-probe output

point_probe_r_nrm : normalized radial coordinates for point-probe output

point_probe_theta_nrm : normalized theta coordinates for point-probe output

point_probe_phi_nrm : normalized phi coordinates for point-probe output

point_probe_cache_size : number of time-samples to save before accessing the disk

Point-probes allow us to sample a simulation at an arbitrary set of points. This output type serves two purposes: 1. It provides an analog to laboratory measurements where slicing and averaging are difficult, but taking high-time-cadence using (for example) thermistors is common-practice. 2. It provides an alternative method of slicing a model ( for when equatorial, meridional, or shell slices do yield the desired result).

## IX.1 Specifying Point-Probe Locations#

Point-probe locations are indicated by specifying a grid. The user does not supply a set of ordered coordinates (r,theta,phi). Instead, the user specifies nodes on the grid using the namelist variables described above. Examples follow.

**Example 1: 4-point Coarse Grid**

This example would produce point probes at the four coordinates { (0.25, 0.5, 0.2), (0.25, 0.5, 0.8), (0.5, 0.5, 0.2), (0.5,0.5,0.8) } (r,theta,phi; normalized coordinates).

**Example 2: “Ring” in Phi**

This example describes a ring in longitude, sampled at mid-shell, in the equatorial plane. We have made use of the positional range feature here by indicating normalized phi coordinates of 0.0, -1.0. Rayleigh intreprets this as an instruction to sample all phi coordinates.

** Example 3: 2-D Surface in (r,phi) **

This example uses the positional range feature along with normalized coordinates to generate a 2-D slice in r-phi at theta = 45 degrees (theta_nrm = 0.25). Using the syntax 0,-1.0 instructs *Rayleigh* to grab all r and phi coordinates.

** Example 4: 3-D Meridional “Wedges” **

This example generates two 3-D wedges described by all r,theta points and all longitudes in the ranges [72 deg, 108 deg] and [252 deg, 288 deg].

## IX.2 Point-Probe Caching#

When performing sparse spatial sampling using point-probes, it may be desireable to output with a high-time cadence. As this may cause disk-access patterns characterized by frequent, small writes, the point-probes are programmed with a caching feature. This feature is activated by specifing the **point_probe_cache_size** variable in the output namelist.

This variable determines how many time-samples are saved in memory before a write is performed. Its default value is 1, which means that the disk is accessed with a frequency of **point_probe_frequency**. If the cache size is set to 10 (say), then samples are still peformed at **point_probe_frequency** but they are only written to disk after 10 have been collected in memory.

**NOTE:** Be sure that **point_probe_cache_size** divides evenly into **point_probe_nrec**.

## IX.3 Example: Force-Balance with Point Probes#

Our example input file specifies a coarse, six-point grid. Examining the *main_input* file, we see that all variables necessary to examine the force balance in each direction have been specified. (see *rayleigh_output_variables.pdf* for mathematical formulae):

Menu Code |
Description |
---|---|

1 |
Radial Velocity |

2 |
Theta Velocity |

3 |
Phi Velocity |

1201 |
Radial Advection (v dot grad v) |

1202 |
Theta Advection |

1203 |
Phi Advection |

1216 |
Buoyancy Force (ell=0 component subtracted) |

1219 |
Radial Coriolis Force |

1220 |
Theta Coriolis Force |

1221 |
Phi Coriolis Force |

1228 |
Radial Viscous Force |

1229 |
Theta Viscous Force |

1230 |
Phi Viscous Force |

**Note that the pressure force appears to be missing.** This is not an oversight. The diagnostic nature of the Pressure equation in incompressible/anelastic models, coupled with the second-order Crank-Nicolson time-stepping scheme, means that the pressure field can exhibit an even/odd sawtoothing in time. The *effective* pressure force (as implemented through the Crank-Nicolson scheme) is always a weighted average over two time steps **and is always well-resolved in time**.

When sampling at regular intervals as we have here, if we directly sample the pressure force, we will sample either the high or low end of the sawtooth envelope, and the force balance will be off by a large factor. The easiest fix is to output the velocity field and compute its time derivative. This, in tandem with the sum of all other forces, can be used to calculate the effective pressure as a post-processing step. The (undesireable) alternative is to output once every time step and compute the effective pressure using the Crank-Nicolson weighting.

We demonstrate how to compute the effective pressure force via post-processing in the example below.

```
[ ]:
```

```
from rayleigh_diagnostics import Point_Probes, build_file_list
import numpy
from matplotlib import pyplot as plt
#Decide which direction you want to look at (set direction = {radial,theta, or phi})
#This is used to determine the correct quantity codes below
radial = 0
theta = 1
phi = 2
direction=radial
# Build a list of all files ranging from iteration 0 million to 1 million
files = build_file_list(0,1000000,path='Point_Probes')
nfiles = len(files)-1
for i in range(nfiles):
pp = Point_Probes(files[i],path='')
if (i == 0):
nphi = pp.nphi
ntheta = pp.ntheta
nr = pp.nr
nq = pp.nq
niter = pp.niter
vals=numpy.zeros( (nphi,ntheta,nr,nq,niter*nfiles),dtype='float64')
time=numpy.zeros(niter*nfiles,dtype='float64')
vals[:,:,:,:, i*niter:(i+1)*niter] = pp.vals
time[i*niter:(i+1)*niter]=pp.time
istring='00040000' # iteration to examine
help(pp)
##################################################
# We choose the coordinate indices **within**
# the Point-Probe array that we want to examine
# These indices start at zero and run to n_i-1
# where n_i is the number of points sampled in
# the ith direction
# Use help(pp) after loading the Point-Probe file
# to see the Point-Probe class structure
pind = 0 # phi-index to examine
rind = 0 # r-index to examine
tind = 0 # theta-index to examine
pp = Point_Probes(istring)
lut = pp.lut
nt = pp.niter
#######################################################################
# Grab velocity from the point probe data
u = vals[pind,0,rind,pp.lut[1+direction],:]
dt=time[1]-time[0]
###########################################################################
# Use numpy to compute time-derivative of u
# (necessary to compute a smooth effective pressure without outputing every timestep)
#Depending on Numpy version, gradient function takes either time (array) or dt (scalar)
try:
dudt = numpy.gradient(u,time)
except:
dt = time[1]-time[0] # Assumed to be constant...
dudt = numpy.gradient(u,dt)
################################################################
# Forces (modulo pressure)
# Note the minus sign for advection. Advective terms are output as u dot grad u, not -u dot grad u
advec = -vals[ pind, tind, rind, lut[1201 + direction], :]
cor = vals[ pind, tind, rind, lut[1219 + direction], :]
visc = vals[ pind, tind, rind, lut[1228 + direction], :]
forces = visc+cor+advec
if (direction == radial):
buoy = vals[ pind, tind, rind, lut[1216], :]
forces = forces+buoy
############################################3
# Construct effective pressure force
pres = dudt-forces
forces = forces+pres
############################################################
# Set up the plot
yfsize='xx-large' # size of y-axis label
ustrings = [r'u_r', r'u_\theta', r'u_\phi']
ustring=ustrings[direction]
dstring = r'$\frac{\partial '+ustring+'}{\partial t}$'
fstrings = [r'$\Sigma\,F_r$' , r'$\Sigma\,F_\theta$' , r'$\Sigma\,F_\phi$' ]
fstring = fstrings[direction]
diff_string = dstring+' - '+fstring
pstring = 'pressure'
cstring = 'coriolis'
vstring = 'viscous'
bstring = 'buoyancy'
fig, axes = plt.subplots(nrows=2, figsize=(7*2.54, 9.6))
ax0 = axes[0]
ax1 = axes[1]
########################################
# Upper: dur/dt and F_total
#mpl.rc('xtick', labelsize=20) --- still trying to understand xtick label size etc.
#mpl.rc('ytick', labelsize=20)
ax0.plot(time,forces, label = fstring)
ax0.plot(time,pres,label=pstring)
ax0.plot(time,cor,label=cstring)
ax0.plot(time,visc,label=vstring)
if (direction == radial):
ax0.plot(time,buoy,label=bstring)
ax0.set_xlabel('Time', size=yfsize)
ax0.set_ylabel('Acceleration', size=yfsize)
ax0.set_title('Equilibration Phase',size=yfsize)
ax0.set_xlim([0,0.1])
leg0 = ax0.legend(loc='upper right', shadow=True, ncol = 1, fontsize=yfsize)
##########################################
# Lower: Numpy Gradient Approach
ax1.plot(time,forces,label=fstring)
ax1.plot(time,pres,label=pstring)
ax1.plot(time,cor,label=cstring)
ax1.plot(time,visc,label=vstring)
if (direction == radial):
ax1.plot(time,buoy,label=bstring)
ax1.set_title('Late Evolution',size=yfsize)
ax1.set_xlabel('Time',size=yfsize)
ax1.set_ylabel('Acceleration', size =yfsize)
ax1.set_xlim([0.2,4])
leg1 = ax1.legend(loc='upper right', shadow=True, ncol = 1, fontsize=yfsize)
plt.tight_layout()
plt.show()
```

# X. Modal Outputs#

**Summary:** Spherical Harmonic Spectral Coefficients sampled at discrete radii and degree ell.

**Subdirectory:** SPH_Modes

**main_input prefix:** sph_mode

**Python Classes:** SPH_Modes

**Additional Namelist Variables:**

sph_mode_levels (indicial) : indices along radial grid at which to output spectral coefficients.

sph_mode_levels_nrm (normalized) : normalized radial grid coordinates at which to output spectral coefficients.

sph_mode_ell : Comma-separated list of spherical harmonic degree ell to output.

The Modal output type allows us to output a restricted set of complex spherical harmonic coefficients at discrete radii. For each specified ell-value, all associated azimuthal wavenumbers are output.

This output can be useful for storing high-time-cadence spectral data for a few select modes. In the example below, we illustrate how to read in this output type, and we plot the temporal variation of the real and complex components of radial velocity for mode ell = 4, m = 4.

Examining the *main_input* file, we see that the following output values have been denoted for the Shell Spectra (see *rayleigh_output_variables.pdf* for mathematical formulae):

Menu Code |
Description |
---|---|

1 |
Radial Velocity |

2 |
Theta Velocity |

3 |
Phi Velocity |

We also see that ell=2,4,8 have been selected in the *main_input* file, leading to power at the following modes:

ell-value |
m-values |
---|---|

2 |
0,1,2 |

4 |
0,1,2,3,4 |

8 |
0,1,2,3,4,5,6,7,8 |

```
[ ]:
```

```
from rayleigh_diagnostics import SPH_Modes, build_file_list
import matplotlib.pyplot as plt
import numpy
qind = 1 # Radial velocity
rind = 0 # First radius stored in file
files = build_file_list(0,1000000,path='SPH_Modes')
nfiles = len(files)
for i in range(nfiles):
spm = SPH_Modes(files[i],path='')
if (i == 0):
nell = spm.nell
nr = spm.nr
nq = spm.nq
niter = spm.niter
lvals = spm.lvals
max_ell = numpy.max(lvals)
nt = niter*nfiles
vr = spm.lut[qind]
vals=numpy.zeros( (max_ell+1,nell,nr,nq,nt),dtype='complex64')
time=numpy.zeros(nt,dtype='float64')
vals[:,:,:,:, i*niter:(i+1)*niter] = spm.vals
time[i*niter:(i+1)*niter]=spm.time
help(spm)
#####################################################3
# Print some information regarding the bookkeeping
print('...........')
print(' Contents')
print(' nr = ', nr)
print(' nq = ', nq)
print(' nt = ', nt)
for i in range(nell):
lstring=str(lvals[i])
estring = 'Ell='+lstring+' Complex Amplitude : vals[0:'+lstring+','+str(i)+',0:nr-1,0:nq-1,0:nt-1]'
print(estring)
print(' First dimension is m-value.')
print('...........')
######################################
# Create a plot of the ell=4, m=4 real and imaginary amplitudes
radius = spm.radius[rind]
lfour_mfour = vals[4,1,rind,vr,:]
fig, ax = plt.subplots()
ax.plot(time,numpy.real(lfour_mfour), label='real part')
ax.plot(time,numpy.imag(lfour_mfour), label='complex part')
ax.set_xlabel('Time')
ax.set_ylabel('Amplitude')
rstring = "{0:4.2f}".format(radius)
ax.set_title(r'Radial Velocity ( $\ell=4$ , m=4, radius='+rstring+' ) ')
ax.legend(shadow=True)
ax.set_xlim([0.5,4.0])
plt.show()
```

```
[ ]:
```

```
```