Metadata-Version: 2.4
Name: planetmodel
Version: 1.2.0
Summary: Spherically layered planetary models: reference bodies, fields, mappings and meshes
License-Expression: BSD-3-Clause
License-File: LICENSE
Keywords: geophysics,seismology,planetary models,PREM,spectral elements,meshing
Author: David Al-Attar
Author-email: da380@cam.ac.uk
Requires-Python: >=3.12
Classifier: Development Status :: 3 - Alpha
Classifier: Intended Audience :: Science/Research
Classifier: Programming Language :: Python :: 3
Classifier: Programming Language :: Python :: 3.12
Classifier: Programming Language :: Python :: 3.13
Classifier: Topic :: Scientific/Engineering :: Physics
Provides-Extra: harmonics
Provides-Extra: meshing
Provides-Extra: mfem
Provides-Extra: netcdf
Provides-Extra: notebook
Provides-Extra: plot
Requires-Dist: ducc0 (>=0.41) ; extra == "harmonics"
Requires-Dist: gmsh (>=4.13) ; extra == "meshing"
Requires-Dist: ipykernel (>=6.29) ; extra == "notebook"
Requires-Dist: matplotlib (>=3.11,<4) ; extra == "plot"
Requires-Dist: mfem (>=4.8) ; extra == "mfem"
Requires-Dist: netCDF4 (>=1.7) ; extra == "netcdf"
Requires-Dist: numpy (>=2.0,<3)
Requires-Dist: pyshtools (>=4.14) ; extra == "harmonics"
Requires-Dist: scipy (>=1.15,<2)
Project-URL: Homepage, https://github.com/da380/planetmodel
Project-URL: Repository, https://github.com/da380/planetmodel
Description-Content-Type: text/markdown

# 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](https://github.com/da380/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

```python
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

