An Explicit Mapped Tent Pitching Scheme for Maxwell Equations

TL;DR

This paper introduces a novel explicit numerical scheme for Maxwell equations using unstructured spacetime tents, enabling high-order accuracy, variable time steps, and efficient parallel computation on modern architectures.

Contribution

It develops a new mapped tent pitching method that overcomes limitations of standard explicit schemes for hyperbolic PDEs, allowing flexible, high-order, and parallelizable solutions.

Findings

Enables explicit high-order solutions on unstructured meshes.

Supports variable time steps and local mesh refinement.

Achieves efficient parallelization suitable for modern architectures.

Abstract

We present a new numerical method for solving time dependent Maxwell equations, which is also suitable for general linear hyperbolic equations. It is based on an unstructured partitioning of the spacetime domain into tent-shaped regions that respect causality. Provided that an approximate solution is available at the tent bottom, the equation can be locally evolved up to the top of the tent. By mapping tents to a domain which is a tensor product of a spatial domain with a time interval, it is possible to construct a fully explicit scheme that advances the solution through unstructured meshes. This work highlights a difficulty that arises when standard explicit Runge Kutta schemes are used in this context and proposes an alternative structure-aware Taylor time-stepping technique. Thus explicit methods are constructed that allow variable time steps and local refinements without compromising high order accuracy in space and time. These Mapped Tent Pitching (MTP) schemes lead to highly parallel algorithms, which utilize modern computer architectures extremely well.

Tables

Table 1. Table 1: Number of degrees of freedom and simulation times for spatial polynomial orders p = 2 , 3 𝑝 2 3 p=2,3 . This data was generated using a shared memory server with 4 E7-8867 CPUs with 16 cores each.

	$p=2$	$p=3$
number of spatial dof	$2.938\times {10}^7$	$5.875\times {10}^7$
number of spacetime dof per slab	$1.908\times {10}^9$	$7.632\times {10}^9$
simulation time per slab	4.6 s	49.2 s
total simulation time	20 min	3 h 33 min

Equations

$${\partial }_t\epsilon E=\nabla \times H,{\partial }_t\mu H=-\nabla \times E,$$

$$\phi (x,\hat{t})=(1-\hat{t}){\phi }_b(x)+\hat{t}{\phi }_t(x).$$

$$\operatorname{skew}E=\left(\begin{array}[]{ccc}0&E_{z}&-E_{y}\\ -E_{z}&0&E_{x}\\ E_{y}&-E_{x}&0\end{array}\right)\;,$$

$$g(u)=\begin{bmatrix}\varepsilon E\\ \mu H\end{bmatrix}\;,\qquad f(u)=\left[\begin{array}[]{r}-\operatorname{skew}H\\ \operatorname{skew}E\end{array}\right]\;.$$

$$F(u)=\left[f(u)\;\;g(u)\right]=\left[\begin{array}[]{rc}-\operatorname{skew}H&\;\varepsilon E\\ \operatorname{skew}E&\;\mu H\end{array}\right]\;,$$

$${\mathrm{div}}_{x,t}F(u)=0.$$

$${\Phi }^{\prime }=\left[\begin{array}{cc}I& 0\\ \nabla {\phi }^T& \delta \end{array}\right]\text{and}({\Phi }^{\prime }{)}^{-T}=\left[\begin{array}{cc}I& -{\delta }^{-1}\nabla \phi \\ 0& {\delta }^{-1}\end{array}\right].$$

$${\partial }_{\hat{t}}(g(\hat{u})-f(\hat{u})\nabla \phi )+{\mathrm{div}}_x(\delta f(\hat{u}))=0,$$

$$\int_{\omega_{V}}\partial_{\hat{t}}\big{[}g(\hat{u})-f(\hat{u})\nabla\varphi\big{]}\,v_{h}-\sum_{T\subset\omega_{V}}\int_{T}\delta f(\hat{u})\nabla v_{h}+\sum_{F\subset\omega_{V}}\int_{F}\delta f_{n}(\hat{u}^{+},\hat{u}^{-})\llbracket v\rrbracket=0$$

$$f_{n}(\hat{u}^{+},\hat{u}^{-})=\left[\begin{array}[]{r}\{\hat{H}\}\times n+\llbracket\hat{E}_{t}\rrbracket\\ -\{\hat{E}\}\times n+\llbracket\hat{H}_{t}\rrbracket\end{array}\right]\;,$$

$$\frac{d}{d\hat{t}}\left[MU\right](\hat{t})-AU(\hat{t})=0,\hat{t}\in (0,1),$$

$$\frac{d}{d\hat{t}}Y-A{M}^{-1}Y=0,$$

$$\frac{d}{d\hat{t}}Y=AU,Y=MU.$$

$$M(\hat{t})=M_i+(\hat{t}-{\hat{t}}_i)M^{\prime },\hat{t}\in ({\hat{t}}_i,{\hat{t}}_{i+1})$$

$$Y_i(\hat{t})=\sum _{n=0}^q\frac{(\hat{t}-{\hat{t}}_i{)}^n}{n!}{Y}_{i,n}{U}_i(\hat{t})=\sum _{n=0}^{q-1}\frac{(\hat{t}-{\hat{t}}_i{)}^n}{n!}{U}_{i,n},\hat{t}\in ({\hat{t}}_i,{\hat{t}}_{i+1}),$$

$$Y^{(n+1)}(\hat{t})$$

$$Y^{(n)}(\hat{t})$$

$$\begin{array}[]{r c l @{\qquad}l}Y_{i,n}&=&AU_{i,n-1}\;,&1\leq n\leq q\;,\\ M_{i}U_{i,n}&=&Y_{i,n}-nM^{\prime}U_{i,n-1}\;,&1\leq n\leq q-1\;,\end{array}$$

$$Y_{i,0}$$

$${\int }_{{\omega }_V}[g(\hat{u})-f(\hat{u})\nabla \phi ]\cdot \hat{v}={\int }_{{\omega }_V}({\hat{w}}_E,{\hat{w}}_H)\cdot \hat{v},\text{ for all }\hat{v}\in V_h.$$

$$a(\hat{u},\hat{u})$$

$$\hat{E}$$

$$\hat{H}$$

$${\partial }_t\epsilon E_z={\partial }_x{H}_y-{\partial }_y{H}_x,{\partial }_t\mu H_x=-{\partial }_y{E}_z,{\partial }_t\mu H_y={\partial }_x{E}_z,$$

$$E_z$$

$$H_x$$

$$H_y$$

$$\sum _{i=1}^{N_{\mathrm{tents}}}{N}_{\mathrm{dof},i}mq=\left(\sum _{i=1}^{N_{\mathrm{tents}}}{N}_{\mathrm{dof},i}\right)2p(p+1).$$

Full text

1 Introduction

Electromagnetic waves propagate at the speed of light. Thus, the field at a certain point in space and time depends only on field values within a dependency cone. A tent pitching method introduces a special “causal” spacetime mesh that respects this finite speed of propagation. It is not limited to Maxwell equations, but can be applied to general hyperbolic equations. A tent pitching method requires a numerical scheme to discretize the equation on that mesh. Discontinuous Galerkin (DG) methods are of particular interest since they offer a systematic avenue to build high order methods. For a given initial condition at the bottom of a tent, the discrete equations may be solved within each individual tent, up to the tent top. The computed solution at the tent top provides initial conditions for the tents that follow later in time. This method is highly parallel, since many tents can be solved independently. Methods using such tent-pitched meshes may be traced back to LowriRoeLeer95 ; Richt94 . More recent works AbediPetraHaber06 ; MR05 ; ASHT00 develop Spacetime DG (SDG) methods within tents by formulating local variational problems, for which linear systems are set up and solved. Although these systems are local, the matrix size can grow rapidly with the polynomial order, especially in four-dimensional spacetime tents. In this context it is natural to ask if one can develop explicit schemes (which usually perform well under low memory bandwidth) that take advantage of tents.

A key ingredient to answer this question was presented in GSW17 , where Mapped Tent Pitching (MTP) schemes were introduced. The MTP discretization, which proceeds by mapping tents to a spacetime cylinder, allows one to evolve the solution either implicitly or explicitly within tents. The memory requirements of the explicit MTP scheme are limited to what is needed for storing the spatial mesh, the solution coefficients at one time step, and the topology of the tents.

In this work, we show that notwithstanding the above-mentioned advantages of the explicit MTP scheme, one may lose higher order convergence if a naive time stepping strategy (involving a standard explicit Runge-Kutta scheme) is used. We then develop a new Taylor time-stepping for the local problems within tents. Despite its simplicity, our numerical experiments show that it delivers optimal order of convergence.

2 Mesh generation by tent pitching

We start with a conforming spatial mesh consisting of elements ${\mathcal{T}}=\{T\}$ and vertices ${\mathcal{V}}=\{V\}$. We progress in time by defining a sequence of advancing fronts $\tau_{i}$. A front $\tau_{i}$ is given as a standard nodal finite element function on this mesh. It is defined by storing the current time for every vertex of the mesh. We move from $\tau_{i}$ to the next front $\tau_{i+1}$ by moving one vertex forward in time, while keeping all other vertices fixed. The spacetime domain between $\tau_{i}$ and $\tau_{i+1}$ we call a tent. In Fig. 1, the red domain is the tent between $\tau_{i}$ and $\tau_{i+1}$.

Its projection to the spatial domain is exactly the vertex patch $\omega_{V}$ around $V$ of the original mesh. The data to be stored for one tent are the bottom and top-times of the central vertex, plus the times for all neighboring vertices.

Note that although the algorithm is described sequentially, it is highly parallel. Vertices with graph-distance of at least two can be moved forward independently. For example, in Fig. 1, all blue tents can be built and processed in parallel.

The distance for advancing a vertex is limited by the speed of light, a constraint often referred to in the literature as the causality condition. Under this condition, the Maxwell problem inside the tent is solvable using the initial conditions at the tent bottom. Thus, the top boundary is an outgoing boundary and no boundary conditions are needed there.

Note that the spatial mesh is refined towards the right boundary, which leads to smaller tent heights at the right boundary. Hence, smaller time steps in locally refined regions is a very natural feature of tent pitching methods.

3 The MTP discretization

Now, we consider the discretization method for one tent domain $K=\{(x,t):x\in\omega_{V},\varphi_{b}(x)\leq t\leq\varphi_{t}(x))\},$ where $\omega_{V}$ is the union of elements containing the vertex $V$, and $\varphi_{b}$ and $\varphi_{t}$ are the bottom and top fronts, respectively, restricted to $\omega_{V}$. Our aim is to numerically solve the Maxwell system on $K$, namely

[TABLE]

where boundary values for both fields are given at the tent bottom and $\nabla=\nabla_{x}$ denotes the spatial gradient.

The approach of MTP schemes is to map the tent domain to a spacetime cylinder $\omega_{V}\times(0,1)$ and solve the transformed equation there. The transformation from the cylinder to the tent is denoted by $\Phi:\omega_{V}\times(0,1)\rightarrow K$ and is defined by $\Phi(x,\hat{t})=(x,\varphi(x,\hat{t}))$ where

[TABLE]

It is similar to the Duffy transformation mapping a square to a triangle.

With the notation

[TABLE]

we can rephrase the curl operator as $\nabla\times E=\operatorname{div}\operatorname{skew}E,$ where the divergence of the matrix function is taken row-wise. To simplify notation further, we define $u:K\rightarrow{{\mathbb{R}}}^{6}$ by $u=(E,H),$ and set $g:K\rightarrow{\mathbb{R}}^{6}$ and $f:K\rightarrow{\mathbb{R}}^{6\times 3}$ by

[TABLE]

Then (1) may be rewritten as the conservation law $\partial_{t}g(u)+\operatorname{div}_{x}f(u)=0.$ Furthermore, we define $F(u)\in{\mathbb{R}}^{6\times 4}$ as

[TABLE]

which allows us to write Maxwell’s system (1) as the spacetime conservation law

[TABLE]

For each row of $F$, the spacetime divergence $\operatorname{div}_{x,t}$ sums the spatial divergence of the first three components with the time-derivative of the last component.

Now, we apply the Piola transformation to pull back $F$ from the tent $K$ to the cylinder using the mapping $\Phi$. The derivative of $\Phi$ and its transposed inverse are

[TABLE]

The Piola transform of $F$ is $\widehat{F}(\hat{u})={\mathcal{P}}\{F\}=(\det\Phi^{\prime})(F\circ\Phi)(\Phi^{\prime})^{-T}$ with $\hat{u}=u\circ\Phi$. Since the Piola transform provides an algebraic transformation of the divergence, equation (3) is simply transformed to $\operatorname{div}_{x,\hat{t}}\widehat{F}(\hat{u})=0$ on the spacetime cylinder. Then, inserting the Jacobian of $\Phi$ leads us to the transformed equation

[TABLE]

where $\delta(x)=\varphi_{t}(x)-\varphi_{b}(x)$ is the local height of the tent. Note that $\nabla\varphi$ is an affine-linear function in quasi-time $\hat{t}$. Equation (4) describes the evolution of $\hat{u}$ along quasi-time from $\hat{t}=0$ to $\hat{t}=1$. Details of the calculations are given in GSW17 .

The next step is the space discretization of (4) by a standard discontinuous Galerkin method. Let $V_{h}\subset[L_{2}]^{6}$ be the DG finite element space of degree $p$ on $\mathcal{T}$. On each tent we search for $\hat{u}:[0,1]\rightarrow V_{h}$ such that

[TABLE]

holds for all $v_{h}\in V_{h}$ and all $\hat{t}\in[0,1]$. Only the restriction of $V_{h}$ on the patch $\omega_{V}$ is used in this equation. The numerical flux $f_{n}(\hat{u}^{+},\hat{u}^{-})$ depends on the positive trace $\lim_{s\rightarrow 0^{+}}\hat{u}(x+sn)$ and negative trace $\lim_{s\rightarrow 0^{+}}\hat{u}(x-sn)$, where $n$ is a unit normal vector of arbitrary orientation to the face. The jump is defined as usual by $\llbracket\hat{u}\rrbracket:=\hat{u}^{+}-\hat{u}^{-}$ and the mean value by $\{\hat{u}\}:=\tfrac{1}{2}(\hat{u}^{+}+\hat{u}^{-})$. One example is the upwind flux (HW08, , p. 434)

[TABLE]

with the tangential components $\hat{E}_{t}=-(\hat{E}\times n)\times n$ and $\hat{H}_{t}=-(\hat{H}\times n)\times n$ of $\hat{E}=E\circ\Phi$ and $\hat{H}=H\circ\Phi$. Note that the local tent height $\delta$ enters the boundary integrals as a multiplicative factor. At the outer boundary of the vertex patch we have $\delta=0$, so the facet integrals on the outer boundary disappear. For the above semidiscrete system, initial values for the tent problem are given finite element functions at the tent bottom. The finite element solution on the tent top provides the initial conditions for the next level tent. Therefore, no projection of initial values is needed when propagating from one tent to the next.

After the semi-discretization, as usual, we are left to solve a system of $N=\dim V_{h}(\omega_{V})$ ordinary differential equations for $U:[0,1]\rightarrow{\mathbb{R}}^{N}$,

[TABLE]

given $U(0)$. The non-standard feature of (5) is that $M$ is an affine-linear function of the quasi-time $\hat{t}$ (since our mapping enters the mass matrix $M$ through $\nabla\varphi$). The matrix $A$ is independent of $\hat{t}$. A straightforward approach is to substitute $Y=MU$ and solve

[TABLE]

instead of (5). Although first order convergence was observed with this strategy, further numerical studies showed reduced order of convergence if the stage-order of the Runge Kutta (RK) method is not high enough – see Fig. 3 (right). While the implicit MTP schemes discussed in GSW17 do not show this problem, the issue remains critical for explicit schemes. Thus, we propose to use a new type of explicit time-stepping for time discretization, discussed next.

4 Structure-aware Taylor time-stepping

Returning to the ordinary differential equation (5) and continuing to make the substitution $Y=MU$, we now reconsider the previous equation as the following differential-algebraic system:

[TABLE]

We begin by subdividing the interval $(0,1)$ into $m\in\mathbb{N}$ smaller intervals of size $\frac{1}{m}$, defined by $(\hat{t}_{i},\hat{t}_{i+1})=(\frac{i}{m},\frac{i+1}{m})$, for $i\in\mathbb{N}$ and $0\leq i\leq m-1$. Recall that $A$ is independent of quasi-time $\hat{t}$, and $M$ is an affine function of $\hat{t}$, i.e.,

[TABLE]

where $M_{i}=M(\hat{t}_{i})$ and the derivative $M^{\prime}$ is a constant matrix. We want to design a time-stepping scheme that is aware of this structure.

Consider the approximations to $Y,U$ on $(\hat{t}_{i},\hat{t}_{i+1})$ in the form of Taylor polynomials $Y_{i},U_{i}$ of degree $q$, defined by

[TABLE]

where $Y_{i,n}=Y_{i}^{(n)}(\hat{t}_{i})$ and $U_{i,n}=U_{i}^{(n)}(\hat{t}_{i}).$ To find these derivatives, we differentiate both equations of (6) $n$ times to get

[TABLE]

For the second equation we used Leibnitz’ formula $(fg)^{(n)}=\sum_{i=0}^{n}{n\choose i}f^{(i)}g^{(n-i)}$, and the fact that $M$ is affine-linear. Evaluating these equations for the Taylor polynomials $Y_{i},U_{i}$ at $\hat{t}=\hat{t}_{i}$, we obtain a recursive formula for $Y_{i,n}$ and $U_{i,n}$ in terms of $U_{i,n-1}$, namely

[TABLE]

for all $0\leq i\leq m-1$. Given $Y_{0,0}=Y(\hat{t}_{0})$, $M_{0}U_{0,0}=Y_{0,0}$, applying (8) with $i=0$ gives the approximate functions $Y_{0}(\hat{t}),U_{0}(\hat{t})$ in the first subinterval $(\hat{t}_{0},\hat{t}_{1})$. The recursive formulas are initiated for later subintervals at $n=0$ by

[TABLE]

After the final subinterval, we get $Y_{m-1}(t_{m})$, our approximation to $Y(1)$. We shall refer to the new time-stepping scheme generated by (8) as the $q$-stage SAT (structure-aware Taylor) time-stepping.

Note that $Y_{m-1}(t_{m})$ is our approximation to $Y=MU$ at the top of the tent. This value is then passed to the next tent in time. The time dependence of $M$ arises from the time dependence of $\nabla\varphi$. This gradient is continuous along spacetime lines of constant spatial coordinates. Therefore, when passing from one element of a tent to the same element within the next tent in time, $Y$ is continuous (since the solution $U$ is continuous). Of course, on flat fronts $\nabla\varphi=\nabla\tau=0$, so there $M$ is just a diagonal matrix containing the material parameters.

To briefly remark on the expected convergence rate of a $q$-stage SAT time-stepping, recall that due to the mapping of the MTP method we solve for $\hat{u}=u\circ\Phi$, which satisfies $\partial^{n}_{\hat{t}}\hat{u}=\delta^{n}(\partial^{n}_{t}u)\circ\Phi$. The causality condition implies that $\delta\rightarrow 0$ if the mesh size $h\rightarrow 0$. Thus we may expect the $n^{th}$ temporal derivative of $\hat{u}$, and correspondingly $U^{(n)},$ to go to zero at the rate $\mathcal{O}(h^{n})$. By using a $q$-stage SAT time-stepping, we approximate the first $q-1$ terms of the exact Taylor expansion of $U$. Thus we expect the convergence rate to be $O(h^{q})$, the size of the remainder term involving $U^{(q)}$. The next section provides numerical evidence for this.

Before concluding this section, we should note that in (8) and (9), we tacitly assumed that $M_{i}$ is invertible. Let us show that this is indeed the case whenever the causality condition (see §2) $|\nabla\varphi|<\sqrt{\varepsilon\mu}$ is fulfilled. At any quasi-time $\hat{t}$, given a $\hat{w}=(\hat{w}_{E},\hat{w}_{H})\in V_{h}$ whose coefficient vector in the basis expansion is $W\in\mathbb{R}^{N}$, consider the equation $M(\hat{t})U=W$ for the coefficient vector $U$ of $\hat{u}\in V_{h}$. This equation, in variational form, is

[TABLE]

Let $a(\hat{u},\hat{v})$ denote the left hand side of (10). To prove solvability of (10), it suffices to prove that $a(\cdot,\cdot)$ is a coercive bilinear form on $[L_{2}]^{6}$ for any $\hat{t}$. By inserting $g(\hat{u})=[\varepsilon\hat{E},\mu\hat{H}]^{T}$ and $f(\hat{u})=[-\operatorname{skew}\hat{H},\operatorname{skew}\hat{E}]^{T}$ into $a(\hat{u},\hat{u})$,

[TABLE]

where we used the Cauchy-Schwarz inequality and inserted $\sqrt{\varepsilon}$ and $\sqrt{\mu}$ to achieve the desired scaling. By applying Young’s inequality and $|\nabla\varphi|<\sqrt{\varepsilon\mu}$,

[TABLE]

form some constant $C>0$. Thus $M_{i}$ is invertible and the SAT time-stepping is well defined on all tents respecting the causality condition.

One may exploit the specific details of the Maxwell problem to avoid the assembly and the inversion of matrices $M_{i}$ (as we have done in our implementation). In fact, instead of (10), we can explicitly solve the corresponding exact undiscretized equation obtained by replacing $V_{h}$ by $[L_{2}]^{6}$ in (10). The solution $\hat{u}=(\hat{E},\hat{H})$ in closed form reads

[TABLE]

We then perform a projection of these into $V_{h}$ to obtain the coefficients $U(\hat{t}_{i})$. For uncurved elements, this just involves the inversion of a diagonal mass matrix. For the small number of curved elements, we use a highly optimized algorithm which uses an approximation instead of the exact inverse mass matrix.

5 Numerical Results

The MTP discretization in combination with the SAT time-stepping on tents is implemented within the Netgen/NGSolve finite element library. In this section numerical results concerning accuracy as well as performance are reported.

5.1 Convergence studies in two space dimensions

We consider the model problem in two space dimensions

[TABLE]

on the spacetime cube $[0,\pi]^{2}\times[0,\sqrt{2}\pi]$. Parameters are set $\varepsilon=\mu=1$ such that speed of light is $c=1$. Initial and boundary values are set such that the exact solution is given by

[TABLE]

Based on a spatial mesh with mesh size $h$, we generate a tent pitched mesh such that the maximal slope $|\nabla\varphi|$ is bounded by $(2c)^{-1}$ and apply a discontinuous Galerkin method in space using polynomials of order $p$, with $1\leq p\leq 4$. On each cylinder we perform a $(p+1)$-stage SAT time-stepping with $m=2p$ intervals. The spatial $L_{2}$ error of all field components at the final time is reported in the left plot of Fig. 3. We observe that the error goes to zero at the optimal rate of $\mathcal{O}(h^{p+1})$ until we are close to machine precision.

In contrast, the right plot in Fig. 3 illustrates the previously mentioned loss of convergence rates when the classical Runge-Kutta method is used. The convergence rates stagnate at first order no matter what $p$ is used. A similar behavior was also observed for other explicit Runge-Kutta methods.

5.2 Large scale problem in three space dimensions

As a second example we present a simulation on a domain similar to the resonator shown in HPSTW15 . The geometry is given as body of revolution of smooth B-spline curves. The mesh consisting of 489593 curved tetrahedral elements is shown in Fig. 4. Due to higher curvature the mesh is refined along the inner roundings, where the ratio of the largest to the smallest element is approximately 5:1. We used a Gaussian peak (located at the axis of revolution and the position of the fifth inner rounding) for the electric field as initial data. The explicit MTP scheme with SAT time-stepping then computed the solution at $t=260$ using time slabs of height 1, with each slab composed of $N_{\mathrm{tents}}=149072$ tents. On each tent we used a $(p+1)$-stage SAT time-stepping with $m=2p$ intervals, where $p$ denotes the spatial polynomial order. With the spatial degrees of freedom $N_{\mathrm{dof},i}$ of the $i^{th}$ tent and the number of stages $q=p+1$, we obtain the total spacetime degrees of freedom per time slab

[TABLE]

The corresponding numbers of degrees of freedom and the simulation times are shown in Table 1. In HPSTW15 a similar problem is solved using a discontinuous Galerkin method with quadratic elements, combined with a polynomial Krylov subspace method in time. Using 96 cores it took them 7:10 hours to reach the final time. Our simulation with polynomial order $p=3$, which has a comparable number of unknowns, took 3:33 hours on 64 cores. This significant speed up is an illustration of the capability of the new method. The $H_{y}$ component of the obtained solution at $t=260$, using third order polynomials in space, is shown in Fig. 4.

Acknowledgements.

This work was supported in part by the National Science Foundation.