Metadata-Version: 2.4
Name: planetmodel
Version: 1.2.4
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 planet described as concentric
layers with physical properties on each — and the meshes that hand such
a model to a numerical solver.

With this library you can build a model (PREM out of the box, simple
layered models, any mineos deck file, or a model type of your own), ask
it for its properties with the discontinuities treated honestly, deform
it with ellipticity or surface topography, and mesh it: a
spectral-element mesh along the radius, or a 2D or 3D finite-element
mesh for [MFEM](https://mfem.org), refined where you need it. A
sub-package, `planetmodel.randomfield`, draws Matérn random fields on
balls, annuli and layers.

## Design

**A skeleton.** The list of boundary radii, increasing, possibly
starting above zero for a hollow shell. It answers geometric questions
(which layer a radius lies in, the intervals between boundaries) and
supports surgery: refine, truncate, hollow, extend, coarsen.

**A geometry.** A skeleton placed in the physical world by one
continuous mapping from the reference (spherical) body to the physical
one, plus names for the layers and interfaces. Ellipticity or surface
topography is therefore a property of the whole model, stated once,
rather than something each field has to know about. The shipped mapping
is the radial stretch `m(X) = (r + h) e_r` driven by any callable
`h(r, theta, phi)`; any object with `__call__`,
`deformation_gradient` and `jacobian` serves as a mapping, and a
geometry checks the invariants (orientation preserved, continuous,
kinked only on boundaries) when it is built.

**A field on one layer.** An interval, a name, a character (the tensor
rank and weight that say how the values transform), and
`evaluate(r, theta, phi, *, frame)` giving components in the local
spherical frame or in Cartesian ones. A discontinuity is simply two
layers: ask each side and get each side's answer. Radial fields sit on
an algebra that is exact on polynomials — PREM's moduli `rho v^2` are
exact polynomials, not fits — and fields may be complex-valued, which
is what a model frozen at a frequency holds.

**One model class.** A model is a geometry with named fields on every
layer; what a name like `"rho"` means comes from the shipped vocabulary
or from the specs the model is given. A model type (`PREM`,
`LayeredIsotropicElastic`, `MineosModel`, or your own) derives from
`Model` alone — there is no hierarchy of model types — and the mixins
of `planetmodel.behaviours` add the standard derivations: 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
of length, mass and time is in SI, and `converted` re-expresses the
whole model, exactly for polynomials. Nothing else names a unit: radii
are numbers, tolerances are 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 boundary an element boundary,
evaluates a model's fields and gravity on the nodes, and samples a
model on an angular grid.
`planetmodel.mesh3d` meshes a geometry, full or hollow, in 2D or 3D,
with buffer shells outside it when a far-field condition needs room;
element sizes are controlled per interface, and refinements
(`NearPoints`, `near_curve`, or any size function of your own) make the
mesh finer towards a coastline, a source region or a station network
without coarsening anything elsewhere. Every mesh comes with a JSON
manifest saying what each attribute means, and exports to MFEM together
with the mapping's displacement and the model's fields.

**Executable contracts.** `planetmodel.testing` holds one `check_*`
function per protocol, and the shipped implementations and your own 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`.

## A first look

```python
from planetmodel import PREM, RadialMesh, flattening, gravity

m = PREM()                              # exact polynomials, SI units

# Fields live on layers, so a discontinuity gives two honest answers,
# one from each side.
r_cmb = m.geometry.interface("cmb").radius
print(m.layer("outer_core")["rho"](r_cmb))         # density, core side
print(m.layer("lowermost_mantle")["rho"](r_cmb))   # density, mantle side
print(gravity(m, [r_cmb, 6371e3]))                 # gravity at the CMB and surface

# Reshape the model: non-dimensionalise it, flatten it by 1/300, and
# lay a radial spectral-element mesh over it.
nd = m.nondimensionalised().stretched(flattening(1 / 300, rmax=1.0))
mesh = RadialMesh(nd, ngll=5, drmax=0.05)
print(mesh.nodal(nd, "rho").shape)      # density at every mesh node
```

## 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

