Two arbitrary-order constraint-preserving schemes for the Yang–Mills equations on polyhedral meshes

1.   Introduction

In this paper we investigate two arbitrary-order numerical methods for the Yang–Mills equations on general polyhedral meshes, based on the fully discrete serendipity Discrete de Rham (SDDR) complex ^[1]. The first method was proposed (but not tested) in ^[2] for the non-serendipity version of the Discrete de Rham (DDR) sequence, while the second method is novel to this paper. The two discretisations differ in the treatment of one of the nonlinearities present in these equations. In contrast to conforming methods, the discrete structure of the SDDR spaces means that there is no obvious construction of the nonlinear terms, and this can be problematic when specific algebraic manipulations need to be reproduced, for instance to prove consistency estimates. The implementation cost is another important factor to the viability of each approach, which is explored with accompanying numerical results on the convergence and discrete conservation properties of each scheme.

The classical Yang–Mills equations come from a class of non-abelian gauge theories, generalising the abelian $U(1)$ group of electromagnetism to certain non-abelian gauge groups. Once quantised, this theory forms the foundation of the current Standard Model of particle physics. In the classical setting, the non-commutativity of the group manifests as the appearance of nonlinear quantities in addition to the linear Maxwell terms. Analogous to Maxwell, the Yang–Mills equations can be formulated as a set of evolution equations preserving particular constraints (e.g., the conservation of charge), given that the initial data satisfies these constraints. In the linear case, the preservation of these constraints is a direct consequence of the calculus formula div curl = 0, which is linked to the complex property of the de Rham sequence. Designing numerical methods that replicate this property is essential to maintaining constraint preservation at the discrete level, and thus to obtaining stable schemes. Much work has been done in the Finite Element framework to design discrete versions of the de Rham complex, see, e.g., ^{[3,4,5,6,7,8]} and references therein. Finite Element methods are, however, limited to meshes made of specific elements (mostly tetrahedra and hexahedra in 3D), which limits their flexibility in terms of mesh refinement or agglomeration. Recently, discrete polytopal complexes – discrete versions of continuous complexes, that are applicable on meshes made of generic polyhedra – have been introduced, see, e.g., ^[9,10,11,12]. The discrete complex property enabled the design of stable and robust schemes, in particular for magnetostatics ^[10,13], plate problems ^{[14,15,16,17]}, and the Stokes equations ^[18,19]).

Given the importance, for the stability of schemes, of preserving constraints at the discrete level, similar techniques have been explored for the Yang–Mills equations, using either Finite Element or polytopal approaches ^[2,20,21]. For these equations, however, the nonlinearity has proven to be troublesome, and required additional techniques (e.g., the introduction of Lagrange multipliers) beyond a discrete version of the formula div curl = 0. The interest in developing our understanding of such methods is in the application to numerical schemes for Einstein's equations, where the absence of this constraint propagation can cause disastrous error growth ^[22,23,24] in the numerical simulations. Current techniques to control this error include constraint damping ^[22,23,25], where specific terms are added to the evolution equations to suppress the growth of the constraint violations, but methods for exact preservation remain limited. The link with the Yang–Mills equations is that in certain formulations of General Relativity (GR), such as the Einstein-Bianchi system ^[26,27], these equations can resemble greatly those of electromagnetism with additional nonlinear terms. Therefore it is natural to expect that these ideas will aid in designing a constraint preserving scheme for GR based on the framework of discrete polytopal complexes.

An equally important aspect of the design of numerical methods is the feasibility of the implementation and testing under real world conditions. We find more commonly, for the Yang–Mills equations, numerical tests run in only low-order 2D settings ^[20,21]. Any increase in the dimension or the order of the approximation generally leads to schemes that are vastly more expensive to run, and this is compounded by the nonlinearity of the model. Hence working with spaces that are smaller and more refined is an effective way to cut the cost of the simulations. The SDDR complex, introduced in ^[1], is a variant of the DDR complex ^[11,12], where the spaces have undergone a serendipity reduction, eliminating many unknowns, while retaining the complex and consistency properties of the original sequence. This enables the seamless transfer of any DDR scheme and results to the SDDR version, with all the flexibility of the general-order polytopal method at a lower cost. Additionally, this can be combined with other reduction techniques such as static condensation to further increase the efficiency.

The issue with the nonlinearity in the practical implementation is the computations involving 'high dimensional' arrays that are required to deal with all the coefficients. This number grows exponentially with the degree of the multilinearity, and thus takes up majority of the time in the assembly phase of the runtime. For matrices (2-dimensional arrays), there exists many specialised algorithms to speed up calculations, as well as efficient storage structures in the case that it is sparse. The libraries for higher dimensional arrays are less advanced, and often incomplete in their features; as a consequence, operations need to be done manually, introducing another source of possible inefficiencies. Simply rearranging the order of calculations can lead to sizeable differences in the space and time complexities, therefore finding the optimal trade-off is key to measuring the actual performance of the scheme.

The paper is organised as follows: In Section 2, we give a presentation of the SDDR complex and its Lie algebra extension, that is independent of the DDR framework. Section 3 starts with the constrained formulation of the continuous Yang–Mills equations. Based on that, we introduce the two schemes that are considered in the paper and the differing approaches on the nonlinear terms. This is followed by a proof of the preservation of a discrete constraint functional, as well as energy estimates. Section 4 covers the major steps in the implementation of the scheme, showing the impact of Lie algebra tensorisation on the physical data structures, and also highlighting how the trilinear and quadrilinear forms and sums are managed in the tensorised situation. Numerical results for these implementations are found in Section 5, where we test both the convergence and the discrete constraint preservation on three different mesh families in the 3D setting. Expected convergence rates of $k+1$ are mostly seen, as well as the preservation of the initial constraint up to machine precision. We also report on the differences in the results and runtimes, which turned out to be very minor between the two methods. A brief conclusion is provided in Section 6.

2.   Lie Algebra-valued serendipity Discrete de Rham complex

We present here the serendipity version of the arbitrary-order Lie Algebra-valued DDR complex, originally sketched in [2, Section 6]. This complex consists in tensorising the SDDR complex of ^[1], which is built in this reference from the (regular) DDR complex; such a presentation relies on a complete description of the latter complex, together with "extension" and "reduction" maps that link the two complexes. In the following, we adopt a stand-alone description of the SDDR complex, directly translating the formulas resulting from the links with the DDR complex. For this reason, the notations adopted below differ slightly from ^[1]: there, the serendipity spaces and operators are denoted using a hat (the non-hat version referring to the regular DDR spaces and operators, which are not needed here).

2.1.   Mesh notations

We use the same mesh and polynomial space notations as in ^[12]. Let $U$ be a polygonal domain of $\mathbb{R}^3$. A mesh $\mathcal{M}_{h} = \mathcal{T}_{h}\cup\mathcal{F}_{h}\cup\mathcal{E}_{h}\cup\mathcal{V}_{h}$ is a collection of polyhedral elements (gathered in $\mathcal{T}_{h}$, and partitioning $U$), of polygonal faces (gathered in $\mathcal{F}_{h}$), of edges (gathered in $\mathcal{E}_{h}$) and vertices (gathered in $\mathcal{V}_{h}$). Each ${\mathsf{P}}\in \mathcal{M}_{h}$ is assumed to be topologically trivial (simply connected with connected boundary), and we denote by $h_{ {\mathsf{P}}}$ the diameter of ${\mathsf{P}}$; we set $h = \max_{T\in\mathcal{T}_{h}}h_T$. When applicable, the sets $\mathcal{F}_{ {\mathsf{P}}}$ (resp. $\mathcal{E}_{ {\mathsf{P}}}$, resp. $\mathcal{V}_{ {\mathsf{P}}}$) gather the faces (resp. edges, resp. vertices) of ${\mathsf{P}}$. Each face $F\in\mathcal{F}_{h}$ is oriented by the choice of a unit normal $\boldsymbol{n}_F$, and each edge $E\in\mathcal{E}_{h}$ is oriented by the choice of a unit tangent $\boldsymbol{t}_E$. If $T\in\mathcal{T}_{h}$ and $F\in \mathcal{F}_{{T}}$, $\omega_{TF}$ is the relative orientation of $F$ with respect to $T$: $\omega_{TF} = +1$ if $\boldsymbol{n}_F$ points outside $T$, $\omega_{TF} = -1$ otherwise. For $F\in\mathcal{F}_{h}$ and $E\in \mathcal{E}_{{F}}$, $\omega_{FE}$ denotes the relative orientation of $E$ with respect to $F$: $\omega_{FE} = +1$ if, along $\partial F$, $\boldsymbol{t}_E$ points counter-clockwise with respect to the orientation of $F$ induced by $\boldsymbol{n}_F$, and $\omega_{FE} = -1$ otherwise; we also denote by $\boldsymbol{n}_{FE}$ the unit vector such that $(\boldsymbol{t}_E, \boldsymbol{n}_{FE}, \boldsymbol{n}_F)$ defines a right-handed system in $\mathbb{R}^3$. The final orientation is that of the vertices of each edge: for $E\in\mathcal{E}_{h}$ and $V\in \mathcal{V}_{{E}}$, $\omega_{EV} = +1$ if $\boldsymbol{t}_E$ points towards $V$ on $E$, and $\omega_{EV} = -1$ otherwise.

We assume that $\mathcal{T}_{h}\cup\mathcal{F}_{h}$ satisfies the regularity assumption of [28, Definition 1.9] with regularity parameter $\varrho$, and we write $A\lesssim B$ when $A\le CB$ for some $C$ depending only on $U$, $\varrho$ and the possible polynomial degrees involved in $A, B$.

For any mesh face $F\in\mathcal{F}_{h}$ and smooth enough function $r:F\to \mathbb{R}$, ${\bf{grad}}_F r$ is the gradient of $r$ on $F$ and ${\bf{rot}}_F r$ its 2-dimensional vector curl (rotation of ${\bf{grad}}_F r$ by $-\pi/2$ in the plane spanned by $F$). For a smooth function ${z}$ on $F$ with values in the tangent plane of $F$, the divergence of ${z}$ on $F$ is ${\bf{div}}_F {z}$, and its scalar curl (divergence of the rotated by $-\pi/2$ of ${z}$) is ${\rm{rot}}_F{z}$.

If ${\mathsf{P}}\in\mathcal{M}_{h}$ and $\ell\ge 0$ is an integer, $\mathcal{P}^{{\ell}}({\mathsf{P}})$ denotes the space of restrictions to ${\mathsf{P}}$ of three-variate polynomials on $\mathbb{R}^3$ of total degree $\le\ell$, and $\mathcal{P}^{{0, \ell}}({\mathsf{P}})$ is its subspace of polynomials with vanishing integral over ${\mathsf{P}}$. We adopt the convention $\mathcal{P}^{{\ell}}({\mathsf{P}}) = \{0\}$ if $\ell < 0$. If $T\in\mathcal{T}_{h}$, we set ${\mathcal{{P}}}^{{\ell}}(T) = \mathcal{P}^{{\ell}}(T)^3$ and, for $F\in\mathcal{F}_{h}$, ${\mathcal{{P}}}^{{\ell}}(F)$ is the subspace of $\mathcal{P}^{{\ell}}(F)^3$ which take value in the tangent space of $F$. The $L^2$-orthogonal projector on $\mathcal{P}^{{\ell}}({\mathsf{P}})$ is denoted by $\pi_{\mathcal{P}, { {\mathsf{P}}}}^{{\ell}}$. Selecting, for each ${\mathsf{P}}\in\mathcal{T}_{h}\cup\mathcal{F}_{h}$, a point ${x}_ {\mathsf{P}}\in {\mathsf{P}}$ such that ${\mathsf{P}}$ contains a ball centered at ${x}_ {\mathsf{P}}$ and of radius $\gtrsim h_ {\mathsf{P}}$, we recall the following decompositions of vector-valued polynomial spaces: For all $F\in\mathcal{F}_{h}$,

and, for $T\in\mathcal{T}_{h}$,

(here and in the following, when used between two vectors or a vector and a space, × denotes the cross product in $\mathbb{R}^3$). The $L^2$-orthogonal projectors on these spaces are, with obvious notations, ${\boldsymbol\pi}_{{\mathcal{{R}}}, {F}}^{{\ell}}$, ${\boldsymbol\pi}_{{\mathcal{{R}}}, {F}}^{ {{\rm{c}}}, {\ell}}$, ${\boldsymbol\pi}_{{\mathcal{{R}}}, {T}}^{{\ell}}$, ${\boldsymbol\pi}_{{\mathcal{{R}}}, {T}}^{ {{\rm{c}}}, {\ell}}$, ${\boldsymbol\pi}_{{\mathcal{{G}}}, {T}}^{{\ell}}$ and ${\boldsymbol\pi}_{{\mathcal{{G}}}, {T}}^{ {{\rm{c}}}, {\ell}}$.

2.2.   Serendipity DDR complex

For each element or face ${\mathsf{P}}\in\mathcal{T}_{h}\cup\mathcal{F}_{h}$, we select on the boundary of the element (resp. face) a set $\mathcal{B}_{{ {\mathsf{P}}}}$ of $\eta_ {\mathsf{P}}\ge 2$ faces (resp. edges) that are not pairwise coplanar (resp. aligned) and such that, for each ${\mathsf{b}}\in\mathcal{B}_{{ {\mathsf{P}}}}$, ${\mathsf{P}}$ lies entirely on one side of the affine space spanned by ${\mathsf{b}}$. From here on, we fix a polynomial degree $k\ge 0$, measuring the accuracy of the discrete complex, and we set

2.2.1.   Spaces and serendipity operators

The SDDR versions of the $H^1(U)$, ${\boldsymbol H}({\bf{curl}}; {U})$, ${\boldsymbol H}({\bf{div}}; {U})$ and $L^2(U)$ spaces appearing in the continuous de Rham complex are the following spaces.

The interpolators on these spaces consist in projecting continuous scalar/vector fields (or some of their traces) onto the polynomial components of the spaces. Specifically,

and

are defined as:

where $\boldsymbol{{v}}_{{\rm t}, F} = \boldsymbol{n}_F\times(\boldsymbol{{v}}_{|F}\times \boldsymbol{n}_F)$ is the tangential trace of $\boldsymbol{{v}}$ on $F$.

As usual in fully discrete complexes, we adopt the underlined notation for vectors of polynomial components, and we replace the index $h$ with ${\mathsf{P}}$ to denote the restriction of these spaces (and the operators defined on them) to a mesh entity ${\mathsf{P}}\in\mathcal{T}_{h}\cup\mathcal{F}_{h}\cup\mathcal{E}_{h}$ and its boundary entities. So, for example, a vector $\underline{\boldsymbol{{{v}}}}_F\in\underline{\boldsymbol{{{X}}}}_{{\bf{curl}}, {F}}^{{k}}$ corresponds to $\underline{\boldsymbol{{{v}}}}_F = (\boldsymbol{{v}}_{\boldsymbol{{\mathcal{{R}}}}, F}, \boldsymbol{{v}}_{\boldsymbol{{\mathcal{{R}}}}, F}^ {{\rm{c}}}, (v_E)_{E\in \mathcal{E}_{{F}}})$.

The DDR spaces correspond to the spaces above with the choice $\ell_F = \ell_T = k-1$ (that is, $\eta_F = \eta_P = 2$). This implies in particular that, for $k = 0$ (which forces $\ell_ {\mathsf{P}} < 0$ for all ${\mathsf{P}}\in\mathcal{T}_{h}\cup\mathcal{F}_{h}$), the standard and serendipity DDR spaces are identical. However, as soon as $k\ge 1$, the SDDR spaces have lower dimensions, while still encoding the same level of polynomial consistency as the DDR spaces. This is due to the existence of two families of key operators, the serendipity gradient and curl operators. Specifically, for ${\mathsf{P}}\in\mathcal{T}_{h}\cup\mathcal{F}_{h}$, the role of the gradient serendipity operator $\boldsymbol{{S}}_{{\bf{grad}}, { {\mathsf{P}}}}^k:\underline{X}_{{\bf{grad}}, { {\mathsf{P}}}}^{{k}}\to\boldsymbol{{\mathcal{{P}}}}^{{k}}({\mathsf{P}})$ is to reconstruct a consistent gradient, while the curl serendipity operator $\boldsymbol{{S}}_{{\bf{curl}}, { {\mathsf{P}}}}^k:\underline{\boldsymbol{{{X}}}}_{{\bf{curl}}, { {\mathsf{P}}}}^{{k}}\to\boldsymbol{{\mathcal{{P}}}}^{{k}}({\mathsf{P}})$ reconstructs a consistent vector potential. The consistencies in questions are expressed by the following relations (see [1, Proposition 18]):

We do not present the precise definitions of these operators, which are not essential to describe the SDDR complex, and refer the reader to ^[1].

In the next three sections we define operators acting on these spaces, with values in full polynomial spaces, mimicking the gradient, curl, and divergence. It should be noted that, in the original presentation of the SDDR complex in ^[1], these operators were not explicitly defined – only the discrete operators (projections on the complex spaces, see Section 2.2) and discrete inner products were detailed, based on those of the DDR complex. The polynomial operators below correspond to those of the DDR complex composed with the extension operators linking the DDR and SDDR complex; for ease of reference, we indicate which formulas from ^[1] yield the definitions presented here.

2.2.2.   Operators on the gradient space

For each edge $E\in\mathcal{E}_{h}$ we define the edge gradient $G_E^{{k}}:\underline{X}_{{\bf{grad}}, {E}}^{{k}}\to\mathcal{P}^{{k}}(E)$ and potential reconstruction $\gamma_E^{{k+1}}:\underline{X}_{{\bf{grad}}, {E}}^{{k}}\to\mathcal{P}^{{k+1}}(E)$ by: For all $\underline{q}_E = (q_E, (q_V)_{V\in \mathcal{V}_{{E}}})\in\underline{X}_{{\bf{grad}}, {E}}^{{k}}$,

The definition of $G_E^{{k}}\underline{q}_E$, in which the derivative $r_E'$ is taken in the direction $\boldsymbol{t}_E$, mimics an integration-by-parts formula; it can be checked that

For each $F\in\mathcal{F}_{h}$, combining [1, Eqs (4.2), (5.18) and (6.6)] together with ${\rm{div}}_F {\bf{rot}}_F = 0$ yields the following definition of the face gradient $\boldsymbol{\mathsf{G}}_F^{{k}}:\underline{X}_{{\bf{grad}}, {F}}^{{k}}\to\boldsymbol{{\mathcal{{P}}}}^{{k}}(F)$: For all $\underline{q}_F\in\underline{X}_{{\bf{grad}}, {F}}^{{k}}$,

Using this face gradient and [1, Eqs (4.3) and (5.18)], the scalar potential reconstruction on $F\in\mathcal{F}_{h}$ is then $\gamma_F^{{k+1}}:\underline{X}_{{\bf{grad}}, {F}}^{{k}}\to\mathcal{P}^{{k+1}}(F)$ defined by: For all $\underline{q}_F\in\underline{X}_{{\bf{grad}}, {F}}^{{k}}$,

Finally, for $T\in\mathcal{T}_{h}$, we use [1, Eqs (4.4),(5.32) and (6.6)] to write the element gradient $\boldsymbol{\mathsf{G}}_T^{{k}}:\underline{X}_{{\bf{grad}}, {T}}^{{k}}\to\boldsymbol{{\mathcal{{P}}}}^{{k}}(T)$ as: For all $\underline{q}_T\in\underline{X}_{{\bf{grad}}, {T}}^{{k}}$,

The potential reconstruction $P_{{\bf{grad}}, T}^{{k+1}}:\underline{X}_{{\bf{grad}}, {T}}^{{k}}\to\mathcal{P}^{{k+1}}(T)$ is such that: For all $\underline{q}_T\in\underline{X}_{{\bf{grad}}, {T}}^{{k}}$,

Remark 1 (Approximation properties of the potential reconstructions in $\underline{X}_{{\bf{grad}}, {T}}^{{k}}$). As demonstrated by [12, Theorem 6], the potential reconstructions on the space $\underline{X}_{{\bf{grad}}, {T}}^{{k}}$ have optimal approximation properties of degree $k+1$. This is however an exception to the rule of spaces and potential reconstructions in the DDR complex; the reasons for this exception are better understood when translating this complex in the language of differential forms (see ^[29], especially Remarks 7 and 18 therein).

2.2.3.   Operators on the curl space

For $F\in\mathcal{F}_{h}$, using [1, Eqs (4.6), (5.19) and (6.7)] we define the face curl $C_F^{{k}}:\underline{\boldsymbol{{{X}}}}_{{\bf{curl}}, {F}}^{{k}}\to\mathcal{P}^{{k}}(F)$ by: For all $\underline{\boldsymbol{{{v}}}}_F\in\underline{\boldsymbol{{{X}}}}_{{\bf{curl}}, {F}}^{{k}}$,

This definition is actually identical to the face curl in the DDR complex and does not invoke $\boldsymbol{{S}}_{{\bf{curl}}, {F}}^k$, as the extension operators between the SDDR and DDR curl face space do not modify the components on $\boldsymbol{{\mathcal{{R}}}}^{{k-1}}(F)$ and on $\times_{E\in \mathcal{E}_{{F}}}\mathcal{P}^{{k}}(E)$. The curl serendipity operator is however involved in the definition of the component on $\boldsymbol{{\mathcal{{R}}}}^{ {{\rm{c}}}, {k}}(F)$ of the extension (see [1, Eq (5.19)]), and therefore in the vector potential reconstruction $\boldsymbol{{\gamma}}_{{\rm t}, F}^{{k}}:\underline{\boldsymbol{{{X}}}}_{{\bf{curl}}, {F}}^{{k}}\to\boldsymbol{{\mathcal{{P}}}}^{{k}}(F)$ on $F$, defined as: For all $\underline{\boldsymbol{{{v}}}}_F\in\underline{\boldsymbol{{{X}}}}_{{\bf{curl}}, {F}}^{{k}}$,

Similar considerations apply to the element curl and vector potential. For $T\in\mathcal{T}_{h}$, the element curl $\boldsymbol{\mathsf{C}}_T^{{k}}:\underline{\boldsymbol{{{X}}}}_{{\bf{curl}}, {T}}^{{k}}\to\boldsymbol{{\mathcal{{P}}}}^{{k}}(T)$ is defined by: For all $\underline{\boldsymbol{{{v}}}}_T\in\underline{\boldsymbol{{{X}}}}_{{\bf{curl}}, {T}}^{{k}}$,

The vector potential $\boldsymbol{{P}}_{{\bf{curl}}, T}^{{k}}:\underline{\boldsymbol{{{X}}}}_{{\bf{curl}}, {T}}^{{k}}\to\boldsymbol{{\mathcal{{P}}}}^{{k}}(T)$ is given by: For all $\underline{\boldsymbol{{{v}}}}_T\in\underline{\boldsymbol{{{X}}}}_{{\bf{curl}}, {T}}^{{k}}$,

2.2.4.   Operators on the divergence space

The discrete divergence and potential on $\underline{\boldsymbol{{{X}}}}_{{\rm{div}}, {T}}^{{k}}$, for $T\in\mathcal{T}_{h}$, are identical to those of the DDR complex since no serendipity reduction is actually possible on this space (see [1, Section 6.5]): $D_T^{{k}}:\underline{\boldsymbol{{{X}}}}_{{\rm{div}}, {T}}^{{k}}\to\mathcal{P}^{{k}}(T)$ and $\boldsymbol{{P}}_{{\rm{div}}, T}^{{k}}:\underline{\boldsymbol{{{X}}}}_{{\rm{div}}, {T}}^{{k}}\to\boldsymbol{{\mathcal{{P}}}}^{{k}}(T)$ are such that, for all $\underline{\boldsymbol{{{w}}}}_T\in\underline{\boldsymbol{{{X}}}}_{{\rm{div}}, {T}}^{{k}}$,

2.2.5.   Serendipity DDR complex

The serendipity DDR complex is

where the discrete differential operators $\underline{\boldsymbol{{{G}}}}_h^{{k}}$, $\underline{\boldsymbol{{{C}}}}_h^{{k}}$ and $D_h^{{k}}$ are obtained projecting the edge/face/element operators onto the proper spaces (dictated by the co-domains):

It was proved that this sequence is indeed a complex ^[11,12], and has the same cohomology as the de Rham complex ^[30].

2.2.6.   Discrete $L^2$-inner products

To design numerical schemes based on the SDDR complex, an essential ingredient, besides the discrete differential operators, are consistent $L^2$-inner products on the spaces of the complex. A scheme can then be designed by replacing, in the weak formulation of the PDE, the continuous differential operators and $L^2$-products by the discrete operators of the complex and the $L^2$-inner products on its spaces.

The design of these discrete $L^2$-inner products rely on the element potential reconstructions defined in the previous sections. Specifically, if $\underline{X}_{\bullet, {h}}^{{k}}$ is one of the space $\underline{X}_{{\bf{grad}}, {h}}^{{k}}$, $\underline{\boldsymbol{{{X}}}}_{{\bf{curl}}, {h}}^{{k}}$ or $\underline{\boldsymbol{{{X}}}}_{{\rm{div}}, {h}}^{{k}}$ and $P_{\bullet, T}^{{k}}$ is the associated potential in the element $T$, the discrete $L^2$-product on $\underline{X}_{\bullet, {h}}^{{k}}$ is defined by

where the dot product in the integral is replaced by a multiplication if $\bullet = {\bf{grad}}\ $, and the stabilisation term $\mathrm{s}_{\bullet, T}$ penalises the difference between traces of the element potential and potential reconstructions on the face/edges (where relevant). The precise definition of the stabilisation term therefore depends on the space, and the available traces:

An important property of the potential reconstruction on each mesh entity is their polynomial consistency: Applied to interpolates of polynomials of the correct degree $\ell$ ($\ell = k+1$ for the gradient space, $\ell = k$ for the curl and divergence spaces), they return the polynomial itself. This translates into the following polynomial consistency of the $L^2$-inner products:

2.3.   Lie algebra-valued serendipity DDR complex

Since the Yang–Mills equations involve Lie algebra-valued functions, a Lie algebra-valued complex is required to discretise them. This complex is simply obtained by tensorisation of the real-valued complex, as in ^[2]: the spaces are made of Lie algebra-valued polynomials, and the operators of the complex act component by component on the Lie algebra.

In the following, we consider a Lie algebra $\mathfrak{g}$, that is, a finite-dimensional vector space endowed with a bilinear bracket $[\cdot, \cdot]: \mathfrak{g}\times \mathfrak{g}\to \mathfrak{g}$ and an inner product $\langle\cdot, \cdot\rangle: \mathfrak{g}\times \mathfrak{g}\to \mathbb{R}$ which satisfy the Jacobi identity

and the Ad-invariance property, which implies

We denote the Lie algebra-valued SDDR spaces by appending an exponent $\mathfrak{g}$ after the degree $k$. So, for example, the gradient space in the LASDDR (Lie algebra SDDR) complex is

We note that, selecting a basis $(e_I)_I$ of $\mathfrak{g}$, for any ${\mathsf{P}}\in\mathcal{M}_{h}$ we have

Here and in the following we use the implicit summation convention so, for example, $\phi^Ie_I = \sum_I \phi^Ie_I$. For a general space $X$, an element $v\in X\otimes \mathfrak{g}$ can be uniquely decomposed as $v = v^I\otimes e_I$. Any linear operator $L:X\to Y$ acting between two SDDR spaces $X, Y$ (or an SDDR space and a polynomial space) – such as a discrete differential operator, a potential reconstruction, etc. – then gives rise to the corresponding LASDDR operator $L^ \mathfrak{g}:X\otimes \mathfrak{g}\to Y\otimes \mathfrak{g}$ defined as $L^ \mathfrak{g}(v) = (L(v^I))\otimes e_I\in Y\otimes \mathfrak{g}$; this definition is independent of the choice of the basis in $\mathfrak{g}$. With these notations, the LASDDR complex is

In each space an inner product is obtained by tensorising the inner product of the corresponding SDDR space and of the Lie algebra. So, if $\bullet\in\{ \bf grad, {\bf{curl}}, div\}$,

Practical implementations of the LASDDR complex and related schemes can be easily done, in principle, by tensorising the operators and inner products of an SDDR implementation. Early tensorisation can however lead to unduly expensive calculations, especially when nonlinear terms are involved. We discuss in Section 4 the main considerations that must be taken into account to limit the assembly cost in implementations of LASDDR-based schemes.

3.   Two DDR-based schemes for the Yang–Mills equations

We propose two schemes for the Yang–Mills equations, which only differ in the handling of the nonlinear terms appearing in the equations. The first was introduced at the lowest order in ^[2], in which a discrete 'bracket' was constructed to approximate the value in the $\underline{\boldsymbol{{{X}}}}_{{\rm{div}}, {{h}}}^{{{k}, \mathfrak{g}}}$ space. This term is used in the discrete $L^2$-products, and the exact preservation of a discrete constraint, as well as energy estimates, are proven. Numerical tests in ^[2] however only considered the lowest-order $k = 0$ of the method.

The second method we present here is new, and leverages instead the continuous $L^2$-product and nonlinear bracket, as well as the elemental potential reconstructions, to achieve the same goal.

3.1.   Weak constrained form of the equations

Deriving from the bracket on the Lie algebra, the following two bilinear maps are defined:

We use in the discretisation a weak constrained formulation of the Yang–Mills equations, appearing previously in ^[20]: Find $(\boldsymbol{{A}}, \boldsymbol{{E}}, \lambda):[0, T]\to(\boldsymbol{{H}}({\bf{curl}}; {U})\otimes \mathfrak{g})^2\times (H^1({U})\otimes \mathfrak{g})$ such that

Note that the right-hand side of (3.3b) is equal to ${\int}_U\langle\boldsymbol{{B}}, {\bf{curl}}\ \boldsymbol{{v}}+ {\star}[{\boldsymbol{{A}}}, {\boldsymbol{{v}}}]^{}\rangle$, where the magnetic field is defined as $\boldsymbol{{B}}:= {\bf{curl}}\ \boldsymbol{{A}}+\frac12 {\star}[{\boldsymbol{{A}}}, {\boldsymbol{{A}}}]^{}$. We have developed this expression as it will drive different choices of discretisations. Solutions to these equations preserve the quantity

The use of this particular constrained form facilitates the preservation of a discrete counterpart in the numerical scheme; discussions on the derivation and implications of these continuous equations can be found in ^[2,20].

3.2.   Schemes

We consider a time discretisation $0 = t^0 < t^1 < \ldots < t^N = T$ of $[0, T]$ and denote the step in time between $n$ and $n+1$ as $\delta t^{n+\frac{1}{2}}:= t^{n+1}-t^n$. Then define for a family $v = (v^n)_n$, $\delta_t^{n+1}v = \frac{v^{n+1}-v^n}{ \delta t^{n+\frac{1}{2}}}.$

Starting from initial conditions $(\underline{\boldsymbol{{{A}}}}_h^0, \underline{\boldsymbol{{{E}}}}_h^0)\in (\underline{\boldsymbol{{{X}}}}_{{\bf{curl}}, {{h}}}^{{{k}, \mathfrak{g}}})^2$, the constrained scheme based on (3.3) is: Find families $(\underline{\boldsymbol{{{A}}}}_h^n)_n$, $(\underline{\boldsymbol{{E}}}_h^n)_n$, $(\underline{\lambda}_h^n)_n$ such that for all $n$, $(\underline{\boldsymbol{{{A}}}}_h^n, \underline{\boldsymbol{{E}}}_h^n, \underline{\lambda}_h^n)\in(\underline{\boldsymbol{{{X}}}}_{{\bf{curl}}, {{h}}}^{{{k}, \mathfrak{g}}})^2\times(\underline{X}_{{\bf{grad}}, {{h}}}^{{{k}, \mathfrak{g}}})$ and

We refer the reader to [2, Section 4] for a discussion on the choice of the initial conditions $(\underline{\boldsymbol{{{A}}}}_h^0, \underline{\boldsymbol{{{E}}}}_h^0)$, and also for alternative choices to the fully implicit time stepping selected here.

In (3.5b), $\mathfrak{N}\in\{\mathfrak{N}_1, \mathfrak{N}_2\}$ is one of the following two discretisations of the nonlinear terms in the right-hand side of (3.3b):

Above, we have set $\underline{\boldsymbol{{{A}}}}_h^{n+\frac{1}{2}} = \frac12(\underline{\boldsymbol{{{A}}}}_h^n+\underline{\boldsymbol{{{A}}}}_h^{n+1})$. Moreover, in $\mathfrak{N}_1$, we have made use of the discrete version ${\star}[{{\cdot}}, {{\cdot}}]^{{\rm{div}}, k, h}:(\underline{\boldsymbol{{{X}}}}_{{\bf{curl}}, {{h}}}^{{{k}, \mathfrak{g}}})^2\to\underline{\boldsymbol{{{X}}}}_{{\rm{div}}, {{h}}}^{{{k}, \mathfrak{g}}}$ of the map (3.2), defined for all $\underline{\boldsymbol{{{v}}}}_h, \underline{\boldsymbol{{{w}}}}_h\in\underline{\boldsymbol{{{X}}}}_{{\bf{curl}}, {{h}}}^{{{k}, \mathfrak{g}}}$ through its components by:

In $\mathfrak N_2$ we have used the global piecewise polynomial curl $\boldsymbol{\mathsf{C}}_h^{{k, \mathfrak{g}}}$ defined by patching the element curls: $(\boldsymbol{\mathsf{C}}_h^{{k, \mathfrak{g}}}\underline{\boldsymbol{{{v}}}}_h)_{|T} = \boldsymbol{\mathsf{C}}_T^{{k, \mathfrak{g}}}\underline{\boldsymbol{{{v}}}}_T$ for all $T\in\mathcal{T}_{h}$ and $\underline{\boldsymbol{{{v}}}}_h\in\underline{\boldsymbol{{{X}}}}_{{\bf{curl}}, {{h}}}^{{{k}, \mathfrak{g}}}$.

Remark 2 (Motivation for the discretisation of the nonlinear terms). The nonlinear terms in the right-hand side of (3.3b) are

When discretising these terms, the continuous fields $\boldsymbol{{A}}, \boldsymbol{{v}}\in\boldsymbol{{H}}({\bf{curl}}; {U})\otimes \mathfrak{g}$ are replaced by fully discrete objects $\underline{\boldsymbol{{{A}}}}_h, \underline{\boldsymbol{{{v}}}}_h\in\underline{\boldsymbol{{{X}}}}_{{\bf{curl}}, {{h}}}^{{{k}, \mathfrak{g}}}$, and we have to give meaning to the terms in (3.9) after this substitution – which is not straightforward since $\underline{\boldsymbol{{{X}}}}_{{\bf{curl}}, {{h}}}^{{{k}, \mathfrak{g}}}$ is not a subspace of $\boldsymbol{{H}}({\bf{curl}}; {U})\otimes \mathfrak{g}$.

Applying the standard DDR procedure on (3.9), we build these terms by replacing the inner product ${\int}_U\langle\cdot, \cdot\rangle$ and differential ${\bf{curl}}\ $ by the corresponding discrete notions found in the LASDDR complex. The only missing element is a discrete version of the bracket ${\star}[{\cdot}, {\cdot}]^{}$ on $\underline{\boldsymbol{{{X}}}}_{{\bf{curl}}, {{h}}}^{{{k}, \mathfrak{g}}}\times \underline{\boldsymbol{{{X}}}}_{{\bf{curl}}, {{h}}}^{{{k}, \mathfrak{g}}}$ which produces consistent discrete approximations in $\underline{\boldsymbol{{{X}}}}_{{\rm{div}}, {{h}}}^{{{k}, \mathfrak{g}}}$. This is what (3.8) provides, and this approach leads to $\mathfrak N_1$.

Another approach to discretising (3.9) is in a sense more straightforward (but only works because none of these terms, in the weak formulation, comes from integrating-by-parts the strong form of the model): Since we can reconstruct piecewise polynomial reconstructions and curls from elements in $\underline{\boldsymbol{{{X}}}}_{{\bf{curl}}, {{h}}}^{{{k}, \mathfrak{g}}}$, we can decide to simply substitute all the terms $\boldsymbol{{A}}, \boldsymbol{{v}}$ by these polynomial reconstruction based on $\underline{\boldsymbol{{{A}}}}_h, \underline{\boldsymbol{{{v}}}}_h$ and keep the other elements (integrals, brackets) exactly the same. This idea leads to $\mathfrak N_2$.

Remark 3 (Discretisation of the linear terms). The same way we used, in $\mathfrak N_2$, the piecewise polynomial potentials and element curl, we could consider replacing, in (3.5b), the term $(\underline{\boldsymbol{{{C}}}}_h^{{{k}, \mathfrak{g}}}\underline{\boldsymbol{{{A}}}}_h^{n+1}, \underline{\boldsymbol{{{C}}}}_h^{{{k}, \mathfrak{g}}}\underline{\boldsymbol{{{v}}}}_h)_{{\rm{div}}, \mathfrak{g}, h}$ with ${\int}_U \langle\boldsymbol{\mathsf{C}}_h^{{k, \mathfrak{g}}}\underline{\boldsymbol{{{A}}}}_h^{n+1}, \boldsymbol{\mathsf{C}}_h^{{k, \mathfrak{g}}}\underline{\boldsymbol{{{v}}}}_h\rangle$. This would however not lead to a suitable scheme, for the following reason.

Consider the pure Maxwell model, discretised using a linear unconstrained scheme (that is, (3.5a)–(3.5b) without the terms involving $\underline{\lambda}_h$, without the nonlinear terms and with $\mathfrak{g} = \mathbb{R}$ with the trivial Lie bracket):

Using the results in [12, Section 6] it can easily be shown that, for a smooth enough potential $\boldsymbol{{A}}$, solution of the continuous model, the consistency error (as defined in ^[31]) of the scheme satisfies

where $\|{{\cdot}}\|_{{{\bf{curl}}, \mathfrak{g}, h}}$ and $\|{{\cdot}}\|_{{{\rm{div}}, \mathfrak{g}, h}}$ denote the norms respectively associated with the inner products $(\cdot, \cdot)_{{\bf{curl}}, \mathfrak{g}, h}$ and $(\cdot, \cdot)_{{\rm{div}}, \mathfrak{g}, h}$. The scheme (3.10) is stable for the norm $\|{{\cdot}}\|_{{{\bf{curl}}, \mathfrak{g}, h}}+\|{\underline{\boldsymbol{{{C}}}}_h^{{k}}{\cdot}}\|_{{{\rm{div}}, \mathfrak{g}, h}}$, so (3.11) and the 3rd Strang Lemma ^[31] provide an $\mathcal O(\delta t +h^{k+1})$ error estimate.

However, replacing $(\underline{\boldsymbol{{{C}}}}_h^{{{k}, \mathfrak{g}}}\underline{\boldsymbol{{{A}}}}_h^{n+1}, \underline{\boldsymbol{{{C}}}}_h^{{{k}, \mathfrak{g}}}\underline{\boldsymbol{{{v}}}}_h)_{{\rm{div}}, \mathfrak{g}, h}$ with ${\int}_U \langle\boldsymbol{\mathsf{C}}_h^{{k, \mathfrak{g}}}\underline{\boldsymbol{{{A}}}}_h^{n+1}, \boldsymbol{\mathsf{C}}_h^{{k, \mathfrak{g}}}\underline{\boldsymbol{{{v}}}}_h\rangle$ in (3.10b) results in a scheme that is stable for the weaker norm $\|{{\cdot}}\|_{{{\bf{curl}}, \mathfrak{g}, h}}+\|{\boldsymbol{\mathsf{C}}_h^{{k, \mathfrak{g}}}\cdot}\|_{{L^2(U)}}$ (which does not control, in particular, the face curls). On the other hand, the consistency estimate remains (3.11), in the stronger norm. As a consequence, this estimate and the weaker stability cannot be combined together to obtain error estimates on the scheme. As a matter of fact, numerical tests (not reported in this paper) show that, on some mesh families, this alternative scheme does not converge as the mesh size and time step are refined.

3.3.   Discrete energy and constraint preservation

We define the discrete conserved quantity through the constraint functional $\mathfrak C^n:\underline{X}_{{\bf{grad}}, {{h}}}^{{{k}, \mathfrak{g}}}\to \mathbb{R}$:

Proposition 4 (Constraint preservation). For any choice of $\mathfrak{N}$, if $(\underline{\boldsymbol{{{A}}}}_h^n, \underline{\boldsymbol{{{E}}}}_h^n, \underline{\lambda}_h^n)$ solve (3.5) then, for all $\underline{q}_h\in\underline{X}_{{\bf{grad}}, {{h}}}^{{{k}, \mathfrak{g}}}$, the quantity $\mathfrak C^n(\underline{q}_h)$ is independent of $n$.

Proof. We note that the proof for the preservation of constraint in [2, Proposition 7] is in fact independent of the discretisation of the second equation (3.5b), and also applies in the case for general $k$ (see [2, Section 6.2]). □

To state the energy dissipation property, we introduce the discrete magnetic fields based on $\boldsymbol{{B}}$. Their nature depends on the chosen discretisation of the nonlinear terms in (3.3b). If $\mathfrak N = \mathfrak N_1$, exploiting the discrete bracket ${\star}[{{\cdot}}, {{\cdot}}]^{{\rm{div}}, k, h}$ we can define the discrete magnetic field as an element of $\underline{\boldsymbol{{{X}}}}_{{\bf{curl}}, {{h}}}^{{{k}, \mathfrak{g}}}$:

If $\mathfrak N = \mathfrak N_2$, the nonlinear terms being discretised as piecewise polynomial functions, the discrete magnetic field has the same nature:

Proposition 5 (Energy dissipation). If the initial conditions $(\underline{\boldsymbol{{{A}}}}_h^0, \underline{\boldsymbol{{{E}}}}_h^0)$ are such that $\mathfrak C^0\equiv 0$, then we have the decay of energy in the sense that, for all $n$,

where

where $\mathrm{s}_{{\rm{div}}, h}^{ \mathfrak{g}}$ is the stabilisation form involved in the definition of the inner product $(\cdot, \cdot)_{{\rm{div}}, \mathfrak{g}, h}$ (that is, the tensorisation of $\mathrm{s}_{{\rm{div}}, h}$ defined by (2.1)).

Proof. The proof for $\mathfrak N = \mathfrak N_1$ is identical to the one for $k = 0$ done in [2, Proposition 8] (see also [2, Section 6.2]).

Let us consider the case $\mathfrak N = \mathfrak N_2$. Choosing $\underline{\boldsymbol{{{v}}}}_h = \underline{\boldsymbol{{{E}}}}_h^{n+1}$ in (3.5b) and multiplying by $\delta t^{n+\frac12}$, the first term in the LHS is

while the remaining two form the constraint $\mathfrak C^{n+1}(\underline{\lambda}_h^{n+1})$, that vanishes by Proposition 4 and the assumption that $\mathfrak C^0\equiv 0$.

On the RHS, we use (3.5a) to substitute instead $\underline{\boldsymbol{{{v}}}}_h = -\delta_t^{n+1}\underline{\boldsymbol{{{A}}}}_h$, noting the cancellation of $\delta t^{n+\frac12}$ after multiplying. Expanding the discrete $L^2$-product by its definition and invoking [12, Proposition 7] (which can easily be extended to the SDDR complex using [1, Eq (2.2)]) to write $\boldsymbol{{P}}_{{\rm{div}}, h}^{{{k}, \mathfrak{g}}}\underline{\boldsymbol{{{C}}}}_h^{{{k}, \mathfrak{g}}} = \boldsymbol{\mathsf{C}}_h^{{k, \mathfrak{g}}}$, the first term is

Combining with the integrals in $\mathfrak N_2$ (see (3.7)), expanding $\underline{\boldsymbol{{{A}}}}_h^{n+\frac12} = \frac12(\underline{\boldsymbol{{{A}}}}_h^n+\underline{\boldsymbol{{{A}}}}_h^{n+1})$, recalling the definition (3.13) of $\boldsymbol{{B}}_{h, 2}$, then using the symmetry and bilinearity of the bracket (3.2), the RHS becomes

where the conclusion follows by applying the relation $b(x, y-x) = \frac12 b(y, y)-\frac12 b(x, x)-\frac12 b(y-x, y-x)$ to the symmetric bilinear forms $b = (\cdot, \cdot)_{L^2(U)\otimes \mathfrak{g}}$ and $b = \mathrm{s}_{{\rm{div}}, h}^ \mathfrak{g}$.

Finally arranging both side of the equation, moving the pure $n+1$ terms to the left and the rest to the right, we get

proving the statement, since the norm and $\mathrm{s}_{{\rm{div}}, h}^{ \mathfrak{g}}$ are both positive semidefinite. □

4.   Implementation

We cover in this section the broad mechanisms of how the schemes are implemented. For our numerical simulations, this implementation was done in the $\texttt{HArDCore3D}$ library (see https://github.com/jdroniou/HArDCore) starting from the serendipity DDR spaces and operators described in Section 2. This library contains a fully automated construction of these objects, including the computation of the degree depletions $\ell_ {\mathsf{P}}$ defined at the start of Section 2.2.

We first eliminate $\underline{\boldsymbol{{{A}}}}_h^{n+1}$ by using the first equation (3.5a) to write $\underline{\boldsymbol{{{A}}}}_h^{n+1} = \underline{\boldsymbol{{{A}}}}_h^n- \delta t^{n+\frac12}\underline{\boldsymbol{{{E}}}}_h^{n+1}$. Denoting the resulting equation by $F(\boldsymbol{{X}}^{n+1}) = \boldsymbol{{b}}$, where $\boldsymbol{{X}}^{n+1}$ represents the combined vector of $\underline{\boldsymbol{{{E}}}}_h^{n+1}, \lambda_h^{n+1}$, we then employ the iterative Newton method to find a solution up to an accuracy of $\epsilon$. The quantity $\underline{\boldsymbol{{{A}}}}_h^{n+1}$, which is required at the next time step, is finally recovered via back substitution.

At a single time step $n$, the most costly part of this process lies in the repeated assembly and resolution of the linear Newton problem: Find vectors $(\boldsymbol{{X}}^{n+1, i+1})_{i\in\mathbb{N}}$ such that

Although the derivative matrix $DF_{\boldsymbol{{X}}^{n+1, i}}$ is simple to determine because $F$ is multilinear, it must be rebuilt at every iteration, with terms stemming from bilinear, trilinear, and even quadrilinear forms (see the product of bilinear brackets in the last terms of (3.6) and (3.7)). In the rest of this section, we discuss how to perform these calculations without it becoming too expensive in either memory space or computational time.

The other major expense is resolving the linear system, for which we use the $\texttt{Intel}\ \texttt{MKL}\ \texttt{PARADISO}$ library (see https://software.intel.com/en-us/mkl), which provides a multi-threaded direct solver. An efficient technique to reduce this solver cost is to apply to the linear systems the static condensation process, which eliminates all elemental unknowns of the system prior to solving. This is made possible by the specific stencil resulting from a hybrid method like (S)DDR, which couples the unknowns inside one element only with the unknowns on the faces, edges and vertices of that element (inter-element unknowns are never directly coupled). We emphasize here the difference between static condensation and the serendipity DDR process. Both reduce the number of degrees of freedom, but the SDDR spaces and operators are leaner from the start; as a result, any scheme built from it will see an improved performance in every aspect of the implementation at no additional cost. In contrast, static condensation reduces the number of unknowns only after all contributions are assembled, and thus its effect is more limited as it only reduces the cost of solving the global system, not the cost of assembling that system. Additionally, it must be repeated every iteration, with some overhead (solving a smaller linear system in each element) each time. Finally, we should highlight that, since static condensation only eliminates unknowns in the elements while serendipity also eliminates degrees of freedom on the edges/faces, the linear systems resulting from a statically condensed DDR scheme remain larger in general than the linear systems resulting from a statically condensed SDDR scheme.

4.1.   LASDDR tensorisation

The numerical construction of the LASDDR complex primarily consists of wrapping a layer of matrix tensorisation around the existing SDDR code. For a code like $\texttt{HArDCore3D}$ based on the $\texttt{Eigen3}$ library (see http://eigen.tuxfamily.org), this is easily achieved using the $\texttt{KroneckerProduct}$ functionality.

Fixing a basis $(e_I)_I$ of the $d$-dimensional Lie algebra $\mathfrak{g}$, each Lie algebra-valued degree of freedom can by expressed by $d$ real values. In other words, we can think of an element of an LASDDR space as being made up of $d$ SDDR vectors, one for each basis $e_I$; i.e., $\underline{\boldsymbol{{{v}}}}_h\equiv\underline{\boldsymbol{{{v}}}}_h^I\otimes e_I$. Fixing an ordering that combines everything into a single vector fixes the physical interpretation of all the remaining operators. We choose to store the values associated to each mesh entity (vertex, edge, face, element) sequentially, but for ease of distinguishing the significance of each entry, they are doubly indexed: $(\underline{\boldsymbol{{{v}}}}_h^i)^I$. The lowercase letter numbers the mesh entity it originates from, and the capital letter labels the Lie algebra basis it is attached to. As an example, if $d = 3$, we have the following structures:

Then, considering a linear LASDDR operator $L^ \mathfrak{g}$, the matrix representation $\boldsymbol{\sf L}^{ \mathfrak{g}}$ must by definition act as $\boldsymbol{\sf L}^{ \mathfrak{g}}(\underline{\boldsymbol{{{v}}}}_h)\equiv\boldsymbol{\sf L}(\underline{\boldsymbol{{{v}}}}_h^I)\otimes e_I$, where $\boldsymbol{\sf L}$ is the corresponding SDDR matrix operator. From some basic arithmetic, we can conclude that $\boldsymbol{\sf L}^{ \mathfrak{g}}$ has the form

where $\boldsymbol{\sf I}_{{d}}$ is the $d\times d$ identity matrix. For bilinear operators such as the discrete inner products $(\cdot, \cdot)_{\bullet, \mathfrak{g}, h}$, and the bracket terms which we deal with in the next section, the idea is very much the same; the difference lies in the usage of a more general $d\times d$ matrix $\boldsymbol{\sf M}$ in place of $\boldsymbol{\sf I}_{{d}}$, indicating an interaction of the Lie algebra bases. For example, denoting the LASDDR (resp. SDDR) product matrix by $\boldsymbol{\sf B}^ \mathfrak{g}$ (resp. $\boldsymbol{\sf B}$), with action previously defined as $(\underline{\boldsymbol{{{v}}}}_h)^T\boldsymbol{\sf B}^ \mathfrak{g}(\underline{\boldsymbol{{{w}}}}_h) = ((\underline{\boldsymbol{{{v}}}}_h^I)^T\boldsymbol{\sf B}(\underline{\boldsymbol{{{w}}}}_h^J))\langle e_I, e_J\rangle$, the $\boldsymbol{\sf M}$ becomes evidently the mass matrix of the Lie algebra:

4.2.   Bracket terms

The nonlinear terms in both schemes create operators which can not be encoded by a single matrix, but rather (local or global) 3 or 4-dimensional arrays:

The tools available in $\texttt{Eigen3}$ for dealing with these objects are not nearly as developed as the ones for matrices. In most cases, the multidimensional storage is done using the $\texttt{Boost.MultiArray}$ library (see https://www.boost.org/doc/libs/1_61_0/libs/multi_array/doc/index.html), but the $\texttt{Eigen::Map}$ function is used to interpret the data, so that we can still perform the usual matrix operations.

These generic objects are independent of time, so seemingly the most time efficient method would be to pre-compute them once and for all, and recall them when necessary. Unfortunately, the extra dimensionalities mean these operators are much larger than the (bi)linear ones appearing in LASDDR; in fact the memory usage to store these terms grows exponentially with the number of entries. The Lie algebra tensorisation only exacerbates this problem, even accounting for the many symmetries that could be exploited. This memory issue is particularly sensitive in an implementation – such as ours (which follows the HArDCore general strategy) – that assumes that each element can have its own geometry, which forces the local arrays to be computed/stored independently for each element (in a situation where the mesh elements can be classified using a few reference elements, all memory issues disappear as only local multilinear maps in reference elements need to be stored). In this context, the immense amount of memory required to store just a single one of these global maps means that there is no choice but to recalculate them at each time step and locally (mesh entity by mesh entity) as needed.

Another important effect on the runtime lies in the order in which some tensorisation-related computations are performed. Taking the sum of smaller matrices multiple times is sometimes preferable to doing it once with larger matrices (which often have lots of zeros). Therefore, delaying the tensorisation until after the vectors have been evaluated can improve both the memory usage and the speed, even though some calculations and the tensorisation have to be repeated.

We furnish these ideas with an example for the nonlinear terms of the form

in which $\underline{\boldsymbol{{{v}}}}_T$ is a given vector in $\underline{\boldsymbol{{{X}}}}_{{\bf{curl}}, {{T}}}^{{{k}, \mathfrak{g}}}$. This term is represented by a coefficient matrix with entries (doubly indexed by $((i, I), (k, K))$) given by:

where the rows and columns range over the basis vectors, that are defined as $(\phi_i)_I\equiv\phi_i\otimes e_I$ (resp. $(\psi_k)_K)$ where $(\phi_i)_i$ is a basis for $\underline{\boldsymbol{{{X}}}}_{{\bf{curl}}, {T}}^{{k}}$ (resp. $(\psi_k)_k$ a basis of $\underline{X}_{{\bf{grad}}, {T}}^{{k}}$). Pulling the vector out (recall that we use implicit summation), and using the definition of the $L^2$-product, this can be viewed as

where $M$ is indeed a trilinear form (indices $((i, I), (j, J), (k, K))$), represented by a 3-dimension block of $d^3\times(\dim(\underline{\boldsymbol{{{X}}}}_{{\bf{curl}}, {T}}^{{k}}))\times(\dim(\underline{\boldsymbol{{{X}}}}_{{\bf{curl}}, {T}}^{{k}}))\times(\dim(\underline{X}_{{\bf{grad}}, {T}}^{{k}}))$ coefficients. As mentioned, although $M$ would be useful to calculate in itself, since it can then be used to find other terms in the scheme, performing the contraction with each vector $\underline{\boldsymbol{{{v}}}}_T$ is actually quite slow because of the large matrix sums.

Instead we do something less intuitive, by combining the vector with the integral portion of the product first:

These integral coefficients represent the smaller SDDR trilinear form $\displaystyle {\int}_T(\boldsymbol{{P}}_{{\bf{curl}}, T}^{{k}}\cdot)\cdot(\boldsymbol{{P}}_{{\bf{curl}}, T}^{{k}}\cdot)(P_{{\bf{grad}}, T}^{{k+1}}\cdot)$; this term is first combined (through the sum over $j$) with the vector $(\underline{\boldsymbol{{{v}}}}_T)^J$, before performing the tensorisation with $N_{I, J, K}$ (on $(i, I), (k, K)$), and finally taking the much shorter matrix sum over $J$. This delay in combining the Lie algebra indices initially seems like extra work, as this $M^J_{i, k}$ sum is unique to each vector $(\underline{\boldsymbol{{{v}}}}_T)^J$, and so the tensorisation must be re-done each time, but testing showed that the tradeoff for smaller matrix sums is worth it in this case.

The ${\star}[{\cdot}, {\cdot}]^{}$ bracket is dealt with differently because it can appear twice in a single product. In this double bracket case, if we use the same summation methods as in (4.2), then we have the appearance of quadrilinear forms instead of trilinear forms, that would require the computation of a 4-dimensional array. To avoid the extra dimension, we treat the bracket terms as independent objects, that can be manipulated separately to the product matrix. For example, in the first discretisation $\mathfrak N_1$ described in (3.6), two terms involve this bracket, which we compute the following way (underbrace denotes the dimension of the array representation with the implied transposition, and we write formal products to show how the calculation could be decomposed in the code):

We note that the map ${\star}[{{\cdot}}, {{\cdot}}]^{{\rm{div}}, k, h}:(\underline{\boldsymbol{{{X}}}}_{{\bf{curl}}, {{h}}}^{{{k}, \mathfrak{g}}})^2\to\underline{\boldsymbol{{{X}}}}_{{\rm{div}}, {{h}}}^{{{k}, \mathfrak{g}}}$ is bilinear, but requires an extra dimension in the corresponding array to represent the output in $\underline{\boldsymbol{{{X}}}}_{{\rm{div}}, {{h}}}^{{{k}, \mathfrak{g}}}$. Thus the terms ${\star}[{{\underline{\boldsymbol{{{v}}}}_h}}, {{\cdot}}]^{{\rm{div}}, k, h}$ and ${\star}[{{\underline{\boldsymbol{{{w}}}}_h}}, {{\cdot}}]^{{\rm{div}}, k, h}$ are indeed represented by matrices, calculated using the same principle (described in (4.2)) of avoiding the tensorisation until the last step. The order of operations is also crucial in (4.4); the vector-matrix multiplication must come first, to ensure that the 3-dimensional array is only contracted with a vector. We stress again however, that the full set of ${\star}[{{\cdot}}, {{\cdot}}]^{{\rm{div}}, k, h}$ coefficients is never constructed, and all calculations are done in the way of (4.2).

The same idea is implemented in the second discretisation $\mathfrak N_2$ (see (3.7)), by introducing a basis for $\boldsymbol{{\mathcal{{P}}}}^{{2k}}(T)$, and working with the map

The numerical decomposition of the term analogous to (4.3) is

where the integral is realised by the tensorisation of the mass matrix of the basis on $\boldsymbol{{\mathcal{{P}}}}^{{2k}}(T)$, and the mass matrix of the Lie algebra (see (4.1)). With orthonormal choices of bases of $\boldsymbol{{\mathcal{{P}}}}^{{2k}}(T)$ and $\mathfrak g$ (which is the default in the $\texttt{HArDCore}$ library), the calculations here are greatly simplified; the components of ${\star}[{\boldsymbol{{P}}_{{\bf{curl}}, T}^{{{k}, \mathfrak{g}}}\cdot}, {\boldsymbol{{P}}_{{\bf{curl}}, T}^{{{k}, \mathfrak{g}}}\cdot}]^{}$ can be found using only the triple integrals of the bases of $\mathcal{P}^{{k}}(T)$ and $\boldsymbol{{\mathcal{{P}}}}^{{2k}}(T)$ (without requiring to solve a linear system afterwards), and the integral in (4.5) is just given by the identity matrix

5.   Numerical tests

We present a numerical comparison of the convergence of the two schemes, as well as the exact constraint preservation that is expected. The tests were performed on a Dell Precision 5820 desktop with a 14-core Intel Xeon processor (W-2275) clocked at 3.3 GHz and equipped with 128 GB of DDR4 RAM, running Ubuntu 22.04.1 LTS. The discretisation setting is identical to that of [2, Section 5], which only contained tests for $k = 0$ of the scheme pertaining to $\mathfrak N_1$. Here we expand on these results for higher orders ($k = 0, 1, 2$), and also consider the performance in relation to the second discretisation.

Let us recall the setting of these tests. The Lie algebra is $\mathfrak{g} = \mathfrak{su}(2)$, with basis

The time interval is [0, 1] and the space domain is the unit cube $(0, 1)^3$; the spatial discretisation is based on three families of Voronoi, tetrahedral, and cubic cell meshes. For each mesh of size $h$, the time interval is uniformly{\rm{div}}ided into $\max\{10, \left \lceil{5/h}^{k+1}\right \rceil\}$ time steps; given that we use an implicit time discretisation, the expected rate of convergence is in $\delta t+h^{k+1}$, and the choice of time step is designed so that, when plotted against $h$, the errors should decay as $h^{k+1}$. To set non-zero initial conditions, and to assess the convergence properties of the schemes, we select a manufactured solution based on

from which $\boldsymbol{{E}}(t)$ is calculated using (3.3a). Then for all tests in this section, the initial conditions are assumed to be $\underline{\boldsymbol{{{A}}}}_h^0 = \underline{\boldsymbol{{{I}}}}_{{\bf{curl}}, {{h}}}^{{{k}, \mathfrak{g}}}\boldsymbol{{A}}(0), \underline{\boldsymbol{{{E}}}}_h^0 = \underline{\boldsymbol{{{I}}}}_{{\bf{curl}}, {{h}}}^{{{k}, \mathfrak{g}}}\boldsymbol{{E}}(0)$.

Remark 6 (Convergence results for $\underline{\lambda}_h$). Although the fields $\boldsymbol{{A}}$, $\boldsymbol{{E}}$ of a solution to the constrained formulation (3.3) solve the Yang–Mills equations, there is no proof that the $\lambda$ obtained is unique. Testing performed in [2, Section 5] suggest indeed that there are infinitely many solutions; the values obtained for $\underline{\lambda}_h^n$ are therefore not very instructive, and have been omitted from the graphs. Discussion around the solvability of the linear system deriving from such a scheme can also be found in the cited section; we experienced a similar success, with a worst residual for the linear solver of the order 1e-09.

5.1.   Convergence tests

In addition, appropriate boundary conditions and forcing terms are introduced to balance the equations (see [2, Section 5.1] for details). The errors for both schemes are measured by calculating the difference $\|{\underline{\boldsymbol{{{E}}}}_h^{n+1}-\underline{\boldsymbol{{{I}}}}_{{\bf{curl}}, {{h}}}^{{{k}, \mathfrak{g}}}\boldsymbol{{E}}(1)}\|_{{{\bf{curl}}, \mathfrak{g}, h}}$ (resp. $\underline{\boldsymbol{{{A}}}}_h^{n+1}, \boldsymbol{{A}}$) at the final time, and dividing by the norm of $\underline{\boldsymbol{{{E}}}}_h^{n+1}$ (resp. $\underline{\boldsymbol{{{A}}}}_h^{n+1}$). These relative errors are plotted in Figure 1 for $\boldsymbol{{E}}$, and Figure 2 for $\boldsymbol{{A}}$.

We remark immediately that for $\boldsymbol{{A}}$, the errors are indistinguishable from the figure alone. The precise difference for the Voronoi and cubic sequences are calculated in Table 1, with the trend applying identically to the tetrahedral family. As either $k$ increases or $h$ decreases, the difference between the errors shrink accordingly, and this is seen also in the errors for $\boldsymbol{{E}}$ for $k = 1, 2$. However for the electric field, there is a more visible separation when $k = 0$, with the $\mathfrak N_1$-based discretisation performing slightly better. These tests seem to indicate that, overall, both choices $\mathfrak N_1, \mathfrak N_2$ lead to acceptable and similar results.

The schemes converge on the Voronoi and cubic meshes at the expected rate of $k+1$ for both $\boldsymbol{{E}}$ and $\boldsymbol{{A}}$, but this behaviour was not as stable for $\boldsymbol{{E}}$ on the tetrahedral line, where we see a rate that jumps around $3$ for every $k$. This might be due to the asymptotic regime not been reached yet on these meshes (we note that, for $k = 0$ for example, the simulation on the finest mesh seem to indicate that the convergence rate slows down). The magnitude of these errors are still ordered in the expected way, except for the coarsest cubic mesh in Figure 1, where it is corrected after the first refinement.

5.2.   Constraint preservation

The tests for the preservation of constraint are run with the same initial conditions as the convergence tests. For proper solutions to the Yang–Mills equations, and for Proposition 5, it is expected that these discrete fields prescribe a small or vanishing initial constraint $\mathfrak C^0$. This can be achieved by projecting the initial conditions (see ^[2]), however for the purpose of testing Proposition 4, which does not depend on the initial values, it suffices to measure the maximum change in the functional $\mathfrak C^n-\mathfrak C^0$ over all times. The differences are presented in Table 2 for selected meshes from each sequence. We see for both methods that the constraint is stationary up to machine precision, with the small drift coming from rounding errors present at each iteration.

For reference, we also provide in Table 3 a comparison of the runtimes for the two different choices of discretisation of the nonlinear tests. The performances of both schemes are very comparable on a variety of meshes and degrees $k$, with the exception of a $10\%$ difference in favour of $\mathfrak N_1$ for $k = 2$ on the Tetrahedral meshes, that shows up consistently through our tests. This difference is likely due to the nontrivial calculation attached to each face component (3.8a) in the definition of the discrete bracket of $\mathfrak N_1$. These calculations are necessary because the face values also contribute to the $L^2$-product, but they result in a larger dependency of the runtime on the number of faces in a particular mesh. In comparison, the $\mathfrak N_2$ discretisation is less affected by the shape of the elements; numerically, an extra face only represents an increase in size of the SDDR operators (which equally affects $\mathfrak N_1$) used in the integrals. This is supported by the results for the Voronoi sequence, that has a higher face to element ratio than the tetrahedral sequence, where $\mathfrak N_2$ starts to slightly outperform $\mathfrak N_1$.

6.   Conclusions

We designed two schemes for the Yang–Mills equations based on the DDR method, both displaying arbitrary orders of accuracy and applications on generic polyhedral meshes. Thanks to the complex property of DDR and to the usage of a Lagrange multiplier, both schemes also preserve a discrete nonlinear constraint deriving from the Yang–Mills equations, and satisfy energy bounds. The schemes only differ in their treatment of the nonlinearity akin to a cross product combined with the Lie bracket for Lie algebra-valued vector functions. The first scheme reconstructs a discrete version of the continuous product bracket, that can then be used in the discrete $L^2$-products of the DDR complex. The second scheme uses the DDR potential reconstructions to get polynomials in each element, on which the continuous product bracket can be applied. We show how a clever ordering of the algebraic operations in the assembly of the schemes can help keep the computational cost at a reasonable level, despite needing to deal with multidimensional arrays and high system sizes due to the Lie algebra components. Numerical results are presented which show a good behaviour and an expected rate of convergence, with respect to the mesh size, in $h^{k+1}$.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

Funded by the European Union (ERC, NEMESIS, No. 101115663). Views and opinions expressed are however those of the authors only and do not necessarily reflect those of the European Union or the European Research Council Executive Agency. Neither the European Union nor the granting authority can be held responsible for them.

Conflict of interest

The authors declare no conflicts of interest.