Gplates

gplates

ScalarFieldConnector(source, output, *, mesh=None, interpolation=None, gc_collect_frequency=10)

Composition of a Source and an OutputStrategy into a single time- varying scalar field on a mesh.

Construction validates output.requires <= source.provides; this catches the obvious mis-pairings (e.g. PolygonSource + GeothermERFOutput, which needs "age") at the point of wiring rather than at the first get_indicator call.

Two consumers (e.g. an indicator and a geotherm) that share the same Source instance see a single, coherent advance of the source's underlying state per geological age — the source's per-age cache enforces that, so the order of get_indicator calls between consumers is immaterial.

gc_collect_frequency controls how often a full gc.collect() runs in the update loop (every Nth get_indicator). The default of 10 matches gtrack's own internal GC cadence and periodically breaks the pygplates C++ reference cycles without paying a collection on every call. Set it to None to disable the connector-level collect entirely (relying on gtrack's internal collect plus Python's automatic generational GC) when GC is documented as hot and confirmed memory stays bounded; set it to 1 for a lithosphere spin-up or very-long adjoint run where the connector is driven for thousands of ages and the tightest bound on C++ cycle accumulation is wanted (the per-call cost is negligible there against the per-age step_to).

Source code in g-adopt/gadopt/gplates/connectors.py

def __init__(
    self,
    source: Source,
    output: OutputStrategy,
    *,
    mesh: MeshConfig | None = None,
    interpolation: InterpolationConfig | None = None,
    gc_collect_frequency: int | None = 10,
):
    # Validate the pairing of source and output
    if not output.requires <= source.provides:
        missing = output.requires - source.provides
        raise ValueError(
            f"{type(output).__name__} requires {sorted(missing)} which "
            f"{type(source).__name__} does not provide "
            f"(provides={sorted(source.provides)})."
        )
    # Validate the GC collect frequency
    if gc_collect_frequency is not None and gc_collect_frequency < 1:
        raise ValueError(
            f"gc_collect_frequency must be >= 1 or None, "
            f"got {gc_collect_frequency}"
        )

    self.source = source
    self.output = output
    self.mesh = mesh or MeshConfig()
    self.interpolation = interpolation or InterpolationConfig()
    self.gc_collect_frequency = gc_collect_frequency

    # Result cache: keyed on (age, identity of the target_coords buffer).
    # Distinct from the source's per-age cache, which only sees the age
    # axis. GplatesScalarFunction.mesh_coords is allocated once and held
    # for the SF lifetime, so a weakref to that buffer is a sound O(1)
    # key — no need to hash ~24 MB of coordinates on every call.
    self.reconstruction_age: float | None = None
    self._cached_result: np.ndarray | None = None
    self._cached_coords_ref: weakref.ref | None = None
    self._gc_call_counter = 0

get_indicator(target_coords, ndtime)

Evaluate the scalar field at target_coords for time ndtime.

Source code in g-adopt/gadopt/gplates/connectors.py

def get_indicator(
    self,
    target_coords: np.ndarray,
    ndtime: float,
) -> np.ndarray:
    """Evaluate the scalar field at ``target_coords`` for time ``ndtime``."""
    age = self.source.ndtime2age(ndtime)
    self.source.validate_age(age)

    use_cache = self._check_cache(age, target_coords)
    # All ranks must agree, else a rank that misses will enter the
    # collective broadcast inside source.prepare while a rank that hits
    # returns early and hangs the collective.
    use_cache = self.comm.allreduce(use_cache, op=MPI.MIN)
    if use_cache:
        return self._cached_result

    # If the cache is not suitable, prepare the source and compute the result
    sources_dict = self.source.prepare(age)
    result = self._compute(sources_dict, target_coords)
    self._update_cache(age, target_coords, result)

    self._gc_call_counter += 1
    if (self.gc_collect_frequency is not None
            and self._gc_call_counter % self.gc_collect_frequency == 0):
        gc.collect()

    return result

InterpolationConfig(kernel='idw', k_neighbors=50, distance_threshold=0.1, gaussian_sigma=0.04) dataclass

Spherical kNN-interpolation knobs.

The interpolation step normalises both source and target points to the unit sphere, builds a cKDTree over sources, queries k_neighbors nearest neighbours per target node, and computes a weighted average using either inverse-distance ("idw") or Gaussian ("gaussian") weights.

Target nodes whose nearest source seed is farther than distance_threshold (in radians on the unit sphere) are flagged as too_far and handled by the OutputStrategy (which decides whether to fill with a default thickness, switch to mantle temperature, etc).

Defaults: 0.1 rad ≈ 640 km — generous; tighten to e.g. 0.02 rad (~127 km) for sharper boundaries. For polygon sources where the boundary is encoded by zero-thickness halo seeds, the threshold degenerates into a pathological-query guard and the actual roll-off length is controlled by gaussian_sigma instead.

GplatesScalarFunction(function_space, *args, indicator_connector=None, name=None, **kwargs)

Bases: Function

Firedrake Function for scalar indicator/geotherm fields from plate reconstructions.

Wraps any ScalarFieldConnector (built either directly via composition of a Source and an OutputStrategy, or via one of the convenience factories lithosphere_indicator, lithosphere_geotherm, polygon_indicator, polygon_geotherm). Calling update_plate_reconstruction(ndtime) recomputes the field and assigns it onto the underlying Firedrake DoFs.

Only scalar function spaces are supported — vector / tensor elements raise at construction time.

Source code in g-adopt/gadopt/gplates/gplates.py

def __init__(
    self,
    function_space,
    *args,
    indicator_connector: ScalarFieldConnector | None = None,
    name: str | None = None,
    **kwargs
):
    super().__init__(function_space, *args, name=name, **kwargs)

    if indicator_connector is None:
        return

    if not isinstance(indicator_connector, ScalarFieldConnector):
        raise TypeError(
            f"indicator_connector must be a ScalarFieldConnector, "
            f"got {type(indicator_connector).__name__}"
        )

    scalar_element = function_space.ufl_element()
    # Tensor/vector spaces silently produce ill-shaped coords arrays and
    # later fail deep inside Function construction with a cryptic
    # NotImplementedError. Catch it up front.
    if isinstance(scalar_element, (fd.VectorElement, fd.TensorElement)):
        raise TypeError(
            "GplatesScalarFunction requires a scalar function space; "
            f"got {type(scalar_element).__name__}. Use one component or "
            "switch to a scalar element."
        )

    self.indicator_connector = indicator_connector

    mesh = extract_unique_domain(self)
    vector_element = fd.VectorElement(scalar_element)
    coords_space = fd.FunctionSpace(mesh, vector_element)
    coords_func = fd.Function(coords_space)
    coords_func.interpolate(fd.SpatialCoordinate(mesh))
    self.mesh_coords = coords_func.dat.data_ro_with_halos.copy()

update_plate_reconstruction(ndtime)

Recompute and assign the field for the given non-dimensional time.

Skipped silently when the requested age is within delta_t of the last computed age (the connector's own cache would short-circuit anyway; this just avoids paying for the collective MPI dance).

Source code in g-adopt/gadopt/gplates/gplates.py

def update_plate_reconstruction(self, ndtime: float):
    """Recompute and assign the field for the given non-dimensional time.

    Skipped silently when the requested age is within ``delta_t`` of the
    last computed age (the connector's own cache would short-circuit
    anyway; this just avoids paying for the collective MPI dance).
    """
    connector = self.indicator_connector
    age = connector.ndtime2age(ndtime)

    if connector.reconstruction_age is not None:
        if abs(age - connector.reconstruction_age) < connector.delta_t:
            return

    with stop_annotating():
        values = connector.get_indicator(self.mesh_coords, ndtime)
        self.dat.data_with_halos[:] = values

    if annotate_tape():
        self.create_block_variable()

GplatesVelocityFunction(function_space, *args, gplates_connector=None, top_boundary_marker='top', quad_degree=None, solver_parameters=None, **kwargs)

Bases: GPlatesFunctionalityMixin, SolverConfigurationMixin, Function

Extends firedrake.Function to incorporate velocities calculated by Gplates, coming from plate tectonics reconstruction.

GplatesVelocityFunction is designed to associate a Firedrake function with a GPlates connector, allowing the integration of plate tectonics reconstructions. This is particularly useful when setting "top" boundary condition for the Stokes systems when performing data assimilation (sequential or adjoint). Note that we subtract the radial component of the velocity field to ensure that the velocity field is tangential to the surface in a FEM sense.

Attributes:

Name	Type	Description
dbc	DirichletBC	A Dirichlet boundary condition that applies the function only to the top boundary.
boundary_coords	ndarray	The coordinates of the function located at the "top" boundary, normalised, so that it is meaningful for PyGPlates.
gplates_connector		The GPlates connector instance used for fetching plate tectonics data.

Parameters:

Name	Type	Description	Default
function_space		The function space on which the GplatesVelocityFunction is defined.	required
gplates_connector		An instance of a pyGplatesConnector, used to integrate GPlates functionality or data. See Documentation for pyGplatesConnector.	None
top_boundary_marker	defaults to "top"	marker for the top boundary.	'top'
val	optional	Initial values for the function. Defaults to None.	required
name	str	Name for the function. Defaults to None.	required
dtype	data type	Data type for the function. Defaults to ScalarType.	required
quad_degree	int	Quadrature degree for the surface measure. If None, defaults to 2 * element_degree + 1.	None
solver_parameters	dict	Solver parameters for the tangential projection solver. If provided, these will be merged with the default parameters. If None, uses default parameters.	None
kwargs	dict	Additional keyword arguments passed to the Firedrake Function.	{}

Methods:

Name	Description
update_plate_reconstruction	Updates the function values based on plate velocities from GPlates for a given model time. Note model time is non-dimensionalised

Examples:

>>> gplates_function = GplatesVelocityFunction(V,
...                                    gplates_connector=pl_rec_model,
...                                    name="GplateVelocity")
>>> gplates_function.update_plate_reconstruction(ndtime=0.0)

Source code in g-adopt/gadopt/gplates/gplates.py

def __init__(self, function_space, *args, gplates_connector=None, top_boundary_marker="top", quad_degree=None, solver_parameters=None, **kwargs):
    # Initialise as a Firedrake Function
    super().__init__(function_space, *args, **kwargs)

    # When gplates_connector is None (e.g. during adjoint perturbations via
    # _ad_mul/_ad_add, or when Firedrake internally creates subfunctions),
    # skip all the plate reconstruction setup.
    if gplates_connector is None:
        return

    # Set the class name required by SolverConfigurationMixin
    self._class_name = self.__class__.__name__

    # Cache all the necessary information that will be used to assign surface velocities
    # the marker for surface boundary. This is typically "top" in extruded mesh.
    self.top_boundary_marker = top_boundary_marker

    # establishing the DirichletBC that will be used to find surface nodes
    # Note: If one day we are moving towards adaptive mesh, we need to be dynamic with this.
    self.dbc = fd.DirichletBC(
        self.function_space(),
        self,
        sub_domain=self.top_boundary_marker)

    # coordinates of surface points
    self.boundary_coords = fd.Function(
        self.function_space(),
        name="coordinates").interpolate(
            fd.SpatialCoordinate(
                extract_unique_domain(self))).dat.data_ro_with_halos[self.dbc.nodes]
    self.boundary_coords /= np.linalg.norm(self.boundary_coords, axis=1)[:, np.newaxis]

    # Store the GPlates connector
    self.gplates_connector = gplates_connector

    # Store the kwargs and the quad_degree for surface measure
    self.kwargs = kwargs
    self.quad_degree = quad_degree

    # Initialise solver configuration
    self.reset_solver_config(
        default_config=self.tangential_project_solver_parameters,
        extra_config=solver_parameters
    )
    self.register_update_callback(self.setup_solver)

    # Set up the solver for tangential projection
    self.setup_solver()

setup_solver()

Set up the linear variational solver for removing radial component.

Source code in g-adopt/gadopt/gplates/gplates.py

def setup_solver(self):
    """Set up the linear variational solver for removing radial component."""
    # Define the solution in the velocity function space
    # If velocity is discontinuous, we need to use a continuous equivalent
    if not is_continuous(self):
        raise ValueError(
            "GplatesVelocityFunction: Velocity field is discontinuous."
            "Removing the radial component of discontinuous velocity fields is not supported."
        )
    else:
        V = self.function_space()

    self.tangential_velocity = fd.Function(V, name="tangential_velocity")

    # Normal vector (radial direction for spherical geometry)
    r = fd.FacetNormal(V.mesh())

    # Define the test and trial functions for a projection solve
    phi = fd.TestFunction(V)
    v = fd.TrialFunction(V)

    # Define the tangential projection: v_tangential = v - (v·r)r
    tangential_expr = self - fd.dot(self, r) * r

    # Getting ds measure: this can be set by the user through quad_degree parameter
    # We are only looking at the surface measure at the top boundary
    if self.quad_degree is not None:
        degree = self.quad_degree
    else:
        degree = 2 * V.ufl_element().degree()[0] + 1

    domain = self.function_space().mesh()

    if domain.extruded:
        self.ds = fd.ds_t(domain=domain, degree=degree)
    else:
        self.ds = fd.ds(domain=domain, degree=degree, subdomain_id=self.top_boundary_marker)

    # Setting up a manual projection
    # Project onto the tangential space: solve for tangential component
    a = fd.inner(phi, v) * self.ds
    L = fd.inner(phi, tangential_expr) * self.ds

    # Setting up boundary condition, problem and solver
    # The field is only meaningful on the boundary, so set zero everywhere else
    self.interior_null_bc = InteriorBC(V, 0., [self.top_boundary_marker])

    self.problem = fd.LinearVariationalProblem(a, L, self.tangential_velocity,
                                               bcs=self.interior_null_bc,
                                               constant_jacobian=True)
    self.solver = fd.LinearVariationalSolver(
        self.problem,
        solver_parameters=self.solver_parameters,
        options_prefix="gplates_projection"
    )
    self._solver_is_set_up = True

remove_radial_component()

Project velocity to tangent plane by removing radial component using linear variational solver.

This method uses a pre-built linear variational solver for efficiency, avoiding the overhead of rebuilding the projection operator each time.

Source code in g-adopt/gadopt/gplates/gplates.py

def remove_radial_component(self):
    """
    Project velocity to tangent plane by removing radial component using linear variational solver.

    This method uses a pre-built linear variational solver for efficiency, avoiding the overhead
    of rebuilding the projection operator each time.
    """
    # Use the a linear variational solver for removing the radial component
    if not getattr(self, "_solver_is_set_up", False):
        self.setup_solver()

    # Project the velocity to the tangential plane
    self.solver.solve()

    # Copy the tangential velocity back to self
    self.dat.data_wo_with_halos[self.dbc.nodes, :] = (
        self.tangential_velocity.dat.data_ro_with_halos[self.dbc.nodes, :]
    )

PlateModelFiles(continental_polygons=None, static_polygons=None) dataclass

Plate-model polygon file paths for the indicator/geotherm Sources.

Owned by the user and passed to LithosphereSource / PolygonSource. These files drive the continental filtering and plate-id assignment that the Sources need; they are not used by the velocity reconstruction, so they live here rather than on pyGplatesConnector.

pyGplatesConnector(rotation_filenames, topology_filenames, oldest_age, delta_t=1.0, scaling_factor=1.0, nseeds=100000.0, nneighbours=4, kappa=1e-06)

Bases: object

An interface to PyGPlates, used for updating top Dirichlet boundary conditions using plate tectonic reconstructions.

This class provides functionality to assign plate velocities at different geological ages to the boundary conditions specified with dbc. Due to potential challenges in identifying the plate id for a given point with PyGPlates, especially in high-resolution simulations, this interface employs a method of calculating velocities at a number of equidistant points on a sphere. It then interpolates these velocities for points assigned a plate id. A warning is raised for any point not assigned a plate id.

Parameters:

Name	Type	Description	Default
rotation_filenames	Union[str, List[str]]	Collection of rotation file names for PyGPlates.	required
topology_filenames	Union[str, List[str]]	Collection of topology file names for PyGPlates.	required
oldest_age	float	The oldest age present in the plate reconstruction model.	required
delta_t	Optional[float]	The t window range outside which plate velocities are updated.	1.0
scaling_factor	Optional[float]	Scaling factor for surface velocities.	1.0
nseeds	Optional[int]	Number of seed points used in the Fibonacci sphere generation. Higher seed point numbers result in finer representation of boundaries in pygpaltes. Notice that the finer velocity variations will then result in more challenging Stokes solves.	100000.0
nneighbours	Optional[int]	Number of neighbouring points when interpolating velocity for each grid point. Default is 4.	4
kappa		(Optional[float]): Diffusion constant used for don-dimensionalising time. Default is 1e-6.	1e-06

Examples:

>>> connector = pyGplatesConnector(rotation_filenames, topology_filenames, oldest_age)
>>> connector.get_plate_velocities(ndtime=100)

Source code in g-adopt/gadopt/gplates/gplates.py

def __init__(self,
             rotation_filenames,
             topology_filenames,
             oldest_age,
             delta_t=1.,
             scaling_factor=1.0,
             nseeds=1e5,
             nneighbours=4,
             kappa=1e-6):
    """An interface to PyGPlates, used for updating top Dirichlet boundary conditions
    using plate tectonic reconstructions.

    This class provides functionality to assign plate velocities at different geological
    ages to the boundary conditions specified with dbc. Due to potential challenges in
    identifying the plate id for a given point with PyGPlates, especially in high-resolution
    simulations, this interface employs a method of calculating velocities at a number
    of equidistant points on a sphere. It then interpolates these velocities for points
    assigned a plate id. A warning is raised for any point not assigned a plate id.

    Arguments:
        rotation_filenames (Union[str, List[str]]): Collection of rotation file names for PyGPlates.
        topology_filenames (Union[str, List[str]]): Collection of topology file names for PyGPlates.
        oldest_age (float): The oldest age present in the plate reconstruction model.
        delta_t (Optional[float]): The t window range outside which plate velocities are updated.
        scaling_factor (Optional[float]): Scaling factor for surface velocities.
        nseeds (Optional[int]): Number of seed points used in the Fibonacci sphere generation. Higher
                seed point numbers result in finer representation of boundaries in pygpaltes. Notice that
                the finer velocity variations will then result in more challenging Stokes solves.
        nneighbours (Optional[int]): Number of neighbouring points when interpolating velocity for each grid point. Default is 4.
        kappa: (Optional[float]): Diffusion constant used for don-dimensionalising time. Default is 1e-6.

    Examples:
        >>> connector = pyGplatesConnector(rotation_filenames, topology_filenames, oldest_age)
        >>> connector.get_plate_velocities(ndtime=100)
    """

    # Store original filenames for downstream use (e.g., the Sources, which
    # share the rotation/topology files with the velocity reconstruction).
    self.rotation_filenames = rotation_filenames
    self.topology_filenames = topology_filenames

    # Rotation model(s)
    self.rotation_model = pygplates.RotationModel(rotation_filenames)

    # Topological plate polygon feature(s).
    self.topology_features = []
    for fname in topology_filenames:
        for f in pygplates.FeatureCollection(fname):
            self.topology_features.append(f)

    # oldest_age coincides with non-dimensionalised time ndtime=0
    self.oldest_age = oldest_age

    # time window for recalculating velocities
    self.delta_t = delta_t

    # seeds are equidistantial points generate on a sphere
    self.seeds = self._fibonacci_sphere(samples=nseeds)

    # number of neighbouring points that will be used for interpolation
    self.nneighbours = nneighbours

    # PyGPlates velocity features
    self.velocity_domain_features = (
        self._make_GPML_velocity_feature(
            self.seeds
        )
    )

    # last reconstruction time
    self.reconstruction_age = None

    # Factor to scale plate velocities to RMS velocity of model,
    # the definition: RMS_Earth / RMS_Model is used
    # This is specially used in low-Rayleigh-number simulations
    self.scaling_factor = scaling_factor

    # Factor to non-dimensionalise time
    self.kappa = kappa

velocity_non_dim_factor property

u [metre/sec] * velocity_non_dim_factor = u []

time_dim_factor property

d^2/kappa

t [] X time_dim_factor = t [sec]

time_dimDmyrs2sec property

Converts the time dimension factor from Myrs (million years) to seconds.

This method calculates the conversion factor to transform time values from Myrs to seconds using the predefined time_dim_factor and the constant myrs2sec from the pyGplatesConnector module.

Returns:

Name	Type	Description
float		The time dimension factor in seconds.

velocity_dimDcmyr property

Converts velocity from cm/year to non-dimensionalised units.

Returns:

Name	Type	Description
float		The conversion factor for velocity.

get_plate_velocities(target_coords, ndtime)

Returns plate velocities for the specified target coordinates at the top boundary of a sphere, for a given non-dimensional time, by integrating plate tectonic reconstructions from PyGPlates.

This method calculates new plate velocities. It utilizes the cKDTree data structure for efficient nearest neighbor searches to interpolate velocities onto the mesh nodes.

Parameters:

Name	Type	Description	Default
target_coords	array - like	Coordinates of the points at the top of the sphere.	required
ndtime	float	The non-dimensional time for which plate velocities are to be calculated and assigned. This time is converted to geological age and used to extract relevant plate motions from the PyGPlates model.	required

Raises:

Type	Description
Exception	If the requested ndt ime is a negative age (in the future), indicating an issue with the time conversion.

Notes

• The method uses conversions between non-dimensionalised time and geologic age.
• Velocities are non-dimensionalised and scaled for the simulation before being applied.

Source code in g-adopt/gadopt/gplates/gplates.py

def get_plate_velocities(self, target_coords, ndtime):
    """Returns plate velocities for the specified target coordinates at the top boundary of a sphere,
    for a given non-dimensional time, by integrating plate tectonic reconstructions from PyGPlates.

    This method calculates new plate velocities.
    It utilizes the cKDTree data structure for efficient nearest neighbor
    searches to interpolate velocities onto the mesh nodes.

    Args:
        target_coords (array-like): Coordinates of the points at the top of the sphere.
        ndtime (float): The non-dimensional time for which plate velocities are to be calculated and assigned.
            This time is converted to geological age and used to extract relevant plate motions
            from the PyGPlates model.

    Raises:
        Exception: If the requested ndt ime is a negative age (in the future), indicating an issue with
            the time conversion.

    Notes:
        - The method uses conversions between non-dimensionalised time and geologic age.
        - Velocities are non-dimensionalised and scaled for the simulation before being applied.
    """

    # Raising an error if the user is asking for invalid time
    if self.ndtime2age(ndtime=ndtime) < 0:
        max_time = self.oldest_age / (self.time_dimDmyrs2sec / self.scaling_factor)
        raise ValueError(
            "Input non-dimensionalised time corresponds to negative age (it is in the future)!"
            f" Maximum non-dimensionalised time: {max_time}"
        )

    # cache the reconstruction age
    self.reconstruction_age = self.ndtime2age(ndtime)
    # interpolate velicities onto our grid
    self.interpolated_u = self._interpolate_seeds_u(target_coords)
    return self.interpolated_u

ndtime2age(ndtime)

Converts non-dimensionalised time to age (Myrs before present day).

Parameters:

Name	Type	Description	Default
ndtime	float	The non-dimensionalised time to be converted.	required

Returns:

Name	Type	Description
float		The converted geologic age in millions of years before present day(Ma).

Source code in g-adopt/gadopt/gplates/gplates.py

def ndtime2age(self, ndtime):
    """Converts non-dimensionalised time to age (Myrs before present day).

    Args:
        ndtime (float): The non-dimensionalised time to be converted.

    Returns:
        float: The converted geologic age in millions of years before present day(Ma).
    """

    return self.oldest_age - float(ndtime) * self.time_dimDmyrs2sec / self.scaling_factor

age2ndtime(age)

Converts geologic age (years before present day in Myrs (Ma) to non-dimensionalised time.

    Args:
        age (float): geologic age (before present day in Myrs)

    Returns:

• float: non-dimensionalised time

Source code in g-adopt/gadopt/gplates/gplates.py

    def age2ndtime(self, age):
        """Converts geologic age (years before present day in Myrs (Ma) to non-dimensionalised time.

        Args:
            age (float): geologic age (before present day in Myrs)

        Returns:
-            float: non-dimensionalised time
        """
        return (self.oldest_age - age) * (self.scaling_factor / self.time_dimDmyrs2sec)

MeshConfig(r_outer=2.208, depth_scale=2890.0) dataclass

Non-dimensional mesh geometry shared by every output that converts radius to physical depth.

r_outer is the radial coordinate of Earth's surface in mesh units (default 2.208 for an Earth-like spherical shell with r_inner=1.208). depth_scale is the physical depth (km) per non-dimensional unit (default 2890 — the mantle's depth ratio).

OutputStrategy

Bases: ABC

Map interpolated source arrays at target coords to a scalar field.

Subclasses declare requires — the set of source-dict keys (excluding "xyz") they read from. The connector validates output.requires <= source.provides at construction time.

compute(interpolated, r_target, too_far, mesh) abstractmethod

Return the scalar field at the target points.

Parameters:

Name	Type	Description	Default
interpolated	dict[str, ndarray]	dict of arrays, one per key in self.requires, already kNN-interpolated onto target coords.	required
r_target	ndarray	norms of target coords (one value per target node).	required
too_far	ndarray	boolean mask; True where no source seed is within InterpolationConfig.distance_threshold.	required
mesh	MeshConfig	mesh geometry for radius↔depth conversion.	required

Source code in g-adopt/gadopt/gplates/outputs.py

@abstractmethod
def compute(
    self,
    interpolated: dict[str, np.ndarray],
    r_target: np.ndarray,
    too_far: np.ndarray,
    mesh: MeshConfig,
) -> np.ndarray:
    """Return the scalar field at the target points.

    Args:
        interpolated: dict of arrays, one per key in ``self.requires``,
            already kNN-interpolated onto target coords.
        r_target: norms of target coords (one value per target node).
        too_far: boolean mask; True where no source seed is within
            ``InterpolationConfig.distance_threshold``.
        mesh: mesh geometry for radius↔depth conversion.
    """

Source

Bases: ABC

Time-dependent source of (xyz, properties) on a sphere.

Subclasses expose a provides set listing the property keys (excluding "xyz") that prepare(age) returns. provides is declared here as an abstract property, which a subclass may satisfy either with a plain class constant (the concrete sources do this, keeping class-level access like LithosphereSource.provides) or with an instance property (handy for test doubles and dynamically-configured sources). The ABC implements a per-age cache so that consumers sharing this Source can call prepare in any order and the underlying gtrack state advances at most once per age.

prepare is collective across self.comm. Subclasses implement _compute_sources(age), which runs on rank 0 only; the ABC handles broadcast and caching.

provides abstractmethod property

The set of property keys (excluding 'xyz') that prepare(age) returns.

prepare(age)

Return source arrays at age (collective across self.comm).

Source code in g-adopt/gadopt/gplates/sources.py

def prepare(self, age: float) -> dict[str, np.ndarray]:
    """Return source arrays at ``age`` (collective across self.comm)."""
    # Cache hit: deterministic on age alone, so the decision is identical
    # on every rank — no allreduce needed.
    if (self._cached_age is not None
            and abs(self._cached_age - age) < self.delta_t):
        return self._cached_dict

    self._ensure_loaded()
    sources = self._compute_sources(age) if self._is_root else None
    sources = self.comm.bcast(sources, root=0)

    self._cached_age = age
    self._cached_dict = sources
    # The source cloud just changed, so any cached interpolation geometry
    # (cKDTree indices/weights over the previous cloud) is stale. Clear it
    # in lockstep with _cached_dict so the two caches never disagree on age.
    if self._interp_geometry_cache is not None:
        self._interp_geometry_cache.clear()
    return sources

get_or_build_geometry(key, build_fn)

Return the cached interpolation-geometry bundle for key, building it via build_fn() on a miss.

The cache is rank-local and keyed on (target_coords content hash, id(interpolation_config)); the connector owns the geometry math (build_fn) while the source owns the cache so siblings sharing this source reuse one build. Cleared by prepare whenever the source cloud advances to a new age.

Source code in g-adopt/gadopt/gplates/sources.py

def get_or_build_geometry(self, key, build_fn):
    """Return the cached interpolation-geometry bundle for ``key``, building
    it via ``build_fn()`` on a miss.

    The cache is rank-local and keyed on
    ``(target_coords content hash, id(interpolation_config))``; the
    connector owns the geometry math (``build_fn``) while the source owns
    the cache so siblings sharing this source reuse one build. Cleared by
    ``prepare`` whenever the source cloud advances to a new age.
    """
    if self._interp_geometry_cache is None:
        self._interp_geometry_cache = {}
    bundle = self._interp_geometry_cache.get(key)
    if bundle is None:
        bundle = build_fn()
        self._interp_geometry_cache[key] = bundle
    return bundle

validate_age(age)

Raise if age is outside the plate model's range.

Subclasses with extra state (e.g. forward-only ocean trackers) may override to add further checks. Called by ScalarFieldConnector before prepare().

Source code in g-adopt/gadopt/gplates/sources.py

def validate_age(self, age: float) -> None:
    """Raise if ``age`` is outside the plate model's range.

    Subclasses with extra state (e.g. forward-only ocean trackers) may
    override to add further checks. Called by ScalarFieldConnector before
    prepare().
    """
    if age > self.oldest_age:
        raise ValueError(
            f"Requested age {age:.2f} Ma is older than the plate model's "
            f"oldest age ({self.oldest_age:.2f} Ma)."
        )
    if age < 0:
        raise ValueError(
            f"Requested age {age:.2f} Ma is negative (in the future). "
            f"Ages must be >= 0 (present day)."
        )