Meshes

AbstractMesh{dim}

Abstract supertype of meshes, which describe how a domain is discretized into elements. Subtypes: IntervalMesh, RectilinearMesh, and the AbstractCubedSphere meshes.

A mesh is lightweight and exists on all MPI ranks; a mesh stored in a file, for example, holds only the filename.

Notes

Face and vertex numbering. In 1D, faces and vertices coincide, and both are numbered [1, 2]. In 2D, a face is a line segment between two vertices, and both are numbered [1, 2, 3, 4] in counter-clockwise order.

 v4        f3        v3
   o-----------------o
   |                 |	    face    vertices
   |                 |	      f1 =>  v1 v2
f4 |                 | f2     f2 =>  v2 v3
   |                 |	      f3 =>  v3 v4
   |                 |        f4 =>  v4 v1
   |                 |
   o-----------------o
  v1       f1        v2

Interface. A subtype of AbstractMesh defines the following methods:

• domain(mesh)
• elements(mesh)
• is_boundary_face(mesh, elem, face)
• boundary_face_name(mesh, elem, face)
• opposing_face(mesh, elem, face)
• coordinates(mesh, elem, vert)
• containing_element (optional)

The following types and methods are provided for every AbstractMesh:

• SharedVertices(mesh, elem, vert)
• face_connectivity_matrix(mesh[,elemorder])
• vertex_connectivity_matrix(mesh[,elemorder])

AbstractMesh1D

Alias for AbstractMesh{1} (see AbstractMesh): the supertype of one-dimensional meshes, i.e. IntervalMesh, whose elements are intervals.

AbstractMesh2D

Alias for AbstractMesh{2} (see AbstractMesh): the supertype of two-dimensional meshes, i.e. RectilinearMesh and the AbstractCubedSphere meshes, whose elements are quadrilaterals.

IntervalMesh <: AbstractMesh

One-dimensional mesh on an IntervalDomain. Elements are numbered 1:nelems.

Fields

• stretch: The stretching rule used to place the faces.
• domain: The IntervalDomain.
• faces: Vector of face coordinates, of length nelems + 1.
• meta: Stretching-specific metadata (e.g. the solved stretching parameter), or nothing.
• reverse_mode: Whether the smallest element is at the top of the domain.

Constructor

IntervalMesh(domain::IntervalDomain, faces::AbstractVector)

Construct a 1D mesh with face locations at faces.

IntervalMesh(domain::IntervalDomain, stretching = Uniform(); nelems, reverse_mode = false)

Construct a 1D mesh on domain with nelems elements, using stretching. Possible values of stretching are:

• Uniform()
• ExponentialStretching(H)
• GeneralizedExponentialStretching(dz_bottom, dz_top)
• HyperbolicTangentStretching(dz_surface)
• BoundaryLayerHyperbolicTangentStretching(dz_surface, dz_top, z_pbl)

If reverse_mode is true, the smallest element is at the top, and the largest at the bottom (this is typical for land model configurations).

Examples

using ClimaCore: Domains, Geometry, Meshes
domain = Domains.IntervalDomain(
    Geometry.ZPoint(0.0),
    Geometry.ZPoint(30000.0);
    boundary_names = (:bottom, :top),
)
uniform_mesh = Meshes.IntervalMesh(domain; nelems = 30)
face_mesh = Meshes.IntervalMesh(domain, Geometry.ZPoint.([0.0, 100.0, 1000.0, 30000.0]))

RectilinearMesh <: AbstractMesh2D

Rectilinear mesh of a RectangleDomain: the product of two IntervalMeshes. Elements are indexed by CartesianIndex{2}.

Constructor

RectilinearMesh(domain::RectangleDomain, n1, n2)

Construct a RectilinearMesh of n1 by n2 equally spaced elements on domain.

RectilinearMesh(intervalmesh1::IntervalMesh, intervalmesh2::IntervalMesh)

Construct the product mesh of intervalmesh1 and intervalmesh2.

AbstractCubedSphere <: AbstractMesh2D

Abstract supertype of cubed-sphere meshes on SphereDomains. Subtypes: EquiangularCubedSphere, EquidistantCubedSphere, and ConformalCubedSphere. A cubed-sphere mesh has 6 panels, laid out as follows:

                                          :   Panel 1   :
                            +-------------+-------------+
                            |     +x1     |     +x1     |
                            |             |             |
                            |    Panel    |    Panel    |
                            |+x3   5   -x3|-x2   6   +x2|
                            |     -x2     |     -x3     |
                            |             |             |
                            |     -x1     |     -x1     |
              +-------------+-------------+-------------+
              |     -x2     |     -x2     |
              |             |             |
              |    Panel    |    Panel    |
              |+x1   3   -x1|+x3   4   -x3|
              |     +x3     |     -x1     |
              |             |             |
              |     +x2     |     +x2     |
+-------------+-------------+-------------+
|     +x3     |     +x3     |
|             |             |
|    Panel    |    Panel    |
|-x2   1   +x2|+x1   2   -x1|
|     +x1     |     +x2     |
|             |             |
|     -x3     |     -x3     |
+-------------+-------------+
:   Panel 6   :

This is the panel ordering of the S2 Geometry library, with 1-based instead of 0-based numbering.

Elements are indexed by a CartesianIndex{3}, whose components are:

• the horizontal element index (left to right) within each panel,
• the vertical element index (bottom to top) within each panel,
• the panel number.

Fields

Subtypes have the following fields:

• domain: The SphereDomain.
• ne: Number of elements across each panel.
• localelementmap: The LocalElementMap from the reference element to the sphere.

Notes

External references:

• S2Geometry library.
• MIT GCM exch2.

EquiangularCubedSphere <: AbstractCubedSphere

An equiangular gnomonic mesh proposed by [6]. Uses the element indexing convention of AbstractCubedSphere.

Constructor

EquiangularCubedSphere(
    domain::Domains.SphereDomain,
    ne::Integer,
    localelementmap = NormalizedBilinearMap(),
)

Construct an EquiangularCubedSphere on domain with ne elements across each panel.

EquidistantCubedSphere <: AbstractCubedSphere

An equidistant gnomonic mesh outlined in [49] and [50]. Uses the element indexing convention of AbstractCubedSphere.

Constructor

EquidistantCubedSphere(
    domain::Domains.SphereDomain,
    ne::Integer,
    localelementmap = NormalizedBilinearMap(),
)

Construct an EquidistantCubedSphere on domain with ne elements across each panel.

ConformalCubedSphere <: AbstractCubedSphere

A conformal mesh outlined in [49]. Uses the element indexing convention of AbstractCubedSphere.

Constructor

ConformalCubedSphere(
    domain::Domains.SphereDomain,
    ne::Integer,
    localelementmap = NormalizedBilinearMap(),
)

Construct a ConformalCubedSphere on domain with ne elements across each panel.

LocalElementMap

Abstract supertype of mappings from the reference element to the physical domain. Subtypes: IntrinsicMap and NormalizedBilinearMap.

IntrinsicMap()

LocalElementMap that uses the intrinsic mapping of the cubed sphere to map the reference element to the physical domain.

NormalizedBilinearMap()

LocalElementMap for meshes on spherical domains of [51]. It uses bilinear interpolation between the Cartesian coordinates of the element vertices, then normalizes the result to lie on the sphere.

StretchingRule

Abstract type of the rules that place the faces of an IntervalMesh constructed from a number of elements. Each subtype S has a method IntervalMesh(domain, rule::S; nelems, kwargs...) that builds the faces, and, if it has parameters, a method Meshes._convert_stretch(FT, rule::S) that converts them to the float type of the domain. Subtypes: Uniform, ExponentialStretching, GeneralizedExponentialStretching, HyperbolicTangentStretching, and BoundaryLayerHyperbolicTangentStretching.

Uniform()

Use uniformly-sized elements.

ExponentialStretching(H::FT)

Apply exponential stretching to the domain when constructing elements. H is the scale height [m]; a typical atmospheric scale height is H ≈ 7.5e3.

For an interval $[z_0, z_1]$, the elements are uniformly spaced in $ζ$, where

\[\zeta = \frac{1 - e^{-\eta/h}}{1-e^{-1/h}},\]

where $\eta = \frac{z - z_0}{z_1-z_0}$, and $h = \frac{H}{z_1-z_0}$ is the nondimensional scale height. The resulting faces are reference heights without terrain warping.

Fields

• H: Scale height [m].

Keyword Arguments

Keyword arguments of IntervalMesh(domain, ExponentialStretching(H); ...):

• nelems: Number of elements.
• reverse_mode = false: Whether the smallest element is at the top and the largest at the bottom, as in land model configurations.

Examples

using ClimaCore: Domains, Geometry, Meshes
domain = Domains.IntervalDomain(
    Geometry.ZPoint(0.0),
    Geometry.ZPoint(30000.0);
    boundary_names = (:bottom, :top),
)
mesh = Meshes.IntervalMesh(domain, Meshes.ExponentialStretching(7500.0); nelems = 30)

GeneralizedExponentialStretching(dz_bottom::FT, dz_top::FT)

Apply exponential stretching with a scale height that varies with height, so that the elements have spacing dz_bottom at the bottom and dz_top at the top of the domain [m].

For an interval $[z_0, z_1]$ with $N$ elements, the elements are uniformly spaced in $\zeta$, where

\[\eta = -h(\zeta) \ln\left[1 - \left(1 - e^{-1/h(\zeta)}\right)\zeta\right]\]

and $\eta = \frac{z - z_0}{z_1-z_0}$. The nondimensional scale height $h(\zeta)$ varies linearly from $h_\text{bottom}$ at $\zeta = 1/N$ to $h_\text{top}$ at $\zeta = 1 - 1/N$; both are solved so that the first and last elements have the target spacings. With constant $h$, this is ExponentialStretching. In typical atmosphere configurations, dz_bottom is the smallest spacing and dz_top the largest. In typical land configurations, dz_bottom is the largest and dz_top the smallest, and reverse_mode = true mirrors the profile so that its fine end is at the top. The resulting faces are reference heights without terrain warping.

Fields

• dz_bottom: Target element spacing at the bottom of the domain [m].
• dz_top: Target element spacing at the top of the domain [m].

Keyword Arguments

Keyword arguments of IntervalMesh(domain, GeneralizedExponentialStretching(dz_bottom, dz_top); ...):

• nelems: Number of elements, at least 2.
• FT_solve = Float64: Float type of the root solves for $h_\text{bottom}$ and $h_\text{top}$.
• tol = sqrt(eps(FT_solve)) * min(dz_bottom, dz_top): Residual tolerance of the root solves, a spacing [m].
• reverse_mode = false: Whether dz_bottom ≥ dz_top, with the fine end of the profile at the top, as in land model configurations.

Examples

using ClimaCore: Domains, Geometry, Meshes
domain = Domains.IntervalDomain(
    Geometry.ZPoint(-10.0),
    Geometry.ZPoint(0.0);
    boundary_names = (:bottom, :top),
)
stretch = Meshes.GeneralizedExponentialStretching(1.5, 0.1)
mesh = Meshes.IntervalMesh(domain, stretch; nelems = 20, reverse_mode = true)

HyperbolicTangentStretching(dz_surface::FT)

Apply hyperbolic tangent stretching with element spacing dz_surface at the surface [m]. The surface is the bottom of the domain, or the top with reverse_mode = true, as in land model configurations.

For an interval $[z_0, z_1]$, the elements are uniformly spaced in $\zeta$, where

\[\eta = 1 - \frac{\tanh[\gamma(1-\zeta)]}{\tanh(\gamma)}\]

and $\eta = \frac{z - z_0}{z_1-z_0}$. The stretching parameter $\gamma$ is solved so that the first element has spacing dz_surface, which requires dz_surface < (z_1 - z_0) / nelems. The resulting faces are reference heights without terrain warping.

Fields

• dz_surface: Target element spacing at the surface [m].

Keyword Arguments

Keyword arguments of IntervalMesh(domain, HyperbolicTangentStretching(dz_surface); ...):

• nelems: Number of elements, at least 2.
• FT_solve = Float64: Float type of the root solve for $\gamma$.
• tol = 1e-6 * dz_surface: Residual tolerance of the root solve, a spacing [m].
• reverse_mode = false: Whether the surface, and with it the smallest element, is at the top, as in land model configurations.

Examples

using ClimaCore: Domains, Geometry, Meshes
domain = Domains.IntervalDomain(
    Geometry.ZPoint(0.0),
    Geometry.ZPoint(30000.0);
    boundary_names = (:bottom, :top),
)
mesh = Meshes.IntervalMesh(domain, Meshes.HyperbolicTangentStretching(30.0); nelems = 43)

BoundaryLayerHyperbolicTangentStretching(dz_surface::FT, dz_top::FT, z_pbl::FT)

Apply two-scale hyperbolic tangent stretching that refines the boundary layer with scale height z_pbl [m] while matching element spacing dz_surface at the surface [m] and dz_top at the top of the domain [m].

For an interval $[z_0, z_1]$ of height $H = z_1 - z_0$ with $N$ elements and height above the surface $u = z - z_0 \in [0, H]$ with $\eta = u / H$, the cumulative level index $k(u) \in [0, N]$ is the sum of a background hyperbolic tangent profile and a boundary-layer hyperbolic tangent profile:

\[k(u; \gamma) = n_{\text{bg}}(\gamma) \left[1 - \frac{\operatorname{atanh}[(1 - \eta)\tanh(\gamma)]} {\gamma}\right] + \big(N - n_{\text{bg}}(\gamma)\big)\, \frac{\tanh(u / z_{\text{pbl}})}{\tanh(H / z_{\text{pbl}})},\]

where $n_{\text{bg}}(\gamma) = \gamma / \operatorname{atanh}[(\text{dz\_top} / H)\tanh(\gamma)]$ sets the top element spacing of the background profile to dz_top, and $\gamma > 0$ is solved so that $k(\text{dz\_surface}; \gamma) = 1$. The faces are the heights at which $k$ takes integer values, and the element spacings increase monotonically from dz_surface to dz_top. The resulting faces are reference heights without terrain warping.

nelems has a lower and an upper bound. Outside them, the constructor throws an ArgumentError that names the admissible range.

Fields

• dz_surface: Target element spacing at the surface [m].
• dz_top: Target element spacing at the far end of the domain [m].
• z_pbl: Scale height of the boundary-layer refinement [m].

Keyword Arguments

Keyword arguments of IntervalMesh(domain, BoundaryLayerHyperbolicTangentStretching(dz_surface, dz_top, z_pbl); ...):

• nelems: Number of elements, within the admissible range.
• FT_solve = Float64: Float type of the root solves.
• tol = sqrt(eps(FT_solve)): Residual tolerance of the root solves for $\gamma$ and the faces, in level-index units [-].
• reverse_mode = false: Whether the profile is mirrored so that dz_surface is the spacing at the top of the domain and dz_top the spacing at the bottom, as in land model configurations.

Examples

using ClimaCore: Domains, Geometry, Meshes
domain = Domains.IntervalDomain(
    Geometry.ZPoint(0.0),
    Geometry.ZPoint(60000.0);
    boundary_names = (:bottom, :top),
)
stretch = Meshes.BoundaryLayerHyperbolicTangentStretching(28.0, 3000.0, 2200.0)
mesh = Meshes.IntervalMesh(domain, stretch; nelems = 75)

Extended help

The lower bound on nelems is the number of levels of the single-scale HyperbolicTangentStretching profile that matches both dz_surface and dz_top; at it, the boundary-layer weight $N - n_{\text{bg}}$ is zero. The upper bound is the count beyond which the boundary-layer profile, which receives the levels the background profile does not, places the first face below dz_surface for every $\gamma$. In between, up to two values of $\gamma$ satisfy the surface condition. The constructor takes the larger, which assigns the most levels to the background profile.

Since $n_{\text{bg}}$ increases with $\gamma$, $N \ge n_{\text{bg}}(\gamma)$ holds for every admissible $\gamma$, so $k'(u) > 0$ and $k''(u) < 0$ on $(0, H)$ and the spacings increase monotonically. The boundary-layer profile is exponentially small near the top, so the top element spacing matches dz_top to within a relative error of order $e^{-2H / z_{\text{pbl}}}$.

truncate_mesh(parent_mesh::IntervalMesh, trunc_domain::IntervalDomain)

Construct an IntervalMesh on trunc_domain by truncating parent_mesh. The truncated parent faces run from the last parent face at or below the bottom of trunc_domain to the first at or above its top. The result has as many elements as these faces span, and uses a GeneralizedExponentialStretching with their bottom and top spacings; reverse_mode is set when the bottom spacing is the larger, as in land model meshes.

Meshes.domain(mesh::AbstractMesh)

Return the domain (a subtype of Domains.AbstractDomain) on which mesh is defined.

Meshes.elements(mesh::AbstractMesh)

Return an iterator over the elements of mesh. Elements can be of any type.

nelements(mesh::AbstractMesh)

Return the number of elements in mesh.

n_elements_per_panel_direction(::AbstractCubedSphere)

Return the number of elements per panel along each direction.

Meshes.coordinates(mesh, elem, vert::Int)
Meshes.coordinates(mesh, elem, ξ::SVector)

Return the physical coordinates of a point in an element elem of mesh. The position of the point can either be a vertex number vert or the coordinates ξ in the reference element.

Elements of a 2D mesh are quadrilaterals, so vert must be in 1:4; elements of a 1D mesh have two vertices, so vert must be in 1:2.

Meshes.containing_element(mesh::AbstractMesh, coord)

Return the element of mesh containing the coordinate coord. If the coordinate falls on the boundary between two or more elements, one of them is chosen arbitrarily.

Meshes.reference_coordinates(mesh::AbstractMesh, elem, coord)

Return an SVector ξ of coordinates in the reference element such that Meshes.coordinates(mesh, elem, ξ) == coord. Used for interpolation to a specific point.

Meshes.linearindices(elemorder)

Given a data structure elemorder with elemorder[i] = elem that orders elements, return the inverse map orderindex such that orderindex[elem] = i.

The result is a LinearIndices for CartesianIndices, a dense Int array for a vector of CartesianIndex, and a Dict otherwise.

Meshes.element_horizontal_length_scale(mesh::AbstractMesh)

Return the approximate length scale of the elements of mesh, in the units of the domain coordinates.