Python codes for processing and plotting the output of the NormalModes library.
- Dependencies
- How to install the code
- Files required for post-processing
- How to run a 'quick projection'
- How to run automatic mode identification
- How to run a 'full projection'
- How to do projections in parallel
- Description of the processing steps
- Known issues
- Contributors and version history
- How to contribute
- How to acknowledge this code
It is necessary to install Python and the following Python libraries:
- NumPy and SciPy;
- MatPlotLib and Cartopy for plotting;
- SHTns, for projection into vector spherical harmonics;
- ssrfpy, for interpolation on a sphere.
Make a copy of this repository. It should then be ready to use.
The PlanetaryModels mesh files used by NormalModes are required for post-processing, in a single directory. Post-processing reads the .ele, .node, and .neigh files, although currently only the .node file is actually used.
You should also place a file in this directory called radii.txt. The first line is the state of the outermost spherical shell (solid or fluid). The remaining lines are the radial coordinates of the solid-fluid discontinuities in decreasing order, preceded by the radius of the planet (in km). If the planet has no fluid-solid discontinuities, this file only has two lines (state and radius). An example is shown for an Earth model:
solid
6371.0
3480.0
1221.5
A successful NormalModes run will produce the files listed below, in a single directory. Items 1–4 are necessary for post-processing. These files are detected by NMPostProcess by matching regular expressions, so, to avoid errors, you should not add other files to this directory.
- The mode files. In the non-rotating case, the modes are real and there is one file per mode (e.g.
mod.1_JOB1_pod1_np48_0.1000000_1.000000_1.dat). In the rotating case, the modes are complex and there are two files per mode containing the real and imaginary parts (e.g.model.1_JOB5_pod1_np48_0.1000000_1.000000_imag_1.datandmodel.1_JOB5_pod1_np48_0.1000000_1.000000_real_1.dat). In either case, the modes are listed in order of increasing frequency, as indicated by the final integer in the file name. Each file contains displacement vector field for one mode: the three components of the mode's displacement at each solution point (including second-order points). - The vertex index file (e.g.
mod.1_pod1_np48_vlist.dat). The vertex indices are integers (one per solution point) which map between the solution points and the mesh node points. This is necessary because some node points have displacement specified twice, on either side of a boundary. For more information, see the docstring forprocess.read_eigenvector(). - The vertex attribute file (e.g.
mod.1_pod1_np48_stat.dat). The vertex attributes are integers (one per mesh node) indicating whether the grid point is solid, fluid, or solid-fluid, and whether it is a first- or second-order point*. For more information, see the docstring forprocess.py. - The
sbatchstandard output log (e.g.LLSVP_6222017.txt). This contains detailed progress information as well as printing out the relative errors and the mode frequencies. Alternatively, some versions of NormalModes produce a file*eigs.txtwhich is an increasing list of the squares of the angular frequencies in radians per second. Only one of these two files can be present (becauseprocess.pysearches for a file ending in.txt). - The log file (e.g.
mod.1_0.log) contains brief progress information output by the code. Sometimes it is useful for error messages or timing information. - The
sbatcherror log (e.g.LLSVP_6222017.err). Hopefully this file will be empty or only contain warnings, but if something goes wrong, it will be reported here.
*Note that second-order points are discarded immediately. Typically, this does not degrade the analysis of the displacement patterns, because the number of nodes required to accurately resolve the frequencies is much larger than the number of nodes is required to accurately resolve the displacement patterns.
All the .dat files (items 1, 2 and 3 above) are stored in a Fortran binary format. They can be read in Python using numpy.fromfile(), e.g.
vert_idxs = numpy.fromfile('mod.1_pod1_np48_vlist.dat', dtype = '<i') The file input_NMPostProcess.txt controls the post-processing. The first three lines are
- The path to the directory containing the PlanetaryModels input files.
- The path to the directory containing the NormalModes input files.
- An string specifying the processing options described in the following sections:
quickFind the vector spherical harmonic expansion for a single radial coordinate.
An example of the first three lines of the NMPostProcess input file is here:
/example/model_dir/ # Directory containing model files.
/example/eigvec_dir/ # Directory containing eigenvector files.
quick # Option ('quick' or 'full').
...
To select 'quick projection', set the third line input_NMPostProcess.txt to quick. The remaining two lines should be
- The maximum angular order of the expansion.
- The number of the mode.
For example, to project mode number 15 using a maximum -value of 10, call the command
python3 project.py
with the following input file (comments should be removed)
/example/model_dir/ # Directory containing model files.
/example/eigvec_dir/ # Directory containing eigenvector files.
quick # Option ('quick' or 'full').
10 # l_max (maximum angular order, must be integer).
15 # i_mode (mode ID: integer, 'all', or 'parallel').
Notice that the mode number can also be set as all to loop over all modes (this can be slow for large runs) or parallel to loop over all modes using all available cores.
Choosing the maximum -value should be based on the expected spatial wavelength of the modes, which depends on the upper frequency limit of the NormalModes run. For example, on Earth, we see that modes have
-values less than approximately
, where
where
is the frequency in mHz. In this example, the run had an upper frequency limit of 1.0 mHz, so a maximum
-value of 10 was appropriate. Similar relationships can deduced for other planets using spherically-symmetric codes such as Mineos. Choosing a maximum
-value that is much too high will slow down the processing and may increase numerical noise. If it is too low, the expansion of the spatial patterns will be truncated, causing aliasing and mode misidentification.
The maximum -value also controls the spacing of the regular longitude-latitude grid used in the interpolation step (see following section). The number of latitude grid points is chosen to be
, and the number of longitude grid points is twice this many, which should be sufficient to avoid aliasing.
The projection will create a number of output files in the NormalModes output directory, in a new subdirectory called processed/. Some of them are simply the NormalModes output re-written in a more convenient form:
eigenvalue_list.txtA list of mode number and frequency.indices_discon_*.txtEach of these files contains a list of the indices of the sample points belonging to one of the solid-fluid discontinuities. They are labelled0,1, ... starting at the innermost discontinuity. The last one is the outer surface (this will be the only file if there are no internal discontinuities).indices_interior_*.txtSimilar to the discontinuity files, listing the indices of sample points belonging to each interior region, starting with the innermost region.nodes.txtThe coordinates of the nodes of the mesh (including second-order nodes).
The projected output is a series of SHTns complex vector spherical harmonics, which are written to NumPy binary files quick_spectral_*.npy in the processed/spectral/ subdirectory. The mode number is indicated in the file name. There are three kinds of coefficients (radial , consoidal
and toroidal
) and one header line, so each file stores array of shape
where
is the number of coefficients. The number of coefficients is
. The coefficients are stored in the order used by SHTns. In the rotating case, the real and imaginary parts of the vector field are projected separately, so the output arrays have shape
, with the real coefficients stored first.
The VSH projection can be plotted with the command
python3 plot/plot_displacement.py
which is controlled by the input file input_plotting.txt. The first line of this file controls the type of plot: spatial or spectral, which are described below. The plots will appear in the processed/plots/ subdirectory.
A spatial plot displays a map representation of the three components (radial, consoidal and spheroidal) of the vector field in three panels. Here we have an example of a low-frequency mode of a spherically-symmetric Earth model, neglecting gravity:
We can see this mode is a spheroidal mode (no toroidal component). Note that the consoidal and toroidal components (lower two panels) are surface vector fields with both a direction (indicated by arrows) and a magnitude. Their magnitude is shown using only the positive part of the colour bar. The radial component (top panel) has a fixed direction can be positive (outwards) or negative (inward).
A spatial plot input file looks like this
quick # Option: 'quick' for quick-mode plots.
spatial 90 # 'spatial n_lat' or 'spectral'.
15 # i_mode: The mode ID (integer or 'all').
png # Figure output format ('png' or 'pdf').
real # Which component to plot ('real', 'complex', 'complex_real_only',
# 'complex_imag_only')
where the first line is the type of processing (quick or full), the second line is the plot type (spatial or spectral) and the number of latitude grid points (only used for spatial; see discussion in 'Using the code: Quick projection' for finding the appropriate number of points), the third line is the mode number (can also be all), and fourth line is the output figure format (currently supports png and pdf). The spatial field is calculated by re-projecting the VSH coefficients.
In the rotating case, the displacement field is complex. By setting the last line to complex, both real and imaginary parts will be plotted. You must also specify the rotation period in hours (e.g. complex 23.9345; see https://en.wikipedia.org/wiki/Sidereal_time to see why Earth's rotation period is not 24 hours). If you only wish to plot either the real or imaginary parts of the complex field, instead use complex_real_only or complex_imag_only.
The figure title gives the radial coordinate of the projection (the radius at which the maximum displacement occurs). An optional file shell_names.txt can be placed in the processed/ directory to give more descriptive names to the regions between solid-fluid discontinuities. For example, for Earth (with no crust), an appropriate shell-name file is
mantle
outer core
inner core
For the same example as above, if we change the plot input file to
quick # Option: 'quick' for quick-mode plots.
spectral # 'spatial n_lat' or 'spectral'.
15 # i_mode: The mode ID (integer or 'all').
png # Figure output format ('png' or 'pdf').
real # Which component to plot ('real', 'complex', 'complex_real_only',
# 'complex_imag_only')
we get a spectral plot:
This shows the VSH coefficients (radial , consoidal
and toroidal
). They have been converted from complex to real format (so they range from
to
). (Note that the last line of the plot input file refers to the real/complex nature of the eigenvectors, not the coefficients.) We can see that this mode is dominantly spheroidal with
. We also know the frequency of the mode (0.415 mHz) from the
eigenvalue_list.txt file. Comparison with calculations of a spherically-symmetric Earth then allows us to identify this mode as one of the seven nearly-degenerate modes of Earth's multiplet.
quick # Option: 'quick' for quick-mode plots.
animation 90 # The number is n_lat, same as for 'spatial' plots.
15 # The mode ID (integer; 'all' not allowed for animation).
none # Figure output format line (ignored for animation).
real # What type of mode is being plotted ('real' or 'complex').
The dominant -value and mode type (spheroidal or toroidal) can be determined automatically from the power distribution VSH coefficients. To do this, process each mode as described above using
process.py. Then run
python3 characterise.py
This first creates a file characterisation_quick.npy in the processed output directory. This stores an array of shape where
is the number of modes. The rows of the array are:
- 0, 1, 2: The spectral 'power',
, in the
,
and
components, for example
.
- 3, ..., 3 +
: The spectral 'power' in each
-band, defined as
.
- 4 +
- 5 +
This information is then used to identify the modes. The mode identification is saved in the file mode_ids.txt. The file has one row per mode. Each row contains four integers, giving
- The mode number.
- The
- The mode type (0 for radial, 1 for spheroidal, and 2 for toroidal).
- The shell within which the maximum displacement occurs. This is used for distinguishing toroidal modes in different solid regions.
For example, the mode plotted in the previous examples has the following line:
15 3 1 1
meaning mode 15 has , is spheroidal, with maximum displacement in the outer core.
To plot the mode identifications, run
python3 plot/plot_mode_diagram.py
which produces a plot such as
This plot shows all of the spheroidal, toroidal and radial modes, as indicated in the legend. Note that 'T0' indicates 'toroidal modes from shell 0' where shell 0 is the outermost shell, in this case the mantle. At higher frequencies, we would also expect T2 modes (inner-core toroidal modes).
For spherically-symmetrical models, modes are expected to occur in 'multiplets' containing modes of identical frequency. More generally, the multiplets will show small frequency splittings. The plotting code attempts to group the modes into clusters and will label multiplets with too many modes (in red) or too few modes (in green) for the given
-value. In this example, we see that the multiplet
is missing 3 modes. This is because the mode has a frequency of almost exactly 1.0 mHz so falls on the boundary of the eigenvalue filter. All of the other modes are recovered (note that toroidal modes with
or
and the modes
and
cannot be generated by an internal source and are not returned by NormalModes).
The mode clustering information can be helpful for identifying numerical errors (usually caused by a too-coarse mesh) or coupled modes in non-spherical models. However, in more complicated models at higher frequencies it becomes difficult to automate the clustering, and the clustering information may not be meaningful. Users can try changing the clustering tolerance f_tol in the script characterise.py.
'Full projection' calculates VSH coefficients on spherical surfaces with various radii from the centre of the planet to the surface, instead of the one radius used in quick projection. This means that full projection is slower, but it provides a full description of the mode displacement field throughout the planet. This can give more robust mode identifications.
Using 'full projection' is very similar to 'quick projection'. The only differences are in the input file (input_NMPostProcess.txt). The third line must be changed to full n_radii where n_radii is the number of radial sampling points. For example:
/example/model_dir/ # Directory containing model files.
/example/eigvec_dir/ # Directory containing eigenvector files.
full 20 # Option and number of radii (integer).
10 # l_max (maximum angular order, must be integer).
15 # i_mode (mode ID: integer, 'all', or 'parallel').
will run a full projection at 20 radii for mode 15. The code will choose the radii to be as evenly-distributed as possible, with at least three radii in each shell (one at the outer boundary, one at the inner boundary, and one within the shell).
Plotting displacement maps and coefficients is very similar to 'quick projection', except that you must also specify the integer ID of the radial 'slice' you wish to plot (or use 'all' to plot all the slices). So, for example, the following input file
full all # Option and radius ID (integer or all).
spatial 90 # 'spatial n_lat' or 'spectral'.
15 # i_mode: The mode ID (integer or 'all').
png # Figure output format ('png' or 'pdf').
would produce a full-projection spatial plot for mode 15 at all of the radial slices defined during the processing step. Note that the mode ID can also be 'all', but this will take a long time. As before, changing 'spatial' to 'spectral' will yield a spectral plot for each specified depth slice. The animation below shows the displacement pattern for low-frequency mode of a large planet:
The indiviudal PNG files were stitched into a GIF using the following ImageMagick command:
convert -dispose previous -delay 20 -loop 0 displacement_full_00015_0* displacement_full_00015.gif
where -dipose previous is necessary with the transparent background.
To plot the variation in the U, V and W components as a function of radius, run the command
python3 plot/plot_eigenfunction.py
As before, the mode is specific in the input_plotting.txt file (the other lines are ignored). This will create a plot such as this:
The calculation of the VSH coefficients of one mode is independent of the calculation of the coefficients of any other mode. Therefore, it is straightforward to process all of the modes in parallel. This speeds up the calculation, which is especially useful for 'full projection'.
As mentioned before, to use parallel processing, change the fifth line of input_NMPostProcess.txt from all to parallel.
Parallel processing is often more effective if you have access to a cluster with many cores. An example job script is given in example_slurm.sbatch for a Slurm cluster. This script would be submitted to the cluster with the command sbatch example_slurm.sbatch. The script is very simple, but you must also make sure that the NMPostProcesss file is correct, the required Python modules are available on the cluster, and any local cluster rules are obeyed.
Vector spherical harmonics (VSHs) are the natural basis to describe the modes of spherically-symmetrical planets (see Dahlen and Tromp, 1998, chapter 8). They are also a good starting point for non-spherically-symmetric planets if the deviations from spherical symmetry are small. The main task of NMPostProcess is to project the vector displacement fields (eigenfunctions) of the modes calculated by NormalModes into VSHs to allow the modes to be identified. In particular, from the VSH expansion it is trivial to separate toroidal and spheroidal modes, and determine the -value of the mode. For these tasks it is usually sufficient to use the VSH expansion onto a single spherical surface. This is called 'quick projection'. Sometimes it is helpful to calculate the VSH expansion at various radial distances spanning the interior of the planet; we call this 'full projection'.
The first step is to process information about the input model, creating convenience files with sample coordinates and indices identifying each region and boundary. Then the projection is performed as follows:
- Load the eigenvector samples and the coordinates of these samples.
- Discard samples from second-order nodes1.
- Find the radial coordinate of the sample with the largest displacement. This is chosen as the radius of the sphere onto which the projection is done2.
- Interpolate the 3D displacement field onto a regular grid over the sphere with the chosen radius (discontinuities are treated carefully).
- Convert the Cartesian components of the displacement field to polar components (north, east, radial).
- Transform the polar components to VSH coefficients.
- Save the spectral displacement fields of the chosen radius.
1This step is not necessary, but it decreases the computational cost and in most circumstances will not impede the mode identification (a large number of nodes is required to accurately compute the mode frequencies, but the displacement patterns are usually long-wavelength and can be calculated and described with fewer nodes).
2The displacement patterns of a mode of a spherically-symmetrical planet is separable into a function of radius and a function of angular position. Therefore the results of the projection will differ only in amplitude if a different radius is used. However, it is still wise to choose the radius of maximum displacement, to avoid radial nodes where the displacement is zero.
Full projection uses the same processing steps as quick projection (steps 4–7), except there is an additional loop over the different values of radius.
Drawing filled contour maps with Cartopy can cause errors such as
AttributeError: 'PreparedGeometry' object has no attribute 'is_valid'
and other times it will fail to draw some contour levels without raising an error. Usually these issues go away if you change the map projection or the number of contour levels.
- Version 1: 27th August 2020: Created by Harry Matchette-Downes. Includes 'quick' processing, and plotting of maps, spectra and mode diagrams.
- Version 2: 24th September 2020: Added 'full' projection.
We intend to add the following features:
- Support for spheroidal (flattened) planets.
Please feel free to contribute in the form of bug reports, bug fixes, feature requests, or new code.
If you use this code as part of published work, please include an acknowledgement such as 'the NormalModes output was projected into vector spherical harmonics using the NMPostProcess software (https://github.com/harrymd/NMPostProcess)'.