planetmodel.mesh1d package

Submodules

planetmodel.mesh1d.gll module

Gauss-Lobatto-Legendre nodes, weights and Lagrange machinery.

Everything here lives on the reference element [-1, 1] and knows nothing of models or physics: the GLL quadrature rule, the differentiation matrix of the associated Lagrange basis, and evaluation of that basis at arbitrary points.

planetmodel.mesh1d.gll.gll_points_weights(ngll: int) → tuple[ndarray, ndarray]

Gauss-Lobatto-Legendre points and weights on [-1, 1].

The ngll >= 2 nodes are the endpoints plus the ngll - 2 zeros of P’_N with N = ngll - 1, computed stably as the zeros of the Jacobi polynomial P^(1,1)_(N-1); the weights are 2 / (N (N+1) P_N(x)^2). The quadrature integrates polynomials of degree <= 2N - 1 exactly.

planetmodel.mesh1d.gll.lagrange_basis(nodes: ArrayLike, x: ArrayLike) → ndarray

Values L[j, i] = l_i(x_j) of the Lagrange basis at points x.

Uses the second barycentric form for stability; points coinciding with a node, or so near one that the barycentric quotient would overflow, take that node’s basis exactly. Useful for interpolating a spectral-element solution off its nodes.

planetmodel.mesh1d.gll.lagrange_derivative_matrix(x: ArrayLike) → ndarray

Derivative matrix D[k, i] = l_i’(x_k) on the GLL nodes x.

Closed form for GLL nodes (e.g. Canuto et al. 1988): off the diagonal D[k, i] = P_N(x_k) / (P_N(x_i) (x_k - x_i)); the diagonal vanishes except for D[0, 0] = -N(N+1)/4 and D[N, N] = +N(N+1)/4. The index convention (derivative of basis function i evaluated at node k) matches the ‘hp(knode, inode)’ arrays of standard SEM codes, so quadrature sums read sum_k w_k f(x_k) D[k, i] D[k, j].

planetmodel.mesh1d.gravity module

Gravity and mass of a spherically symmetric model, through its fields.

For a spherically symmetric density the acceleration is g(r) = G M(r) / r^2 with M(r) = 4 pi int_0^r rho(s) s^2 ds, the mass inside r accumulated layer by layer. Each layer’s integral is asked of the field algebra: rho * s^2 is a field whose layer function integrates itself, exactly when rho is a polynomial and by quadrature otherwise. This is the gravity of the reference body: rho must be a radial field on every layer, a density that depends on direction is refused by name, and the geometry’s mapping does not enter, so a model placed in the physical world by a non-trivial mapping is answered with the gravity of its spherical reference. Nothing here is a number with a unit: G is the model’s, the radii are in the model’s length, and a hollow model has no mass inside its inner boundary.

planetmodel.mesh1d.gravity.gravity(model: Model, radii: ArrayLike) → np.ndarray

g(r) = G M(r) / r^2 at radii, in the model’s units.

radii broadcast to any shape and must lie in the model, to the geometry’s tolerance; g is zero at the centre and on the inner boundary of a hollow model. M(r) is the antiderivative of each layer’s rho s^2, as gravity_fields has it.

planetmodel.mesh1d.gravity.gravity_fields(model: Model) → tuple[RadialField, ...]

g(r) = G M(r) / r^2 as one radial field per layer, named g.

Each field is exact in value and in its first derivative where the layer’s density is polynomial: M(r) is the polynomial antiderivative of 4 pi rho r^2 plus the mass below the layer, and dg/dr = G (4 pi rho - 2 M / r^3). The fields are continuous across every boundary and vanish at the centre; a model attaches them by with_field(i, "g", field).

planetmodel.mesh1d.gravity.mass(model: Model, *, radius: float | None = None) → float

The mass inside radius, the outer boundary by default.

planetmodel.mesh1d.mesh module

Conforming GLL meshes over an interval and over a skeleton.

Mesh1D is pure geometry on an interval: pinned breakpoints, each span between them cut uniformly into elements no wider than drmax, ngll Gauss-Lobatto-Legendre nodes per element, a global node numbering, the reference-element data a weak form assembles with, and the radial rules between nodes: polynomial, the Lagrange polynomial through each element’s GLL nodes, and linear, straight lines between consecutive nodes, which agree when an element has two nodes. Each rule is offered two ways that agree, as an exact piecewise-polynomial view of nodal values (to_ppoly) and as the weights of an element’s nodes at any point (interpolation_weights). Nodal arrays are stored per element, shape (nspec, ngll), so a shared node carries one value from each side, and at a shared boundary both views take the element above.

RadialMesh lays such a mesh over a Skeleton so that every skeleton boundary is an element boundary, and records the layer each element lies in. Given a model, nodal evaluates the field of each element’s own layer at the element’s nodes, so a shared node carries both one-sided values, and nodal_gravity the model’s gravity there. Everything is numbers: the mesh knows nothing about units.

class planetmodel.mesh1d.mesh.Mesh1D(breakpoints: ArrayLike, *, ngll: int = 5, drmax: float | None = None)

Bases: object

A conforming GLL mesh over an interval with pinned breakpoints.

Each span between consecutive breakpoints is subdivided uniformly into elements no wider than drmax (one element per span when drmax is None). Attributes: nspec elements with left, right and jac (half-widths) of shape (nspec,); node coordinates r of shape (nspec, ngll); the global numbering gmap of shape (nspec, ngll) into nglob = nspec (ngll - 1) + 1 shared nodes with coordinates rglob; and the reference-element data xi, w, deriv (see lagrange_derivative_matrix for the index convention).

element_of(x: float) → int

The element whose interval contains x.

A coordinate on a shared boundary resolves to the element above it; a coordinate outside the mesh clips to the first or last element.

interpolation_weights(x: ArrayLike, *, method: str = 'polynomial', nu: int = 0) → tuple[ndarray, ndarray]

The element of each point, and the weights of that element’s nodes.

Returns (element, weights), element of the shape of x and weights of that shape followed by (ngll,), so that nodal values u of shape (nspec, ngll) take the value sum_i weights[..., i] * u[element, i] at x, the value to_ppoly(u, method=method) gives. The element is element_of’s, the upper one at a shared boundary and the end one beyond the mesh, whose piece continues there. Under "polynomial" the weights are the Lagrange basis of the element’s GLL nodes; under "linear" the two hat functions of the nodes bracketing x, zero elsewhere. With nu=1 they are the weights of the first derivative in x instead: the derivatives of the Lagrange basis, a polynomial of lower degree written in the same basis, or the slope of the bracketing segment.

to_ppoly(nodal: ArrayLike, *, elements: tuple[int, int] | None = None, method: str = 'polynomial') → PPoly

The exact piecewise-polynomial (scipy PPoly) view of nodal values.

nodal has shape (nspec, ngll), or that followed by any trailing shape, in which case the values of the result have the trailing shape; elements, when given, is a half-open (start, stop) element range restricting the view. Under method="polynomial" the result within each element is the interpolating polynomial through the GLL nodes, so its derivative and integral are those of the spectral-element function; under "linear" it is the straight line between each pair of consecutive nodes. Discontinuous data is representable; at a shared breakpoint the upper element’s piece is evaluated, and beyond the ends the end pieces continue. Complex values stay complex. The pieces are in the local monomial form PPoly holds, which loses digits as the order grows on data that is rough on the element (about 1e-11 of the data at ngll 9, 1e-8 at ngll 13), where interpolation_weights stays at rounding; smooth data is held to rounding at every order used.

planetmodel.mesh1d.mesh.RADIAL_METHODS = ('polynomial', 'linear')

the Lagrange polynomial through its GLL nodes, or straight lines between consecutive nodes.

Type:

The radial rules between the nodes of an element

class planetmodel.mesh1d.mesh.RadialMesh(skeleton: Skeleton | Geometry, *, ngll: int = 5, drmax: float | None = None, lmax: int | None = None, rmin: float | None = None, rmax: float | None = None, edges: ArrayLike | None = None, rtol: float = 1e-09)

Bases: Mesh1D

A GLL mesh over [rmin, rmax] of a skeleton.

Exactly one of drmax, lmax and edges sizes the mesh. With drmax or lmax the element boundaries are every skeleton boundary strictly inside the range, and every knot a Geometry’s mapping declares there, plus a uniform subdivision no coarser than drmax, where lmax sets drmax = 0.1 rmax / (lmax + 1); rmin and rmax default to the skeleton’s ends. With edges the element boundaries are taken as given, must lie inside the skeleton, and no element may straddle a skeleton boundary (to rtol of the range). A pin within that tolerance of a skeleton boundary is the boundary. layer records the skeleton layer of every element. A Geometry or a Model is accepted in place of a skeleton, and the knots of its mapping are pinned.

element_at(r: float) → int

The element whose interval contains r; a boundary resolves upward.

nodal(model: Model, name: str, *, nu: int = 0, missing: str = 'refuse') → np.ndarray

A radial field of a model at the nodes, per element: (nspec, ngll) for rank 0, with the stored components appended for higher rank.

Each element takes the field of its own layer, so a node on a skeleton boundary carries both one-sided values; nu asks for the nu-th radial derivative. An element whose layer lacks the name is refused by name, or filled with NaN when missing="nan". Complex where any layer’s field is complex. A new array each call.

nodal_gravity(model: Model) → np.ndarray

The model’s gravity at the nodes, per element, over the whole model.

start_element(l: int, *, eps: float = 1e-08) → int

The first element of the sub-mesh for a degree-l solve.

truncation_radius(l: int, *, eps: float = 1e-08) → float

The radius below which a degree-l solution is negligible.

With the interior decay (r / rmax)^(l + 1), the returned radius satisfies (r / rmax)^(l + 1) = eps.

Module contents

Radial spectral-element meshes: GLL nodes per element over a skeleton.

Mesh1D is the mesh of an interval; RadialMesh lays one over a skeleton with every boundary an element boundary and evaluates a model’s fields and gravity on its nodes. Both read nodal values between the nodes by one of the RADIAL_METHODS. The GLL quadrature, basis and differentiation matrix of the reference element are exported beside them.