planetmodel

Spherically layered planetary models: a skeleton of boundary radii, a geometry that places it in the physical world through one continuous mapping, fields on each layer, and the meshes that hand a model to a solver.

The library covers the skeleton, the geometry and its mappings; fields on one interval with an exact polynomial algebra; the model with its units, named model types (PREM from its polynomials, simple layered models, any mineos deck) and the mixins that complete them with moduli, velocities, gravity and linear rheologies; a radial spectral-element mesh; 2D and 3D meshes via gmsh with a manifest and export to MFEM; and a sub-package that consumes the radial mesh, planetmodel.randomfield, which draws Matern random fields on balls, annuli and layers. The Love numbers of a model, the loading and tidal problem solved on the radial mesh, are computed by pyslfp, which builds on this library. A netCDF file for 3D models is next.

The ideas

A skeleton. A strictly increasing list of boundary radii, possibly starting above zero for a shell. It answers geometric questions (intervals, where a radius lies, with the two sides of a boundary left to the caller) and supports surgery: refine, truncate, hollow, extend, coarsen.

A geometry. A skeleton, one mapping from the reference ball to the physical body, and the names of layers and interfaces. The mapping must be orientation-preserving, continuous, and kinked only on skeleton boundaries; a geometry checks those invariants when it is built. The shipped mapping is the radial stretch m(X) = (r + h) e_r driven by any callable h(r, theta, phi), with closed forms for its deformation gradient, Jacobian, validity, inverse and linearisation; any object with __call__, deformation_gradient and jacobian is a mapping.

A field on one interval. Data on one layer: an interval, a character (the tensor rank and weight that say how it transforms), a name, and evaluate(r, theta, phi, *, frame) giving components in the local spherical frame or in Cartesian ones. A discontinuity is two layers asked separately. Radial fields sit on layer functions whose algebra is exact on polynomials, so PREM's moduli rho v^2 are exact polynomials; analytic formulas, pointwise compositions, and fields pushed forward through a mapping are fields too. A field may be complex-valued, which is what a model frozen at a frequency holds.

One model class. A geometry with a bag of fields on every layer. What a name means comes from the shipped vocabulary or the specs a model is given. A model type (PREM, LayeredIsotropicElastic, MineosModel, or your own) is a class derived from Model alone; the mixins of planetmodel.behaviours add the shared derivations as methods, written once as free functions: the Love moduli beside the velocities, the elastic tensor and its Voigt average, gravity, and the constant-Q and Maxwell rheologies frozen at a frequency.

Units in one place. The model's Scales say what one stored unit is in SI, and converted re-expresses the whole model by name, exactly for polynomials; G is read in the model's units. Nothing else names a unit: radii are numbers, every tolerance is relative, and the meshers hand the geometry's numbers to gmsh unchanged.

Two meshers. planetmodel.mesh1d lays Gauss-Lobatto-Legendre elements along the radius with every skeleton boundary an element boundary, evaluates a model's fields and gravity on its nodes, and samples a model on an angular grid. planetmodel.mesh3d meshes a geometry, full or hollow, in 2D or 3D, with shells outside it, writes a JSON manifest saying what every attribute means, and exports the mesh, the mapping's displacement and the model's fields to MFEM.

Executable contracts. planetmodel.testing holds one check_* function per protocol; the shipped implementations and yours are held to the same call.

Installing

pip install planetmodel                     # numpy and scipy only
pip install 'planetmodel[meshing]'          # 2D and 3D meshes via gmsh
pip install 'planetmodel[mfem]'             # export to MFEM (PyMFEM)
pip install 'planetmodel[harmonics]'        # grid transforms via pyshtools and ducc0
pip install 'planetmodel[plot]'             # matplotlib, for the figures
pip install 'planetmodel[notebook]'         # ipykernel, to run tutorials cell by cell

Python 3.12 or later. Nothing optional is imported by import planetmodel.

Twelve lines

import numpy as np
from planetmodel import PREM, RadialMesh, flattening, gravity

m = PREM()                                                  # exact polynomials, SI
oc, mantle = m.layer("outer_core"), m.layer("lowermost_mantle")
cmb = m.geometry.interface("cmb").radius
print(oc["rho"](cmb), mantle["rho"](cmb), m.is_fluid("outer_core"))  # both sides
print(m.elastic_moduli("lowermost_mantle")(cmb, 0.3, 0.0)[:3, :3] / 1e9)  # Voigt, GPa
print(gravity(m, [cmb, 6371e3]), m.moduli_at("lower_mantle", 2 * np.pi / 3600)["L"](5e6))

nd = m.nondimensionalised().stretched(flattening(1 / 300, rmax=1.0))
mesh = RadialMesh(nd, ngll=5, drmax=0.05)
print(nd.G, mesh.nodal(nd, "rho").shape, nd.geometry.validity())

Where to go next

• examples/tutorials/: eleven walkthroughs, from a skeleton to random fields and deck files, each a # %% script that runs headless.
• src/planetmodel/mesh3d/manifest.py: the manifest beside every mesh, its schema described from the consumer's side.
• CONTRIBUTING.md: the development setup, the hooks, the test selections and how a release is made.

Tests

poetry run pytest                                   # the fast suite
poetry run pytest -m "not slow"                     # with gmsh and MFEM
poetry run ruff check .

Licence

BSD-3