Globally constraint-preserving FR/DG scheme for Maxwell’s equations at all orders

The numerical solution of Maxwell’s equations plays an incredibly important role in the computational solution of many problems in science and engineering. The Finite Difference Time Domain (FDTD) method, originally proposed by Yee [57] and greatly developed since the seminal papers by Taflove and Brodwin [49], [48], has been the mainstay for computational electrodynamics (CED) chiefly because it globally preserves the involution constraints that are inherent in Maxwell’s equations. Several modern texts and reviews document the development of FDTD (Taflove and Hagness [50], Taflove, Oskooi and Johnson 2013, Gedney [34]). As a result, given the enormously well-cited works by Yee and Taflove, it is very desirable to include global preservation of involution constraints into more modern schemes for CED.

FDTD had been formulated well before the explosion in modern higher order Godunov schemes, which started with the pioneering work of van Leer [51], [52]. Since that work, there has been a strong drive to include the physics of wave propagation into the numerical solution of hyperbolic systems; CED being a case in point. Two strong strains of early effort to develop higher order Godunov schemes for CED include finite-volume time-domain (FVTD) methods (Munz et al. [41], Ismagilov [37], Barbas and Velarde [20]) and discontinuous Galerkin time-domain (DGTD) methods (Cockburn & Shu [26], [27]; Cockburn et al. [25], Hesthaven and Warburton [35]). There has also been a strong effort in the engineering CED community to design DGTD schemes for CED (Chen & Liu [24], Ren et al. [42], Angulo et al. 2015) and some of those methods indeed use locally constraint-preserving bases. However, none of those DGTD methods incorporated the very beneficial globally constraint-preserving aspect of FDTD. For that next phase of evolution, one had to wait for developments that emerged in the field of numerical MHD and are now rapidly finding their way into CED.

In MHD, one evolves Faraday’s law in addition to the equations of computational fluid dynamics (Brecht et al. [21], Evans and Hawley [33], DeVore [30], Dai and Woodward [28], Ryu et al. [43], Balsara and Spicer [18]). While studying adaptive mesh refinement and numerical schemes for MHD, advances were made in constraint-preserving reconstruction of magnetic fields (Balsara [1], [2], [3], Balsara and Dumbser [9], Xu et al. [56]). This made it possible to start with the face-centered magnetic induction fields in the Yee-type mesh and specify them at all locations within a computational zone. The edge-centered electric fields that are inherent to a constraint-preserving update of Faraday’s law would then have to be multidimensionally upwinded. This multidimensional upwinding was achieved by using a newly-designed multidimensional Riemann solver (Balsara [4], [5], [6], [7], Balsara, Dumbser and Abgrall [11], Balsara and Dumbser [10], Balsara and Nkonga [16]). These twin innovations, consisting of constraint-preserving reconstruction of vector fields, and multidimensional Riemann solvers, permitted a logically complete description of numerical MHD. Along the way, a third innovation in ADER (Arbitrary DERivatives in space and time) timestepping schemes was added which, while not essential, greatly simplified the accurate temporal evolution of MHD variables (Dumbser et al. [31], [32], Balsara et al. [17], [15]). DG schemes for the induction equation that were globally constraint-preserving were also devised by Balsara and Käppeli [13]. The stage was now set for migrating these innovations back again to CED.

FVTD schemes for CED that were based on the above-mentioned three innovations were developed in Balsara et al. ([8], [19], [12]). The constraint-preservation was accomplished by making a constraint-preserving reconstruction of the magnetic induction and the electric displacement. Unlike FDTD that operates on a pair of staggered control volumes, the present methods operate on the same control volume, see Fig. 1 from Balsara et al. [19]. To ensure the mimetic preservation of the constraints, the primal variables were taken to be the facially collocated normal components of the electric displacement and the magnetic induction; where both vector fields were collocated on the same faces. The facially collocated magnetic induction evolves in response to the edge collocated electric fields, yielding a discrete representation of Faraday’s law. The facially collocated electric displacement evolves in response to the edge collocated magnetic fields, yielding a discrete representation of Ampere’s law. The resulting methods were indeed globally constraint-preserving in that the magnetic induction always remains divergence-free at all locations on the mesh and the electric displacement always satisfies the constraint imposed by Gauss’ Law at any location on the mesh. All these above-mentioned advances that were made to mimetic FVTD schemes have been recently embedded into the DGTD schemes that we report on below.

Globally constraint-preserving DGTD schemes were also explored in Balsara and Käpelli (2018, BK henceforth). BK found that higher order DGTD schemes were almost totally free of dispersion error, having a dispersion error that was almost 75 times smaller than FDTD. Since the methods were based on Riemann solvers, some dissipation is inevitable. Even so, BK found that the higher order DGTD schemes were almost free of dissipation even when electromagnetic waves spanned only a few zones. The first paper by BK was strongly focused on von Neumann stability analysis of DGTD schemes for CED. The present paper extends this study in several ways, which we list in the rest of this paragraph. First, BK focused on DGTD schemes up to fourth order, whereas the present paper presents flux reconstruction (FR) based DGTD schemes¹¹ 1 We will use DGTD to refer to the current scheme though strictly it is a combination of FR and DG schemes. that go up to fifth order, and the reconstruction schemes at fourth and fifth order presented in this work are new. This drive to higher order enables us to generalize an observation from BK who found that from fourth order and upwards one has to include some volumetrically-evolved modes in addition to the facially-evolved modes. Second, BK did not touch on the topic of limiting DGTD schemes, specifically because they realized that the non-linear hybridization of DGTD schemes for CED should be done with the utmost carefulness, if even it is needed. The limiting of any DG scheme always diminishes the optimal accuracy of a DG scheme. This reduction occurs to a greater or lesser extent based on whether the limiter is applied more or less aggressively, respectively. BK did find that DGTD schemes did not need any limiting up to fourth order but they only tried situations where the permittivity and permeability were constant. In this paper we report the very favorable finding that DGTD schemes don’t seem to require any limiting even when the permittivity and permeability vary by almost an order of magnitude. This result is quite useful; however, it is predicated on the assumption that the conductivity is zero. Dealing with non-zero conductivity will be the topic of a subsequent study. Third, BK realized that when conductivity is zero the Maxwell’s equations conserve a quadratic energy. While this energy conservation was not built into the scheme, BK found the very desirable result that quadratic energy was effectively conserved by fourth order DGTD schemes when the waves spanned only a few zones. This paper extends this study to fifth order DGTD schemes. Furthermore, we study the energy conservation of high order DGTD when the permittivity and permeability vary with space. A fourth offering in this paper is a proof that DGTD schemes for CED that are constructed according to the principles in BK and this paper are indeed energy stable.

The rest of the paper is organized as follows. Section (2) introduces the Maxwell’s equations and the simplified 2-D model that we consider in this paper. Section (3) explains the polynomial spaces used to approximate the solution variables and section (4) explains the divergence-free reconstruction scheme at fourth order of accuracy. The numerical descritization of the Maxwell’s equations using flux reconstruction and DG method is shown in section (5) together with constraint preserving property. In section (6), we perform the stability analysis of the semi-discrete scheme at first order and show the dissipative character coming from the 1-D and 2-D Riemann solvers. Section (7) presents many test cases to demonstrate the performance of the scheme. Finally, the Appendix explains the reconstruction scheme at other orders and also briefly discusses the Riemann solvers.

The Maxwell’s equations are a system of linear partial differential equations that model the wave propagation behaviour of electric and magnetic fields in free space and material media. They can be written in vector form as

$$\frac{\partial\bm{B}}{\partial t}+\nabla\times\bm{E}=0,\qquad\frac{\partial\bm{D}}{\partial t}-\nabla\times\bm{H}=-\bm{J}$$

where

$\bm{B}$ = magnetic flux density	$\bm{D}$ = electric flux density
$\bm{E}$ = electric field	$\bm{H}$ = magnetic field
	$\bm{J}$ = electric current density

The fields are related to one another by constitutive laws

$$\bm{B}=\mu\bm{H},\qquad\bm{D}=\varepsilon\bm{E},\qquad\bm{J}=\sigma\bm{E}\qquad\mu,\varepsilon\in\mathbb{R}^{3\times 3}\textrm{ symmetric}$$

where the coefficients

	$\displaystyle\varepsilon$	$\displaystyle=\textrm{permittivity tensor}$
	$\displaystyle\mu$	$\displaystyle=\textrm{magnetic permeability tensor}$
	$\displaystyle\sigma$	$\displaystyle=\textrm{conductivity}$

are material properties and are in general tensorial functions of spatial coordinates. In free space $\varepsilon=\varepsilon_{0}=8.85\times 10^{-12}$ F/m and $\mu=\mu_{0}=4\pi\times 10^{-7}$ . The divergence of the electric flux gives the electric charge density $\rho$, which itself obeys a conservation law

$$\nabla\cdot\bm{D}=\rho,\qquad\frac{\partial\rho}{\partial t}+\nabla\cdot\bm{J}=0$$

Moreover, since magnetic monopoles have never been observed in nature, the divergence of the magetic flux must be zero $\nabla\cdot\bm{B}=0$, which is an additional constraint that must be satisfied by the solution. Note that if the initial condition is divergence-free, then under the time evolution induced by the Maxwell’s equations, the divergence remains zero at future times also.

In the present work, we will consider a 2-D model of the Maxwell’s equations (TE polarization) for which the equations can be written in Cartesian coordinates as

	$\displaystyle\frac{\partial D_{x}}{\partial t}-\frac{\partial H_{z}}{\partial y}$	$\displaystyle=$	$\displaystyle 0$		(1)
	$\displaystyle\frac{\partial D_{y}}{\partial t}+\frac{\partial H_{z}}{\partial x}$	$\displaystyle=$	$\displaystyle 0$		(2)
	$\displaystyle\frac{\partial B_{z}}{\partial t}+\frac{\partial E_{y}}{\partial x}-\frac{\partial E_{x}}{\partial y}$	$\displaystyle=$	$\displaystyle 0$		(3)

where

$$(E_{x},E_{y})=\frac{1}{\varepsilon}(D_{x},D_{y}),\qquad H_{z}=\frac{1}{\mu}B_{z}$$

and $\mu$, $\varepsilon$ are scalars which may depend on spatial coordinates. We will consider a constraint on the divergence of the electric flux density $\bm{D}$ instead of the magnetic field, which has only one component in the above model. The above system of three PDE can be written in conservation form

$$\frac{\partial\bm{U}}{\partial t}+\frac{\partial\bm{F}}{\partial x}+\frac{\partial\bm{G}}{\partial y}=0$$

(4)

where

$$\bm{U}=\begin{bmatrix}D_{x}\\ D_{y}\\ B_{z}\end{bmatrix},\qquad\bm{F}=\begin{bmatrix}0\\ H_{z}\\ E_{y}\end{bmatrix}=\begin{bmatrix}0\\ \frac{1}{\mu}B_{z}\\ \frac{1}{\varepsilon}D_{y}\end{bmatrix},\qquad\bm{G}=\begin{bmatrix}-H_{z}\\ 0\\ -E_{x}\end{bmatrix}=\begin{bmatrix}-\frac{1}{\mu}B_{z}\\ 0\\ -\frac{1}{\varepsilon}D_{x}\end{bmatrix}$$

This is a system of hyperbolic conservation laws for which a Riemann problem can be solved exactly to determine the fluxes required in the numerical schemes like finite volume and DG method. This is explained in the Appendix. While for simplicity, we consider scalar material properties, the constraint preserving nature of the scheme holds for general material properties. In fact, everything we describe in this paper holds for the general case, and the only additional change required is to use the Riemann solvers for the general case, which are explained e.g., in [19].

In the absence of currents, the Maxwell’s equations conserve the total energy provided there is no net gain of energy at the boundaries or if we have periodic boundaries. For the 2-D model that we consider in this work, we can first show that the following additional conservation law holds

$$\frac{\partial}{\partial t}\left[\frac{1}{2\varepsilon}(D_{x}^{2}+D_{y}^{2})+\frac{1}{2\mu}B_{z}^{2}\right]+\frac{\partial}{\partial x}(H_{z}E_{y})-\frac{\partial}{\partial y}(H_{z}E_{x})=0$$

which implies that the quantity

$$\mathcal{E}(t)=\int_{\Omega}\left[\frac{1}{2\varepsilon}(D_{x}^{2}+D_{y}^{2})+\frac{1}{2\mu}B_{z}^{2}\right]\mbox{d}x\mbox{d}y$$

which is the total energy, is conserved under periodic boundary conditions or if the net flux at the boundaries is zero.

We would like to approximate the vector field $\bm{D}$ such that its divergence is zero inside the cell if there is no charge density. If there is some electric charge, then the divergence of $\bm{D}$ must match this charge density, but we do not deal with this case in the present work. The approach we take to ensure divergence-free property is to use the divergence-free reconstruction ideas of Balsara [1], [2], [3] which makes use of known values of the normal component of $\bm{D}$ on the faces of the cell and then reconstruct the vector field inside the cell by enforcing appropriate constraints on its divergence.

We will approximate the normal components of $\bm{D}$ on the faces by one dimensional polynomials of degree $k\geq 0$. We map each cell to the reference cell $[-{\frac{1}{2}},+{\frac{1}{2}}]\times[-{\frac{1}{2}},+{\frac{1}{2}}]$ with coordinates $(\xi,\eta)$. Let $\mathbb{P}_{k}(\xi)$ denote one dimensional polynomials of degree at most $k$ in the variable $\xi$, and $\mathbb{P}_{k}(\xi,\eta)$ denote two dimensional polynomials of degree at most $k$. On the two vertical faces of a cell, the normal component is given by

$$D_{x}^{\pm}(\eta)=\sum_{j=0}^{k}a_{j}^{\pm}\phi_{j}(\eta)\in\mathbb{P}_{k}(\eta)$$

while on the two horizontal faces, the corresponding normal component is given by

$$D_{y}^{\pm}(\xi)=\sum_{j=0}^{k}b_{j}^{\pm}\phi_{j}(\xi)\in\mathbb{P}_{k}(\xi)$$

The location of these polynomials is illustrated in Figure (1). The basis functions $\phi_{j}$ are mutually orthogonal polynomials given by

$$\phi_{0}(\xi)=1,\quad\phi_{1}(\xi)=\xi,\quad\phi_{2}(\xi)=\xi^{2}-\tfrac{1}{12},\quad\phi_{3}(\xi)=\xi^{3}-\tfrac{3}{20}\xi,\quad\phi_{4}(\xi)=\xi^{4}-\tfrac{3}{14}\xi^{2}+\tfrac{3}{560}$$

$$\phi_{5}(\xi)=\xi^{5}-\frac{5}{18}\xi^{3}+\frac{5}{336}\xi,\quad\textnormal{etc.}$$

The magnetic field $B_{z}$ will be approximated inside each cell by two dimensional polynomials of degree $k$ given by

$$B_{z}(\xi,\eta)=\sum_{i=0}^{N(k)-1}\alpha_{i}\Phi_{i}(\xi,\eta)\in\mathbb{P}_{k}(\xi,\eta),\qquad N(k)={\frac{1}{2}}(k+1)(k+2)$$

where the two dimensional basis functions are given by

	$\displaystyle\Phi_{i}\in\{$	$\displaystyle 1,\ \phi_{1}(\xi),\ \phi_{1}(\eta),$	$\displaystyle\textrm{at second order}$		(5)
		$\displaystyle\phi_{2}(\xi),\ \phi_{1}(\xi)\phi_{1}(\eta),\ \phi_{2}(\eta),$	$\displaystyle\textrm{at third order}$
		$\displaystyle\phi_{3}(\xi),\ \phi_{2}(\xi)\phi_{1}(\eta),\ \phi_{1}(\xi)\phi_{2}(\eta),\ \phi_{3}(\eta),$	$\displaystyle\textrm{at fourth order}$
		$\displaystyle\phi_{4}(\xi),\ \phi_{3}(\xi)\phi_{1}(\eta),\ \phi_{2}(\xi)\phi_{2}(\eta),\ \phi_{1}(\xi)\phi_{3}(\eta),\ \phi_{4}(\eta)\}$	$\displaystyle\textrm{at fifth order}$

Figure (1) shows the location of the above solution polynomials. By Gauss Theorem

$$0=\int_{C}\nabla\cdot\bm{D}\mbox{d}x\mbox{d}y=\int_{\partial C}\bm{D}\cdot\bm{n}$$

which implies that

$$(a_{0}^{+}-a_{0}^{-})\Delta y+(b_{0}^{+}-b_{0}^{-})\Delta x=0$$

(6)

The above constraint will be satisfied by the initial condition, and the update scheme we devise will ensure that it is satisfied at future times also. Note that the above constraint depends only on the face averages of the solution variables stored on the faces.

The previous paragraphs have shown us how the modes that are the primal variables in our DG scheme are collocated (for the most part) at the faces of the mesh. While this is needed in order to formulate a globally constraint-preserving scheme, we should realize that we are actually interested in the solution of a PDE. Because of the Cauchy problem, the time-evolution of Maxwell’s equations, just like the time-evolution of any hyperboloic PDE, relies on having all the spatial gradients. The reconstruction strategy that we describe below ensures that we can start with the modes at the skeleton (facial) mesh and obtain from it the variation of the electric displacement and magnetic induction at all locations on the mesh in a manner that is consistent with the involution constraint. Using the information of $D_{x}^{\pm}$, $D_{y}^{\pm}$ on the faces which are one dimensional polynomials of degree $k$, we have to reconstruct $\bm{D}=\bm{D}(\xi,\eta)\in\mathbb{V}_{k}(\xi,\eta)$ inside each cell such that the following conditions are satisfied.

1. 1.

   The normal components of $\bm{D}(\xi,\eta)$ match the known values on the faces

   $$D_{x}(\pm{\tfrac{1}{2}},\eta)=D_{x}^{\pm}(\eta),\quad\forall\eta\in[-{\tfrac{1}{2}},+{\tfrac{1}{2}}],\qquad D_{y}(\xi,\pm{\tfrac{1}{2}})=D_{y}^{\pm}(\xi),\quad\forall\xi\in[-{\tfrac{1}{2}},+{\tfrac{1}{2}}]$$

2. 2.

   The divergence of $\bm{D}$ is zero everywhere inside the cell

   $$\nabla\cdot\bm{D}(\xi,\eta)=0,\qquad\forall\xi,\eta\in[-{\tfrac{1}{2}},+{\tfrac{1}{2}}]$$

The polynomial space $\mathbb{V}_{k}$ will be chosen so that the reconstruction problem is uniquely solvable. Note that we would like to have $(k+1)$’th order accurate approximations inside the cell which implies that $\mathbb{P}_{k}(\xi,\eta)\subset\mathbb{V}_{k}(\xi,\eta)$ must be satisfied. However the space $\mathbb{V}_{k}$ must be bigger than $\mathbb{P}_{k}$ in order to be able to satisfy the matching conditions on the faces and the divergence-free condition inside the cells. The precise form of the polynomial $\bm{D}(\xi,\eta)$ and the solution of the above reconstruction problem at various orders will be explained in the next section and in Appendix. For $k=0,1,2$ (upto third order accuracy), the reconstruction problem can be solved using the information of normal components on the faces but for $k\geq 3$, we require additional information from inside the cells to solve the reconstruction problem. For $k=3$ we specify an additional cell moment $\omega_{1}$ while for $k=4$, we specify three additional cell moments, $\omega_{1},\omega_{2},\omega_{3}$, see Figure (1).

In this section, we explain how to reconstruct the field $\bm{D}$ inside the cell given the values of the normal component on the faces of the cell, and in such a way that $\nabla\cdot\bm{D}=0$. We explain the procedure for the case $k=3$ which leads to a fourth order approximation. The lower orders can be obtained from the fourth order solution and the reader can consult the Appendix. The fifth order case $(k=4)$ is also detailed in the Appendix. For solving the reconstruction problem, it is useful to note down the following results related to the 1-D orthogonal polynomials.

$$\phi_{1}(\pm{\tfrac{1}{2}})=\pm{\frac{1}{2}},\quad\phi_{2}(\pm{\tfrac{1}{2}})=\frac{1}{6},\quad\phi_{3}(\pm{\tfrac{1}{2}})=\pm\frac{1}{20},\qquad\phi_{4}(\pm{\tfrac{1}{2}})=\frac{1}{70}$$

$$\phi_{1}^{\prime}(\xi)=1,\quad\phi_{2}^{\prime}(\xi)=2\phi_{1}(\xi),\quad\phi_{3}^{\prime}(\xi)=3\phi_{2}(\xi)+\frac{1}{10},\quad\phi_{4}^{\prime}(\xi)=4\phi_{3}(\xi)+\frac{6}{35}\phi_{1}(\xi)$$

$$\phi_{5}^{\prime}(\xi)=5\phi_{4}(\xi)+\frac{5}{21}\phi_{2}(\xi)+\frac{1}{126}$$

We will assume the following polynomial form for the vector field $\bm{D}$ inside the cell

	$\displaystyle D_{x}(\xi,\eta)=$	$\displaystyle\ a_{00}+a_{10}\phi_{1}(\xi)+a_{01}\phi_{1}(\eta)+a_{20}\phi_{2}(\xi)+a_{11}\phi_{1}(\xi)\phi_{1}(\eta)+a_{02}\phi_{2}(\eta)+$
		$\displaystyle\ a_{30}\phi_{3}(\xi)+a_{21}\phi_{2}(\xi)\phi_{1}(\eta)+a_{12}\phi_{1}(\xi)\phi_{2}(\eta)+a_{03}\phi_{3}(\eta)+a_{40}\phi_{4}(\xi)+$
		$\displaystyle\ a_{31}\phi_{3}(\xi)\phi_{1}(\eta)+a_{22}\phi_{2}(\xi)\phi_{2}(\eta)+a_{13}\phi_{1}(\xi)\phi_{3}(\eta)$
	$\displaystyle D_{y}(\xi,\eta)=$	$\displaystyle\ b_{00}+b_{10}\phi_{1}(\xi)+b_{01}\phi_{1}(\eta)+b_{20}\phi_{2}(\xi)+b_{11}\phi_{1}(\xi)\phi_{1}(\eta)+b_{02}\phi_{2}(\eta)+$
		$\displaystyle\ b_{30}\phi_{3}(\xi)+b_{21}\phi_{2}(\xi)\phi_{1}(\eta)+b_{12}\phi_{1}(\xi)\phi_{2}(\eta)+b_{03}\phi_{3}(\eta)+b_{31}\phi_{3}(\xi)\phi_{1}(\eta)+$
		$\displaystyle\ b_{22}\phi_{2}(\xi)\phi_{2}(\eta)+b_{13}\phi_{1}(\xi)\phi_{3}(\eta)+b_{04}\phi_{4}(\eta)$

Note that $D_{x}$ has the form of a polynomial $\mathbb{P}_{4}(\xi,\eta)$ except that the terms corresponding to $\eta^{4}$ are not included. Similarly, $D_{y}$ belongs to $\mathbb{P}_{4}(\xi,\eta)$ except for the term $\xi^{4}$ which is not included. Such polynomial spaces to approximate vector fields in a divergence conforming manner were introduced in [22] and are called BDFM polynomials. Both the components completely include $\mathbb{P}_{3}(\xi,\eta)$ polynomials. Matching the cell solution to the face solution, we get the following 16 equations

	$\displaystyle a_{00}\pm{\tfrac{1}{2}}a_{10}+\tfrac{1}{6}a_{20}\pm\tfrac{1}{20}a_{30}+\tfrac{1}{70}a_{40}\;$	$\displaystyle=a_{0}^{\pm}$
	$\displaystyle a_{01}\pm{\tfrac{1}{2}}a_{11}+\tfrac{1}{6}a_{21}\pm\tfrac{1}{20}a_{31}\;$	$\displaystyle=a_{1}^{\pm}$
	$\displaystyle a_{02}\pm{\tfrac{1}{2}}a_{12}+\tfrac{1}{6}a_{22}\;$	$\displaystyle=a_{2}^{\pm}$
	$\displaystyle a_{03}\pm{\tfrac{1}{2}}a_{13}\;$	$\displaystyle=a_{3}^{\pm}$
	$\displaystyle b_{00}\pm{\tfrac{1}{2}}b_{01}+\tfrac{1}{6}b_{02}\pm\tfrac{1}{20}b_{03}+\tfrac{1}{70}b_{04}\;$	$\displaystyle=b_{0}^{\pm}$
	$\displaystyle b_{10}\pm{\tfrac{1}{2}}b_{11}+\tfrac{1}{6}b_{12}\pm\tfrac{1}{20}b_{13}\;$	$\displaystyle=b_{1}^{\pm}$
	$\displaystyle b_{20}\pm{\tfrac{1}{2}}b_{21}+\tfrac{1}{6}b_{22}\;$	$\displaystyle=b_{2}^{\pm}$
	$\displaystyle b_{30}\pm{\tfrac{1}{2}}b_{31}\;$	$\displaystyle=b_{3}^{\pm}$

The divergence of the vector field $\bm{D}$ is a polynomial of degree 3 and making it zero inside the cell yields the following set of 10 equations

	$\displaystyle(a_{10}+\tfrac{1}{10}a_{30})\Delta y+(b_{01}+\tfrac{1}{10}b_{03})\Delta x$	$\displaystyle=0$
	$\displaystyle(2a_{20}+\tfrac{6}{35}a_{40})\Delta y+(b_{11}+b_{13}/10)\Delta x$	$\displaystyle=0$
	$\displaystyle(a_{11}+a_{31}/10)\Delta y+(2b_{02}+\tfrac{6}{35}b_{04})\Delta x$	$\displaystyle=0$
	$\displaystyle 3a_{30}\Delta y+b_{21}\Delta x$	$\displaystyle=0$
	$\displaystyle 2a_{21}\Delta y+2b_{12}\Delta x$	$\displaystyle=0$
	$\displaystyle a_{12}\Delta y+3b_{03}\Delta x$	$\displaystyle=0$
	$\displaystyle 4a_{40}\Delta y+b_{31}\Delta x$	$\displaystyle=0$
	$\displaystyle 3a_{31}\Delta y+2b_{22}\Delta x$	$\displaystyle=0$
	$\displaystyle 2a_{22}\Delta y+3b_{13}\Delta x$	$\displaystyle=0$
	$\displaystyle a_{13}\Delta y+4b_{04}\Delta x$	$\displaystyle=0$

The first equation in the above set is redundant since it is contained in the other equations due to the constraint (6). Ignoring this equation, we can solve for some of the coefficients $a_{ij}$, $b_{ij}$ in terms of the face solution as follows:

$\displaystyle a_{00}=$ $\displaystyle{\tfrac{1}{2}}(a_{0}^{-}+a_{0}^{+})+\tfrac{1}{12}(b_{1}^{+}-b_{1}^{-})\tfrac{\Delta x}{\Delta y}$ $\displaystyle a_{10}=$ $\displaystyle a_{0}^{+}-a_{0}^{-}+\tfrac{1}{30}(b_{2}^{+}-b_{2}^{-})\tfrac{\Delta x}{\Delta y}$ $\displaystyle a_{20}=$ $\displaystyle-{\tfrac{1}{2}}(b_{1}^{+}-b_{1}^{-})\tfrac{\Delta x}{\Delta y}+\tfrac{3}{140}(b_{3}^{+}-b_{3}^{-})\tfrac{\Delta x}{\Delta y}$ $\displaystyle a_{30}=$ $\displaystyle-\tfrac{1}{3}(b_{2}^{+}-b_{2}^{-})\tfrac{\Delta x}{\Delta y}$ $\displaystyle a_{03}=$ $\displaystyle{\tfrac{1}{2}}(a_{3}^{-}+a_{3}^{+})$ $\displaystyle a_{12}=$ $\displaystyle a_{2}^{+}-a_{2}^{-}$ $\displaystyle a_{13}=$ $\displaystyle a_{3}^{+}-a_{3}^{-}$ $\displaystyle a_{40}=$ $\displaystyle-\tfrac{1}{4}(b_{3}^{+}-b_{3}^{-})\tfrac{\Delta x}{\Delta y}$   $\displaystyle b_{00}=$ $\displaystyle{\tfrac{1}{2}}(b_{0}^{-}+b_{0}^{+})+\tfrac{1}{12}(a_{1}^{+}-a_{1}^{-})\tfrac{\Delta y}{\Delta x}$ $\displaystyle b_{01}=$ $\displaystyle b_{0}^{+}-b_{0}^{-}+\tfrac{1}{30}(a_{2}^{+}-a_{2}^{-})\tfrac{\Delta y}{\Delta x}$ $\displaystyle b_{02}=$ $\displaystyle-{\tfrac{1}{2}}(a_{1}^{+}-a_{1}^{-})\tfrac{\Delta y}{\Delta x}+\tfrac{3}{140}(a_{3}^{+}-a_{3}^{-})\tfrac{\Delta y}{\Delta x}$ $\displaystyle b_{30}=$ $\displaystyle{\tfrac{1}{2}}(b_{3}^{-}+b_{3}^{+})$ $\displaystyle b_{03}=$ $\displaystyle-\tfrac{1}{3}(a_{2}^{+}-a_{2}^{-})\tfrac{\Delta y}{\Delta x}$ $\displaystyle b_{21}=$ $\displaystyle b_{2}^{+}-b_{2}^{-}$ $\displaystyle b_{31}=$ $\displaystyle b_{3}^{+}-b_{3}^{-}$ $\displaystyle b_{04}=$ $\displaystyle-\tfrac{1}{4}(a_{3}^{+}-a_{3}^{-})\tfrac{\Delta y}{\Delta x}$

We are now in a position to explain the constraint preserving scheme. Recall that we have several solution polynomials, some of which are stored on the faces and some are stored inside the cells. The basic solution variables have been illustrated in Figure (1). The solution polynomials $D_{x}(\xi,\eta)$, $D_{y}(\xi,\eta)$ are not independent and are obtained by the divergence-free reconstruction process described in section (4) and in the Appendices.

1. 1.

   The normal component of $\bm{D}$ stored on the faces will be evolved by a flux reconstruction scheme applied on each face.

2. 2.

   At fourth and fifth orders, we have additional quantities $\omega_{i}$ which are located inside the cells and we will devise a DG scheme for these quantities.

3. 3.

   The magnetic flux density has only one component which is stored inside the cells and this will be evolved by a standard DG scheme.

The FR and DG schemes require some numerical fluxes that are obtained from 1-D and 2-D Riemann problems. We give a short summary of these numerical fluxes in the Appendices (C.1)-(C.2).