Meshing Stuff

AbstractTopology{dim}

Represents the connectivity of individual elements, with local dimension dim.

BoxElementTopology{dim, T} <: AbstractTopology{dim}

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

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

• hasboundary

  boolean for whether or not this topology has a boundary

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

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

StackedBrickTopology{dim, T} <: AbstractTopology{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}

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

StackedCubedSphereTopology{3, T} <: AbstractTopology{3}

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

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}

Generate a cubed shell mesh with the number of elements along each dimension of the cubes being Nelem. This topology actual creates a cube mesh, and the warping should be done after the grid is created using the cubedshellwarp 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 elements on each cube, total 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.cubedshellwarp(x1[n], x2[n], x3[n])
end

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

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 cubedshellwarp 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, total 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.cubedshellwarp(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 with 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 flatten 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

cubedshellwarp(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 Ronchi, Iacono, Paolucci (1996) https://doi.org/10.1006/jcph.1996.0047

@article{RonchiIaconoPaolucci1996,
  title={The ``cubed sphere'': a new method for the solution of partial
         differential equations in spherical geometry},
  author={Ronchi, C. and Iacono, R. and Paolucci, P. S.},
  journal={Journal of Computational Physics},
  volume={124},
  number={1},
  pages={93--114},
  year={1996},
  doi={10.1006/jcph.1996.0047}
}

hasboundary(topology::AbstractTopology)

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

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

CutoffFilter(grid, Nc=polynomialorder(grid))

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

ExponentialFilter(grid, Nc=0, s=32, α=-log(eps(eltype(grid))))

Returns the spectral filter with the filter function

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

@article{doi:10.1175/MWR-D-16-0220.1,
  author = {Light, Devin and Durran, Dale},
  title = {Preserving Nonnegativity in Discontinuous Galerkin
           Approximations to Scalar Transport via Truncation and Mass
           Aware Rescaling (TMAR)},
  journal = {Monthly Weather Review},
  volume = {144},
  number = {12},
  pages = {4771-4786},
  year = {2016},
  doi = {10.1175/MWR-D-16-0220.1},
}

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.

InterpolationBrick{
FT <: AbstractFloat,
T <: Int,
FTV <: AbstractVector{FT},
FTVD <: AbstractVector{FT},
TVD <: AbstractVector{T},
FTA2 <: Array{FT, 2},
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

• poly_order

  Polynomial order of spectral element approximation

• 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

• m1_r

  GLL points

• m1_w

  GLL weights

• wb

  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
• xres: Resolution of the interpolation grid in x1, x2 and x3 directions

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

• poly_order

  Polynomial order of spectral element approximation

• 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

• m1_r

  GLL points

• m1_w

  GLL weights

• wb

  Barycentric weights

• 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 rad, lat and long 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