# Meshing Stuff

## Topologies

Topologies encode the connectivity of the elements, spatial domain interval and MPI communication.

### Types

ClimateMachine.Mesh.Topologies.BoxElementTopologyType
BoxElementTopology{dim, T, nb} <: AbstractTopology{dim,T,nb}

The local topology of a larger MPI-distributed topology, represented by dim-dimensional box elements, with nb boundary conditions.

This contains the necessary information for the connectivity elements of the elements on the local process, along with "ghost" elements from neighbouring processes.

Fields

• mpicomm

MPI communicator for communicating with neighbouring processes.

• elems

Range of element indices

• realelems

Range of real (aka nonghost) element indices

• ghostelems

Range of ghost element indices

• ghostfaces

Ghost element to face is received; ghostfaces[f,ge] == true if face f of ghost element ge is received.

• sendelems

Array of send element indices

• sendfaces

Send element to face is sent; sendfaces[f,se] == true if face f of send element se is sent.

• interiorelems

Array of real elements that do not have a ghost element as a neighbor.

• exteriorelems

Array of real elements that have at least on ghost element as a neighbor.

Note that this is different from sendelems because sendelems duplicates elements that need to be sent to multiple neighboring processes.

• elemtocoord

Element to vertex coordinates; elemtocoord[d,i,e] is the dth coordinate of corner i of element e

Note

currently coordinates always are of size 3 for (x1, x2, x3)

• elemtoelem

Element to neighboring element; elemtoelem[f,e] is the number of the element neighboring element e across face f. If it is a boundary face, then it is the boundary element index.

• elemtoface

Element to neighboring element face; elemtoface[f,e] is the face number of the element neighboring element e across face f. If there is no neighboring element then elemtoface[f,e] == f."

• elemtoordr

element to neighboring element order; elemtoordr[f,e] is the ordering number of the element neighboring element e across face f. If there is no neighboring element then elemtoordr[f,e] == 1.

• elemtobndy

Element to boundary number; elemtobndy[f,e] is the boundary number of face f of element e. If there is a neighboring element then elemtobndy[f,e] == 0.

• nabrtorank

List of the MPI ranks for the neighboring processes

• nabrtorecv

Range in ghost elements to receive for each neighbor

• nabrtosend

Range in sendelems to send for each neighbor

• origsendorder

original order in partitioning

• bndytoelem

Tuple of boundary to element. bndytoelem[b][i] is the element which faces the ith boundary element of boundary b.

• bndytoface

Tuple of boundary to element face. bndytoface[b][i] is the face number of the element which faces the ith boundary element of boundary b.

• elemtouvert

Element to locally unique vertex number; elemtouvert[v,e] is the vth vertex of element e

• vtconn

Vertex connectivity information for direct stiffness summation; vtconn[vtconnoff[i]:vtconnoff[i+1]-1] == vtconn[lvt1, lelem1, lvt2, lelem2, ....] for each rank-local unique vertex number, where lvt1 is the element-local vertex number of element lelem1, etc. Vertices, not shared by multiple elements, are ignored.

• vtconnoff

Vertex offset for vertex connectivity information

• fcconn

Face connectivity information for direct stiffness summation; fcconn[ufc,:] = [lfc1, lelem1, lfc2, lelem2, ordr] where ufc is the rank-local unique face number, and lfc1 is the element-local face number of element lelem1, etc., and ordr is the relative orientation. Faces, not shared by multiple elements are ignored.

• edgconn

Edge connectivity information for direct stiffness summation; edgconn[edgconnoff[i]:edgconnoff[i+1]-1] == edgconn[ledg1, orient, lelem1, ledg2, orient, lelem2, ...] for each rank-local unique edge number, where ledg1 is the element-local edge number, orient1 is orientation (forward/reverse) of dof along edge. Edges, not shared by multiple elements, are ignored.

• edgconnoff

Edge offset for edge connectivity information

source
ClimateMachine.Mesh.Topologies.SingleExponentialStretchingType
SingleExponentialStretching(A)

Apply single-exponential stretching: A > 0 will increase the density of points at the lower boundary, A < 0 will increase the density at the upper boundary.

Reference

• "Handbook of Grid Generation" J. F. Thompson, B. K. Soni, N. P. Weatherill (Editors) RCR Press 1999, §3.6.1 Single-Exponential Function
source

### Constructors

ClimateMachine.Mesh.Topologies.BrickTopologyMethod
BrickTopology{dim, T}(mpicomm, elemrange; boundary, periodicity)

Generate a brick mesh topology with coordinates given by the tuple elemrange and the periodic dimensions given by the periodicity tuple.

The elements of the brick are partitioned equally across the MPI ranks based on a space-filling curve.

By default boundary faces will be marked with a one and other faces with a zero. Specific boundary numbers can also be passed for each face of the brick in boundary. This will mark the nonperiodic brick faces with the given boundary number.

Examples

We can build a 3 by 2 element two-dimensional mesh that is periodic in the $x2$-direction with


using ClimateMachine.Topologies
using MPI
MPI.Init()
topology = BrickTopology(MPI.COMM_SELF, (2:5,4:6);
periodicity=(false,true),
boundary=((1,2),(3,4)))

This returns the mesh structure for

         x2

^
|
6-  +-----+-----+-----+
|  |     |     |     |
|  |  3  |  4  |  5  |
|  |     |     |     |
5-  +-----+-----+-----+
|  |     |     |     |
|  |  1  |  2  |  6  |
|  |     |     |     |
4-  +-----+-----+-----+
|
+--|-----|-----|-----|--> x1
2     3     4     5

For example, the (dimension by number of corners by number of elements) array elemtocoord gives the coordinates of the corners of each element.

julia> topology.elemtocoord
2×4×6 Array{Int64,3}:
[:, :, 1] =
2  3  2  3
4  4  5  5

[:, :, 2] =
3  4  3  4
4  4  5  5

[:, :, 3] =
2  3  2  3
5  5  6  6

[:, :, 4] =
3  4  3  4
5  5  6  6

[:, :, 5] =
4  5  4  5
5  5  6  6

[:, :, 6] =
4  5  4  5
4  4  5  5

Note that the corners are listed in Cartesian order.

The (number of faces by number of elements) array elemtobndy gives the boundary number for each face of each element. A zero will be given for connected faces.

julia> topology.elemtobndy
4×6 Array{Int64,2}:
1  0  1  0  0  0
0  0  0  0  2  2
0  0  0  0  0  0
0  0  0  0  0  0

Note that the faces are listed in Cartesian order.

source
ClimateMachine.Mesh.Topologies.StackedBrickTopologyMethod
StackedBrickTopology{dim, T}(mpicomm, elemrange; boundary, periodicity)

Generate a stacked brick mesh topology with coordinates given by the tuple elemrange and the periodic dimensions given by the periodicity tuple.

The elements are stacked such that the elements associated with range elemrange[dim] are contiguous in the element ordering.

The elements of the brick are partitioned equally across the MPI ranks based on a space-filling curve. Further, stacks are not split at MPI boundaries.

By default boundary faces will be marked with a one and other faces with a zero. Specific boundary numbers can also be passed for each face of the brick in boundary. This will mark the nonperiodic brick faces with the given boundary number.

Examples

We can build a 3 by 2 element two-dimensional mesh that is periodic in the $x2$-direction with


using ClimateMachine.Topologies
using MPI
MPI.Init()
topology = StackedBrickTopology(MPI.COMM_SELF, (2:5,4:6);
periodicity=(false,true),
boundary=((1,2),(3,4)))

This returns the mesh structure stacked in the $x2$-direction for

         x2

^
|
6-  +-----+-----+-----+
|  |     |     |     |
|  |  2  |  4  |  6  |
|  |     |     |     |
5-  +-----+-----+-----+
|  |     |     |     |
|  |  1  |  3  |  5  |
|  |     |     |     |
4-  +-----+-----+-----+
|
+--|-----|-----|-----|--> x1
2     3     4     5

For example, the (dimension by number of corners by number of elements) array elemtocoord gives the coordinates of the corners of each element.

julia> topology.elemtocoord
2×4×6 Array{Int64,3}:
[:, :, 1] =
2  3  2  3
4  4  5  5

[:, :, 2] =
2  3  2  3
5  5  6  6

[:, :, 3] =
3  4  3  4
4  4  5  5

[:, :, 4] =
3  4  3  4
5  5  6  6

[:, :, 5] =
4  5  4  5
4  4  5  5

[:, :, 6] =
4  5  4  5
5  5  6  6

Note that the corners are listed in Cartesian order.

The (number of faces by number of elements) array elemtobndy gives the boundary number for each face of each element. A zero will be given for connected faces.

julia> topology.elemtobndy
4×6 Array{Int64,2}:
1  0  1  0  0  0
0  0  0  0  2  2
0  0  0  0  0  0
0  0  0  0  0  0

Note that the faces are listed in Cartesian order.

source
ClimateMachine.Mesh.Topologies.CubedShellTopologyMethod
CubedShellTopology(mpicomm, Nelem, T) <: AbstractTopology{dim,T,nb}

Generate a cubed shell mesh with the number of elements along each dimension of the cubes being Nelem. This topology actually creates a cube mesh, and the warping should be done after the grid is created using the cubed_sphere_warp function. The coordinates of the points will be of type T.

The elements of the shell are partitioned equally across the MPI ranks based on a space-filling curve.

Note that this topology is logically 2-D but embedded in a 3-D space

Examples

We can build a cubed shell mesh with 10*10 elements on each cube face, i.e., the total number of elements is 10 * 10 * 6 = 600, with

using ClimateMachine.Topologies
using MPI
MPI.Init()
topology = CubedShellTopology(MPI.COMM_SELF, 10, Float64)

# Typically the warping would be done after the grid is created, but the cell
# corners could be warped with...

x1, x2, x3 = ntuple(j->topology.elemtocoord[j, :, :], 3)
for n = 1:length(x1)
x1[n], x2[n], x3[n] = Topologies.cubed_sphere_warp(EquiangularCubedSphere(), x1[n], x2[n], x3[n])
end

# in case a unitary equiangular cubed sphere is desired, or

x1, x2, x3 = ntuple(j->topology.elemtocoord[j, :, :], 3)
for n = 1:length(x1)
x1[n], x2[n], x3[n] = Topologies.cubed_sphere_warp(EquidistantCubedSphere(), x1[n], x2[n], x3[n], 10)
end

# in case an equidistant cubed sphere of radius 10 is desired.
source
ClimateMachine.Mesh.Topologies.StackedCubedSphereTopologyMethod

StackedCubedSphereTopology(mpicomm, Nhorz, Rrange; boundary=(1,1)) <: AbstractTopology{3}

Generate a stacked cubed sphere topology with Nhorz by Nhorz cells for each horizontal face and Rrange is the radius edges of the stacked elements. This topology actual creates a cube mesh, and the warping should be done after the grid is created using the cubed_sphere_warp function. The coordinates of the points will be of type eltype(Rrange). The inner boundary condition type is boundary[1] and the outer boundary condition type is boundary[2].

The elements are stacked such that the vertical elements are contiguous in the element ordering.

The elements of the brick are partitioned equally across the MPI ranks based on a space-filling curve. Further, stacks are not split at MPI boundaries.

Examples

We can build a cubed sphere mesh with 10 x 10 x 5 elements on each cube face, i.e., the total number of elements is 10 * 10 * 5 * 6 = 3000, with

using ClimateMachine.Topologies
using MPI
MPI.Init()
Nhorz = 10
Nstack = 5
Rrange = Float64.(accumulate(+,1:Nstack+1))
topology = StackedCubedSphereTopology(MPI.COMM_SELF, Nhorz, Rrange)

x1, x2, x3 = ntuple(j->reshape(topology.elemtocoord[j, :, :],
2, 2, 2, length(topology.elems)), 3)
for n = 1:length(x1)
x1[n], x2[n], x3[n] = Topologies.cubed_sphere_warp(EquiangularCubedSphere(),x1[n], x2[n], x3[n])
end

Note that the faces are listed in Cartesian order.

source

### Functions

ClimateMachine.Mesh.Topologies.cubedshellmeshFunction
cubedshellmesh(T, Ne; part=1, numparts=1)

Generate a cubed mesh where each of the "cubes" has an Ne X Ne grid of elements.

The mesh can optionally be partitioned into numparts and this returns partition part. This is a simple Cartesian partition and further partitioning (e.g, based on a space-filling curve) should be done before the mesh is used for computation.

This mesh returns the cubed spehere in a flattened fashion for the vertex values, and a remapping is needed to embed the mesh in a 3-D space.

The mesh structures for the cubes is as follows:

x2
^
|
4Ne-           +-------+
|           |       |
|           |   6   |
|           |       |
3Ne-           +-------+
|           |       |
|           |   5   |
|           |       |
2Ne-           +-------+
|           |       |
|           |   4   |
|           |       |
Ne-   +-------+-------+-------+
|   |       |       |       |
|   |   1   |   2   |   3   |
|   |       |       |       |
0-   +-------+-------+-------+
|
+---|-------|-------|------|-> x1
0      Ne      2Ne    3Ne
source
ClimateMachine.Mesh.Topologies.cubed_sphere_warpFunction
cubed_sphere_warp(::EquiangularCubedSphere, a, b, c, R = max(abs(a), abs(b), abs(c)))

Given points (a, b, c) on the surface of a cube, warp the points out to a spherical shell of radius R based on the equiangular gnomonic grid proposed by C. Ronchi , R. Iacono , P. S. Paolucci (1996)

source
cubed_sphere_warp(a, b, c, R = max(abs(a), abs(b), abs(c)))

Given points (a, b, c) on the surface of a cube, warp the points out to a spherical shell of radius R based on the equidistant gnomonic grid outlined in M. Ran\v { c } i\' { c } , R. J. Purser , F. Mesinger (1996) and Ramachandran D Nair , Stephen J Thomas , Richard D Loft (2005)

source
cubed_sphere_warp(::ConformalCubedSphere, a, b, c, R = max(abs(a), abs(b), abs(c)))

Given points (a, b, c) on the surface of a cube, warp the points out to a spherical shell of radius R based on the conformal grid proposed by M. Ran\v { c } i\' { c } , R. J. Purser , F. Mesinger (1996)

source
ClimateMachine.Mesh.Topologies.cubed_sphere_unwarpFunction
cubed_sphere_unwarp(x1, x2, x3)

The inverse of cubed_sphere_warp. This function projects a given point (x_1, x_2, x_3) from the surface of a sphere onto a cube

source
cubed_sphere_unwarp(x1, x2, x3)

The inverse of cubed_sphere_warp. This function projects a given point (x_1, x_2, x_3) from the surface of a sphere onto a cube

source
cubed_sphere_unwarp(::ConformalCubedSphere, x1, x2, x3)

The inverse of cubed_sphere_warp. This function projects a given point (x_1, x_2, x_3) from the surface of a sphere onto a cube M. Ran\v { c } i\' { c } , R. J. Purser , F. Mesinger (1996)

source
ClimateMachine.Mesh.Topologies.cubed_sphere_topo_warpFunction
cubed_sphere_topo_warp(a, b, c, R = max(abs(a), abs(b), abs(c));
topography = NoTopography())

Given points (a, b, c) on the surface of a cube, warp the points out to a spherical shell of radius R based on the equiangular gnomonic grid proposed by C. Ronchi , R. Iacono , P. S. Paolucci (1996). Assumes a user specified modified radius using the computeanalyticaltopography function. Defaults to smooth cubed sphere unless otherwise specified via the AnalyticalTopography type.

source
ClimateMachine.Mesh.Topologies.grid1dFunction
grid1d(a, b[, stretch::AbstractGridStretching]; elemsize, nelem)

Discretize the 1D interval [a,b] into elements. Exactly one of the following keyword arguments must be provided:

• elemsize: the average element size, or
• nelem: the number of elements.

The optional stretch argument allows stretching, otherwise the element sizes will be uniform.

Returns either a range object or a vector containing the element boundaries.

source

## Geometry

ClimateMachine.Mesh.Geometry.LocalGeometryType
LocalGeometry

The local geometry at a nodal point.

Constructors

LocalGeometry{Np, N}(vgeo::AbstractArray{T}, n::Integer, e::Integer)

Extracts a LocalGeometry object from the vgeo array at node n in element e with Np being the number of points in the element and N being the polynomial order

Fields

• polyorder

polynomial order of the element

• coord

local degree of freedom Cartesian coordinate

• invJ

Jacobian from Cartesian to element coordinates: invJ[i,j] is $∂ξ_i / ∂x_j$

• vgeo

Global volume geometry array

• n

element local linear node index

• e

process local element index

source
ClimateMachine.Mesh.Geometry.resolutionmetricFunction
resolutionmetric(g::LocalGeometry)

The metric tensor of the discretisation resolution. Given a unit vector u in Cartesian coordinates and M = resolutionmetric(g), sqrt(u'*M*u) is the degree-of-freedom density in the direction of u.

source

## Brick Mesh

ClimateMachine.Mesh.BrickMesh.partitionFunction
partition(comm::MPI.Comm, elemtovert, elemtocoord, elemtobndy,
faceconnections)

This function takes in a mesh (as returned for example by brickmesh) and returns a Hilbert curve based partitioned mesh.

source
ClimateMachine.Mesh.BrickMesh.brickmeshFunction
brickmesh(x, periodic; part=1, numparts=1; boundary)

Generate a brick mesh with coordinates given by the tuple x and the periodic dimensions given by the periodic tuple.

The brick can optionally be partitioned into numparts and this returns partition part. This is a simple Cartesian partition and further partitioning (e.g, based on a space-filling curve) should be done before the mesh is used for computation.

By default boundary faces will be marked with a one and other faces with a zero. Specific boundary numbers can also be passed for each face of the brick in boundary. This will mark the nonperiodic brick faces with the given boundary number.

Examples

We can build a 3 by 2 element two-dimensional mesh that is periodic in the $x_2$-direction with

julia> (elemtovert, elemtocoord, elemtobndy, faceconnections) =
brickmesh((2:5,4:6), (false,true); boundary=((1,2), (3,4)));

This returns the mesh structure for

         x_2

^
|
6-  9----10----11----12
|  |     |     |     |
|  |  4  |  5  |  6  |
|  |     |     |     |
5-  5-----6-----7-----8
|  |     |     |     |
|  |  1  |  2  |  3  |
|  |     |     |     |
4-  1-----2-----3-----4
|
+--|-----|-----|-----|--> x_1
2     3     4     5

The (number of corners by number of elements) array elemtovert gives the global vertex number for the corners of each element.

julia> elemtovert
4×6 Array{Int64,2}:
1  2  3   5   6   7
2  3  4   6   7   8
5  6  7   9  10  11
6  7  8  10  11  12

Note that the vertices are listed in Cartesian order.

The (dimension by number of corners by number of elements) array elemtocoord gives the coordinates of the corners of each element.

julia> elemtocoord
2×4×6 Array{Int64,3}:
[:, :, 1] =
2  3  2  3
4  4  5  5

[:, :, 2] =
3  4  3  4
4  4  5  5

[:, :, 3] =
4  5  4  5
4  4  5  5

[:, :, 4] =
2  3  2  3
5  5  6  6

[:, :, 5] =
3  4  3  4
5  5  6  6

[:, :, 6] =
4  5  4  5
5  5  6  6

The (number of faces by number of elements) array elemtobndy gives the boundary number for each face of each element. A zero will be given for connected faces.

julia> elemtobndy
4×6 Array{Int64,2}:
1  0  0  1  0  0
0  0  2  0  0  2
0  0  0  0  0  0
0  0  0  0  0  0

Note that the faces are listed in Cartesian order.

Finally, the periodic face connections are given in faceconnections which is a list of arrays, one for each connection. Each array in the list is given in the format [e, f, vs...] where

• e is the element number;
• f is the face number; and
• vs is the global vertices that face associated with.

I the example

julia> faceconnections
3-element Array{Array{Int64,1},1}:
[4, 4, 1, 2]
[5, 4, 2, 3]
[6, 4, 3, 4]

we see that face 4 of element 5 is associated with vertices [2 3] (the vertices for face 1 of element 2).

source
ClimateMachine.Mesh.BrickMesh.connectmeshFunction
connectmesh(comm::MPI.Comm, elemtovert, elemtocoord, elemtobndy,
faceconnections)

This function takes in a mesh (as returned for example by brickmesh) and returns a connected mesh. This returns a NamedTuple of:

• elems the range of element indices
• realelems the range of real (aka nonghost) element indices
• ghostelems the range of ghost element indices
• ghostfaces ghost element to face is received; ghostfaces[f,ge] == true if face f of ghost element ge is received.
• sendelems an array of send element indices
• sendfaces send element to face is sent; sendfaces[f,se] == true if face f of send element se is sent.
• elemtocoord element to vertex coordinates; elemtocoord[d,i,e] is the dth coordinate of corner i of element e
• elemtoelem element to neighboring element; elemtoelem[f,e] is the number of the element neighboring element e across face f. If there is no neighboring element then elemtoelem[f,e] == e.
• elemtoface element to neighboring element face; elemtoface[f,e] is the face number of the element neighboring element e across face f. If there is no neighboring element then elemtoface[f,e] == f.
• elemtoordr element to neighboring element order; elemtoordr[f,e] is the ordering number of the element neighboring element e across face f. If there is no neighboring element then elemtoordr[f,e] == 1.
• elemtobndy element to bounday number; elemtobndy[f,e] is the boundary number of face f of element e. If there is a neighboring element then elemtobndy[f,e] == 0.
• nabrtorank a list of the MPI ranks for the neighboring processes
• nabrtorecv a range in ghost elements to receive for each neighbor
• nabrtosend a range in sendelems to send for each neighbor
source
ClimateMachine.Mesh.BrickMesh.connectmeshfullFunction
connectmeshfull(comm::MPI.Comm, elemtovert, elemtocoord, elemtobndy,
faceconnections)

This function takes in a mesh (as returned for example by brickmesh) and returns a corner connected mesh. It is similar to connectmesh but returns a corner connected mesh, i.e. ghostelems will also include remote elements that are only connected by vertices. This returns a NamedTuple of:

• elems the range of element indices
• realelems the range of real (aka nonghost) element indices
• ghostelems the range of ghost element indices
• ghostfaces ghost element to face is received; ghostfaces[f,ge] == true if face f of ghost element ge is received.
• sendelems an array of send element indices
• sendfaces send element to face is sent; sendfaces[f,se] == true if face f of send element se is sent.
• elemtocoord element to vertex coordinates; elemtocoord[d,i,e] is the dth coordinate of corner i of element e
• elemtoelem element to neighboring element; elemtoelem[f,e] is the number of the element neighboring element e across face f. If there is no neighboring element then elemtoelem[f,e] == e.
• elemtoface element to neighboring element face; elemtoface[f,e] is the face number of the element neighboring element e across face f. If there is no neighboring element then elemtoface[f,e] == f.
• elemtoordr element to neighboring element order; elemtoordr[f,e] is the ordering number of the element neighboring element e across face f. If there is no neighboring element then elemtoordr[f,e] == 1.
• elemtobndy element to bounday number; elemtobndy[f,e] is the boundary number of face f of element e. If there is a neighboring element then elemtobndy[f,e] == 0.
• nabrtorank a list of the MPI ranks for the neighboring processes
• nabrtorecv a range in ghost elements to receive for each neighbor
• nabrtosend a range in sendelems to send for each neighbor

The algorithm assumes that the 2-D mesh can be gathered on single process, which is reasonable for the 2-D mesh.

source
ClimateMachine.Mesh.BrickMesh.centroidtocodeFunction
centroidtocode(comm::MPI.Comm, elemtocorner; coortocode, CT)

Returns a code for each element based on its centroid.

These element codes can be used to determine a linear ordering for the partition function.

The communicator comm is used to calculate the bounding box for representing the centroids in coordinates of type CT, defaulting to CT=UInt64. These integer coordinates are converted to a code using the function coortocode, which defaults to hilbertcode.

The array containing the element corner coordinates, elemtocorner, is used to compute the centroids. elemtocorner is a dimension by number of corners by number of elements array.

source

## GeometricFactors

GeometricFactors groups data structures that collect geometric terms data needed at each quadrature point, in each element.

### Types

ClimateMachine.Mesh.GeometricFactors.VolumeGeometryType
VolumeGeometry{Nq, AA <: AbstractArray, A <: AbstractArray}

A struct that collects VolumeGeometry fields:

• array: Array contatining the data stored in a VolumeGeometry struct (the following fields are views into this array)
• ∂ξk/∂xi: Derivative of the Cartesian reference element coordinate ξ_k with respect to the Cartesian physical coordinate x_i
• ωJ: Mass matrix. This is the physical mass matrix, and thus contains the Jacobian determinant, J .* (ωᵢ ⊗ ωⱼ ⊗ ωₖ), where ωᵢ are the quadrature weights and J is the Jacobian determinant, det(∂x/∂ξ)
• ωJI: Inverse mass matrix: 1 ./ ωJ
• ωJH: Horizontal mass matrix (used in diagnostics), J .* norm(∂ξ3/∂x) * (ωᵢ ⊗ ωⱼ); for integrating over a plane (in 2-D ξ2 is used, not ξ3)
• xi: Nodal degrees of freedom locations in Cartesian physical space
• JcV: Metric terms for vertical line integrals norm(∂x/∂ξ3) (in 2-D ξ2 is used, not ξ3)
• ∂xk/∂ξi: Inverse of matrix ∂ξk/∂xi that represents the derivative of Cartesian physical coordinate x_i with respect to Cartesian reference element coordinate ξ_k
source
ClimateMachine.Mesh.GeometricFactors.SurfaceGeometryType
SurfaceGeometry{Nq, AA, A <: AbstractArray}

A struct that collects VolumeGeometry fields:

• array: Array contatining the data stored in a SurfaceGeometry struct (the following fields are views into this array)
• ni: Outward pointing unit normal in physical space
• sωJ: Surface mass matrix. This is the physical mass matrix, and thus contains the surface Jacobian determinant, sJ .* (ωⱼ ⊗ ωₖ), where ωᵢ are the quadrature weights and sJ is the surface Jacobian determinant
• vωJI: Volume mass matrix at the surface nodes (needed in the lift operation, i.e., the projection of a face field back to the volume). Since DGSEM is used only collocated, volume mass matrices are required.
source

## Metrics

Metrics encode the computation of metric terms defined at each quadrature point, in each element.

### Functions

ClimateMachine.Mesh.Metrics.creategrid!Function
creategrid!(vgeo, elemtocoord, ξ)

Create a 1-D grid using elemtocoord (see brickmesh) using the 1-D (-1, 1) reference coordinates ξ (in 1D, ξ = ξ1). The element grids are filled using linear interpolation of the element coordinates.

If Nq = length(ξ) and nelem = size(elemtocoord, 3) then the preallocated array vgeo.x1 should be Nq * nelem == length(x1).

source
creategrid!(vgeo, elemtocoord, ξ)

Create a 2-D tensor product grid using elemtocoord (see brickmesh) using the tuple ξ = (ξ1, ξ2), composed by the 1D reference coordinates ξ1 and ξ2 in (-1, 1)^2. The element grids are filled using bilinear interpolation of the element coordinates.

If Nq = (length(ξ1), length(ξ2)) and nelem = size(elemtocoord, 3) then the preallocated arrays vgeo.x1 and vgeo.x2 should be prod(Nq) * nelem == size(vgeo.x1) == size(vgeo.x2).

source
creategrid!(vgeo, elemtocoord, ξ)

Create a 3-D tensor product grid using elemtocoord (see brickmesh) using the tuple ξ = (ξ1, ξ2, ξ3), composed by the 1D reference coordinates ξ1, ξ2, ξ3 in (-1, 1)^3. The element grids are filled using trilinear interpolation of the element coordinates.

If Nq = (length(ξ1), length(ξ2), length(ξ3)) and nelem = size(elemtocoord, 3) then the preallocated arrays vgeo.x1, vgeo.x2, and vgeo.x3 should be prod(Nq) * nelem == size(vgeo.x1) == size(vgeo.x2) == size(vgeo.x3).

source
ClimateMachine.Mesh.Metrics.compute_reference_to_physical_coord_jacobian!Function
compute_reference_to_physical_coord_jacobian!(vgeo, nelem, D)

Input arguments:

• vgeo::VolumeGeometry, a struct containing the volumetric geometric factors
• D::NTuple{2,Int}, a tuple of derivative matrices, i.e., D = (D1,), where:
• D1::DAT2, 1-D derivative operator on the device in the first dimension

Compute the Jacobian matrix, ∂x / ∂ξ, of the transformation from reference coordinates, ξ1, to physical coordinates, vgeo.x1, for each quadrature point in element e.

source
compute_reference_to_physical_coord_jacobian!(vgeo, nelem, D)

Input arguments:

• vgeo::VolumeGeometry, a struct containing the volumetric geometric factors
• D::NTuple{2,Int}, a tuple of derivative matrices, i.e., D = (D1, D2), where:
• D1::DAT2, 1-D derivative operator on the device in the first dimension
• D2::DAT2, 1-D derivative operator on the device in the second dimension

Compute the Jacobian matrix, ∂x / ∂ξ, of the transformation from reference coordinates, ξ1, ξ2, to physical coordinates, vgeo.x1, vgeo.x2, for each quadrature point in element e.

source
compute_reference_to_physical_coord_jacobian!(vgeo, nelem, D)

Input arguments:

• vgeo::VolumeGeometry, a struct containing the volumetric geometric factors
• D::NTuple{3,Int}, a tuple of derivative matrices, i.e., D = (D1, D2, D3), where:
• D1::DAT2, 1-D derivative operator on the device in the first dimension
• D2::DAT2, 1-D derivative operator on the device in the second dimension
• D3::DAT2, 1-D derivative operator on the device in the third dimension

Compute the Jacobian matrix, ∂x / ∂ξ, of the transformation from reference coordinates, ξ1, ξ2, ξ3 to physical coordinates, vgeo.x1, vgeo.x2, vgeo.x3 for each quadrature point in element e.

source
ClimateMachine.Mesh.Metrics.computemetric!Function
computemetric!(vgeo, sgeo, D)

Input arguments:

• vgeo::VolumeGeometry, a struct containing the volumetric geometric factors
• sgeo::SurfaceGeometry, a struct containing the surface geometric factors
• D::NTuple{1,Int}, tuple with 1-D derivative operator on the device

Compute the 1-D metric terms from the element grid arrays vgeo.x1. All the arrays are preallocated by the user and the (square) derivative matrix D should be consistent with the reference grid ξ1 used in creategrid!.

If Nq = size(D, 1) and nelem = div(length(x1), Nq) then the volume arrays x1, J, and ξ1x1 should all have length Nq * nelem. Similarly, the face arrays sJ and n1 should be of length nface * nelem with nface = 2.

source
computemetric!(vgeo, sgeo, D)

Input arguments:

• vgeo::VolumeGeometry, a struct containing the volumetric geometric factors
• sgeo::SurfaceGeometry, a struct containing the surface geometric factors
• D::NTuple{2,Int}, a tuple of derivative matrices, i.e., D = (D1, D2), where:
• D1::DAT2, 1-D derivative operator on the device in the first dimension
• D2::DAT2, 1-D derivative operator on the device in the second dimension

Compute the 2-D metric terms from the element grid arrays vgeo.x1 and vgeo.x2. All the arrays are preallocated by the user and the (square) derivative matrice D1 and D2 should be consistent with the reference grid ξ1 and ξ2 used in creategrid!.

If Nq = (size(D1, 1), size(D2, 1)) and nelem = div(length(vgeo.x1), prod(Nq)) then the volume arrays vgeo.x1, vgeo.x2, vgeo.ωJ, vgeo.ξ1x1, vgeo.ξ2x1, vgeo.ξ1x2, and vgeo.ξ2x2 should all be of size (Nq..., nelem). Similarly, the face arrays sgeo.sωJ, sgeo.n1, and sgeo.n2 should be of size (maximum(Nq), nface, nelem) with nface = 4

source
computemetric!(vgeo, sgeo, D)

Input arguments:

• vgeo::VolumeGeometry, a struct containing the volumetric geometric factors
• sgeo::SurfaceGeometry, a struct containing the surface geometric factors
• D::NTuple{3,Int}, a tuple of derivative matrices, i.e., D = (D1, D2, D3), where:
• D1::DAT2, 1-D derivative operator on the device in the first dimension
• D2::DAT2, 1-D derivative operator on the device in the second dimension
• D3::DAT2, 1-D derivative operator on the device in the third dimension

Compute the 3-D metric terms from the element grid arrays vgeo.x1, vgeo.x2, and vgeo.x3. All the arrays are preallocated by the user and the (square) derivative matrice D1, D2, and D3 should be consistent with the reference grid ξ1, ξ2, and ξ3 used in creategrid!.

If Nq = size(D1, 1) and nelem = div(length(vgeo.x1), Nq^3) then the volume arrays vgeo.x1, vgeo.x2, vgeo.x3, vgeo.ωJ, vgeo.ξ1x1, vgeo.ξ2x1, vgeo.ξ3x1, vgeo.ξ1x2, vgeo.ξ2x2, vgeo.ξ3x2, vgeo.ξ1x3,vgeo.ξ2x3, and vgeo.ξ3x3 should all be of length Nq^3 * nelem. Similarly, the face arrays sgeo.sωJ, sgeo.n1, sgeo.n2, and sgeo.n3 should be of size Nq^2 * nface * nelem with nface = 6.

The curl invariant formulation of Kopriva (2006), equation 37, is used.

Reference:

source

## Grids

Grids specify the approximation within each element, and any necessary warping.

### Functions

ClimateMachine.Mesh.Grids.get_zFunction
get_z(grid; z_scale = 1, rm_dupes = false)

Get the Gauss-Lobatto points along the Z-coordinate.

• grid: DG grid
• z_scale: multiplies z-coordinate
• rm_dupes: removes duplicate Gauss-Lobatto points
source
ClimateMachine.Mesh.Grids.min_node_distanceFunction
min_node_distance(::AbstractGrid, direction::Direction=EveryDirection() )

Returns an approximation of the minimum node distance in physical space.

source
min_node_distance(
::DiscontinuousSpectralElementGrid,
direction::Direction=EveryDirection())
)

Returns an approximation of the minimum node distance in physical space along the reference coordinate directions. The direction controls which reference directions are considered.

source
ClimateMachine.Mesh.Grids.DiscontinuousSpectralElementGridType
DiscontinuousSpectralElementGrid(topology; FloatType, DeviceArray,
polynomialorder,
meshwarp = (x...)->identity(x))

Generate a discontinuous spectral element (tensor product, Legendre-Gauss-Lobatto) grid/mesh from a topology, where the order of the elements is given by polynomialorder. DeviceArray gives the array type used to store the data (CuArray or Array), and the coordinate points will be of FloatType.

The polynomial order can be different in each direction (specified as a NTuple). If only a single integer is specified, then each dimension will use the same order. If the topology dimension is 3 and the polynomialorder has dimension 2, then the first value will be used for horizontal and the second for the vertical.

The optional meshwarp function allows the coordinate points to be warped after the mesh is created; the mesh degrees of freedom are orginally assigned using a trilinear blend of the element corner locations.

source
ClimateMachine.Mesh.Grids.computegeometryFunction
computegeometry(elemtocoord, D, ξ, ω, meshwarp)

Compute the geometric factors data needed to define metric terms at each quadrature point. First, compute the so called "topology coordinates" from reference coordinates ξ. Then map these topology coordinate to physical coordinates. Then compute the Jacobian of the mapping from reference coordinates to physical coordinates, i.e., ∂x/∂ξ, by calling compute_reference_to_physical_coord_jacobian!. Finally, compute the metric terms by calling the function computemetric!.

source

## DSS

Computes the direct stiffness summation of fields in the MPIStateArray.

ClimateMachine.Mesh.DSS.dss!Function
dss3d(Q::MPIStateArray,
grid::DiscontinuousSpectralElementGrid)

This function computes the 3D direct stiffness summation for all variables in the MPIStateArray.

Fields

• Q: MPIStateArray
• grid: Discontinuous Spectral Element Grid
source
ClimateMachine.Mesh.DSS.dss_vertex!Function
dss_vertex!(
vtconn,
vtconnoff,
vertmap,
data,
ivar,
vx,
::Type{FT},

) where {FT}

This function computes the direct stiffness summation for the vertex vx.

Fields

• vtconn: vertex connectivity array
• vtconnoff: offsets for vertex connectivity array
• vertmap: map to vertex degrees of freedom: vertmap[vx] contains the degree of freedom located at vertex vx.
• data: data field of MPIStateArray
• ivar: variable # in the MPIStateArray
• vx: unique edge number
• ::Type{FT}: Floating point type
source
ClimateMachine.Mesh.DSS.dss_edge!Function
dss_edge!(
edgconn,
edgconnoff,
edgemap,
Nemax,
data,
ivar,
ex,
::Type{FT},

) where {FT}

This function computes the direct stiffness summation for all degrees of freedom corresponding to edge ex. dss_edge! applies only to interior (non-vertex) edge nodes.

Fields

• edgconn: edge connectivity array
• edgconnoff: offsets for edge connectivity array
• edgemap: map to edge degrees of freedom: edgemap[i, edgno, orient] contains the element node index of the ith interior node on edge edgno, under orientation orient.
• Nemax: # of relevant degrees of freedom per edge (other dof are marked as -1)
• data: data field of MPIStateArray
• ivar: variable # in the MPIStateArray
• ex: unique edge number
• ::Type{FT}: Floating point type
source
ClimateMachine.Mesh.DSS.dss_face!Function
dss_face!(fcconn, facemap, Nfmax, data, ivar, fx)

This function computes the direct stiffness summation for all degrees of freedom corresponding to face fx. dss_face! applies only to interior (non-vertex and non-edge) face nodes.

Fields

• fcconn: face connectivity array
• facemap: map to face degrees of freedom: facemap[ij, fcno, orient] contains the element node index of the ijth interior node on face fcno under orientation orient
• Nfmax: # of relevant degrees of freedom per face (other dof are marked as -1)
• data: data field of MPIStateArray
• ivar: variable # in the MPIStateArray
• fx: unique face number
source

## Filters

There are methods used to cleanup state vectors.

ClimateMachine.Mesh.Filters.BoydVandevenFilterType
BoydVandevenFilter(grid, Nc=0, s=32)

Returns the spectral filter using the logarithmic error function of the form:

$$$σ(η) = 1/2 erfc(2*sqrt(s)*χ(η)*(abs(η)-0.5))$$$

whenever s ≤ i ≤ N, and 1 otherwise. The function χ(η) is defined as

$$$χ(η) = sqrt(-log(1-4*(abs(η)-0.5)^2)/(4*(abs(η)-0.5)^2))$$$

if x != 0.5 and 1 otherwise. Here, s is the filter order, the filter starts with polynomial order Nc, and alpha is a parameter controlling the smallest value of the filter function.

References

source
ClimateMachine.Mesh.Filters.ExponentialFilterType
ExponentialFilter(grid, Nc=0, s=32, α=-log(eps(eltype(grid))))

Returns the spectral filter with the filter function

$$$σ(η) = xp(-α η^s)$$$

where s is the filter order (must be even), the filter starts with polynomial order Nc, and alpha is a parameter controlling the smallest value of the filter function.

source
ClimateMachine.Mesh.Filters.TMARFilterType
TMARFilter()

Returns the truncation and mass aware rescaling nonnegativity preservation filter. The details of this filter are described in Devin Light , Dale Durran (2016)

Note this needs to be used with a restrictive time step or a flux correction to ensure that grid integral is conserved.

Examples

This filter can be applied to the 3rd and 4th fields of an MPIStateArray Q with the code

Filters.apply!(Q, (3, 4), grid, TMARFilter())

where grid is the associated DiscontinuousSpectralElementGrid.

source
ClimateMachine.Mesh.Filters.apply!Function
Filters.apply!(Q::MPIStateArray,
target,
grid::DiscontinuousSpectralElementGrid,
filter::AbstractSpectralFilter;
kwargs...)

Applies filter to Q given a grid and a custom target.

A target can be any of the following:

• a tuple or range of indices
• a tuple of symbols or strings of variable names
• a colon (:) to apply to all variables
• a custom [AbstractFilterTarget]

The following keyword arguments are supported for some filters:

• direction: for AbstractSpectralFilter controls if the filter is applied in the horizontal and/or vertical directions. It is assumed that the trailing dimension on the reference element is the vertical dimension and the rest are horizontal.
• state_auxiliary: if target requires auxiliary state to compute its argument or results.

Examples

Specifying the target via indices:

Filters.apply!(Q, :, grid, TMARFilter())
Filters.apply!(Q, (1, 3), grid, CutoffFilter(grid); direction=VerticalDirection())

Speciying target via symbols or strings:

Filters.apply!(Q, (:ρ, "energy.ρe"), grid, TMARFilter())
Filters.apply!(Q, ("moisture.ρq_tot",), grid, CutoffFilter(grid);
direction=VerticalDirection())
source
ClimateMachine.Mesh.Filters.apply_async!Function
Filters.apply_async!(Q, target, grid::DiscontinuousSpectralElementGrid,
filter::AbstractFilter;
dependencies,
kwargs...)

An asynchronous version of Filters.apply!, returning an Event object. dependencies should be an Event or tuple of Events which need to finish before applying the filter.

compstream = Filters.apply_async!(Q, :, grid, CutoffFilter(grid); dependencies=compstream)
wait(compstream)
source

## Interpolation

### Types

ClimateMachine.Mesh.Interpolation.InterpolationBrickType
InterpolationBrick{
FT <: AbstractFloat,CuArrays
UI16VD <: AbstractVector{UInt16},
I32V <: AbstractVector{Int32},
} <: InterpolationTopology

This interpolation data structure and the corresponding functions works for a brick, where stretching/compression happens only along the x1, x2 & x3 axis. Here x1 = X1(ξ1), x2 = X2(ξ2) and x3 = X3(ξ3).

Fields

• Nel

Number of elements

• Np

Total number of interpolation points

• Npl

Total number of interpolation points on local process

• xbnd

Domain bounds in x1, x2 and x3 directions

• x1g

Interpolation grid in x1 direction

• x2g

Interpolation grid in x2 direction

• x3g

Interpolation grid in x3 direction

• ξ1

Unique ξ1 coordinates of interpolation points within each spectral element

• ξ2

Unique ξ2 coordinates of interpolation points within each spectral element

• ξ3

Unique ξ3 coordinates of interpolation points within each spectral element

• flg

Flags when ξ1/ξ2/ξ3 interpolation point matches with a GLL point

• fac

Normalization factor

• x1i

x1 interpolation grid index of interpolation points within each element on the local process

• x2i

x2 interpolation grid index of interpolation points within each element on the local process

• x3i

x3 interpolation grid index of interpolation points within each element on the local process

• offset

Offsets for each element

• m_ξ1

GLL points in ξ1 direction

• m_ξ2

GLL points in ξ2 direction

• m_ξ3

GLL points in ξ3 direction

• wb1

Barycentric weights

• wb2

Barycentric weights

• wb3

Barycentric weights

• Np_all

Number of interpolation points on each of the processes

• x1i_all

x1 interpolation grid index of interpolation points within each element on all processes stored only on proc 0

• x2i_all

x2 interpolation grid index of interpolation points within each element on all processes stored only on proc 0

• x3i_all

x3 interpolation grid index of interpolation points within each element on all processes stored only on proc 0

Usage

InterpolationBrick(
grid::DiscontinuousSpectralElementGrid{FT},
xbnd::Array{FT,2},
xres,
) where FT <: AbstractFloat

This interpolation structure and the corresponding functions works for a brick, where stretching/compression happens only along the x1, x2 & x3 axis. Here x1 = X1(ξ1), x2 = X2(ξ2) and x3 = X3(ξ3).

Arguments for the inner constructor

• grid: DiscontinousSpectralElementGrid
• xbnd: Domain boundaries in x1, x2 and x3 directions
• x1g: Interpolation grid in x1 direction
• x2g: Interpolation grid in x2 direction
• x3g: Interpolation grid in x3 direction
source
ClimateMachine.Mesh.Interpolation.InterpolationCubedSphereType
InterpolationCubedSphere{
FT <: AbstractFloat,
T <: Int,
FTV <: AbstractVector{FT},
FTVD <: AbstractVector{FT},
TVD <: AbstractVector{T},
UI16VD <: AbstractVector{UInt16},
I32V <: AbstractVector{Int32},
} <: InterpolationTopology

This interpolation structure and the corresponding functions works for a cubed sphere topology. The data is interpolated along a lat/long/rad grid.

-90⁰ ≤ lat ≤ 90⁰

-180⁰ ≤ long ≤ 180⁰

Rᵢ ≤ r ≤ Rₒ

Fields

• Nel

Number of elements

• Np

Number of interpolation points

• Npl

Number of interpolation points on local process

• n_rad

Number of interpolation points in radial direction

• n_lat

Number of interpolation points in lat direction

• n_long

Number of interpolation points in long direction

• rad_grd

• lat_grd

Interpolation grid in lat direction

• long_grd

Interpolation grid in long direction

• ξ1

Device array containing ξ1 coordinates of interpolation points within each element

• ξ2

Device array containing ξ2 coordinates of interpolation points within each element

• ξ3

Device array containing ξ3 coordinates of interpolation points within each element

• flg

flags when ξ1/ξ2/ξ3 interpolation point matches with a GLL point

• fac

Normalization factor

• radi

Radial coordinates of interpolation points withing each element

• lati

Latitude coordinates of interpolation points withing each element

• longi

Longitude coordinates of interpolation points withing each element

• offset

Offsets for each element

• m_ξ1

GLL points in ξ1 direction

• m_ξ2

GLL points in ξ2 direction

• m_ξ3

GLL points in ξ3 direction

• wb1

Barycentric weights in ξ1 direction

• wb2

Barycentric weights in ξ2 direction

• wb3

Barycentric weights in ξ3 direction

• Np_all

Number of interpolation points on each of the processes

• radi_all

Radial interpolation grid index of interpolation points within each element on all processes stored only on proc 0

• lati_all

Latitude interpolation grid index of interpolation points within each element on all processes stored only on proc 0

• longi_all

Longitude interpolation grid index of interpolation points within each element on all processes stored only on proc 0

Usage

InterpolationCubedSphere(grid::DiscontinuousSpectralElementGrid, vert_range::AbstractArray{FT}, nhor::Int, lat_res::FT, long_res::FT, rad_res::FT) where {FT <: AbstractFloat}

This interpolation structure and the corresponding functions works for a cubed sphere topology. The data is interpolated along a lat/long/rad grid.

-90⁰ ≤ lat ≤ 90⁰

-180⁰ ≤ long ≤ 180⁰

Rᵢ ≤ r ≤ Rₒ

Arguments for the inner constructor

• grid: DiscontinousSpectralElementGrid
• vert_range: Vertex range along the radial coordinate
• lat_res: Resolution of the interpolation grid along the latitude coordinate in radians
• long_res: Resolution of the interpolation grid along the longitude coordinate in radians
• rad_res: Resolution of the interpolation grid along the radial coordinate
source

### Functions

ClimateMachine.Mesh.Interpolation.interpolate_local!Function
interpolate_local!(
intrp_brck::InterpolationBrick{FT},
sv::AbstractArray{FT},
v::AbstractArray{FT},
) where {FT <: AbstractFloat}

This interpolation function works for a brick, where stretching/compression happens only along the x1, x2 & x3 axis. Here x1 = X1(ξ1), x2 = X2(ξ2) and x3 = X3(ξ3)

Arguments

• intrp_brck: Initialized InterpolationBrick structure
• sv: State Array consisting of various variables on the discontinuous Galerkin grid
• v: Interpolated variables
source
interpolate_local!(intrp_cs::InterpolationCubedSphere{FT},
sv::AbstractArray{FT},
v::AbstractArray{FT}) where {FT <: AbstractFloat}

This interpolation function works for cubed spherical shell geometry.

Arguments

• intrp_cs: Initialized cubed sphere structure
• sv: Array consisting of various variables on the discontinuous Galerkin grid
• v: Array consisting of variables on the interpolated grid
source
ClimateMachine.Mesh.Interpolation.project_cubed_sphere!Function
project_cubed_sphere!(intrp_cs::InterpolationCubedSphere{FT},
v::AbstractArray{FT},
uvwi::Tuple{Int,Int,Int}) where {FT <: AbstractFloat}

This function projects the velocity field along unit vectors in radial, lat and long directions for cubed spherical shell geometry.

Fields

• intrp_cs: Initialized cubed sphere structure
• v: Array consisting of x1, x2 and x3 components of the vector field
• uvwi: Tuple providing the column numbers for x1, x2 and x3 components of vector field in the array. These columns will be replaced with projected vector fields along unit vectors in long, lat and rad directions.
source
ClimateMachine.Mesh.Interpolation.accumulate_interpolated_data!Function
accumulate_interpolated_data!(intrp::InterpolationTopology,
iv::AbstractArray{FT,2},
fiv::AbstractArray{FT,4}) where {FT <: AbstractFloat}

This interpolation function gathers interpolated data onto process # 0.

Fields

• intrp: Initialized interpolation topology structure
• iv: Interpolated variables on local process
• fiv: Full interpolated variables accumulated on process # 0
source