.ipynb .pdf

repository

Binder

Contents

1.c Coordinate Geometry

Quagmire mesh variables are numerically differentiable and the various operators defined in the quagmire functions module also support limited (symbolic) differentiation

This is to allow the construction of more complicated operators that are equivalent on the flat plane and the surface of the spehre.

This functionality is supported through the quagmire.function.geometry package.

import numpy as np
from quagmire.function import display

from quagmire import QuagMesh 
from quagmire import tools as meshtools
from quagmire import function as fn
from quagmire.function import display
from mpi4py import MPI

# import lavavu
import stripy
comm = MPI.COMM_WORLD

import matplotlib.pyplot as plt
from matplotlib import cm
%matplotlib inline

Naked coordinate systems

Quagmire supports flat 2D meshes and logically 2D meshes on the surface of a sphere. When a mesh is defined, it attaches a geometry to itself which is the starting point for computing any differential operators that include mesh variables. It is also possible to access the coordinate system definitions directly to see how the various operators are defined. Let us first access a simple, 2D, Cartesian Coordinate Geometry

CartesianCoordinates = fn.geometry.CartesianCoordinates2D()
SphericalCoordinates = fn.geometry.SphericalSurfaceLonLat2D()

# These are the coordinate directions is symbolic form. They are also functions that 
# mask the appropriate member of a coordinate tuple

x = CartesianCoordinates.xi0
y = CartesianCoordinates.xi1

print(x.evaluate(1.0,2.0), y.evaluate(1.0,2.0))

# Another example:

points = np.zeros((20,2))
points[:,0] = np.linspace(0.0,2.0*np.pi,20)

print(x.evaluate(points))

S = fn.math.sin(x)

# These are equivalent

print(S.evaluate(points) - np.sin(points[:,0]))

1.0 2.0
[0.         0.33069396 0.66138793 0.99208189 1.32277585 1.65346982
 1.98416378 2.31485774 2.64555171 2.97624567 3.30693964 3.6376336
 3.96832756 4.29902153 4.62971549 4.96040945 5.29110342 5.62179738
 5.95249134 6.28318531]
[0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]

xi0 = SphericalCoordinates.xi0
xi1 = SphericalCoordinates.xi1

Vector operators

Partial differential equations are often the balancing of gradients of different quantities. The physical gradients may be very different when expressed in different coordinate systems and therefore we have to be careful to develop coordinate-independent expression for gradients that appear in PDEs. We do this by constructing operators such as div, grad, curl and the Laplacian that are expressed independently of the geometry but understand the underlying coordinate system.

When we calculate a derivative, they should be symbolically equivalent but the gradients are not.

for CoordinateSystem in [CartesianCoordinates, SphericalCoordinates]:
    
    xi0 = CoordinateSystem.xi0
    xi1 = CoordinateSystem.xi1
    
    ## Derivatives:
    
    A = xi0 * fn.math.sin(xi0**2) +  xi1 * fn.math.sin(xi1**2) +  xi0 * xi1
    
    print("Derivatives")
    ddx0 = A.derivative(0)
    ddx1 = A.derivative(1)
    display(ddx0)
    display(ddx1)

    gradA = CoordinateSystem.grad(A)
    print("Grad 0")
    gradA[0].display()
    print("Grad 1")
    gradA[1].display()
    
    print("div.grad")
    CoordinateSystem.div(gradA, expand=True).display()
    
    print("Laplacian")
    CoordinateSystem.laplacian(A, 1, expand=True).display()  # Note this is written for variable coefficient problems

    CoordinateSystem.laplacian(A, xi0*xi1, expand=False).display()  ## Neater for display purposes only

Derivatives

\[\displaystyle \sin\left(x^{2}\right) \, + \, x \;2 \;x \;\cos\left(x^{2}\right) \, + \, y\]

\[\displaystyle \sin\left(y^{2}\right) \, + \, y \;2 \;y \;\cos\left(y^{2}\right) \, + \, x\]

Grad 0

\[\displaystyle \sin\left(x^{2}\right) \, + \, x \;2 \;x \;\cos\left(x^{2}\right) \, + \, y\]

Grad 1

\[\displaystyle \sin\left(y^{2}\right) \, + \, y \;2 \;y \;\cos\left(y^{2}\right) \, + \, x\]

div.grad

\[\displaystyle 2 \;x \;\cos\left(x^{2}\right) \, + \, 2 \;x \;\cos\left(x^{2}\right) \, + \, x \;(2 \;\cos\left(x^{2}\right) \, + \, 2 \;x \;-2 \;x \;\sin\left(x^{2}\right)) \, + \, 2 \;y \;\cos\left(y^{2}\right) \, + \, 2 \;y \;\cos\left(y^{2}\right) \, + \, y \;(2 \;\cos\left(y^{2}\right) \, + \, 2 \;y \;-2 \;y \;\sin\left(y^{2}\right))\]

Laplacian

\[\displaystyle 2 \;x \;\cos\left(x^{2}\right) \, + \, 2 \;x \;\cos\left(x^{2}\right) \, + \, x \;(2 \;\cos\left(x^{2}\right) \, + \, 2 \;x \;-2 \;x \;\sin\left(x^{2}\right)) \, + \, 2 \;y \;\cos\left(y^{2}\right) \, + \, 2 \;y \;\cos\left(y^{2}\right) \, + \, y \;(2 \;\cos\left(y^{2}\right) \, + \, 2 \;y \;-2 \;y \;\sin\left(y^{2}\right))\]

\[\displaystyle \nabla \left( x \;y \cdot \nabla \left( x \;\sin\left(x^{2}\right) \, + \, y \;\sin\left(y^{2}\right) \, + \, x \;y \right)\right)\]

Derivatives

\[\displaystyle \sin\left(\theta^{2}\right) \, + \, \theta \;2 \;\theta \;\cos\left(\theta^{2}\right) \, + \, \phi' \]

\[\displaystyle \sin\left(\phi' ^{2}\right) \, + \, \phi' \;2 \;\phi' \;\cos\left(\phi' ^{2}\right) \, + \, \theta\]

Grad 0

\[\displaystyle \frac{ \sin\left(\theta^{2}\right) \, + \, \theta \;2 \;\theta \;\cos\left(\theta^{2}\right) \, + \, \phi' }{ \cos\left(\phi' \right) }\]

Grad 1

\[\displaystyle \sin\left(\phi' ^{2}\right) \, + \, \phi' \;2 \;\phi' \;\cos\left(\phi' ^{2}\right) \, + \, \theta\]

div.grad

\[\displaystyle \frac{ \frac{ 2 \;\theta \;\cos\left(\theta^{2}\right) \, + \, 2 \;\theta \;\cos\left(\theta^{2}\right) \, + \, \theta \;(2 \;\cos\left(\theta^{2}\right) \, + \, 2 \;\theta \;-2 \;\theta \;\sin\left(\theta^{2}\right)) }{ \cos\left(\phi' \right) } }{ \cos\left(\theta\right) } \, + \, \frac{ (2 \;\phi' \;\cos\left(\phi' ^{2}\right) \, + \, 2 \;\phi' \;\cos\left(\phi' ^{2}\right) \, + \, \phi' \;(2 \;\cos\left(\phi' ^{2}\right) \, + \, 2 \;\phi' \;-2 \;\phi' \;\sin\left(\phi' ^{2}\right))) \;\cos\left(\theta\right) }{ \cos\left(\theta\right) }\]

Laplacian

\[\displaystyle \frac{ \frac{ 2 \;\theta \;\cos\left(\theta^{2}\right) \, + \, 2 \;\theta \;\cos\left(\theta^{2}\right) \, + \, \theta \;(2 \;\cos\left(\theta^{2}\right) \, + \, 2 \;\theta \;-2 \;\theta \;\sin\left(\theta^{2}\right)) }{ \cos\left(\phi' \right) } }{ \cos\left(\theta\right) } \, + \, \frac{ (2 \;\phi' \;\cos\left(\phi' ^{2}\right) \, + \, 2 \;\phi' \;\cos\left(\phi' ^{2}\right) \, + \, \phi' \;(2 \;\cos\left(\phi' ^{2}\right) \, + \, 2 \;\phi' \;-2 \;\phi' \;\sin\left(\phi' ^{2}\right))) \;\cos\left(\theta\right) }{ \cos\left(\theta\right) }\]

\[\displaystyle \nabla \left( \theta \;\phi' \cdot \nabla \left( \theta \;\sin\left(\theta^{2}\right) \, + \, \phi' \;\sin\left(\phi' ^{2}\right) \, + \, \theta \;\phi' \right)\right)\]

Symbols

These are placeholders for developing a function that you can then substitute later. This is useful as an aid to debugging functions or for providing an equation template

k0 = fn.symbol(name="k_0", lname="k_0")
k1 = fn.symbol(name="k_1", lname="k_1")
display((k0, k1))

H = x * fn.math.sin(k0*x + k0/(k1+x))
display(H)

try:
    H.evaluate((0.1,1.0))
except:
    print("Cannot evaluate H ... try making a subsitution")

k0.substitute(fn.parameter(4))
k1.substitute(fn.parameter(0.0001))

display(k1)
display(H)

print("H evaluated at (0.1,1.0): ",H.evaluate((0.1,1.0)))


k0.substitute(x)
k1.substitute(fn.parameter(0.01))

display(k1)
display(H)

\[\displaystyle k_0\]

\[\displaystyle k_1\]

\[\displaystyle x \;\sin\left(k_0 \;x \, + \, \frac{ k_0 }{ k_1 \, + \, x }\right)\]

Cannot evaluate H ... try making a subsitution

\[\displaystyle \left\{ k_1 \leftarrow 0.0001\right\}\]

\[\displaystyle x \;\sin\left(\left\{ k_0 \leftarrow 4\right\} \;x \, + \, \frac{ \left\{ k_0 \leftarrow 4\right\} }{ \left\{ k_1 \leftarrow 0.0001\right\} \, + \, x }\right)\]

H evaluated at (0.1,1.0):  [0.04623685]

\[\displaystyle \left\{ k_1 \leftarrow 0.01\right\}\]

\[\displaystyle x \;\sin\left(\left\{ k_0 \leftarrow x\right\} \;x \, + \, \frac{ \left\{ k_0 \leftarrow x\right\} }{ \left\{ k_1 \leftarrow 0.01\right\} \, + \, x }\right)\]

You can “finalise” a symbol which means that it is no longer possible to make a substitution and it will generally behave like a normal quagmire function (including the fact that it no longer displays the substitution). The machinery for differentiation and evaluation is still a little more complicated than that of a regular function, so there may be some performance issues.

display(k0)

\[\displaystyle \left\{ k_0 \leftarrow x\right\}\]

k0.finalise()
display(k0)

\[\displaystyle x\]

The same sort of thing but with some mesh variables

Mesh variables have an associated geometry that they obtain from their mesh and they have some numerical differentiation routines that are used when a derivative or gradient needs to be evaluated.

## Add a mesh

st0 = stripy.spherical_meshes.icosahedral_mesh(refinement_levels=6, include_face_points=True)
dm = meshtools.create_spherical_DMPlex(np.degrees(st0.lons), np.degrees(st0.lats), st0.simplices)
mesh = QuagMesh(dm, downhill_neighbours=1, permute=True)

lon = mesh.geometry.xi0
lat = mesh.geometry.xi1

Underlying Mesh type: sTriMesh
0 - Delaunay triangulation 0.21779666599999814s
0 - Calculate node weights and area 0.03310762499999953s
0 - Find boundaries 0.014876750000002659s
0 - cKDTree 0.01549533399999703s
0 - Construct neighbour cloud arrays 0.7254632500000007s, (0.5542094169999991s + 0.1712229160000014s)
0 - Construct rbf weights 0.045508500000000396s

F = mesh.add_variable("F", lname="\cal{F}")

l0 = fn.parameter(5.0)
l1 = fn.parameter(2.0)

G = fn.math.cos(l0 * fn.math.radians(lon)) * fn.math.sin(l1 * fn.math.radians(lat))
F.data = G.evaluate(mesh)

F.display()
G.display()

\[\displaystyle \cal{F}\left(\theta,\phi' \right)\]

\[\displaystyle \cos\left(5 \;\frac{\pi}{180}\theta\right) \;\sin\left(2 \;\frac{\pi}{180}\phi' \right)\]

display(G.fn_gradient(0))
display(G.derivative(0))
mesh.geometry.grad(G)[0]

\[\displaystyle \frac{ -0.0873 \;\sin\left(5 \;\frac{\pi}{180}\theta\right) \;\sin\left(2 \;\frac{\pi}{180}\phi' \right) }{ \cos\left(\phi' \right) }\]

\[\displaystyle -0.0873 \;\sin\left(5 \;\frac{\pi}{180}\theta\right) \;\sin\left(2 \;\frac{\pi}{180}\phi' \right)\]

\[\displaystyle \frac{ -0.0873 \;\sin\left(5 \;\frac{\pi}{180}\theta\right) \;\sin\left(2 \;\frac{\pi}{180}\phi' \right) }{ \cos\left(\phi' \right) }\]

## Something with symbols we can substitute

m0 = fn.symbol(name="m_0", lname="m_0")
m1 = fn.symbol(name="m_1", lname="m_1")

H = fn.math.cos(m0 * fn.math.radians(lon)) * fn.math.sin(m1 * fn.math.radians(lat))
H.display()

m0.substitute(F**2)
m1.substitute(F**4)

H.display()

m0.substitute(F+1)
m1.substitute(F-1)

H.display()

# G.fn_gradient(1)

\[\displaystyle \cos\left(m_0 \;\frac{\pi}{180}\theta\right) \;\sin\left(m_1 \;\frac{\pi}{180}\phi' \right)\]

\[\displaystyle \cos\left(\left\{ m_0 \leftarrow \cal{F}\left(\theta,\phi' \right)^{2}\right\} \;\frac{\pi}{180}\theta\right) \;\sin\left(\left\{ m_1 \leftarrow \cal{F}\left(\theta,\phi' \right)^{4}\right\} \;\frac{\pi}{180}\phi' \right)\]

\[\displaystyle \cos\left(\left\{ m_0 \leftarrow \cal{F}\left(\theta,\phi' \right) \, + \, 1\right\} \;\frac{\pi}{180}\theta\right) \;\sin\left(\left\{ m_1 \leftarrow \cal{F}\left(\theta,\phi' \right) - 1\right\} \;\frac{\pi}{180}\phi' \right)\]

H.evaluate(mesh)

array([-1.        , -0.09259669,  0.03942575, ...,  0.91069208,
        0.03873612,  0.03113295])

H.derivative(1)

\[\displaystyle -\left. \nabla(\cal{F}\left(\theta,\phi' \right))\right|_1 \;\frac{\pi}{180}\theta \;\sin\left(\left\{ m_0 \leftarrow \cal{F}\left(\theta,\phi' \right) \, + \, 1\right\} \;\frac{\pi}{180}\theta\right) \;\sin\left(\left\{ m_1 \leftarrow \cal{F}\left(\theta,\phi' \right) - 1\right\} \;\frac{\pi}{180}\phi' \right) \, + \, \cos\left(\left\{ m_0 \leftarrow \cal{F}\left(\theta,\phi' \right) \, + \, 1\right\} \;\frac{\pi}{180}\theta\right) \;(\left. \nabla(\cal{F}\left(\theta,\phi' \right))\right|_1 \;\frac{\pi}{180}\phi' \, + \, \left\{ m_1 \leftarrow \cal{F}\left(\theta,\phi' \right) - 1\right\} \;0.0175) \;\cos\left(\left\{ m_1 \leftarrow \cal{F}\left(\theta,\phi' \right) - 1\right\} \;\frac{\pi}{180}\phi' \right)\]

m0.finalise()
m1.finalise()

H.derivative(1)

\[\displaystyle -\left. \nabla(\cal{F}\left(\theta,\phi' \right))\right|_1 \;\frac{\pi}{180}\theta \;\sin\left(\left( \cal{F}\left(\theta,\phi' \right) \, + \, 1 \right) \;\frac{\pi}{180}\theta\right) \;\sin\left(\left( \cal{F}\left(\theta,\phi' \right) - 1 \right) \;\frac{\pi}{180}\phi' \right) \, + \, \cos\left(\left( \cal{F}\left(\theta,\phi' \right) \, + \, 1 \right) \;\frac{\pi}{180}\theta\right) \;(\left. \nabla(\cal{F}\left(\theta,\phi' \right))\right|_1 \;\frac{\pi}{180}\phi' \, + \, (\left( \cal{F}\left(\theta,\phi' \right) - 1 \right)) \;0.0175) \;\cos\left(\left( \cal{F}\left(\theta,\phi' \right) - 1 \right) \;\frac{\pi}{180}\phi' \right)\]

previous

1.b Functions

next

1.d Quagmire cloud