1.b Functions¶

Like Underworld, quagmire provides a function interface that can be used to compose data and operations on the fly in order to construct model equations independent of whatever approach is used for solution.

Noteably, these are lazy function that are only evaluated when needed. More importantly, when evaluated, they also use the current state of any variables in their definition and so can be placed within timestepping loops where they will always use the information of the current timestep.

There are three kinds of lazy functions available in quagmire:

MeshVariable data containers that hold information on the mesh and can return that information at any point by interpolation (or, less reliably by extrapolation) and can also provide the gradient of the data using a cubic spline interpolant (see the documentation for stripy for details).
parameter is a floating point value that can be used for coefficients in an equation. The value of the parameter can be updated.
virtual variables which are operations on MeshVariables and parameters and contain no data record.

These lazy functions are members of the LazyEvaluation class that defines the following methods / data

evaluate(mesh | X, Y) computes a snapshot of the value(s) at the mesh points of mesh or at the points given by X and Y
fn_gradient(dir) is a lazy function that can be evaluated to obtain the gradient in the direction dir=(0|1)
description is a string describing the result returned by evaluate. This is helpful because the function may be a cascade of operations. It very much helps to provide short, useful names for your mesh variables to get back reasonable descriptions.

Note: at present no error checking is done for consistency between the mesh provided to evaluate and the one used to store the original data. This is very bad on our part !

from quagmire.tools import meshtools
from quagmire import QuagMesh
from quagmire.mesh import MeshVariable
from quagmire import function as fn
import numpy as np

Working mesh¶

First we create a basic mesh so that we can define mesh variables and obbtain gradients etc.

minX, maxX = -5.0, 5.0
minY, maxY = -5.0, 5.0,
dx, dy = 0.02, 0.02

from stripy.cartesian_meshes import elliptical_base_mesh_points
epointsx, epointsy, ebmask = elliptical_base_mesh_points(10.0, 7.5, 0.25, remove_artifacts=True)
dm = meshtools.create_DMPlex_from_points(epointsx, epointsy, bmask=ebmask, refinement_levels=1)

mesh = QuagMesh(dm, downhill_neighbours=1)

Basic usage¶

The functions can be demonstrated on the most basic example the parameter which is constant everywhere in the mesh. In fact, these operations work without any reference to the mesh since they are the same at all locations and their gradient is zero everywhere.

A = fn.parameter(10.0)
B = fn.parameter(100.0)

print("Exp({}) = {}".format(A.value, fn.math.exp(A).evaluate(0.0,0.0)))
print("Exp({}) = {}".format(B.value, fn.math.exp(B).evaluate(0.0,0.0)))

# A is a proper lazy variable too so this is required to work

print("Exp({}) = {}".format(A.evaluate(0.0,0.0), fn.math.exp(A).evaluate(0.0,0.0)))

# And this is how to update A

A.value = 100.0
print("Exp({}) = {}".format(A.evaluate(0.0,0.0), fn.math.exp(A).evaluate(0.0,0.0)))

# This works too ... and note the floating point conversion
A(101)
print("Exp({}) = {}".format(A.evaluate(0.0,0.0), fn.math.exp(A).evaluate(0.0,0.0)))

# More complicated examples
print((fn.math.sin(A)**2.0 + fn.math.cos(A)**2.0).evaluate(0.0,0.0))

Descriptions¶

The lazy function carries a description string that tells you approximately what will happen when the function is evaluated. There is also a .math() method that gives a \(\LaTeX\) string which displays nicely in a jupyter notebook. The quagmire.function.display function tries to make this notebook display easy for a variety of cases.

Examples:

print(A.description)
print((fn.math.sin(A)+fn.math.cos(B)).description)
print((fn.math.sin(A)**2.0 + fn.math.cos(A)**2.0).description)

# the description is printed by default if you call print on the function 

print((fn.math.sin(A)**2.0 + fn.math.cos(A)**2.0))

# the latex version is accessed like this:

fn.display(fn.math.sin(A)+fn.math.cos(B))
fn.display(fn.math.sin(A)**2.0 + fn.math.cos(A)**2.0)

Predefined Mesh Functions¶

There are some predefined mesh functions that we can use in building more complicated functions that depend on the mesh geometry. The details are described in the Ex1c-QuagmireCoordinateGeometry.py notebook.

The mesh.coordinates.xi0/1 functions are symbols representing the coordinate directions and can be evaluated to extract the relevant mesh directions. That means we can do this:

X = mesh.coordinates.xi0
Y = mesh.coordinates.xi1

display(fn.math.sin(X))
print(X.evaluate(mesh))

S = fn.math.sin(X+Y) - fn.math.cos(X*Y)
display(S)
print(S.evaluate(mesh))

X.evaluate(mesh)

Mesh Variables as functions¶

Mesh variables (MeshVariables) are also members of the LazyEvaluation class. They can be evaluated exactly as the parameters can, but it is also possible to obtain numerical derivatives. Of course, they also have other properties beyond those of simple functions (see the MeshVariables examples in the previous (Ex1a-QuagmireMeshVariables.py) notebook for details).

Let us first define a mesh variable …

height = mesh.add_variable(name="h(X,Y)", lname="h")
height.data = np.ones(mesh.npoints)
print(height)
display(height)
height.math()

We might introduce a universal scaling for the height variable. This could be useful if, say, the offset is something that we might want to change programmatically within a timestepping loop.

h_scale = fn.parameter(2.0)
h_offset = fn.parameter(1.0)

scaled_height = h_scale * height + h_offset

display(scaled_height)

print(height.evaluate(mesh))
print(scaled_height.evaluate(mesh))

h_offset.value = 10.0
print(scaled_height.evaluate(mesh))

height * h_scale

We might wish to define a rainfall parameter that is a function of height that can be passed in to some existing code. The use of functions is perfect for this.

rainfall_exponent = fn.parameter(2.2)
rainfall = scaled_height**rainfall_exponent
fn.display(rainfall)
print(rainfall.evaluate(mesh))

The rainfall definition is live to any changes in the height but we can also adjust the rainfall parameters on the fly. This allows us to define operators with coefficients that can be supplied as mutable parameters.

height.data = np.sin(mesh.coords[:,0])
print("Height:", height.data)
print("Rainfall Fn evaluation:",rainfall.evaluate(mesh))
print("Rainfall Direct       :",(height.data*2.0+10.0)**2.2)

# change definition of rainfall coefficient but not functional form

rainfall_exponent.value = 1.0
print("Rainfall Fn evaluation:",rainfall.evaluate(mesh))
print("Rainfall Direct       :",(height.data*2.0+10.0)**2.2)

While functions are most useful because they are not computed once and for all, it is also possible to compute their values and assign to a variable. Just be aware that, at this point, numpy has a greater richness of operators than quagmire.function. We can rewrite the above assignment to the height variable using the coord function that extracts values of the x or y ( 0 or 1 ) coordinate direction from the locations given as arguments to evaluate. Note that the rainfall function always used the updated height values.

height.data = fn.math.cos(fn.misc.coord(0)).evaluate(mesh)
print("Height:  ", height.data)
print("Height = ", fn.math.cos(fn.misc.coord(0)).description)

rainfall.evaluate(mesh)

Operator overloading for +, - , *, **, /¶

We define addition / subtraction (negation), multiplication, division, and raising to arbitrary power for mesh variables and parameters and the meaning is carried over from numpy - i.e. generally these are element-by-element operations on the underlying data vector and require the data structures to have compatible sizes.

dhdx, dhdy = mesh.geometry.grad(height)
slope = fn.math.sqrt((dhdx**2 + dhdy**2))

native_slope = mesh.geometry.slope(height)  # This actually just returns the height.slope function

a = fn.parameter(1.3)
k = slope**a
k2 = native_slope ** a

display(slope)
display(native_slope)

display(k)
display(k2)

# Numerical equivalence

print(k.evaluate(mesh))
print(k2.evaluate(mesh))

Gradients¶

Variables associated with a mesh also have the capacity to form spatial derivatives anywhere. This is provided by the stripy gradient routines in the case of triangulations. The gradient can be formed from any lazy function by evaluating it at the mesh points and then obtaining values of derivatives anywhere via stripy. In the case of the spatially invariant parameter objects, the derivatives are identically zero.

gradx = height.fn_gradient(0)
display(gradx)
grady = scaled_height.fn_gradient(1)
display(grady)

Example: It is a common operation to compute a power law of the magnitude of the local slope. In Cartesian geometry, the slope is defined this way

\[ \left| \nabla h \right| \equiv \sqrt{ \frac{\partial h}{\partial x}^2 + \frac{\partial h}{\partial y}^2 } \]

On the sphere, this expression is a little more complicated and this is why the expression is written in terms of components of the gradient operator in the display below

\[ k = \left| \nabla h \right|^a \]

Mesh variables have an optimised numerical shortcut for calculating slopes.

NOTE: The gradient operators are dependent upon the coordinate system itself. This is ususally inherited from a mesh but it can be defined independently of the mesh.

dhdx, dhdy = mesh.geometry.grad(height)
slope = fn.math.sqrt((dhdx**2 + dhdy**2))

native_slope = mesh.geometry.slope(height)  # This actually just returns the height.slope function

a = fn.parameter(1.3)
k = slope**a
k2 = native_slope ** a

display(slope)
display(native_slope)

display(k)
display(k2)

# Numerical equivalence

print(k.evaluate(mesh))
print(k2.evaluate(mesh))

Quagmire.tex