Finite Element, Discontinuous Galerkin, and Finite Difference Evolution Schemes in Spacetime
aa r X i v : . [ g r- q c ] J u l Finite Element, Discontinuous Galerkin, and FiniteDifference Evolution Schemes in Spacetime
G Zumbusch
Institut f¨ur Angewandte Mathematik, Friedrich-Schiller-Universit¨at Jena,07743 Jena, GermanyE-mail: [email protected]
Abstract.
Numerical schemes for Einstein’s vacuum equation are developed.Einstein’s equation in harmonic gauge is second order symmetric hyperbolic. It isdiscretized in four-dimensional spacetime by Finite Differences, Finite Elements, andInterior Penalty Discontinuous Galerkin methods, the latter related to Regge calculus.The schemes are split into space and time and new time-stepping schemes for waveequations are derived. The methods are evaluated for linear and non-linear testproblems of the Apples-with-Apples collection.PACS numbers: 04.25.D, 02.70.Bf, 02.70.Dh, 04.20.Fy
1. Introduction
Numerical methods for the solution of Einstein’s equation in general relativity aremainly based on Finite Differences (FD) and Pseudo-Spectral-Collocation (Bonazzolaet al. 2004, Boyle et al. 2007) schemes space so far. The Finite Element method (FEM),or more generally Galerkin schemes have been used for reduced or auxiliary problemsin numerical relativity (Arnold et al. 1998, Metzger 2004, Sopuerta et al. 2006, Aksoyluet al. 2008, Field et al. 2009). However, Galerkin methods are heavily used for thesolution of wave problems in areas like acoustic and electro-magnetic scattering andelastic waves (Cohen 2002). This is mainly due to their way to deal with heterogeneousmedia and arbitrarily shaped geometric objects, represented by unstructured grids.Furthermore, the convergence theory of Galerkin methods is based on lower regularity(differentiability) requirements than Finite Differences and spectral methods.General relativity is governed by Einstein’s equation, which can be written as asystem of second order partial differential equations in spacetime. In order to define awell-posed initial-value (Cauchy) problem, additional gauge conditions are needed. Forthe numerical solution of the system, spacetime is usually split into space and time 3 + 1and finally a time-stepping scheme is derived. Using a lapse- and a shift-function, asequence of space-like manifolds is constructed, which fixes the gauge freedom. Thereare many improvements of the original ADM (Arnowitt et al. 1962, York 1979) splittinglike BSSN (Shibata & Nakamura 1995, Baumgarte & Shapiro 1999). The equations are
EM, DG, and FD Evolution Schemes in Spacetime ∂ tt u = ∆ u as an example problem is written invariational form as R Ω ( ∂ tt u ) w d x = − R Ω ( ∇ u ) · ( ∇ w ) d x ∀ w with trial functions w , integration over the spatial domain Ω, and zero boundaryconditions. This gives rise to FEM (Dupont 1973, Baker & Bramble 1979) and DG(Ainsworth et al. 2006, Grote et al. 2006, Field et al. 2009) in space schemes, used inconjunction with a standard time integrator like the leapfrog scheme. The first order intime formulation ∂ t v = ∆ u and ∂ t u = v in variational version in time reads as − R T v ( ∂ t w ) dt = R T (∆ u ) w dt ∀ w − R T u ( ∂ t w ) dt = R T vw dt ∀ w on the interval T and without initial value terms. In order to obtain a time-stepping scheme, a time-discontinuous Galerkin method can be constructed (Jamet1978, Eriksson et al. 1985). Note that time continuous functions do not lead to atime-stepping scheme, but a single large equation system for all times. We can combineboth Galerkin schemes to a spacetime FEM like R Ω × T v ( ∂ t w ) dt d x = R Ω × T ( ∇ u ) · ( ∇ w ) dt d x ∀ w − R Ω × T u ( ∂ t w ) dt d x = R Ω × T vw dt d x ∀ w , continuous (French & Peterson 1996, Anderson & Kimn 2007) and discontinuous(Hulbert & Hughes 1990, Monk & Richter 2005) in time. In this paper, however, wewill consider second order in space and time formulations of type R Ω × T ( ∂ t u )( ∂ t w ) dt d x = R Ω × T ( ∇ u ) · ( ∇ w ) dt d x ∀ w , (1) EM, DG, and FD Evolution Schemes in Spacetime R Ω uw d x each time step. The DG method in space iscomputationally more efficient than FEM in general, because the mass matrix is block-diagonal and the equation systems are easier to solve. However, by a special choice ofnumerical quadrature rules (mass-lumping) in FEM, see (Cohen 2002), and a choice oforthogonal ansatz functions in DG, see (Rivi`ere 2008), the mass matrix is diagonal andthe equation systems are trivial to solve.Now we put together the variational formulation of Einstein’s equation and thespacetime Galerkin schemes and we obtain in section 3.6, as the main result, a FEM, asymmetric and a non-symmetric Interior Penalty DG method for Einstein’s full vacuumequation. As an intermediate step we briefly discuss a simpler, linearized version ofEinstein’s equation.Memory requirements for nodal FD and piecewise linear FEM schemes for Einstein’sequation are comparable, namely ten metric component values per grid node. The DGmethods need this storage of 10 values for each element and each ansatz function, i.e.10 · ·
16 for linear or multi-linear functions, thus are more memory intensive.The fields are needed for two previous and the current time-slice in the leapfrogtime-stepping. We put the discrete fields into the variational formulation, which nowtranslates to non-linear equation systems. The matrix entries are computed by numericalquadrature rules. Additional storage may be required for the matrices and solution ofthe equation systems, which depends on the solver.Finally, some numerical experiments inspired by the Apples-with-Apples test suite(Alcubierre et al. 2004, Babiuc et al. 2008) are used to compare both schemes with amore traditional FD scheme in section 4. The Galerkin schemes with piecewise linearfunctions on equidistant, cartesian grids show comparable CFL conditions, comparablesecond order accuracy, similar (sometimes opposite sign) dispersion second order in
EM, DG, and FD Evolution Schemes in Spacetime
2. Einstein’s Vacuum Equation
We start with the standard derivation of Einstein’s equation via the Einstein-Hilbertaction defined by S := Z M R √− gd x in the case of vacuum, in the notation of (Straumann 2004). We consider it as a functionof the metric tensor g αβ and its derivatives. The Ricci tensor R αβ and the Ricci scalar R = g αβ R αβ contain up to second order partial derivatives of g αβ . We are looking foran extremum of S . The variation of S is δS = Z M ( R µν − g µν R )( δg µν ) √− gd x , (2)as long as the variation δg µν vanishes at the boundary of the domain M . Otherwise weobtain an additional boundary term32 Z ∂ M g αβ g µν ( ∂ ν δg µβ − ∂ β δg µν ) n ( g ) α d x (3)which can be used later for boundary conditions using derivatives of g µν . We renamethe variation v µν := δg µν . The variational formulation reads as: We seek a solution g αβ ∈ V a such that δS = 0 forall v αβ ∈ V t with appropriate ansatz and trial spaces. Dirichlet boundary conditions on(parts of) ∂ M can be built into V a and V t : The solution takes the Dirichlet values andthe trial functions vanishes there. Boundary conditions involving derivatives require anadditional boundary term like (3). The variational formulation translates to the strongformulation as R µν − g µν R = 0 or in vacuum R µν = 0 EM, DG, and FD Evolution Schemes in Spacetime α := g ρσ Γ αρσ = 0 , (4)which is a condition on first order derivatives of g αβ . This way, we can modify Einstein’sequation as R ( h ) µν := R µν − g αν ∂ µ Γ α − g αµ ∂ ν Γ α = 0 , (5)with principal part R ( h ) ppµν := − g αβ ∂ α ∂ β g µν . (6)Now, we have a quasi-linear, second order, symmetric hyperbolic differential equation,which we will later discretize by finite differences. Note that this remains true if weswitch to a generalized harmonic gauge. Equation (4) changes to Γ α = H α ( x, g ) with agauge driver H . This driver may depend on coordinates and the metric, but must beindependent of derivatives of g in order to preserve the principal part R ( h ) pp . Galerkin discretizations are based on a variational formulation. We start with thestandard variational formulation (2). By Stokes’ theorem, we can remove the secondorder derivatives. With harmonic gauge (4) we arrive at a variational version of (6) a ( g, v ) := 12 Z M g αβ √− g ( ∂ α g µν )( ∂ β v µν ) d x , (7)which is symmetric in the first order derivatives of g µν and v µν in the special case of afixed background g µν . Again there is an additional boundary term, if the variation v does not vanish on the boundary ∂ M− Z ∂ M g αβ √− g ( ∂ β g µν ) v µν n ( g ) α d x . (8)The remaining terms can be assembled in q ( g, v ) := R M g αβ g ρσ √− g (cid:16) ( ∂ α g ρµ )( ∂ β g σν ) − ( ∂ α g ρµ )( ∂ σ g βν )+ ( ∂ α g ρµ )( ∂ ν g βσ ) + ( ∂ µ g αρ )( ∂ β g σν ) − ( ∂ µ g αρ )( ∂ ν g βσ ) (cid:17) v µν d x , (9)which is quadratic and symmetric in the first order derivatives of g µν , compare also(Fock 1959, App. B). The variational formulation now reads asseek g ∈ V a such that a ( g, v ) + q ( g, v ) = 0 ∀ v ∈ V t and Γ α = 0 . (10) EM, DG, and FD Evolution Schemes in Spacetime g ∈ V a in (10) does not need to have well defined second derivativesas in (5) and may be chosen in an appropriate Sobolev space. In the case of a nonvanishing energy-momentum tensor T µν additional terms of type b ( g, v ) := Z M ( g αν g µβ − g µν g αβ ) √− g T αβ v µν d x appear on the right-hand side of (10).Different types of initial and boundary conditions can be imposed on ∂ M bystandard procedures to define a Cauchy problem: Homogeneous Dirichlet values aredirectly incorporated into all functions in V a and V t . Inhomogeneous Dirichlet conditionsare built into the solution g , either direct in the discrete numerical scheme, or viaan additive splitting into a homogeneous auxiliary solution and a non-homogeneousfunction for the boundary conditions. Neumann boundary conditions and otherconditions based on derivatives of the solution on parts of ∂ M lead to additional termsin a of type (8), where ∂ β g µν is replaced by the given derivatives. The functions in V a and V t do not vanish there. “Natural” boundary conditions can be defined as vanishingterm (8), that is g αβ √− g ( ∂ β g µν ) n ( g ) α = 0. The conditions can be translated back into astrong formulation via (3). In a weak field approximation of Einstein’s equation, we neglect the first order derivativesin (5) and arrive at R ( h ) ppµν = 0 for some background metric ˆ g µν . In the variational version(10), we can neglect q ( g, v ) and solve for a ( g, v ) = 0 instead, again for a fixed backgroundmetric ˆ g µν . a ( g, v ) := 12 Z M ˆ g αβ p − ˆ g ( ∂ α g µν )( ∂ β v µν ) d x (11)The linearized version of the harmonic gauge condition (4) readsˆ g αβ ˆ g µν ( ∂ µ g νβ − ∂ β g µν ) = 0 . (12)Now, we simplify the equations even further and consider a weak field in flat space.The linearization is taken around Minkowski metric ˆ g = η := diag( − , , ,
1) and weobtain the strong formulation − (cid:3) g µν = 0 (13)with ∂ α = η αβ ∂ β and (cid:3) = ∂ α ∂ α . This translates to the variational versionseek g ∈ V a such that a ( g, v ) := 12 Z M η αβ ( ∂ α g µν )( ∂ β v µν ) d x = 0 ∀ v ∈ V t . (14)The harmonic gauge condition (12) reduces to ∂ µ g µν − η µν η αβ ∂ µ g αβ = 0 , which can be further simplified by the substitution h µν := g µν − η αβ g αβ to ∂ µ h µ,ν = 0 . (15) EM, DG, and FD Evolution Schemes in Spacetime − (cid:3) h µν = 0, now with a divergence free h . Thegauge conditions are linear and can be incorporated into the spaces V a and V t .
3. Numerical Schemes
For illustration purposes, the first numerical spacetime scheme will be based on finitedifferences. We consider the discretization of a linear, scalar, second order waveequation − (cid:3) u = 0 with suitable initial and boundary conditions. On a one-dimensional,equidistant grid with grid spacing h , we choose the stencil ( u ( x − h ) − u ( x )+ u ( x + h )) /h ,also abbreviated as [1 − /h , to approximate the second derivative. It is secondorder accurate for u smooth enough. The d’Alembert operator can be obtained byan application of the stencil along each coordinate axis on a cartesian grid. The twodimensional stencil at a grid point ( i, j ) for example is u i − ,j − u i,j + u i +1 ,j h − u i,j − − u i,j + u i,j +1 h = 0 , which gives the explicit time stepping scheme u i +1 ,j = 2 u i,j − u i − ,j + (cid:18) h h (cid:19) ( u i,j − − u i,j + u i,j +1 )using values at time slices i − i to calculate the values at time slice i + 1. This isthe leapfrog scheme in time and can be written as u i +1 = 2 u i − u i − + ( h ) ∆ h u i (16)with a FD approximation of the spatial derivatives ∆. Note that a CFL condition h /h k < k > x = 0 and x = h , the boundaryvalues at x k = 0 and x k = 1. Modifications for other types of initial and boundaryconditions do exist. In order to generalize the FD stencils to mixed first and second order derivatives, weconsider an alternative construction. In the one dimensional case, first derivatives canbe approximated by central stencils u ′ ( x + h/ ≈ ( u ( x + h ) − u ( x )) /h at grid points x + h/
2. The second derivative can be calculated as a central stencil of first derivatives u ′′ ( x ) ≈ ( u ′ ( x + h/ − u ′ ( x − h/ /h which reduces to the one-dimensional FD stencil.However, in two (and more) dimensions the construction differs, if we consider cell-centered first derivatives: We differentiate in one directions and average in the otherdirection(s): ∂ u i +1 / ,j +1 / ≈ (cid:18) u i − ,j − u i,j h + u i − ,j +1 − u i,j +1 h (cid:19) . EM, DG, and FD Evolution Schemes in Spacetime ∂ ∂ ≈ h ) − − − − − − and ∂ ∂ ≈ h h − − . The discretization of the d’Alembert operator again gives a time-stepping scheme fortime slice i + 1. However, the scheme is no more explicit like (16). Let us write thedifference stencil [1 2 1] / M and the stencil [ − − / ( h ) as matrix A . We obtain the scheme M u i +1 = 2 M u i − M u i − − ( h ) Au i . (17)We can compute the values u i +1 at time slice i + 1 by the solution of a linear equationsystem with matrix M using the values u i and u i − at time slices i and i −
1. The matrixis positive definite, symmetric, and of bounded condition number. Hence, the systemis easy to solve numerically for large systems by standard iterative solvers. Again, theCFL condition limits the time step size h . We start with a variational version of the d’Alembert operator (1), a first steptowards (11): 12 Z M η αβ ( ∂ α u )( ∂ β v ) d x = 0 ∀ v (18)Following standard procedures in FEM, we choose a set of global, continuous, piecewisepolynomial ansatz and trial functions ˜ φ i ∈ V a and ˜ ψ j ∈ V t as a basis of finite dimensionalspaces V a and V t , and obtain a finite element method: Find the coefficients ˜ u i of thesolution u = P i ˜ u i φ i ∈ V a , such that (18) holds for all trial functions v ∈ V t . This canbe written in basis functions as12 X i ˜ u i Z M η αβ ( ∂ α ˜ φ i )( ∂ β ˜ ψ j ) d x = 0 ∀ j (19)and in matrix notation ˜ A ˜ u = 0 with solution vector ˜ u and matrix ˜ A = (˜ a ij )˜ a ij = 12 Z M η αβ ( ∂ α ˜ φ i )( ∂ β ˜ ψ j ) d x . This is a spacetime discretization. Introducing a global time step, we split functions˜ φ i ( x ) = φ i ( x ) φ si ( x , x , x ) and ˜ ψ j , and the domain M = T × Ω into time and space.Further, mixed space-time derivatives η a = η b = 0 do not occur with space index a , b .We obtain ˜ a ij = (cid:0)R T η ( ∂ φ i ) ( ∂ ψ j ) dt (cid:1) (cid:0)R Ω φ si ψ sj d x (cid:1) + (cid:0)R T φ i ψ j dt (cid:1) (cid:0)R Ω η ab ( ∂ a φ si ) ( ∂ b ψ sj ) d x (cid:1) . We introduce the mass matrix M and matrix A by m ij = R Ω φ si ψ sj d x and a ij = R Ω η ab ( ∂ a φ si ) ( ∂ b ψ sj ) d x . EM, DG, and FD Evolution Schemes in Spacetime i − i anduse the scheme to calculate the next time slice i + 1. We use piecewise linear functions φ i ( t ) = max(1 − | t − x i | /h ,
0) and ψ j in time for equidistant time-steps h and obtainthe system in time h [ − − M + h [1 4 1] A , which is of leapfrog type (cid:18) M − h A (cid:19) u i +1 = (cid:18) M + 2 h A (cid:19) u i − (cid:18) M − h A (cid:19) u i − . (20)In the 1 + 1 spacetime case, piecewise linear functions on an equidistant space grid, wefurther obtain M = h [1 4 1] and A = h [ − − V a and trial V t spaces: We use piecewise polynomial functions centered at a grid point i at time i and space location ( i , i , i ) for a cartesian grid. The functions are chosenpiecewise linear in time. Let the grid points be in the time domain (0 , k ) with initialconditions at i = 0 and i = 1. We compute the solution for all ansatz functions φ i located at times (2 , k ). However, we choose the trial functions ψ j located at times(1 , k − i + 1 and coincides with a leapfrogscheme for equidistant time steps.The solution of equation systems is the most expensive part of the time steppingprocedure. The advantage of FD schemes for leapfrog (16) is that mass matrix M is the identity and no equation systems have to be solved. However, there is acommon technique in FEM called “mass lumping” to obtain diagonal matrices M , too:Integration in terms of type R Ω φ si φ sj d x and R T φ i φ j dt are approximated by numericalquadrature rules on each finite element. For piecewise (multi-) linear functions,quadrature rules prove to be sufficient, which are based on the function values at theelement vertices only. This is the trapezoidal rule on an edge and its generalizations torectangles, cubes, triangles and tetrahedra. The ansatz functions fulfill φ si ( x j ) = δ ij andthe mixed products ( φ si φ sj )( x k ) = δ ik δ jk = 0 for i = j vanish on all element vertices x k .Hence, off-diagonal entries m ij vanish and M is in fact a diagonal matrix. We arrive atthe computational efficiency of an FD scheme (16), once the integration is done.Note that (19) defines spacetime FEM also for higher order methods in space bypiecewise polynomial functions φ si and ψ sj , by pseudo-spectral Galerkin schemes in space,on unstructured grids in spacetime, and for adaptive grid refinement in spacetime. Theapproach does not easily extend to higher order methods in time due to a lack of stabilityof the respective time-stepping schemes. Again, we start with the variational problem (18). However, we choose piecewisepolynomial ansatz and trial functions φ i and ψ j , which are no longer continuousover element boundaries. This leads to additional terms. Consider a common face EM, DG, and FD Evolution Schemes in Spacetime e ij := ∂E i ∩ ∂E j of two neighbor elements E i and E j and normal unit vector n ij orientedfrom E i to E j . We denote the average by { u } := (( u | E i ) + ( u | E j )) / u ] := ( u | E i ) − ( u | E j ) on the face e ij . Let the volume of the face be | e ij | . We split theintegration over M of (14) into the integration over elements E i and all faces e ij of thegrid. a ( u, v ) := P i R E i η αβ ( ∂ α u )( ∂ β v ) d x − P i 20 2 0 11 − − 40 1 0 2 + c p A , = − c p h + c p h EM, DG, and FD Evolution Schemes in Spacetime i + 1 and position jA , u i +1 ,j = A , − u i,j − + A , u i,j + A , u i,j +1 − A − , u i − ,j . (22)Note that for each element a linear equation system A , needs to be solved. It is ofthe size of number of ansatz functions, which is cheaper to solve than the single largeequation system for the FEM. However, the amount of work can be further reduced: Itis possible to choose the local ansatz functions orthogonal with respect to the bi-linearform such that A , is in fact diagonal or even the identity and no systems need to besolved any more. This way, we obtain an explicit time-stepping scheme like (16).For a second order differential equation in time, we need two initial conditions, like u (0 , x ) and ∂ u (0 , x ). This can be converted into data on two initial time slices i = 0and I 01. However, for the DG schemes, we need an initial spacetime approximation inelements at times slices 0 and 1. For a linear ansatz in time direction, initial data isneeded at least at the beginning and end of both time slices, namely three initial values.These can be computed with a start-up calculation. In order to solve linearized Einstein’s equation (13) resp. (14), we can generalize thescalar schemes for (cid:3) u , apply these to each component g µν , and set the background metricto Minkowski ˆ g = η . The linear gauge condition (15) needs to be fulfilled. Divergence-free initial data guarantees this for all times in the continuous case. However, numericalerrors will lead to a violation of the gauge condition. DG methods easily allow for locallydivergence-free ansatz functions on each element. In contrast, it is difficult to implementglobally divergence-free symmetric tensor fields in FEM analogous to divergence-freevector fields for Maxwell’s equation, see (Nedelec 1980, Nedelec 1986).In the case of a prescribed curved background metric, we have to solve the linear,variable coefficient problem R ( h ) ppµν = 0. The FD stencils are no longer applicable andwe switch to the compact FDM stencils. The FEM implementation is based on thevariational formulation12 X i u i Z M ˆ g αβ p − ˆ g ( ∂ α φ i )( ∂ β ψ j ) d x = 0 ∀ j . (23)The DG method now reads as a ( u, v ) := P i R E i ˆ g αβ √− ˆ g ( ∂ α u )( ∂ β v ) d x − P i We generalize the compact FDM stencils to Einstein’s vacuum equation (5): Thevariable metric g µν and its second order derivatives R ( h ) µν are chosen node centered (atgrid points), but the first order derivatives Γ αβγ are chosen cell centered. The inversemetric g µν is used to calculate Γ αβγ and Γ α and is also cell centered, defined as the inverseof the cell average of the metric g µν . The products of averaged Γ enter the Ricci tensor,as well as the node centered derivatives of Γ. This way, we can use the standard formulasΓ αβγ := g ασ ( ∂ β g γσ + ∂ γ g βσ − ∂ σ g βγ ), R µν := ∂ α Γ αµν − ∂ µ Γ ανα +Γ βαβ Γ αµν − Γ βµα Γ ανβ , (4), and (5)to set up non-linear, discrete Einstein’s equation and derive the time-stepping scheme.Note that no code generated by a symbolic algebra program is needed.The FEM and DG Galerkin schemes can also be generalized to Einstein’s equation.The form a (7) resp. (21) and the quadratic term q (9) define the variational problem(10a). The integration is done numerically. The integral R M is split into integrals overan element P i R E i (and a face e ij in (21)). The integrals over a single element E i andface e ij are approximated by a numerical quadrature rule. The integrands of a and q are evaluated at the quadrature points.The matrices M and A now depend on the current metric g αβ √− g and the equationsystems of type (20) and (22) are non-linear. The time-stepping schemes are implicitand require the solution of a non-linear equation system for each time-slice. The DGmethod leads to a set of easy to solve local equation systems for each element. TheFDM and the FEM have global coupling of the degrees of freedom of a time slice. Inboth cases standard non-linear solvers can be used. Note that the explicit FD methodboth gives an initial guess for a locally fixed background metric g µν and can be used asa preconditioner for the principle part in an iterative solver.The harmonic gauge condition (4) now is a non-linear condition and cannot beincorporated into a linear ansatz space V a . Note that a change of variables leads to aformulation of Einstein’s equation with a new metric g µν := √− gg µν and a linear gaugecondition ∂ µ g µν = 0, which could be built into V a .Note that Regge calculus also discretizes a variational principle in spacetime forsimplicial grids (Sorkin 1975). It can be considered a DG spacetime scheme with piece-wise constant metric tensor g µν . This way, (21) generalizes it to higher order andarbitrary element shapes. However, Regge calculus does not use coordinates and isbased on purely geometric entities like edge lengths and defect angles. Furthermore, thevariation is with respect to the degrees of freedom, which are the squared edge lengthsin Regge calculus and values of the metric in (21). EM, DG, and FD Evolution Schemes in Spacetime −5 −4 −3 −2 −1 10 linear wave in 1+1 dimensionstime m a x e rr o r FDFDMFEMSIPDGNIPDG 0 200 400 600 800 1000 −4 −3 −2 −1 10 linear wave in 1+1 dimensions, FEMtime L2 e rr o r Figure 1. h = 1 / 200 (left) and of the l -error of a FEM solution with n = 1 /h (right). 4. Applications For illustration purposes, we perform some numerical experiments with the schemes ofsection 3. The test cases are adapted from the Apples-with-Apples test suite (Alcubierreet al. 2004, Babiuc et al. 2008). We document and compare convergence and stabilityof the schemes in different settings.We start with a mono-chromatic traveling plane wave for linearized Einstein’sequation (13) and (14) with harmonic gauge (15). We use periodic boundary conditionsand a Courant factor 1 / 2. The one-dimensional (1+1) test case is defined on thespatial unit interval (0 , g = g = − g =sin2 π ( x − x ). We use an equidistant grid and run all schemes of sections 3.1 to 3.4. Notethat the original Apples-with-Apples tests were constructed for non-linear numericalcodes, such that very small wave amplitudes effectively ran the codes in the regime ofthe linearized equations. Standard non-linear solvers like Newton’s method in this casereduce the problem to a linear one. Hence, we directly ran a linear code for the linearproblem. This is why we can use arbitrary amplitudes of the solution rather than verysmall ones (Alcubierre et al. 2004).In figure 1 the time evolution of the spatial maximum error at the grid points fora resolution of h = 1 / 200 is depicted for the FD, FDM, FEM, SIPDG and NIPDG.We use penalty parameters c p = 1 and c p = 2 for SIPDG. Note that continuous errornorms like L (0 , 1) more natural for FEM show a similar behavior with exception ofthe very first time steps, where an additional interpolation error is added to the globalerror. The point-wise divergence is bounded, although we do not take any measures tocontrol it. This does not seem to be necessary. The solution in figure 2 (left) shows thespatial errors at the final time x = 1000. We see mainly dispersion and the phase errorof the different schemes, no errors in the amplitude. This is why the error in fact even EM, DG, and FD Evolution Schemes in Spacetime u12 FDFDMFEMSIPDGNIPDGexact 0 200 400 600 800 1000 −4 −3 −2 −1 10 linear wave in 2+1 dimensionstime m a x e rr o r FDFDMFEMSIPDG Figure 2. Solution of a linear plane wave in 1+1 for h = 1 / 200 at final time x = 1000(left). Evolution of the maximum error of a wave in 2+1 on a cartesian grid with h = 1 / 100 (right). −4 −3 −2 −1 10 linear wave in 2+1 dimensions, FEMtime m a x e rr o r quadtri1tri1*tri2tri2*tri3tri3* quad tri1*tri2* tri3* Figure 3. Linear plane wave in 2+1 on some triangulations with h = 1 / decreases after some time, see figure 1 (right) n = 25 and n = 50. We observe a secondorder convergence of the error, the phase error and the divergence in h for all schemes.The two-dimensional (2+1) test case is defined on the spatial unit square (0 , with periodic boundary conditions. The exact solution and initial data is g = g = g = sin2 π ( x + x − √ x ), g = g = ( √ − g , and g = √ g . We runall schemes on cartesian equidistant grids, see figure 2 (right), except for the NIPDGscheme for a lack of stability. The SIPDG penalty term is chosen as c e = 1 / 2, moreprecisely c p | e ij | ce = 1 /h . The second order convergence is comparable to the 1 + 1 case.In order to test the dependence on the spatial grid, we run the FEM also ona number of triangular grids, both uniform (tri) and randomly distorted (tri*), seefigure 3. Now we obtain a strong dependence of the error on the orientation of theelements. The longest element edges tangential to the direction of the wave leads to alarger approximation error than in normal direction or for quadratic elements. EM, DG, and FD Evolution Schemes in Spacetime −11 −10 −9 −8 −7 −6 10 robust stability in 1+1 dimensionstime m a x e rr o r FDFDMFEMSIPDGNIPDG 0 0.2 0.4 0.6 0.8 1 1.2 −6 −5 −4 −3 −2 −1 10 expanding Gowdy universetime m a x e rr o r FDMFEMSIPDGNIPDG Figure 4. Evolution of the maximum error of the linear robust stability test in 1+1for h = 1 / 200 (left) and of a non-linear Gowdy wave in 1+1 for h = 1 / 100 and h = 1 / 200 (right). Now we consider a stability test for the linear wave equation. The starting point isa random perturbation of the zero solution. We use periodic boundary conditions on(0 , − ǫ, ǫ ] random values for all initial data, with ǫ := 2 . · − ( h ) according to (Alcubierre et al. 2004). In figure 4 (left) we observe stability of all schemeswith oscillatory solutions for NIPDG and compact stencil FDM. The polarized Gowdy spacetime on the Torus T is a model for a gravitational wavein an expanding universe (Gowdy 1971, New et al. 1998). We use periodic boundaryconditions on the spatial unit interval (0 , 1( in x direction. The solution is constantalong x and x direction. Since we use harmonic gauge, time axis x differs from(Alcubierre et al. 2004). We use a Courant factor 1 / 4. The solution g µν is given by g = diag( − e ( λ +3 x ) / , e x + p , e x − p , e ( λ − x ) / ) with p := J (2 πe x )cos(2 πx ) and λ := − πe x J (2 πe x )J (2 πe x )cos (2 πx ) − π J (2 π )J (2 π )+2( πe x ) (J (2 πe x ) + J (2 πe x )) − (2 π ) (J (2 π ) + J (2 π ))We run schemes of section 3.6 with 3rd order Gauss quadrature (two points in eachcoordinate direction) on an element. The SIPDG penalty terms are chosen as c p = . c p = 2. In figure 4 (right) we see the error for spatial resolutions h = 1 / 100 and h = 1 / t due to the exponential growth of some of the solution components(Alcubierre et al. 2004). EM, DG, and FD Evolution Schemes in Spacetime Conclusion We have developed new spacetime Finite Element (FEM) and Interior PenaltyDiscontinuous Galerkin (SIPDG and NIPDG) schemes for second order symmetrichyperbolic wave equations. The Discontinuous Galerkin schemes are computationallymore efficient, but require more memory than FEM and Finite Differences methods.A variational formulation of Einstein’s equation in harmonic gauge was derived, basedon up to first derivatives of solution and trial functions. This led to new Galerkinschemes for numerical relativity. The schemes were presented and tested for secondorder accurate Galerkin schemes with multi-linear functions and global time steps. TheGowdy wave test demonstrated the need for additional numerical stabilization. Thismight be obtained by spatial filtering, artificial viscosity, or streamline diffusion.Extensions to arbitrary spacetime grids or (adaptive) local grid refinement inspacetime are straightforward, but may lead to larger and more expensive to solveequation systems. Higher order polynomials or other more accurate (spectral) functionspaces improve the spatial accuracy of the schemes. However higher order in timeschemes are more difficult to construct. Acknowledgments The author wants to thank G. Sch¨afer for several hints to the literature. Furthermore,helpful comments by S. Husa and the anonymous referees are acknowledged. This workwas partially supported by DFG grant SFB/TR7 “gravitational wave astronomy”. References Ainsworth M, Monk P & Muniz W 2006 J. Scient. Comp. , 5–40.Aksoylu B, Bernstein D, Bond S & Holst M 2008. arXiv:0801.3142v3.Alcubierre M, Allen G, Bona C, Fiske D, Goodale T, Guzman F S, Hawke I, Hawley S H, Husa S,Koppitz M, Lechner C, Pollney D, Rideout D, Salgado M, Schnetter E, Seidel E, Shinkai H,Szil´agyi B, Shoemaker D, Takahashi R & Winicour J 2004 Class. Quant. Grav. , 589.Anderson M & Kimn J H 2007 J. Comp. Phys. , 466–476.Arnold D N, Mukherjee A & Pouly L 1998 in D. F Griffiths, D. J Higham & G. A Watson, eds,‘Numerical Analysis 1997’ Addison Wesley Longman pp. 1–15.Arnowitt R, Deser S & Misner C W 1962 in L Witten, ed., ‘Gravitation: An Introduction to CurrentResearch’ Wiley chapter 7, pp. 227–265.Babiuc M C, Husa S, Alic D, Hinder I, Lechner C, Schnetter E, Szil´agyi B, Zlochower Y, Dorband N,Pollney D & Winicour J 2008 Class. Quant. Grav. , 125012.Baker G A & Bramble J H 1979 RAIRO Anal. Numer. , 75–100.Baker J G, Centrella J, Choi D I, Koppitz M & van Meter J 2006 Phys. Rev. Lett. , 111102.Baumgarte T W & Shapiro S L 1999 Phys. Rev. D , 024007.Bonazzola S, Gourgoulhon E, Grandclement P & Novak J 2004 Phys. Rev. D , 104007.Boyle M, Brown D A, Kidder L E, Mroue A H, Pfeiffer H P, Scheel M A, Cook G B & Teukolsky S A2007 Phys. Rev. D , 124038.Brandt S & Bruegmann B 1997 Phys. Rev. Lett. , 3606–3609.Bruhat Y 1962 in L Witten, ed., ‘Gravitation: An Introduction to Current Research’ Wiley pp. 130–168. EM, DG, and FD Evolution Schemes in Spacetime Campanelli M, Lousto C O, Marronetti P & Zlochower Y 2006 Phys. Rev. Lett. , 111101.Cohen G C 2002 Higher-Order Numerical Methods for Transient Wave Equations. Springer.Dupont T 1973 SIAM J. Numer. Anal. (5), 880–889.Eriksson K, Johnson C & Thom´ee V 1985 RAIRO M.M.A.N. , 611–643.Field S E, Hesthaven J S & Lau S R 2009. arXiv:0902.1287.Fock V 1959 The Theory of Time Space and Gravitation. Pergamon Press.French D A & Peterson T E 1996 Math. Comp. , 491–506.Friedrich H & Rendall A D 2000 Lect. Notes Phys. , 127–224.Gowdy R H 1971 Phys. Rev. Lett. , 826–829.Grote M, Schneebeli A & Sch¨otzau D 2006 SIAM J. Numer. Anal. , 2408–2431.Hulbert G M & Hughes T J R 1990 Comput. Meth. Appl. Mech. Engin. , 327–348.Jamet P 1978 SIAM J. Numer. Anal. , 912–928.Metzger J 2004 Class. Quant. Grav. , 4625–4646.Monk P & Richter G R 2005 J. Scient. Comp. , 443–477.Nedelec J C 1980 Numer. Math. , 315–341.Nedelec J C 1986 Numer. Math. , 57–81.New K C B, Watt K, Misner C W & Centrella J M 1998 Phys. Rev. D , 064022.Pretorius F 2005 Class. Quant. Grav. , 425–452.Regge T 1961 Nuovo Cimento A , 558–571.Reula O A 1998 Living Rev. Relativity (3), 1–40.Rivi`ere B 2008 Discontinuous Galerkin Methods for Solving Elliptic and Parabolic Equations. SIAM.Shibata M & Nakamura T 1995 Phys. Rev. D , 5428–5444.Sopuerta C F, Sun P, Laguna P & Xu J 2006 Class. Quant. Grav. , 251–285.Sorkin R 1975 Phys. Rev. D (2), 385–396.Straumann N 2004 General Relativity. Springer.York, Jr. J W 1979 inin