Curvilinear elements and metric terms

On a curved element the divergence of a constant flux is discretised as $\sum_a \mathbf D_a\,(J\mathbf a^a)\cdot\mathbf F_\text{const}$. It vanishes only if the discrete metric identities

hold. If they do not, a uniform flow generates spurious waves at every curved element. Gyre satisfies the identities to round-off in 2-D and 3-D, including with viscous terms:

The geometry is first interpolated from the geometry nodes to the solution nodes (degree $p$), then differentiated with the same $\mathbf D$ used for the solution.

Two dimensions. The cross-product form is already divergence-free discretely:

Three dimensions. Products such as $y_\eta z_\zeta$ are of degree $2p$ and their derivatives no longer commute after interpolation, so the cross-product form violates the identities. Gyre uses Kopriva’s conservative curl form,

where $I^p$ is nodal interpolation at degree $p$. Being the discrete curl of a polynomial, each contravariant vector is discretely divergence-free by construction.

Outward face normals and surface Jacobians are the face values of the same arrays: $\mathbf n\, s_J = \pm J\mathbf a^a$ on the two faces normal to $\xi_a$, which makes the normals of two elements sharing a face exactly opposite — a requirement for conservation.

Metric terms are not precomputed constants: Discretization.geometry(X) is an ordinary JAX function of the global node coordinates $\mathbf X$. The residual is $\mathbf R(\mathbf U, \mathbf X, \theta)$ and the chain

is differentiated end-to-end by reverse-mode AD. There is no separate “mesh sensitivity” code to keep consistent with the solver. The verification study on mesh deformation checks the derivative of a geometric integral through this chain against finite differences.

The coarsest p-multigrid level uses one node per element. Differentiating a one-point “polynomial” gives no metrics, so this level takes its geometry from the degree-$p_g$ description by exact reduction: the cell Jacobian is the mean $\sum w J / 2^d$, face normals are the face integrals $\int \mathbf n\,\mathrm dS$ divided by the reference face area. For the inviscid terms the resulting scheme is a first-order finite-volume method on the curved cell whose faces close exactly (geometric conservation), verified in test_finite_volume_geometry_reduction. The viscous terms keep only the face-normal two-point difference at this level (penalty parameter).

All built-in generators evaluate a smooth mapping on the LGL lattice, so elements are curved to degree $p_g$:

A detail worth knowing: classical cosine clustering along the chord has zero parametric speed at the trailing edge, which makes the Jacobian vanish at the trailing-edge vertex. The aerofoil generator uses $x = \sin^2(\tfrac{c\pi}{2} v)/\sin^2(\tfrac{c\pi}{2})$ with $c < 1$, which keeps the smooth $\sqrt{x}$ behaviour at the nose but a finite tangent at the trailing edge.

A mesh read from a file (mesh.generator: file) becomes the same global node array and element lattice as a generated one, so everything on this page applies unchanged. For a Gmsh file the reader maps Gmsh’s node ordering of each element onto the tensor-product lattice, re-orients left-handed elements, and interpolates the equispaced Gmsh nodes to the Legendre–Gauss–Lobatto distribution of the same degree (the identity up to degree 2). Lattice nodes are numbered topologically (a vertex by its file node, an edge or face node by the vertices of that edge or face and its position on it), so neighbouring elements share nodes exactly and nodes that merely coincide in space, at a slit or a periodic pair, stay distinct. Elements may meet in any relative orientation. A mesh is refused if the Jacobian determinant is not positive at every geometry node. The imported-mesh study tests these steps.