Mesh warping

[simonbyrne]

Add infrastructure for generic (non-regular, warped) meshes

1. 2D non-regular mesh (e.g. example from Taylor & Fournier, section 3.6)
2. Cubed-sphere mesh, equiangular, S2 geometry

[valeriabarra]

Just putting down here some ideas that @jkozdon suggested and some questions that came up when we chatted about this yesterday, so that we can use these as a starting point for our discussion today.

Questions/ideas:

• We need to define a data structure to store the mesh in (are we interested in reading the mesh from file, generating it internally, or both?)
• Edge connectivity (2 elements/cells are connected to an edge in 2D)
• Corner connectivity (must support arbitrary number of connected elements in 2D or only quads?)

Steps:

• Define a data structure to store the mesh in
• Define an element2vertex map, then from the element2vertex we build the connectivity
• Additional utilities/query functions, eg, for each edge return the pair of vertices v1 and v2 it connects (with orientation, perhaps v1 < v2 --- in Julia we can use tuples for this)
• To test: Start with a simple structured mesh and deform the internal vertices by a sinusoidal, user-defined warping function (leave identity as default warping function)

[simonbyrne]

The required interfaces are specified in https://github.com/CliMA/ClimaCore.jl/blob/main/src/Topologies/Topologies.jl. The three main parts are:

1. Need a way to determine face pairs: map (elem, face) to neighbouring (nbrelem, nbrface, reversed).

   Can do this two ways: a. ClimateMachine.jl way: for each (elem, face), store an arrays of elemtoelem, elemtoface, elemtoordr (which basically matches our reversed): https://github.com/CliMA/ClimateMachine.jl/blob/74964c5f54af2f4d22e9229736f9d353e5ef2c0e/src/Numerics/Mesh/Topologies.jl#L118-L137 b. store all mesh faces as (elem, face, nbrelem, nbrface, reversed)

2. A way to iterate over boundary faces: the easiest option is to store a vector of vectors: for each boundary tag, store a vector of (elem, face) pairs that are on that boundary

3. Iterate over mesh vertices, and the (elem, vert) pairs that are at that point (to support arbitrary vertex valence): a. array of (elem, vert) pairs, grouped by mesh vertex number b. array of integer offsets for each global vertex, corresponding to the index of the first (elem, vert) pair.

[simonbyrne]

Tests for the structured mesh are here: https://github.com/CliMA/ClimaCore.jl/blob/main/test/grid.jl

[jkozdon]

For the unstructured mesh reader I would view it as a sort problem.

For triangles (with DG not CG) this is done with something like this: https://github.com/tcew/nodal-dg/blob/b6864cae74d27b2ff28ad51afec66061ce8336f7/Codes1.1/Codes2D/tiConnect2D.m

[simonbyrne]

To handle the warping, the simplest might be to build some sort of wrapper that does the warps an existing topology, or just make the warping an extra argument to the space constructor. Would need some sort of extra "patch" argument when different sets of warpings are involved. We could have a way to store metadata for each element.

See https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DM/DMGetLabel.html

[valeriabarra]

The required interfaces are specified in https://github.com/CliMA/ClimaCore.jl/blob/main/src/Topologies/Topologies.jl. The three main parts are:

1. Need a way to determine face pairs: map (elem, face) to neighbouring (nbrelem, nbrface, reversed). Can do this two ways: a. ClimateMachine.jl way: for each (elem, face), store an arrays of elemtoelem, elemtoface, elemtoordr (which basically matches our reversed): https://github.com/CliMA/ClimateMachine.jl/blob/74964c5f54af2f4d22e9229736f9d353e5ef2c0e/src/Numerics/Mesh/Topologies.jl#L118-L137 b. store all mesh faces as (elem, face, nbrelem, nbrface, reversed)

As I mentioned in the meeting. Instead of going fully unstructured, I decided to start with tensor-product meshes first (i.e., meshes that still have a fix valence equal to 4 for interior vertices, but where vertices are not necessarily equispaced). I thought this would be an intermediate level of difficulty between a structured grid mesh and a fully unstructured grid mesh. This decision was made not only to break down the tasks in smaller, intermediate tasks, but also based on what Tapio was proposing.

Point 1 is done for tensor-product meshes and fully tested. (I went the 1.b way)

[valeriabarra]

Point 2 is also done for tensor-product meshes and fully tested

[valeriabarra]

Point 3 is also done for tensor-product meshes and fully tested. Note, this is done only for tensor-product meshes for now, so not with arbitrary valence as described in point 3 above.

Next step: test a warped mesh with a sinusoidal distortion

[valeriabarra]

Would the warping function act on the Mesh, Topology or Space object? If it has to be applied similarly to what we had in ClimateMachine, that is, on the physical space variable x (as between L65 and L66 here) and this information does not go back to the Mesh, then there is literally not much difference between a regular rectangular equispaced grid that is distorted afterwards and a tensor-product mesh that is defined to be "irregular" by definition.

[simonbyrne]

We currently have hierarchy of objects

1. domain: e.g. RectangleDomain
2. mesh: e.g. EquispacedRectangleMesh
3. topology: e.g. GridTopology
4. space: e.g. SpectralElementSpace2D

After some discussion with various people, I think we approach this as follows:

Domains

We have two ways to specify this:

1. warp an underlying domain, e.g. apply a coordinate transform w(x') = x, where x' is the coordinate in the "underlying" domain, and x is the coordinate in the new domain. In this case we would have a wrapper type:

   domain = WarpedDomain(underlyingdomain, warpfun)

2. specify the curved domain, eg. Sphere(radius=2.0)

Mesh

Here we define a WarpedMesh type, which takes an underlying mesh (i.e. a mesh on the underlying domain). Correspondingly:

1. Define functions for WarpedDomains

   WarpedEquispacedRectangleMesh(domain::WarpedDomain, args...) = 
   WarpedMesh(domain, EquispacedRectangleMesh(domain.underlying, args...), domain.warpfun)

2. Define aliases for curved domains:

   EquiangularCubedSphere(domain::Sphere, n) =
   WarpedMesh(domain, CubeMesh(n), EquiangularWarp())

Note that in the case of 2., we would need to add some notion of a "patch" (or whatever we want to call it): this is a grouping of elements with a common coordinate system, and the transformation would be applied to all elements in the patch. As a result, it might make sense to define warping as an object with certain interfaces, e.g.

warp(::Warping, x', patch) -> x

for convenience, we could also define inverse (though this might not be unique), and specific methods to transform to different coordinates (Cartesian, Spherical, etc).

Topology

Here the topology needs only be defined on the underlying mesh (e.g. EquispacedRectangleMesh or UnitCubeMesh), but would contain a reference to the warped mesh. Note that we need to define new CubeTopology or something similar.

Space

We would keep the existing SpectralElementSpace2D, but add the following options:

1. placement of quadrature points: This implies some sort of mapping from the reference element to the warped domain. Options to support here are:

   a. bilinear interpolation in the underlying domain, then apply the warp function: i.e. if x1, x2, x3, x4 are the physical coordinates and x1' = w⁻¹(x1), etc are the underlying coordinates, we place nodes by:

   f(ξ1, ξ2) = w(bilinear(x1', x2', x3', x4', ξ1, ξ2))

   where ξ1, ξ2 are the quadrature points in the reference element. b. via bilinear interpolation in physical Cartesian coordinates, followed by a projection to the sphere (as per Guba et al 2014)

   f(ξ1, ξ2) = project(bilinear(x1, x2, x3, x4, ξ1, ξ2)))

2. computation of metric terms (partial derivatives, Jacobians, etc). Options to support are: a. differentiate through the f (i think this matches the "analytic" in the design docs) b. compute the derivative of the interpolating polynomial, which would depend on coordinate system: I think this refers to "numeric" in the design docs. c. curl-invariant metrics (Kopriva), from a or b (not sure if this would work in non-Cartesian coordinates, but can be added later)

[jkozdon]

Seems reasonable to me.

[jakebolewski]

I think we might not want a warped domain object, rather we should just compute a new domain object under some warping passed through to the mesh. That way the 1D vertical (stretching) and 2D horizontal (warping) would have a similar API, and we could get away with a single domain type. We could save the underlying domain similar to how we save the underlying unwarped mesh, and define consistent accessor functions for all forms (1D, 2D, 3D, 1D extruded, 2D extruded, etc.).

https://github.com/CliMA/ClimaCore.jl/pull/82#discussion_r697013564

[simonbyrne]

A few more thoughts:

1. The WarpedMesh object should be lightweight: basically just a wrapper around the underlying mesh:

   struct WarpedMesh
   domain
   underlyingmesh
   warpfun
   end

2. If we just construct the topology on the underlying mesh, we then would need to pass both the WarpedMesh and the topology to the space constructor. This is probably not too bad?

3. Add a facet(topology, elem) function that returns the facet of a given element (can just be an integer). If there is only a single facet, then this is 1. vertex_coordinates would then return coordinate that is local to that facet. e.g. on a CubeTopology, facet would be in 1:6, and each coordinate on a facet be a XYPoint. The orientations should ideally match CubedSpheres.jl.

This should be sufficient to define a CubeTopology

4. The warp function needs to do a few things:
   • transform a (facet, facet_coordinate) to a new global coordinate: for a cubed sphere, this would be from (facet, xypoint) and return a 3D Cartesian point (we need the point in cartesian coordinates to do the bilinear interpolation in 1b above)
   • it should be continuous for each facet, so that we can use AD through it.
   • we then need some mechanism to convert the Cartesian point (and Cartesian vectors arising from the local_geometry terms) to local lat/long and a UVVector: not sure how best to handle this, perhaps some extra "canonicalize" step?
   • it could be helpful to have an unwarp function (but could be done later).

[jakebolewski]

The ECMWF library Atlas calls "facets" -> "zones", this is inherited from the Geodesy community and standards:

https://en.wikipedia.org/wiki/Discrete_global_grid

I don't know if there is more common name for this in the atmospheric sciences community but aligning with the standard name used in Geospatial applications seems not a bad choice.

[tapios]

The ECMWF library Atlas calls "facets" -> "zones", this is inherited from the Geodesy community and standards:

Zones for us are usually latitude bands. Facet isn't bad, I think.

[valeriabarra]

There are several names for this in the literature (facet, chart, patch, panel). I think "facet" is too similar to "face" to type/search in the code.

[valeriabarra]

A few more thoughts:

1. The WarpedMesh object should be lightweight: basically just a wrapper around the underlying mesh:

struct WarpedMesh
   domain
   underlyingmesh
   warpfun
end

Thanks @simonbyrne . I have a couple of questions.

1. Just to make sure I understand, in this intended interface, we would not save the warped mesh, meaning the result of applying the warping function to the underlying mesh. Would this lead to any performance slow-downs if/when we want the warped mesh coordinates outside of the Space constructor? Or we envision that we won't need this for any other purposes, say, visualization?
2. If the warp function is a point-wise function (i.e., it takes a point on the underlying mesh and returns a point on the deformed mesh) we won't have a way to distinguish boundary points. I thought this could have been useful if we wanted to deform only the interior points of a given grid (for instance, like I have done in this warp mesh test that I have now pulled out of the tensor-product mesh PR but needs some reworking). If, instead, the warpfun acted on an entire mesh object, then it would make sense to have a place to store its returning value (ie, the warped mesh).

[simonbyrne]

I personally like "chart": https://en.wikipedia.org/wiki/Atlas_(topology), but patch or panel seems fine.

[valeriabarra]

I personally like "chart": https://en.wikipedia.org/wiki/Atlas_(topology), but patch or panel seems fine.

Yeah, to me chart is more the map, as in the link you shared. I would personally go with panel

[glwagner]

I prefer "patch" between (facet, chart, patch, panel), but I am ok with "panel" as well.