Source-linked AI summary
GridapGeosciences.jl: A Julia finite element package for partial differential equations on general manifolds
Tamara A. Tambyah, Alberto F. Martín, David Lee, Santiago Badia
TL;DR
Many finite element packages rely on separate compiled and interpreted languages, while PDE approximation on general manifolds calls for both intrinsic and extrinsic formulations. GridapGeosciences.jl extends Gridap with a mesh-refined atlas framework, and cubed-sphere applications recover wave dynamics on coarser meshes and demonstrate surface and vertical dynamics in three dimensions.
Problem
Many state-of-the-art finite element packages rely on separate compiled and interpreted languages, creating a two-language barrier.
Source: §1 · p. 1Method
GridapGeosciences.jl extends Gridap with atlas triangulations generated from coarse atlas data through mesh refinement, supporting intrinsic and extrinsic formulations.
Source: §3.1 · p. 4, §5 · p. 16Results
The correct wave number is recovered on a mesh coarser than those used in cited comparison studies, and three-dimensional simulations show oscillating vertical buoyancy and surface pressure waves.
Source: §4.2.2 · p. 13, §4.3.2 · p. 16Takeaways & Limitations
The cubed-sphere applications demonstrate the package’s capacity to represent both vertical buoyancy dynamics and surface pressure dynamics on a three-dimensional manifold.
Source: §4.3.2 · p. 16The study uses conforming meshes whose elements are all at the same refinement level, and leaves extension to evolving manifolds for future consideration.
Source: §5 · pp. 16–17
Abstract
from arXiv · showhide
We present GridapGeosciences$.$jl, a new parallel distributed-memory Julia package for the numerical approximation of partial differential equations on general manifolds. Our abstract approach to generating discrete manifold domains relies on mesh refinement in general, and is implemented within the Gridap finite element library, whose high-level interface enables seamless definition of both intrinsic and extrinsic variational formulations. GridapGeosciences$.$jl features the cubed sphere manifold, which is used in this study to showcase the functionality of GridapGeosciences$.$jl in different applications. The first application involves scalar transport where we develop intrinsic formulations of common stabilisation terms, which are shown to be equivalent to well-known extrinsic stabilisation formalisations. Such stabilisation terms are used to develop an intrinsic formulation of the two-dimensional thermal shallow water equations that includes upwinded numerical fluxes. The final application is a linearisation of the Boussinesq equations where we demonstrate the capacity of GridapGeosciences$.$jl to capture both surface and vertical dynamics on a three-dimensional manifold.
1. Introduction
GridapGeosciences.jl extends Julia’s Gridap finite element ecosystem to PDEs on manifolds, supporting both intrinsic and extrinsic formulations. Its atlas-refinement approach and cubed-sphere applications illustrate the framework’s scope for geophysical simulations.
1. Introduction
- GridapGeosciences.jl extends Gridap’s Julia finite element library to PDEs on manifolds, avoiding the compiled/interpreted two-language split common in FE packages.The package uses Gridap’s high-level interface to define intrinsic variational forms with general metric transformations.Source: §1 · p. 1
- The atlas-triangulation approach starts with one element per chart and uses mesh refinement to produce finer charts for general manifolds.For triangular or quadrilateral charts, affine mappings permit intrinsic variational forms to be integrated in reference space, unlike extrinsic forms computed in ambient space.Source: §1 · p. 1
- Supporting both intrinsic and extrinsic formulations is a key novelty compared with codes that provide only extrinsic functionality.Intrinsic formulations use a smooth geometric map and parametric charts, while extrinsic formulations embed the manifold in ambient space.Source: §1 · pp. 1–2
- The cubed sphere represents Earth in atmospheric applications and avoids the pole problem of latitude-longitude parametrisations.Compatible finite element spaces on this manifold support mimetic and conservation properties in the intrinsic approach, as developed in the companion methods paper [17].Source: §1 · p. 2
- The study demonstrates the framework with scalar transport and high-resolution thermal shallow-water and linearised Boussinesq applications on the cubed sphere.The Boussinesq case includes surface and vertical dynamics; the paper also discusses future extensions to other manifolds and evolving manifolds.Source: §1 · p. 2
2. Mathematical formulation
The formulation represents manifold PDEs in parametric charts, translating geometric operators and integrals through the chart map and metric. It then builds intrinsic finite element spaces and discrete systems on triangulations of the parametric domains.
2. Mathematical formulation
- The motivating Poisson problem seeks a function on a closed, compact, orientable manifold, with a compatible source and solution determined up to a constant.The source must have zero integral over the manifold.Source: §2 · p. 2
2.1. Parametric representation.
- An atlas chart pairs an open, bounded Euclidean parametric domain with a smooth map into the manifold; the framework’s discussion applies from one chart to an atlas.The manifold’s topological dimension n does not exceed its ambient dimension m.Source: §2.1 · p. 2
- The chart Jacobian defines the metric tensor, which relates ambient and parametric vector components and supplies the metric factor for transforming integrals.Parametric Lebesgue measure and manifold surface measure are related through the square root of the metric determinant.Source: §2.1 · p. 3
2.1.1. Surface operators.
- Surface gradient, divergence, and Laplace–Beltrami operators are expressed through Euclidean parametric-space derivatives and the metric.The parametric gradient and divergence are Euclidean operators, while the metric and its determinant encode the manifold geometry.Source: §2.1.1 · p. 3, §2.1 · p. 3
- In two-dimensional parametric space, the skew gradient rotates vectors anticlockwise; intrinsic cross-product and curl-related operators are also discussed for three-dimensional manifolds.The surface curl and vector Laplacian are not required for the atmospheric applications presented in the study.Source: §2.1.1 · p. 3
- Triangulating the Euclidean parametric domains and mapping the triangulation to the manifold yields a conforming manifold triangulation without approximating the parametric elements.The triangulation is assembled across charts for a multi-chart atlas.Source: §2.2 · p. 3
2.2. Parametric finite element methods.
- The intrinsic discretisation uses continuous Lagrangian scalar spaces and compatible de Rham spaces, with vector fields represented through flux proxies and the contravariant Piola map.The paper identifies Raviart–Thomas spaces for V1 and discontinuous polynomial spaces for V2.Source: §2.2 · pp. 3–4
- Integration by parts yields the intrinsic discrete system, whose metric transformations capture the manifold geometry and can be evaluated to machine precision with sufficient quadrature.This construction is described as free of geometric consistency error, in contrast to methods that approximate the curved manifold.Source: §2.2 · pp. 3–4
3. GridapGeosciences.jl
GridapGeosciences.jl extends Gridap to PDEs on manifolds through mesh-refined atlas triangulations, supporting intrinsic and extrinsic formulations. The cubed sphere illustrates the approach, including parallel Poisson convergence examples.
3. GridapGeosciences.jl
- GridapGeosciences.jl extends Gridap with a Julia library for numerical approximation of PDEs on manifolds.Its atlas triangulations use mesh refinement to generate discrete manifold domains.Source: §3 · p. 4, §3.1 · p. 4
3.1. Atlas triangulation.
- Atlas triangulations originate from coarse atlas information and are refined into parametric meshes while retaining each element’s root-chart identity.The coarse information includes a chart-wise geometrical map; chart gluing can be represented as a graph.Source: §3.1 · p. 4
3.1.1. Cubed sphere manifold.
- The cubed sphere uses six panels with local parametric coordinates and a cube-shaped chart topology, without a global parametric coordinate system.Its topology has no edge permutations, simplifying finite-element implementation but not being required in general.Source: §3.1.1 · p. 4
3.1.2. Extension to topological dimension 3.
- The three-dimensional cubed sphere extends the surface radially, with thickness and a topography profile controlling the shell’s radial boundaries.With zero topography, the radial coordinate is R + Tγ; the top shell is flat.Source: §3.1.2 · p. 6
3.2. Implementation.
- Gridap’s high-level API supports variational formulations on the physical domain, while its automatic engine transforms integrals to reference domains.Intrinsic formulations require users to include metric transformations; extrinsic formulations are written in ambient space.Source: §3.2 · p. 6
3.2.1. Example usage.
- The example driver defines an intrinsic Laplace–Beltrami formulation using metric fields within Gridap’s standard triangulation, integration, and finite-element interfaces.The example uses a cubed-sphere mesh with three refinement levels.Source: §3.2.1 · p. 6
3.3. Parallelisation.
- Distributed atlas triangulations use tree-based refinement and support isotropic meshes as well as anisotropic three-dimensional domains.GridapDistributed.jl and GridapP4est.jl provide distributed finite-element support and scalable mesh refinement and partitioning.Source: §3.3 · p. 7
3.3.1. Example usage in parallel.
- The parallel driver shows expected hp-convergence rates on two- and three-dimensional cubed sphere manifolds under h-refinement.For the three-dimensional case, the manufactured right-hand side includes a radial-boundary flux from integration by parts.Source: §3.3.1 · p. 8, Figure 2 · p. 8
4. Application
The application section introduces scalar transport on the cubed sphere and states the mesh-size convention and simulation setting.
4. Application
- The simulations use N = 6 × 2^nℓ spatial elements, with refinement level ℓ tending to L corresponding to h tending to zero.Tests are conducted on the Gadi supercomputer at NCI Australia.Source: §4 · p. 9
- The first application considers transport of a scalar by a known tangent velocity field without sources or sinks.The transport is described by an advection equation.Source: §4 · p. 9
4.1. Scalar transport.
The scalar-transport section develops intrinsic SUPG and upwinding stabilisations on a two-dimensional manifold and compares them with established extrinsic formulations. The tests report expected hp-convergence, equivalent intrinsic and extrinsic solutions, and decreasing mass-conservation error with refinement.
4.1. Scalar transport.
- The transport problem is posed on a two-dimensional manifold, where spurious oscillations can arise in semi-discrete finite-element formulations.The section considers SUPG for divergence-free velocity and upwinding for divergent velocity, developing intrinsic formulations for comparison with extrinsic ones.Source: §4.1 · p. 9
4.1.1. Streamline upwind Petrov–Galerkin method.
- For divergence-free transport, the intrinsic SUPG formulation adds streamline-directed test-function weighting, with stabilisation parameter τ, to the continuous finite-element formulation.The implementation obtains intrinsic velocity through the contravariant mapping and integrates the semi-discrete system with a Runge–Kutta method.Source: §4.1.1 · p. 9, §4.1 · p. 9
4.1.2. Upwinding stabilisation.
- For divergent velocity fields, intrinsic upwinding uses edge conormals and a manifold edge measure to formulate a discontinuous finite-element stabilisation.Across atlas chart interfaces, neighbouring parametric normals need not be oppositely oriented, so the weak formulation uses scalar differences rather than normal-based jump and average operators.Source: §4.1.2 · p. 10
4.1.3. Results.
- Both intrinsic schemes show expected hp-convergence and solutions equivalent to their extrinsic counterparts; the cosine-bell test also recovers its initial condition at t = 5 with mass error decreasing under refinement.For Figure 3, h and the time step reduce simultaneously; the cosine-bell test uses intrinsic SUPG with continuous p = 1 finite elements and a Crank–Nicolson integrator.Source: §4.1.3 · pp. 10–11, Figure 3 · p. 11
4.2. Thermal shallow water equations.
The section formulates thermal shallow water equations intrinsically with upwinded numerical fluxes, then tests the method on a modified Galewsky flow. At day 6, the jet has the correct wave number, comparable to an accepted solution, on a coarser mesh than in cited studies.
4.2. Thermal shallow water equations.
- The prognostic thermal shallow water system evolves velocity, depth, and density-weighted buoyancy, and the paper reformulates a prior compatible finite-element discretisation in an intrinsic framework.The earlier discretisation included upwinded fluxes; the intrinsic formulation follows the systematic derivations in Tambyah et al. [17].Source: §4.2 · pp. 11–12
4.2.1. Semi-discrete system.
- The semi-discrete system separates advection and stabilisation forms, with an upwinding parameter χ(wh) = |wh · n+|/2 for numerical fluxes.The paper also gives a weak representation of vorticity; time integration uses SSPRK3 and mass systems use conjugate gradients.Source: §4.2.1 · pp. 12–13
4.2.2. Results.
- The modified Galewsky test probes nonlinear atmospheric flow whose day-6 solution is sensitive to spurious vorticity waves and grid imprinting.The cited context contrasts the accepted wave number 6 solution with a wave number 4 grid-imprinting pattern.Source: §4.2.2 · p. 13
- On the ℓ = 6, p = 1 cubed-sphere mesh, the modified Galewsky jet reaches the correct wave number by day 6, comparable to an accepted solution and at coarser resolution than in cited studies.The authors attribute the result to the intrinsic framework’s exact, smooth geometry, which is free of geometric consistency errors.Source: §4.2.2 · p. 13, Figure 5 · p. 14
4.3. Linearised Boussinesq equations.
The section formulates linearised Boussinesq dynamics on a three-dimensional atmospheric shell and tests whether the discretisation captures vertical buoyancy and surface pressure waves.
4.3. Linearised Boussinesq equations.
- The linearised system couples fluid velocity, pressure, and buoyancy evolution on an atmospheric shell, with zero radial velocity at the shell boundary.The model variables and radial-boundary condition are specified in the governing equations.Source: §4.3 · p. 14
4.3.1. Semi-discrete system.
- The semi-discrete system defines the outward radial direction from the metric and uses Raviart–Thomas V1 to enforce zero flux at radial boundaries.The radial unit vector is orthogonal to spherical shells at constant γ, and its normalisation uses only the metric.Source: §4.3.1 · p. 15
4.3.2. Results.
- The gravity-wave test uses p = 1 finite elements, four refinement levels, and Crank–Nicolson integration with parameters chosen to match a reduced-Earth grid aspect ratio.The test assesses the semi-discrete formulation’s ability to capture surface and vertical dynamics.Source: §4.3.2 · p. 15
- Buoyancy is zero on the inner and outer shells and evolves vertically between them, while pressure propagates over the surface; both waves oscillate with time.The snapshots show the solution in an x-z cross-section of the three-dimensional manifold.Source: Figure 6 · p. 16, §4.3.2 · p. 16
5. Conclusion and future work
The conclusion presents GridapGeosciences.jl as a general-manifold finite element framework supporting intrinsic and extrinsic formulations. Its demonstrated scope uses conforming meshes and space-only maps, while several extensions remain future work.
5. Conclusion and future work
- GridapGeosciences.jl generates discrete manifolds from coarse atlas information and mesh refinement, supporting intrinsic and extrinsic finite element formulations.The conclusion identifies this dual capability as a novelty relative to studies limited to extrinsic functionality.Source: §5 · p. 16
- High-resolution cubed-sphere simulations in two and three dimensions demonstrate intrinsic and extrinsic frameworks resolving complex flows over long periods.The conclusion describes the library as leveraging Julia JIT compilation, multiple dispatch, and lazy evaluation for computational efficiency.Source: §5 · p. 16
- The study uses conforming meshes at a single refinement level and space-only maps; local adaptation and evolving-manifold methods remain extensions.Non-conforming refinement is noted as advantageous for localized features, while evolving manifolds would require space-time maps and potentially space-time finite element methods.Source: §5 · pp. 16–17
Declaration of competing interest.
The declaration reports no additional data and identifies GridapGeosciences.jl as a free, open-source library.
Declaration of competing interest.
- The article contains no additional data, and GridapGeosciences.jl is a free, open-source library.Source: Data availability · p. 17