Diagnostics

diagnostics

This module provides a class to simplify computing of diagnostics typically encountered in geodynamical simulations. Users instantiate the class by providing relevant parameters and call individual class methods to compute associated diagnostics.

GeodynamicalDiagnostics(z, T=0.0, /, bottom_id=None, top_id=None, *, quad_degree=4)

Typical simulation diagnostics used in geodynamical simulations.

Parameters:

Name	Type	Description	Default
z	Function	Firedrake function for mixed Stokes function space (velocity, pressure)	required
T	Function	Firedrake function for temperature	0.0
bottom_id	int | str | None	Bottom boundary identifier	None
top_id	int | str | None	Top boundary identifier	None
quad_degree	int	Degree of polynomial quadrature approximation	4

Note

All diagnostics are returned as floats.

Methods:

Name	Description
u_rms	Root-mean squared velocity
u_rms_top	Root-mean squared velocity along the top boundary
Nu_top	Nusselt number at the top boundary
Nu_bottom	Nusselt number at the bottom boundary
T_avg	Average temperature in the domain
T_min	Minimum temperature in domain
T_max	Maximum temperature in domain
ux_max	Maximum velocity (first component, optionally over a given boundary)
uv_min	Minimum velocity (vertical component, optionally over a given boundary)

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

def __init__(
    self,
    z: Function,
    T: Function = 0.0,
    /,
    bottom_id: int | str | None = None,
    top_id: int | str | None = None,
    *,
    quad_degree: int = 4,
):
    mesh = extract_unique_domain(z)

    # Allows InternalVariableSolver in stokes_integrators.py
    # (just disp no mixed space) to use same diagnostics.
    # Would be better to separate out into Base classes
    is_mixed_space = isinstance(
        z.function_space().topological, MixedFunctionSpace
    )
    if is_mixed_space:
        self.u, self.p = z.subfunctions[:2]
    else:
        self.u = z

    # vertical component of vel/disp
    self.uv = Function(self.u.function_space().sub(0))

    self.T = T

    self.dx = dx(domain=mesh, degree=quad_degree)
    self.domain_volume = assemble(Constant(1) * self.dx)

    if self.u.function_space().extruded:
        self.ds = CombinedSurfaceMeasure(mesh, quad_degree)
    else:
        self.ds = ds(mesh)

    if bottom_id:
        self.ds_b = self.ds(bottom_id)
        self.bottom_surface = assemble(Constant(1) * self.ds_b)
    if top_id:
        self.ds_t = self.ds(top_id)
        self.top_surface = assemble(Constant(1) * self.ds_t)

    self.n = FacetNormal(mesh)

uv_min(boundary_id=None)

Minimum value of vertical component of velocity/displacement

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

def uv_min(self, boundary_id=None) -> float:
    "Minimum value of vertical component of velocity/displacement"
    self.uv.interpolate(vertical_component(self.u))
    if boundary_id:
        bcu = DirichletBC(self.uv.function_space(), 0, boundary_id)
        uv_data = self.uv.dat.data_ro_with_halos[bcu.nodes]
    else:
        uv_data = self.uv.dat.data_ro[:]

    return self.uv.comm.allreduce(uv_data.min(initial=np.inf), MPI.MIN)