Meshing Stuff

AbstractTopology{dim, T, nb}

Represents the connectivity of individual elements, with local dimension dim with nb boundary conditions types. The element coordinates are of type T.

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

BrickTopology{dim, T, nb} <: AbstractTopology{dim, T, nb}

A simple grid-based topology. This is a convenience wrapper around BoxElementTopology.

StackedBrickTopology{dim, T, nb} <: AbstractStackedTopology{dim}

A simple grid-based topology, where all elements on the trailing dimension are stacked to be contiguous. This is a convenience wrapper around BoxElementTopology.

CubedShellTopology{T} <: AbstractTopology{2, T, 0}

A cube-shell topology. This is a convenience wrapper around BoxElementTopology.

AnalyticalTopography

Abstract type to allow dispatch over different analytical topography prescriptions in experiments.

NoTopography <: AnalyticalTopography

Allows definition of fallback methods in case cubedspheretopo_warp is used with no prescribed topography function.

DCMIPMountain <: AnalyticalTopography

Topography description based on standard DCMIP experiments.

StackedCubedSphereTopology{T, nb} <: AbstractStackedTopology{3, T, nb}

A cube-sphere topology. All elements on the same "vertical" dimension are stacked to be contiguous. This is a convenience wrapper around BoxElementTopology.

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

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.

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.

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...

# Shell radius = 1
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

# Shell radius = 10
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.

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.

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

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)

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)

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)

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

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

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)

equiangular_cubed_sphere_warp(a, b, c, R = max(abs(a), abs(b), abs(c)))

A wrapper function for the cubedspherewarp function, when called with the EquiangularCubedSphere type

equiangular_cubed_sphere_unwarp(x1, x2, x3)

A wrapper function for the cubedsphereunwarp function, when called with the EquiangularCubedSphere type

equidistant_cubed_sphere_warp(a, b, c, R = max(abs(a), abs(b), abs(c)))

A wrapper function for the cubedspherewarp function, when called with the EquidistantCubedSphere type

equidistant_cubed_sphere_unwarp(x1, x2, x3)

A wrapper function for the cubedsphereunwarp function, when called with the EquidistantCubedSphere type

conformal_cubed_sphere_warp(a, b, c, R = max(abs(a), abs(b), abs(c)))

A wrapper function for the cubedspherewarp function, when called with the ConformalCubedSphere type

conformal_cubed_sphere_unwarp(x1, x2, x3)

A wrapper function for the cubedsphereunwarp function, when called with the ConformalCubedSphere type

hasboundary(topology::AbstractTopology)

query function to check whether a topology has a boundary (i.e., not fully periodic)

compute_lat_long(X,Y,δ,faceid)

Helper function to allow computation of latitute and longitude coordinates given the cubed sphere coordinates X, Y, δ, faceid

cubed_sphere_topo_warp(a, b, c, R = max(abs(a), abs(b), abs(c));
                   r_inner = _planet_radius,
                   r_outer = _planet_radius + domain_height,
                   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.

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.

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

lengthscale(g::LocalGeometry)

The effective geometric mean grid resolution at the point.

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.

lengthscale_horizontal(g::LocalGeometry)

The effective horizontal grid resolution at the point.

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.

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

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

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.

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.

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

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.

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

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

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

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.

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.

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.

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.

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

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:

• David A Kopriva (2006)

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

referencepoints(::AbstractGrid)

Returns the points on the reference element.

referencepoints(::DiscontinuousSpectralElementGrid)

Returns the 1D interpolation points used for the reference element.

min_node_distance(::AbstractGrid, direction::Direction=EveryDirection() )

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

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.

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.

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!.

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

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

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

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

CutoffFilter(grid, Nc=polynomialorders(grid))

Returns the spectral filter that zeros out polynomial modes greater than or equal to Nc.

MassPreservingCutoffFilter(grid, Nc=polynomialorders(grid))

Returns the spectral filter that zeros out polynomial modes greater than or equal to Nc while preserving the cell average value. Use this filter if the jacobian is nonconstant.

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

• JP Boyd (1996)

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.

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.

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

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)

InterpolationBrick{
    FT <: AbstractFloat,CuArrays
    UI8AD <: AbstractArray{UInt8, 2},
    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

InterpolationCubedSphere{
FT <: AbstractFloat,
T <: Int,
FTV <: AbstractVector{FT},
FTVD <: AbstractVector{FT},
TVD <: AbstractVector{T},
UI8AD <: AbstractArray{UInt8, 2},
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

  Interpolation grid in radial direction

• 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

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

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

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.

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