Gplates

gplates

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

Bases: GPlatesFunctionalityMixin, Function

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

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).

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

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", **kwargs):
    # Initialize as a Firedrake Function
    super().__init__(function_space, *args, **kwargs)

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

pyGplatesConnector(rotation_filenames, topology_filenames, oldest_age, delta_t=1.0, scaling_factor=1.0, nseeds=100000.0, nneighbours=4)

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 neighboring points when interpolating velocity for each grid point. Default is 4.	4

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):
    """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 neighboring points when interpolating velocity for each grid point. Default is 4.

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

    # 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

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:
        raise ValueError(
            ("Input non-dimensionalised time corresponds to negative age (it is in the future)!"
             f"maximum non-dimensionalised time:"
             f" {self.oldest_age/(pyGplatesConnector.time_dimDmyrs2sec/self.scaling_factor)}")
        )

    # 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) * pyGplatesConnector.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 / pyGplatesConnector.time_dimDmyrs2sec)