Gplates

gplates

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)

    # 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, :]
    )

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

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