Operators: spectral element

grad = Gradient()              # strong form, == Gradient{StrongForm}()
wgrad = Gradient{WeakForm}()   # weak form
grad.(f)

Compute the gradient of f on each element, returning a CovariantVector-field. Gradient() computes the strong form; Gradient{WeakForm}() computes the weak form (the volume-integral contribution obtained after integration by parts), distinguished by the FormType parameter.

Strong form

The $i$th covariant component of the gradient is the partial derivative with respect to the reference element:

\[(\nabla f)_i = \frac{\partial f}{\partial \xi^i}\]

Discretely, this can be written in matrix form as

\[D_i f\]

where $D_i$ is the derivative matrix along the $i$th dimension.

Weak form

The weak gradient is the vector field $\theta \in V_0$ such that for all $\phi \in V_0$

\[\int_\Omega \phi \cdot \theta \, d \Omega = - \int_\Omega (\nabla \cdot \phi) f \, d\Omega\]

where $V_0$ is the space of $f$. It arises from the contribution of the volume integral after applying integration by parts to the weak form expression of the gradient

\[\int_\Omega \phi \cdot (\nabla f) \, d \Omega = - \int_\Omega f (\nabla \cdot \phi) \, d\Omega + \oint_{\partial \Omega} f (\phi \cdot n) \, d \sigma\]

In matrix form, this becomes

\[{\phi^i}^\top W J \theta_i = - ( J^{-1} D_i J \phi^i )^\top W J f\]

which reduces to

\[\theta_i = -W^{-1} D_i^\top W f\]

References

• [22], equation 16.

div = Divergence()              # strong form, == Divergence{StrongForm}()
wdiv = Divergence{WeakForm}()   # weak form
div.(u)

Compute the per-element spectral divergence of a vector field $u$. Divergence() computes the strong form; Divergence{WeakForm}() computes the weak form (the volume-integral contribution obtained after integration by parts), distinguished by the FormType parameter.

Strong form

The divergence of a vector field $u$ is defined as

\[\nabla \cdot u = \sum_i \frac{1}{J} \frac{\partial (J u^i)}{\partial \xi^i}\]

where $J$ is the Jacobian determinant, $u^i$ is the $i$th contravariant component of $u$.

This is discretized by

\[\sum_i I \left\{\frac{1}{J} \frac{\partial (I\{J u^i\})}{\partial \xi^i} \right\}\]

where $I\{x\}$ is the interpolation operator that projects to the unique polynomial interpolating $x$ at the quadrature points. In matrix form, this can be written as

\[J^{-1} \sum_i D_i J u^i\]

where $D_i$ is the derivative matrix along the $i$th dimension.

Weak form

The weak divergence is the scalar field $\theta \in V_0$ such that for all $\phi\in V_0$

\[\int_\Omega \phi \theta \, d \Omega = - \int_\Omega (\nabla \phi) \cdot u \,d \Omega\]

where $V_0$ is the space of $u$. It arises as the contribution of the volume integral after applying integration by parts to the weak form expression of the divergence

\[\int_\Omega \phi (\nabla \cdot u) \, d \Omega = - \int_\Omega (\nabla \phi) \cdot u \,d \Omega + \oint_{\partial \Omega} \phi (u \cdot n) \,d \sigma\]

It can be written in matrix form as

\[ϕ^\top WJ θ = - \sum_i (D_i ϕ)^\top WJ u^i\]

which reduces to

\[θ = -(WJ)^{-1} \sum_i D_i^\top WJ u^i\]

where $W$ is the diagonal matrix of quadrature weights.

References

• [22], equation 15.

split_div = SplitDivergence()
split_div.(ρu, ψ)

Compute the divergence of the product ρu * ψ using a split-form (entropy-stable) discretization.

This operator is designed for the advection of scalar quantities in conservation laws (e.g., thermodynamic variables or tracers). By evaluating the divergence using an average of the conservative and advective forms, this formulation cancels aliasing errors that arise from the product of two spectrally variable fields, which inhibits the growth of quadratic instabilities (such as cold temperature spikes) without requiring hyperviscosity.

Arguments

• ρu: The transport vector field, typically the mass flux. It must be a vector quantity (e.g., Geometry.Contravariant12Vector).
• ψ: The specific scalar quantity to be advected (e.g., specific total energy $e_{tot}$ or specific humidity $q_{tot}$).

Mathematical Formulation

Continuous

The split form of the divergence operator is defined as the arithmetic mean of the conservative and advective forms:

\[\nabla \cdot (\rho \mathbf{u} \psi)|_\textrm{split} = \frac{1}{2} \nabla \cdot (\rho \mathbf{u} \psi) + \frac{1}{2} \left( \psi \nabla \cdot (\rho \mathbf{u}) + \rho \mathbf{u} \cdot \nabla \psi \right)\]

Discrete

The discretized split operator is equivalent to using the strong formulation of the gradient operator and the weak formulation of the divergence operator:

\[\textrm{split_div}(\rho \mathbf{u}, \psi) = \frac{1}{2} \textrm{wdiv}(\rho \mathbf{u} \psi) + \frac{1}{2} \left( \psi \textrm{wdiv}(\rho \mathbf{u}) + \rho \mathbf{u} \cdot \textrm{grad}(\psi) \right)\]

Swapping the weak and strong formulations in the last two terms also results in the same operator. The discrete form of the divergence theorem, which stems from the generalized summation-by-parts (SBP) property, guarantees that the integral of the first term vanishes,

\[\int_\Omega \textrm{wdiv}(\rho \mathbf{u} \psi) dV = 0\]

while the integrals of the other two terms cancel out,

\[\int_\Omega \psi \textrm{wdiv}(\rho \mathbf{u}) dV = -\int_\Omega \rho \mathbf{u} \cdot \textrm{grad}(\psi) dV\]

So, this discretization ensures that the split operator conserves the integral of $\rho \mathbf{u} \psi$.

Two-Point

A more compact representation of the discretized operator can be obtained with the symmetric two-point flux, whose values in one dimension are

\[(F^1)_{ij} = \frac{1}{2} (\rho_i J_i (u^1)_i + \rho_j J_j (u^1)_j) (\psi_i + \psi_j)\]

With $D$ denoting the spectral derivative matrix, the split operator in one dimension can be expressed as

\[\textrm{split_div}(\rho \mathbf{u}, \psi)_i = \frac{1}{J_i} \sum_{j \neq i} D_{ij} (F^1)_{ij}\]

In two dimensions, $F^1$ and the analogous quantity $F^2$ provide a similar expression for the split divergence, with the one-dimensional operator applied sequentially along each dimension.

Properties

1. Conservation: The split operator conserves $\rho \mathbf{u} \psi$.
2. Consistency: If $\psi = 1$, the split operator degenerates to the weak formulation of $\nabla \cdot \rho \mathbf{u}$ (mass continuity).
3. Complexity: The split operator has the same $O(N^2)$ complexity per element as the strong and weak operators, but needs twice as many operations.

References

• Fisher, T. C., & Carpenter, M. H. (2013). High-order entropy stable finite difference schemes for nonlinear conservation laws: Finite domains. Journal of Computational Physics, 252, 518-557. https://doi.org/10.1016/j.jcp.2013.06.014
• Gassner, G. J. (2013). A skew-symmetric discontinuous Galerkin spectral element discretization and its relation to SBP-SAT finite difference methods. SIAM Journal on Scientific Computing, 35, A1233-A1253. https://doi.org/10.1137/120890144

curl = Curl()              # strong form, == Curl{StrongForm}()
wcurl = Curl{WeakForm}()   # weak form
curl.(u)

Compute the per-element spectral curl of a covariant vector field $u$. Curl() computes the strong form; Curl{WeakForm}() computes the weak form (the volume-integral contribution obtained after integration by parts), distinguished by the FormType parameter.

The vector field $u$ must be explicitly converted to a CovariantVector, as then the curl is independent of the local metric tensor.

Strong form

The curl of a vector field $u$ is a vector field with contravariant components

\[(\nabla \times u)^i = \frac{1}{J} \sum_{jk} \epsilon^{ijk} \frac{\partial u_k}{\partial \xi^j}\]

where $J$ is the Jacobian determinant, $u_k$ is the $k$th covariant component of $u$, and $\epsilon^{ijk}$ are the Levi-Civita symbols. In other words

\[\begin{bmatrix} (\nabla \times u)^1 \\ (\nabla \times u)^2 \\ (\nabla \times u)^3 \end{bmatrix} = \frac{1}{J} \begin{bmatrix} \frac{\partial u_3}{\partial \xi^2} - \frac{\partial u_2}{\partial \xi^3} \\ \frac{\partial u_1}{\partial \xi^3} - \frac{\partial u_3}{\partial \xi^1} \\ \frac{\partial u_2}{\partial \xi^1} - \frac{\partial u_1}{\partial \xi^2} \end{bmatrix}\]

In matrix form, this becomes

\[\epsilon^{ijk} J^{-1} D_j u_k\]

Unused dimensions are dropped: e.g. the 2D curl of a Covariant12Vector-field returns a Contravariant3Vector.

Weak form

The weak curl is the vector field $\theta \in V_0$ such that for all $\phi \in V_0$

\[\int_\Omega \phi \cdot \theta \, d \Omega = \int_\Omega (\nabla \times \phi) \cdot u \,d \Omega\]

where $V_0$ is the space of $u$. It arises from the contribution of the volume integral after applying integration by parts to the weak form expression of the curl

\[\int_\Omega \phi \cdot (\nabla \times u) \,d\Omega = \int_\Omega (\nabla \times \phi) \cdot u \,d \Omega - \oint_{\partial \Omega} (\phi \times u) \cdot n \,d\sigma\]

In matrix form, this becomes

\[{\phi_i}^\top W J \theta^i = (J^{-1} \epsilon^{kji} D_j \phi_i)^\top W J u_k\]

which, by using the anti-symmetry of the Levi-Civita symbol, reduces to

\[\theta^i = - \epsilon^{ijk} (WJ)^{-1} D_j^\top W u_k\]

References

• [22], equation 17.

FormType

Supertype of the singleton types StrongForm and WeakForm, which distinguish the variational form of a spectral element operator.

The strong and weak variants of an operator share the same interior computation; they differ only in three form-dependent factors:

• whether the derivative matrix is applied directly or transposed with a sign flip (from integration by parts); see deriv_matrix,
• whether the argument is multiplied by the quadrature weights W or by one of the Jacobian factors J or WJ; see materialize_quadrature_weighted and materialize_jacobian_weighted,
• whether the result is divided by the quadrature weights or by one of the Jacobian factors; see quadrature_unweighted and jacobian_unweighted.

Operators that have strong/weak variants use a FormType parameter; e.g., Divergence() is a Divergence{StrongForm}, and the weak divergence is Divergence{WeakForm}().

StrongForm()

The FormType of an operator that discretizes a derivative directly at the quadrature points (e.g. Divergence, Gradient, and Curl).

WeakForm()

The FormType of an operator that discretizes the volume-integral contribution of the corresponding weak-form expression, obtained after integration by parts (e.g. Divergence{WeakForm}, Gradient{WeakForm}, and Curl{WeakForm}).

scalar_laplacian(χ; weight = nothing)

Return the weak horizontal Laplacian of the scalar χ, wdivₕ(gradₕ(χ)), or wdivₕ(weight * gradₕ(χ)) when a weight field (e.g., a density) is given.

On continuous (CG) spaces this returns a lazy operator expression, so it fuses into consuming broadcasts; the result is element-local, and materialized intermediates must be made continuous with Spaces.weighted_dss! before they are differentiated again (batch several intermediates into one call to share the ghost exchange). On discontinuous (DG) spaces, this returns a materialized Field that already includes the face terms coupling neighbouring elements (see sipg_laplacian_tendency!), and weighted_dss! is a no-op — so the same prep → weighted_dss! → apply sequence is correct for both discretizations.

Each call owns its result on both discretizations, so any number of results may be live at once. The DG result is a fresh Field; in a tendency, where that allocation is not wanted, use scalar_laplacian! to write into a caller-owned field instead.

Examples

Fourth-order hyperdiffusion ∇⁴χ on either discretization:

∇²χ = similar(χ)
Operators.scalar_laplacian!(∇²χ, χ)
Spaces.weighted_dss!(∇²χ)          # no-op on DG
@. dydt -= ν * Operators.scalar_laplacian(∇²χ)

scalar_laplacian!(out, χ; weight = nothing)

Write scalar_laplacian(χ; weight) into the caller-owned field out and return out. Allocation-free on both discretizations, so this is the form to use in a tendency. out may alias χ: the argument is copied into scratch first, since the spectral operators read a whole element while writing it. out may not alias weight.

vector_laplacian(u; divergence_factor = 1)

Return the weak horizontal vector Laplacian of the horizontal covariant vector u, as the grad-div minus curl-curl identity divergence_factor * wgradₕ(divₕ(u)) − C12(wcurlₕ(C3(curlₕ(u)))), with the grad-div part scaled by divergence_factor (used for divergence damping in the second pass of ∇⁴ hyperdiffusion). On spaces with one horizontal dimension the curl-curl part vanishes identically and only the grad-div part is built.

Returns a lazy operator expression on continuous (CG) spaces; the result is element-local and must be made continuous with Spaces.weighted_dss! before it is differentiated again. Throws an error on discontinuous (DG) spaces, which need grad-div/curl-curl face lifting that is not implemented.

i = Interpolate(space)
i.(f)

Interpolate f to space. If space has equal or higher polynomial degree than the space of f, this is exact; otherwise it is lossy.

In matrix form, it is the linear operator

\[I = \bigotimes_i I_i\]

where $I_i$ is the barycentric interpolation matrix in the $i$th dimension.

See also Restrict.

r = Restrict(space)
r.(f)

Compute the projection of a field f on $V_0$ to a lower degree polynomial space space ($V_0^*$). space must be on the same topology as the space of f, but have a lower polynomial degree.

It is defined as the field $\theta \in V_0^*$ such that for all $\phi \in V_0^*$

\[\int_\Omega \phi \theta \,d\Omega = \int_\Omega \phi f \,d\Omega\]

In matrix form, this is

\[\phi^\top W^* J^* \theta = (I \phi)^\top WJ f\]

where $W^*$ and $J^*$ are the quadrature weights and Jacobian determinant of $V_0^*$, and $I$ is the interpolation operator (see Interpolate) from $V_0^*$ to $V_0$. This reduces to

\[\theta = (W^* J^*)^{-1} I^\top WJ f\]

tensor_product!(out, in, M)
tensor_product!(inout, M)

Compute the tensor product out = (M ⊗ M) * in on each element (M * in for 1D elements) and return out; the two-argument form applies it in place to inout.