ClimaCore Spectral Element Extension

ConservativeRegridding.jl provides conservative regridding between finite volume (FV) grids and ClimaCore.jl spectral element (SE) spaces. The extension is loaded automatically when ClimaCore is imported alongside ConservativeRegridding.

Spectral element and finite volume integrals

A finite volume grid carries one mean value $\bar f_i$ per cell, so the integral of a field over cell $i$ is

\[\int_i f \, \mathrm{d}A = \bar f_i \, A_i ,\]

where $A_i$ is the cell area.

A spectral element space carries $N_q^2$ nodal values $f^e_{ij}$ per element $e$ (polynomial degree $N_q - 1$), expanded in a tensor product of one-dimensional Lagrange basis polynomials,

\[f^e(\xi, \eta) = \sum_{i,j} f^e_{ij} \, \phi_i(\xi) \, \phi_j(\eta) ,\]

over the reference square $(\xi, \eta) \in [-1, 1]^2$. The integral over an element is evaluated by Gauss–Legendre–Lobatto (GLL) quadrature,

\[\int_e f \, \mathrm{d}A \approx \sum_{i,j} W^e_{ij} \, f^e_{ij} , \qquad W^e_{ij} = w_i \, w_j \, J^e_{ij} ,\]

where $w_i$ are the quadrature weights, $J^e_{ij}$ the Jacobian of the reference-to-physical map at node $(i,j)$, and $W^e_{ij}$ the weighted Jacobian stored by ClimaCore.Spaces.weighted_jacobian.

Implemented regridding directions

Three directions are conceivable; the extension implements the two that involve a finite volume side:

Direction	Method
SE $\to$ FV	area-weighted average of the SE field over each cell (see SE → FV below)
FV $\to$ SE	per-element $L^2$ projection onto the SE basis (see FV → SE below)
SE $\to$ SE	not supported — compose SE $\to$ FV $\to$ SE

Both implemented directions are built from the same family of weights $B(k,(e,i,j))$, defined and evaluated below.

The intersection weights $B$

For a finite volume cell $k$ and a spectral element $e$, define

\[B(k, (e,i,j)) = \int_{k \cap e} \phi_i(\xi) \, \phi_j(\eta) \, J^e(\xi, \eta) \, \mathrm{d}\xi \, \mathrm{d}\eta ,\]

the integral of the SE basis function $(i,j)$ over the intersection of cell $k$ with element $e$. Because the reference-to-physical map has Jacobian determinant $J^e$, the physical area element is $\mathrm{d}A = J^e \, \mathrm{d}\xi \, \mathrm{d}\eta$, and the weight is equivalently a physical-space integral with no explicit Jacobian,

\[B(k, (e,i,j)) = \int_{k \cap e} \phi_i\big(\xi(x)\big) \, \phi_j\big(\eta(x)\big) \, \mathrm{d}A .\]

The extension evaluates this physical-space form directly (accumulate_principled_b) by (1) fan-triangulating the intersection polygon $k \cap e$ (already computed by ConservativeRegridding.jl via GeometryOps.jl) from its centroid, (2) applying a barycentric Gauss rule on each triangle, and (3) evaluating the basis functions at each quadrature point after inverting the element map to recover $(\xi, \eta)$. On the sphere, the quadrature points are mapped onto each spherical triangle by the flat triangle followed by the radial projection $p \mapsto p/|p|$, whose area density is proportional to $1/|p|^3$; each point is weighted by this density, and the weights are normalized to the triangle's area on the manifold's sphere, so constants are integrated exactly. The physical triangle areas carry the Jacobian, so $J^e$ does not appear explicitly. The integrand is smooth on each triangle, so the integration converges with the order of the rule and is high-order accurate in the element size.

Two identities follow from the definition and are used below. Where the cells $k$ cover element $e$,

\[\sum_k B(k, (e,i,j)) = \int_e \phi_i \, \phi_j \, \mathrm{d}A \approx W^e_{ij} ,\]

and since the elements $e$ tile the domain (and $\sum_{i,j} \phi_i \phi_j = 1$),

\[\sum_{e,i,j} B(k, (e,i,j)) = A_k .\]

SE → FV

The destination cell value is the area-weighted integral of the source SE field over the cell,

\[f^{\mathrm{dst}}_k = \frac{1}{A_k} \sum_{e,i,j} B(k,(e,i,j)) \, f^{\mathrm{src},e}_{ij} .\]

This is assembled as a sparse matrix with entries $B(k,(e,i,j))$ (rows indexed by FV cell, columns by SE node), and regrid! performs the matrix–vector product followed by division by the destination areas. Conservation is exact in the continuous weights: summing over $k$ and using $\sum_k B = \int_e \phi_i \phi_j \, \mathrm{d}A$ recovers $\int f^{\mathrm{src}} \, \mathrm{d}A$, the integral of the SE expansion of the source, which the GLL sum $\sum_{e,i,j} W^e_{ij} f^{\mathrm{src},e}_{ij}$ approximates.

FV → SE

The destination SE field is defined by the weak condition that it agree with the source in every test direction of the SE space: find $f^{\mathrm{dst}}$ such that

\[\int f^{\mathrm{dst}} \, \psi \, \mathrm{d}A = \int f^{\mathrm{src}} \, \psi \, \mathrm{d}A\]

for all basis functions $\psi$ of the SE space. Taking $\psi = 1$ makes this conservative. Expanding $f^{\mathrm{dst}}$ in the SE basis turns the condition into a linear system $\mathbf{M} \, \mathbf{f}^{\mathrm{dst}} = \mathbf{b}$ with the mass matrix $\mathbf{M}_{\alpha\beta} = \int \psi_\alpha \psi_\beta \, \mathrm{d}A$. Both sides are block-diagonal over elements, so the problem reduces to a per-element solve. The local mass matrix is

\[M^e_{(ab),(cd)} = \int_e \phi_a(\xi) \, \phi_b(\eta) \, \phi_c(\xi) \, \phi_d(\eta) \, J^e(\xi, \eta) \, \mathrm{d}\xi \, \mathrm{d}\eta\]

(assembled by compute_local_mass_matrix, with the index pairs flattened by LinearIndices). Its integrand has degree $3(N_q - 1)$ per direction, so it is integrated with an $n$-point GLL rule, which is exact to degree $2n - 3$, using the smallest sufficient $n = \lceil 3N_q/2 \rceil$; the Jacobian $J^e$ is interpolated from its nodal values where needed (element_jacobian_at). The right-hand side reuses the intersection weights,

\[b^e_{ij} = \sum_k B(k,(e,i,j)) \, f^{\mathrm{src}}_k ,\]

and the per-element solve gives the assembled regridding operator,

\[f^{\mathrm{dst},e}_{ij} = \sum_{a,b} \big[(M^e)^{-1}\big]_{(ij),(ab)} \, \sum_k B(k,(e,a,b)) \, f^{\mathrm{src}}_k .\]

Using the full mass-matrix inverse — rather than dividing by the lumped diagonal $W^e_{ij}$ — preserves constants exactly on nodes all of whose elements the source covers and is the optimal $L^2$ fit of the source on the SE basis over each element.

Two corrections appear in the implementation:

Quadrature-consistency rescaling. The mass matrix is integrated by tensor-product GLL on the reference square with the interpolated Jacobian, while $B$ is integrated by triangulation of the physical polygons, so $M^e \mathbf{1}$ and the column sums of $B$ differ slightly. On elements the source covers, each row of $M^e$ is rescaled before solving so that its row sum matches the corresponding column sum of $B$, which makes constants exact.
Partial coverage. Where the source covers an element only partly (for example a tripolar ocean grid that stops short of the South Pole, regridded onto a full-sphere SE space), the uncovered part counts as zero, as in the FV $\to$ FV and SE $\to$ FV transfers. $b^e$ integrates only the covered part, and the unscaled mass matrix gives the $L^2$ projection of the source extended by zero. This is conservative and well conditioned; near the edge of the source grid the result falls off toward zero (though over- and undershoots of an $L^2$ projection of a discontinuous field are to be expected). Elements the source does not reach are set to zero.

SE → SE

Direct SE $\to$ SE regridding is not supported: the Regridder constructor raises an error for two SE arguments and directs the caller to compose the two implemented directions, SE $\to$ FV $\to$ SE, with an intermediate finite volume grid.

Usage

Constructing a regridder

Regridder detects which argument is a ClimaCore.Spaces.AbstractSpectralElementSpace (or a ClimaCore.Fields.Field) and selects the corresponding direction. The first argument is the destination, the second the source.

Note that Regridder, regrid!, and areas are marked public but are not exported, so using ConservativeRegridding does not bring them into scope; import them explicitly (as below) or qualify them as ConservativeRegridding.Regridder, etc.

using ConservativeRegridding, ClimaCore, Oceananigans
using ConservativeRegridding: Regridder, regrid!
using ClimaCore: CommonSpaces

space = CommonSpaces.CubedSphereSpace(; radius = 6.371e6, n_quad_points = 4, h_elem = 16)
grid  = LatitudeLongitudeGrid(size = (360, 180, 1),
                              longitude = (0, 360), latitude = (-90, 90), z = (0, 1))

R_se2fv = Regridder(grid, space)   # SE source → FV destination
R_fv2se = Regridder(space, grid)   # FV source → SE destination

Regridding

regrid! accepts and writes ClimaCore Fields directly on the SE side, flattening to and from the nodal layout internally (and applying a weighted direct stiffness summation when writing an SE destination):

using ClimaCore: Fields

# SE → FV
src_field = Fields.coordinate_field(space).lat
fv_vals   = zeros(360 * 180)
regrid!(fv_vals, R_se2fv, src_field)

# FV → SE
dst_field = Fields.zeros(space)
regrid!(dst_field, R_fv2se, fv_vals)

Converting between ClimaCore fields and flat vectors

const ClimaCoreExt = Base.get_extension(ConservativeRegridding, :ConservativeRegriddingClimaCoreExt)

field = Fields.coordinate_field(space).lat           # an example SE field

src_vec = ClimaCoreExt.se_field_to_vec(field)        # Field → flat vector

dst_field = Fields.zeros(space)
ClimaCoreExt.vec_to_se_field!(dst_field, src_vec)    # flat vector → Field

positions = ClimaCoreExt.se_node_positions(space)    # Vector{UnitSphericalPoint}
weights   = ClimaCoreExt.se_node_weights(space)      # nodal weights W^e_{ij}

Checking conservation

SE → FV conserves the integral of the SE expansion of the source, $\sum_{e,i,j} \big(\sum_k B(k,(e,i,j))\big) f^{\mathrm{src},e}_{ij}$, exactly. ClimaCore.sum(field) computes the GLL sum $\sum_{e,i,j} W^e_{ij} f^e_{ij}$, which approximates it to the GLL quadrature error. A check takes both sides from the regridder, whose areas are on its manifold's sphere:

lat = Fields.coordinate_field(space).lat
src_field = @. 2 + sind(lat)
regrid!(fv_vals, R_se2fv, src_field)
b = vec(sum(R_se2fv.intersections; dims = 1))  # Σ_k B(k,(e,i,j)) = ∫_e ϕ_i ϕ_j dA
@assert sum(fv_vals .* R_se2fv.dst_areas) ≈ sum(b .* ClimaCoreExt.se_field_to_vec(src_field))