Weight-adjusted discontinuous Galerkin methods: matrix-valued weights and elastic wave propagation in heterogeneous media

1.7 MB·

by  Jesse Chan

WEIGHT-ADJUSTED DISCONTINUOUS GALERKIN METHODS: MATRIX-VALUED
WEIGHTS AND ELASTIC WAVE PROPAGATION IN HETEROGENEOUS MEDIA
JESSE CHAN
Abstract. Weight-adjustedinnerproducts[1,2]areeasilyinvertibleapproximationstoweightedL2innerproducts. These
approximations can be paired with a discontinuous Galerkin (DG) discretization to produce a time-domain method for wave
propagation which is low storage, energy stable, and high order accurate for arbitrary heterogeneous media and curvilinear
meshes. In this work, we extend weight-adjusted DG (WADG) methods to the case of matrix-valued weights, with the linear
elasticwaveequationasanapplication. WepresentaDGformulationofthesymmetricformofthelinearelasticwaveequation,
with upwind-like dissipation incorporated through simple penalty fluxes. A semi-discrete convergence analysis is given, and
numericalresultsconfirmthestabilityandhighorderaccuracyofWADGforseveralproblemsinelasticwavepropagation.
7
1
1. Introduction. Efficient and accurate methods for elastic wave propagation form a foundation for
0
a broad range of applications, from seismic and medical imaging to rupture and earthquake simulation.
2
Finitedifferencesarethemostcommonchoiceofmethod[3]; however, finiteelementmethodshavegarnered
n
interestduetotheirlownumericaldispersionandabilitytoaccommodategeometricallyflexibleunstructured
a
J meshes.
1 Typical methods for time-domain wave propagation utilize explicit time stepping, since the hyperbolic
partial differential equations (PDEs) which govern wave propagation admit a reasonable stable time-step

] restriction. However, finite element methods require the inversion of a global mass matrix when paired with
A
explicit time integrators. Spectral element methods (SEM) sidestep this issue on hexahedral meshes by
N
utilizing a mass lumped diagonal mass matrix [4]. The inversion of a globally coupled matrix can also be
. avoided through the use of discontinuous Galerkin (DG) methods, which result in a locally invertible block
h
t diagonal mass matrices. Due to difficulties in extending mass-lumping techniques from hexahedra to tetra-
a
hedra, highorderDGmethodsareoftenemployedforseismicsimulationswhichrequiretheuseofsimplicial
m
meshes [5, 6, 7, 8, 9, 10]. High order DG methods also lend themselves well to efficient implementations
[
using Graphics Processing Units (GPUs) [11, 12, 13, 14].

1 Most high order DG methods on simplicial meshes assume that material coefficients are constant over
v each element, which can result in spurious reflections and the loss of high order accuracy [15, 16, 17].
5
This limitation can be overcome by incorporating sub-element heterogeneities into weighted mass matrices,
1
resultinginaDGmethodwhichisbothhighorderaccurateandenergystable[16,17]. Ontetrahedralmeshes,
2
0 this approach necessitates the pre-computation and storage of factorizations or inverses for each local mass
0 matrix, which greatly increases both storage costs and data transferred at high orders of approximation.
. These costs are especially problematic for accelerator architectures such as GPUs, which possess limited
1
0 memory.
7 Storage costs associated with weighted mass matrices can be avoided by approximating weighted L2
1 inner products using weight-adjusted inner products, which result in easily invertible approximations to
:
v weighted mass matrices [1, 2]. For sufficiently regular weights, high order accuracy is also retained. When
i paired with an energy stable DG formulation, these approximations result in weight-adjusted DG methods
X
(WADG), which preserve energy stability and high order accuracy while retaining a low asymptotic storage
r
a cost. Additionally, unlike mass-lumping techniques, WADG methods do not rely on the use of inexact
quadrature rules, and reduce to the exact inversion of mass matrices for constant weights.
Weight-adjusted DG methods have been applied to acoustic wave propagation in heterogeneous media
andoncurvilinearmeshes[1,2]. Bothofthesepreviousapplicationshaveinvolvedscalarweightingfunctions.
In this work, we extend weight-adjusted inner products to matrix-valued weights, and demonstrate the
usefulness of this generalization to elastic wave propagation in arbitrary media. A new energy stable DG
formulation is introduced based on the symmetric form of the elastic wave equations, with upwind-like
numerical dissipation introduced through simple penalty fluxes [18]. In contrast to the fluxes proposed in
[10], the penalty fluxes used here are independent of material coefficients.
This work proceeds as follows: Sections 2 and 3 present an energy stable DG formulation with simple
penalty fluxes for the symmetric hyperbolic form of the elastic wave equation, and discuss issues related to
storage and inversion of local mass matrices for material coefficients with sub-element variations. Section 4
extends weight-adjusted approximations to weighted L2 inner products and mass matrices to the case of
1
matrix-valued weights, and provides interpolation estimates which account for the regularity of the matrix
weight. Theseresultsareincorporatedintoaweight-adjustedDGmethodforthelinearelasticwaveequations
in Section 5. Finally, numerical results in Section 6 demonstrate the accuracy of this method for several
problems in linear elasticity.
2. Symmetric form of the elastic wave equation. We begin with the linear elastic wave equation
in a domain Ω ∈ Rd. These equations can be written a first order velocity-stress system for velocity v and
symmetric stress tensor S
∂v
ρ =∇·S+f
∂t
∂S = 1C(cid:0)∇v+∇vT(cid:1),
∂t 2
where f is the body force per unit volume, ρ is density, and C is the symmetric constitutive stiffness tensor
relating stress and strain. We rewrite these equations as a symmetric hyperbolic system of PDEs [19] using
Voigt notation
d
∂v (cid:88) ∂σ
ρ = AT +f
∂t i ∂x
i
i=1
d
∂σ (cid:88) ∂v
(1) C−1 = A .
∂t i∂x
i
i=1
where C is the symmetric matrix form of the constitutive tensor C and σ is a vector of length N = d(d+1),
d 2
the number of unique entries of the stress tensor S in d dimensions. We note that the matrices A are
i
spatially constant, while ρ,C, and C−1 can vary spatially. Furthermore, we will assume that ρ and C are
positive-definite and bounded pointwise such that
0<ρ ≤ρ(x)≤ρ <∞,
min max
0<c ≤uTC(x)u≤c <∞
min max
0<c˜ ≤uTC−1(x)u≤c˜ <∞
min max
for all x∈Rd and all u∈RNd.
In two dimensions, v =(v ,v )T and σ =(σ ,σ ,σ )T
1 2 xx yy xy
(cid:18) (cid:19)
σ σ
S = xx xy ,
σ σ
xy yy
while the matrices A are
i
   
1 0 0 0
A1 = 0 0 , A2 = 0 1 .
0 1 1 0
In three dimensions, the velocity is v = (v ,v ,v )T, while σ = (σ ,σ ,σ ,σ ,σ ,σ )T denotes the
1 2 3 xx yy zz yz xz xy
unique entries of the stress tensor S
 
σ σ σ
xx xy xz
S = σxy σyy σyz .
σ σ σ
xz yz zz
The matrices A are then
i
     
1 0 0 0 0 0 0 0 0
 0 0 0   0 1 0   0 0 0 
     
 0 0 0   0 0 0   0 0 1 
A1 = 0 0 0 , A2 = 0 0 1 , A3 = 0 1 0 
     
 0 0 1   0 0 0   1 0 0 
0 1 0 1 0 0 0 0 0
2
Ingeneralanisotropicmedia,C issymmetricandpositive-definite. Fortwo-dimensionalisotropicmedia,
C is given as
 
2µ+λ λ 0
C = λ 2µ+λ 0 ,
0 0 µ
where λ,µ are Lame parameters. For three-dimensional isotropic media, C is given instead by
 
2µ+λ λ λ
 λ 2µ+λ λ 
 .
 λ λ 2µ+λ 
µI3×3
We will consider both spatially varying isotropic and anisotropic media in this work.
3. An energy stable discontinuous Galerkin formulation for elastic wave propagation. The
advantage of the symmetric first order formulation of the elastic wave equations (1) is that it is straight-
forward to derive an energy stable discontinuous Galerkin formulation. We assume that the domain Ω is
Lipschitz and exactly triangulated by a mesh Ω , which consists of elements Dk. We further assume that
h
each element Dk is the image of a reference element D(cid:98) under the local elemental mapping
xk =Φkx,
(cid:98)
where xk = (cid:8)xk,yk,zk(cid:9) denote physical coordinates on Dk and x = {x,y,z} denote coordinates on the
(cid:98) (cid:98) (cid:98) (cid:98)
reference element. We denote the determinant of the Jacobian of Φk as J, and refer to it as the Jacobian
for the remainder of this work.
We will approximate solution components over each element Dk from an approximation space V (cid:0)Dk(cid:1),
h
(cid:16) (cid:17)
which we define as the composition of the mapping Φk and a reference approximation space Vh D(cid:98)
Vh(cid:0)Dk(cid:1)=Φk◦Vh(cid:16)D(cid:98)(cid:17).
TheglobalapproximationspaceV (Ω )isthendefinedasthedirectsumofelementalapproximationspaces
h h
V (Ω )=(cid:77)V (cid:0)Dk(cid:1).
h h h
Dk
(cid:16) (cid:17) (cid:16) (cid:17) (cid:16) (cid:17)
For the remainder of this work, we will take Vh D(cid:98) =PN D(cid:98) , where PN D(cid:98) is the polynomial space of
total degree N on the reference simplex. In two dimensions, PN on a triangle is
PN(cid:16)D(cid:98)(cid:17)=(cid:8)x(cid:98)iy(cid:98)j, 0≤i+j ≤N(cid:9),
and in three dimensions, PN on a tetrahedron is
PN(cid:16)D(cid:98)(cid:17)=(cid:8)x(cid:98)iy(cid:98)jz(cid:98)k, 0≤i+j+k ≤N(cid:9).
We denote the L2 inner product and norm over Dk by (·,·) , such that
L2(Dk)
(cid:90) (cid:90)
(f,g) = f ·gdx= f ·gJ, (cid:107)f(cid:107)2 =(f,f) ,
L2(Dk) L2(Dk) L2(Dk)
Dk D(cid:98)
wheref,garevector-valuedfunctions. GlobalL2innerproductsandnormsareusinglocalL2innerproducts
and norms
(cid:90)
(f,g) = (cid:88) (f,g) , (cid:107)f(cid:107)2 = (cid:88) (cid:107)f(cid:107)2 .
L2(Ω) L2(Dk) L2(Ω) L2(Dk)
Dk
Dk∈Ωh Dk∈Ωh
3
We define also the L2 inner product and norm over the boundary ∂Dk of an element
(cid:90) (cid:90)
(cid:104)u,v(cid:105) = uvdx= (cid:88) uvJfdx, (cid:107)u(cid:107)2 =(cid:104)u,u(cid:105) ,
L2(∂Dk) (cid:98) L2(∂Dk) L2(∂Dk)
∂Dk f∈∂Dk f(cid:98)
where Jf is the Jacobian of the mapping from a reference face f(cid:98)to a physical face f of an element.
Letf beafaceofanelementDk withneighboringelementDk,+ andunitoutwardnormaln. Wedefine
u− and u+ denote interior and exterior values of a discontinuous function u such that
u− = u| , u+ = u| .
f∩∂Dk f∩∂Dk,+
The jump and average of a scalar function u∈V (Ω ) over f are then defined as
h h
u++u−
u =u+−u−, {{u}}= .
(cid:74) (cid:75) 2
Jumps and averages of a vector-valued functions u∈Rm and S ∈Rm×n are then defined component-wise
( u ) = u , 1≤i≤m, ( S ) = S , 1≤i≤m, 1≤j ≤n.
i i ij i
(cid:74) (cid:75) (cid:74) (cid:75) (cid:74) (cid:75) (cid:74) (cid:75)
We can now specify a DG formulation for the linear elastic wave equation (1). Symmetric hyperbolic
systems readily admit a DG formulation based on penalty fluxes [20]. For the linear elastic wave equation
in symmetric first order form, this formulation is given as
(cid:18) (cid:19) (cid:32) d (cid:33) (cid:28) (cid:29) 
D(cid:88)k∈Ωh ρ∂∂vt,w L2(Dk) =D(cid:88)k∈Ωh (cid:88)i=1ATi ∂∂xσi +f,w L2(Dk)+ 12ATn(cid:74)σ(cid:75)+ τ2vATnAn(cid:74)v(cid:75),w L2(∂Dk)
(2)
(cid:18) (cid:19) (cid:32) d (cid:33) (cid:28) (cid:29) 
D(cid:88)k∈Ωh C−1∂∂σt,q L2(Dk) =D(cid:88)k∈Ωh (cid:88)i=1Ai∂∂xvi,q L2(Dk)+ 12An(cid:74)v(cid:75)+ τ2σAnATn(cid:74)σ(cid:75),q L2(∂Dk),
for all w,q ∈ V (Ω ). Here, τ ,τ ≥ 0 are penalty parameters and A is the normal matrix defined on a
h h v σ n
face f as A =(cid:80)d n A . In two dimensions, A is
n i=1 i i n
 
n 0
x
An = 0 ny .
n n
y x
while in three dimensions, A is
n
 
n 0 0
x
 0 ny 0 
 
An = 00 n0z nnyz .
 
 nz 0 nx 
n n 0
y x
We note that, unlike the penalty DG formulation given in [10], material coefficients ρ,C appear only on
the left hand side of (2). Thus, efficient techniques for constant coefficient formulations [21] can be used to
evaluate the right hand side of the formulation, even in the presence of sub-element variations in ρ,C.
3.1. Boundary conditions. In this work, we assume boundary conditions on velocity and traction of
the form
v =v , Sn=t
bc bc
where v and t are given values. Traction boundary conditions where t = 0 are referred to as free-
bc bc bc
surface boundary conditions. We follow [22, 23] and impose boundary conditions on the DG formulation
4
through exterior values and jumps of the solution. Boundary conditions on the normal component of the
stress can be imposed by noting that the numerical flux contains the term ATσ = Sn .
n
For a face which lies on a boundary, velocity boundary conditions are i(cid:74)mpose(cid:75)d by(cid:74) set(cid:75)ting
v =2(cid:0)v −v−(cid:1), ATσ = Sn =0,
bc n
(cid:74) (cid:75) (cid:74) (cid:75) (cid:74) (cid:75)
while traction boundary conditions are enforced through
ATσ = Sn =2(cid:0)t −S−n(cid:1)=2(cid:0)t −ATσ−(cid:1), v =0.
n bc bc n
(cid:74) (cid:75) (cid:74) (cid:75) (cid:74) (cid:75)
For problems which involve the truncation of infinite or large domains, absorbing boundary conditions are
required. For such cases, we impose simple non-reflective boundary conditions [22] through jumps
ATσ = Sn =−S−n=−ATσ−, v =−v−.
n n
(cid:74) (cid:75) (cid:74) (cid:75) (cid:74) (cid:75)
We note that more accurate absorbing conditions can be imposed using, for example, perfectly matched
layers [24] or high order absorbing boundary conditions [25, 26].
In all cases, boundary conditions are imposed by computing numerical fluxes using these modified
jumps. Thisimpositionguaranteesenergystabilityforfreesurface,non-reflective,andhomogeneousvelocity
boundary conditions.
3.2. Energy stability. Integrating the velocity equations of (2) by parts gives
(cid:18)ρ∂v,w(cid:19) =−(cid:32)(cid:88)d σ,A ∂w(cid:33) +(cid:68)AT {{σ}}+ τvATA v ,w(cid:69)
∂t L2(Dk) i=1 i∂xi L2(Dk) n 2 n n(cid:74) (cid:75) L2(∂Dk)
(cid:18) (cid:19) (cid:32) d (cid:33) (cid:28) (cid:29)
C−1∂σ,q = (cid:88)A ∂v ,q + 1A v + τσA AT σ ,q
∂t L2(Dk) i=1 i∂xi L2(Dk) 2 n(cid:74) (cid:75) 2 n n(cid:74) (cid:75) L2(∂Dk)
Taking (w,q)=(v,σ) and adding both equations together yields
(cid:88) 1 ∂ (cid:16)(ρv,v) +(cid:0)C−1σ,σ(cid:1) (cid:17)
2∂t L2(Dk) L2(Dk)
Dk∈Ωh
= (cid:88) (cid:68)AT {{σ}}+ τvATA v ,v(cid:69) +(cid:28)1A v + τσA AT σ ,σ(cid:29)
Dk∈Ωh n 2 n n(cid:74) (cid:75) ∂Dk 2 n(cid:74) (cid:75) 2 n n(cid:74) (cid:75) ∂Dk
= (cid:88) (cid:88) (cid:90) AT {{σ}}v+ τvATA v v+ 1A v σ+ τσA AT σ σ.
f n 2 n n(cid:74) (cid:75) 2 n(cid:74) (cid:75) 2 n n(cid:74) (cid:75)
Dk∈Ωhf∈∂Dk
Let Γ denote the set of unique faces in Ω , and let Γ ,Γ ,Γ denote the parts of the boundary where
h h v σ abc
velocity, traction, and non-reflective boundary conditions are imposed, respectively. We separate surface
terms into contributions from interior shared faces and from boundary faces. On an interior shared face, we
sum contributions from the two adjacent elements to yield
(cid:88) (cid:90) AT {{σ}}v+ τvATA v v+ 1A v σ+ τσA AT σ σ
f n 2 n n(cid:74) (cid:75) 2 n(cid:74) (cid:75) 2 n n(cid:74) (cid:75)
f∈Γh\∂Ω
=− (cid:88) (cid:90)f τ2v |An(cid:74)v(cid:75)|2+ τ2σ (cid:12)(cid:12)ATn(cid:74)σ(cid:75)(cid:12)(cid:12)2.
f∈Γh\∂Ω
For faces which lie on the boundary Γ where velocity boundary conditions are imposed, v = −2v−,
v
ATσ =0, and AT {{σ}}=Sn−, implying that (cid:74) (cid:75)
n n
(cid:74) (cid:75)
(cid:88) (cid:90) AT {{σ}}v+ τvATA v v+ 1A v σ+ τσA AT σ σ
f∈Γv f n 2 n n(cid:74) (cid:75) 2 n(cid:74) (cid:75) 2 n n(cid:74) (cid:75)
(cid:90) (cid:90)
= (cid:88) ATnσ−v−−Anv−σ−−τvATnAn(cid:12)(cid:12)v−(cid:12)(cid:12)2 =− (cid:88) τvATnAn(cid:12)(cid:12)v−(cid:12)(cid:12)2.
f∈Γv f f∈Γv f
5
For faces which lie on Γ , AT σ =−2ATσ−, AT {{σ}}=0, and v =0, yielding a similar contribution
σ n n n
(cid:74) (cid:75) (cid:74) (cid:75)
f(cid:88)∈Γσ(cid:90)fATn {{σ}}v+ τ2vATnAn(cid:74)v(cid:75)v+ 12An(cid:74)v(cid:75)σ+ τ2σAnATn(cid:74)σ(cid:75)σ =−f(cid:88)∈Γσ(cid:90)fτσAnATn (cid:12)(cid:12)σ−(cid:12)(cid:12)2.
Finally, for faces in Γ we have AT {{σ}}= 1ATσ−, AT σ =−ATσ−, and v =−v−, yielding
abc n 2 n n n
(cid:74) (cid:75) (cid:74) (cid:75)
(cid:88) (cid:90) AT {{σ}}v+ τvATA v v+ 1A v σ+ τσA AT σ σ =
f∈Γabc f n 2 n n(cid:74) (cid:75) 2 n(cid:74) (cid:75) 2 n n(cid:74) (cid:75)
− (cid:88) (cid:90) τ2vATnAn(cid:12)(cid:12)v−(cid:12)(cid:12)2+ τ2σAnATn (cid:12)(cid:12)σ−(cid:12)(cid:12)2.
f∈Γabc f
Combining all face contributions together gives the following result:
Theorem 1. The DG formulation (2) is energy stable for τ ,τ ≥0, in the sense that
v σ
(cid:88) 21∂∂t(cid:16)(ρv,v)L2(Dk)+(cid:0)C−1σ,σ(cid:1)L2(Dk)(cid:17)=− (cid:88) (cid:90)f τ2v |An(cid:74)v(cid:75)|2+ τ2σ (cid:12)(cid:12)ATn(cid:74)σ(cid:75)(cid:12)(cid:12)2
Dk∈Ωh f∈Γh\∂Ω
(3) − (cid:88) (cid:90) τv(cid:12)(cid:12)Anv−(cid:12)(cid:12)2− (cid:88) (cid:90) τσ(cid:12)(cid:12)ATnσ−(cid:12)(cid:12)2− (cid:88) (cid:90) τ2v (cid:12)(cid:12)Anv−(cid:12)(cid:12)2+ τ2σ (cid:12)(cid:12)ATnσ−(cid:12)(cid:12)2 ≤0.
f∈Γv f f∈Γσ f f∈Γabc f
Since ρ and C−1 are positive definite, the left hand side of (3) is an L2-equivalent norm on (v,σ), and
Theorem 1 implies that the magnitude of the DG solution (v,σ) is non-increasing in time. This also shows
that dissipation present for positive penalization constants τ ,τ acts on non-conforming components with
v σ
non-zero jumps AT v and A σ . In fact, it was shown in [20] that, in the limit as τ ,τ → ∞, the
n n v σ
eigenspaces of DG d(cid:74)isc(cid:75)retization(cid:74)s s(cid:75)plit into a conforming part consisting of u,σ which satisfy
A u =0, AT σ =0
n n
(cid:74) (cid:75) (cid:74) (cid:75)
and an non-conforming part (defined through the L2 orthogonal complement) corresponding to eigenvalues
containrealpartswhichapproach−∞. Forthelinearelasticwaveequations,theseconditionsareequivalent
to requirements of C0 continuity for u and normal continuity of the stress tensor S n=0.
(cid:74) (cid:75)
3.3. The semi-discrete matrix system for DG. The solution to (2) can be approximated by dis-
cretizinginspaceandusinganexplicittimeintegrator,whichrequiresonlyevaluationsoflocalcontributions
over Dk to the DG formulation
(cid:18) (cid:19) (cid:32) d (cid:33) (cid:28) (cid:29)
ρ∂v,w = (cid:88)AT ∂σ,w + 1AT σ + τvATA v ,w
∂t L2(Dk) i=1 i ∂xi L2(Dk) 2 n(cid:74) (cid:75) 2 n n(cid:74) (cid:75) L2(∂Dk)
(cid:18) (cid:19) (cid:32) d (cid:33) (cid:28) (cid:29)
(4) C−1∂σ,q = (cid:88)A ∂v ,q + 1A v + τσA AT σ ,q .
∂t L2(Dk) i=1 i∂xi L2(Dk) 2 n(cid:74) (cid:75) 2 n n(cid:74) (cid:75) L2(∂Dk)
(cid:16) (cid:17)
Let {φi}Ni=p1 be a basis for PN D(cid:98) . We define the reference mass matrix M(cid:99) and the physical mass
matrix M for an element Dk as
(cid:16) (cid:17) (cid:90) (cid:90) (cid:90)
M(cid:99) = φjφidx(cid:98), (M)ij = φjφidx= φjφiJdx(cid:98).
ij D(cid:98) Dk D(cid:98)
For affine mappings, J is constant and M =JM(cid:99). We also define weak differentiation matrices Sk and face
mass matrix M such that
f
(cid:90) ∂φ (cid:90) (cid:90)
(S ) = jφ dx, (M ) = φ φ dx= φ φ Jfdx,
k ij ∂x i f ij j i j i (cid:98)
Dk k f f(cid:98)
6
where Jf is the Jacobian of the mapping from a reference face f(cid:98)to f. For affinely mapped simplices, Jf is
also constant and Mf =JfM(cid:99)f, where the definition of the reference face mass matrix M(cid:99)f is analogous to
the definition of the reference mass matrix M(cid:99).
Finally, we define weighted mass matrices. Let w(x) ∈ R and W(x) ∈ Rm×n. Then, scalar and
matrix-weighted mass matrices M and M are defined through
w W
 
M ... M
(cid:90) W1,1 W1,n
(Mw)ij = w(x)φj(x)φi(x)dx, MW = ... ... ... ,
Dk
M ... M
Wm,1 Wm,n
where M is the scalar weighted mass matrix weighted by the (i,j) entry of W.
Wi,j
LocalcontributionstotheDGvariationalformmaythenbeevaluatedinaquadrature-freemannerusing
these matrices. Let Σ ,V denote vectors containing degrees of freedom for solution components σ ,v such
i i i i
that
Np
(cid:88)
σ (x,t)= Σ (t)φ (x), 1≤i≤N
i i j d
j=1
Np
(cid:88)
v (x,t)= V (t)φ (x), 1≤i≤d.
i i j
j=1
Then, the local DG formulation can be written as a block system of ordinary differential equations (ODEs)
by concatenating Σ ,V into single vectors Σ,V and using the Kronecker product ⊗
i i
d
M ∂V =(cid:88)(cid:0)AT ⊗S (cid:1)Σ+ (cid:88) (I⊗M )F
ρI ∂t i i f v
i=1 f∈∂Dk
d
∂Σ (cid:88) (cid:88)
(5) M = (A ⊗S )V + (I⊗M )F .
C−1 ∂t i i f σ
i=1 f∈∂Dk
where F and F denote degrees of freedom for velocity and stress numerical fluxes.
v σ
In order to apply standard time integration methods, we must invert M and M to isolate ∂v
ρI C−1 ∂t
and ∂σ on the left hand side. While the inversion of M and M can be parallelized from element to
∂t ρI C−1
element, doing so typically requires either the precomputation and storage of large dense matrix inverses
or the on-the-fly construction and solution of a large dense matrix system at every time-step. The former
option requires a large amount of storage, while the latter option is computationally expensive and difficult
to parallelize. This cost can be overcome for ρ,C which are constant over an element Dk, in which case
M is block diagonal with identical blocks M =ρM, while M reduces to
ρI ρ C−1
 C−1M ... C−1 M 
1,1 1,Nd
MC−1 = ... ... ... =(cid:0)C−1⊗M(cid:1).
C−1 M ... C−1 M
Nd,1 Nd,Nd
(cid:16) (cid:17)
Then, M−1 = 1M−1 = 1 M(cid:99)−1, and M−1 =C⊗M−1 =C⊗ 1M(cid:99)−1 , and each matrix inverse can be
ρ ρ Jρ C−1 J
appliedusingtheinverseofthereferencemassmatrixM(cid:99)−1 andthevaluesofρ,C,andJ overeachelement.
Applying this observation to (5) then yields the following local system of ODEs
d (cid:18) (cid:19) (cid:18) (cid:19)
∂V (cid:88) 1 (cid:88) 1
= AT ⊗D Σ+ I⊗L F
∂t ρ i i ρ f v
i=1 f∈∂Dk
d
∂Σ (cid:88) (cid:88)
= (CA ⊗D )V + (C⊗L )F ,
∂t i i f σ
i=1 f∈∂Dk
7
where we have introduced the differentiation matrix D = M−1S and lift matrix L = M−1M . For
i i f f
affine elements, both derivative and lift matrices may be applied using only geometric factors and reference
derivative and lift matrices.
Unfortunately, if ρ and C vary spatially within an element, the above approach can no longer be used
to invert M and M in an efficient and low-storage manner. For isotropic media, one way to address
ρ C−1
sub-element variations in material parameters is to diagonalize the matrix C through a change of variables
[27]. This results in a local system of ODEs with only scalar weighted mass matrices [16], which can be
treatedusingscalarweight-adjustedapproximations. Wetakeadifferentapproachinaddressingtheseissues
and approximate the matrix-weighted L2 inner product (and corresponding matrix-weighted mass matrix
M ) using a weight-adjusted approximation which is low storage, simple to invert, energy stable, and
C−1
provably high order accurate for spatially varying weights ρ,C with sufficiently regularity.
4. Weight-adjusted inner products for matrix-valued weights. Weight-adjusted inner products
are high order accurate approximations of weighted L2 inner products. These can be interpreted as general-
izationsofmasslumpingtechniques,reducingtomasslumpingwhenintegralsareevaluatedwithappropriate
quadrature rules. These weight-adjusted inner products result in weight-adjusted mass matrices, whose in-
verses approximate the inverse of a weighted L2 mass matrix.
We wish to apply weight-adjusted approximations to avoid the inversion of M and M . Approxi-
ρ C−1
mating the inverse of M can be done using weight-adjusted approximations for scalar weights [1, 2], which
ρ
we review in Section 4.1. We then extend scalar weight-adjusted approximations to matrix-valued weights
in Section 4 to approximate the inverse of M .
C−1
4.1. Scalar weight adjusted inner products. We introduce standard Lebesgue Lp norms and their
associated Lp spaces over a general domain Ω
(cid:18)(cid:90) (cid:19)1/p (cid:110) (cid:111)
(cid:107)u(cid:107) = up , Lp(Ω)= u:Ω→R, (cid:107)u(cid:107) <∞
Lp(Ω) Lp(Ω)
Ω
for 1≤p<∞. For p=∞, these are defined as
(cid:110) (cid:111)
(cid:107)u(cid:107) =inf{C ≥0:|u(x)|≤C ∀x∈Ω}, L∞(Ω)= u:Ω→R, (cid:107)u(cid:107) <∞ .
L∞(Ω) L∞(Ω)
These induce Lp Sobolev seminorms and norms of degree s
 1/p
|u| =(cid:88) (cid:107)Dαu(cid:107)p  , |u| = max(cid:107)Dαu(cid:107)
Ws,p(Ω) Lp(Ω) Ws,∞(Ω) L∞(Ω)
|α|=s
|α|=s
 1/p
(cid:107)u(cid:107) =(cid:88) (cid:107)Dαu(cid:107)p  , (cid:107)u(cid:107) = max(cid:107)Dαu(cid:107) .
Ws,p(Ω) Lp(Ω) Ws,∞(Ω) L∞(Ω)
|α|≤s
|α|≤s
where α={α ,α ,α } is a multi-index such that
1 2 3
∂α1 ∂α2 ∂α3
Dαu= u,
∂xα1 ∂yα2 ∂zα3
We define two operators T :L2(cid:0)Dk(cid:1)→PN(cid:0)Dk(cid:1) and T−1 :L2(cid:0)Dk(cid:1)→PN(cid:0)Dk(cid:1) such that
w w
T u=Π (wu), (cid:0)wT−1u,v(cid:1) =(u,v) .
w N w L2(Dk) L2(Dk)
A weighted L2 inner product (wu,v) can be approximated by a weight-adjusted inner product
L2(Dk)
(cid:16) (cid:17)
(wu,v) =(T u,v) ≈ T−1u,v .
L2(Dk) w L2(Dk) 1/w
L2(Dk)
based on the observation that T−1u ≈ uw. This approximation is made precise by the following estimates
1/w
for approximations of the product uw and weighted moments on affinely mapped elements:
8
Theorem 2 (Theorem 5 in [1]). Let Dk be quasi-regular with representative size h = diam(cid:0)Dk(cid:1). For
N ≥0, w ∈WN+1,∞(cid:0)Dk(cid:1), and u∈WN+1,2(cid:0)Dk(cid:1),
(6) (cid:107)uw−T u(cid:107) ≤C hN+1(cid:107)u(cid:107) ,
w L2(Dk) w WN+1,2(Dk)
(cid:13) (cid:13)
(7) (cid:13)uw−T−1u(cid:13) ≤C hN+1(cid:107)u(cid:107) .
(cid:13) 1/w (cid:13) w WN+1,2(Dk)
L2(Dk)
where Cw =C(cid:107)w(cid:107)L∞(Dk)(cid:13)(cid:13)w1(cid:13)(cid:13)L∞(Dk)(cid:107)w(cid:107)WN+1,∞(Dk).
Thistheorem reliesona scalarweightedinterpolationestimate derivedin[18, 1]fora generalnon-affine
element Dk.
Theorem 3 (Theorem 3.1 in [1]). Let Dk be a quasi-regular element with representative size h =
diam(cid:0)Dk(cid:1). For N ≥0, w ∈WN+1,∞(cid:0)Dk(cid:1), and u∈WN+1,2(cid:0)Dk(cid:1),
(cid:13) (cid:13)
(cid:13) 1 (cid:13)
(cid:13)(cid:13)u− wΠN(wu)(cid:13)(cid:13) ≤
L2(Dk)
ChN+1(cid:13)(cid:13)(cid:13)√J(cid:13)(cid:13)(cid:13)L∞(Dk)(cid:13)(cid:13)(cid:13)(cid:13)√1J(cid:13)(cid:13)(cid:13)(cid:13)L∞(Dk)(cid:13)(cid:13)(cid:13)(cid:13)w1(cid:13)(cid:13)(cid:13)(cid:13)L∞(Dk)(cid:107)w(cid:107)WN+1,∞(Dk)(cid:107)u(cid:107)WN+1,2(Dk).
4.2. Extension to matrix weights. We now generalize weight-adjusted inner products to the case
of matrix-valued weights. We first define appropriate generalizations of norms used in Section 4.1 to vector-
valued functions. Let Dαv denote component-wise differentiation of v with respect to a d-dimensional
multi-index α. Then, vector Lp Sobolev norms for v(x)∈Rm can be defined as
m m
|v|p =(cid:88)|v |p , (cid:107)v(cid:107)p =(cid:88)(cid:107)v (cid:107)p 1≤p<∞,
Wk,p i Wk,p Wk,p i Wk,p
i=1 i=1
|v| =max|v | , (cid:107)v(cid:107) =max(cid:107)v (cid:107) .
Wk,∞ i Wk,∞ Wk,∞ i Wk,∞
i i
The corresponding Sobolev spaces Wk,p and Wk,∞ are defined similarly to the scalar case.
Let W(x) be a matrix-valued weight function which is pointwise symmetric positive-definite
0<wmin ≤(cid:107)W(x)(cid:107)2 ≤wmax <∞, 0<w˜min ≤(cid:13)(cid:13)W−1(x)(cid:13)(cid:13)2 ≤w˜max <∞, ∀x∈Ω.
We define a kth order Sobolev norm for W(x) in terms of the induced p-norm
(cid:107)W(x)(cid:107)p = (cid:88) sup(cid:107)DαW(x)(cid:107)p
k,p,∞ p
x
|α|≤k
where DαW(x) again denotes component-wise differentiation. While this norm is not sub-multiplicative,
the following bound holds
(cid:90)
(cid:107)Wv(cid:107)pWk,p = (cid:88) (cid:107)Dα(Wv)(cid:107)pLp ≤CN (cid:88) (cid:88) (cid:13)(cid:13)(cid:0)DβW(cid:1)(cid:0)Dα−βv(cid:1)(cid:13)(cid:13)pp
|α|≤k |α|≤k|β|≤|α|
 p p
(cid:90)
(cid:88) (cid:88)
≤CN  (cid:107)(DαW)(cid:107)p  (cid:107)Dαv(cid:107)p
|α|≤k |α|≤k
≤C (cid:107)W(cid:107)p (cid:107)v(cid:107)p ,
N k,p,∞ k,p
where we have used Leibniz’s rule, Cauchy-Schwarz, and the arithmetic-geometric mean inequality.
The following theorem extends Theorem 3 to matrix weights by computing weighted interpolation esti-
mates for the quantity W−1Π (Wv).
N
9
Theorem 4. Let Dk be a quasi-regular element with representative size h = diam(cid:0)Dk(cid:1). For N ≥ 0,
W ∈(cid:0)WN+1,∞(cid:0)Dk(cid:1)(cid:1)d×d, and v ∈(cid:0)WN+1,2(cid:0)Dk(cid:1)(cid:1)d,
(cid:13)(cid:13)v−W−1ΠN(Wv)(cid:13)(cid:13)L2(Dk) ≤ChN+1(cid:13)(cid:13)(cid:13)√J(cid:13)(cid:13)(cid:13)L∞(Dk)(cid:13)(cid:13)(cid:13)(cid:13)√1J(cid:13)(cid:13)(cid:13)(cid:13)L∞(Dk)w˜max(cid:107)W(cid:107)N+1,2,∞(cid:107)v(cid:107)WN+1,2(Dk)
Proof. The proof is the similar to the scalar case. Using vector-valued versions of Bramble-Hilbert and
a scaling argument for quasi-regular elements yields
(cid:13)(cid:13)v−W−1ΠN(Wv)(cid:13)(cid:13)L2(Dk) ≤C1(cid:13)(cid:13)(cid:13)√J(cid:13)(cid:13)(cid:13)L∞(Dk)suxp(cid:13)(cid:13)W−1(cid:13)(cid:13)2(cid:107)Wv−ΠN(Wv)(cid:107)L2(Dk)
≤C1(cid:13)(cid:13)(cid:13)√J(cid:13)(cid:13)(cid:13)L∞(Dk)suxp(cid:13)(cid:13)W−1(cid:13)(cid:13)2|Wv|WN+1,2(D(cid:98))
≤C2hN+1(cid:13)(cid:13)(cid:13)√J(cid:13)(cid:13)(cid:13)L∞(Dk)(cid:13)(cid:13)(cid:13)(cid:13)√1J(cid:13)(cid:13)(cid:13)(cid:13)L∞(Dk)suxp(cid:13)(cid:13)W−1(cid:13)(cid:13)2(cid:107)Wv(cid:107)WN+1,2(Dk)
≤C3hN+1(cid:13)(cid:13)(cid:13)√J(cid:13)(cid:13)(cid:13)L∞(Dk)(cid:13)(cid:13)(cid:13)(cid:13)√1J(cid:13)(cid:13)(cid:13)(cid:13)L∞(Dk)w˜max(cid:107)W(cid:107)N+1,2,∞(cid:107)v(cid:107)WN+1,2(Dk).
4.2.1. Weight-adjusted approximations with matrix weights. Let Π u be defined as the L2
N
projection applied to each component of the vector-valued function u. We then define the operator T
W
analogously to the scalar case
T v =Π (Wv).
W N
The inverse operator T−1 is defined implicitly via
W
(cid:0)WT−1v,w(cid:1) =(v,w) .
W L2(Dk) L2(Dk)
This definition is a straightforward generalization of T−1 to matrix-valued weights W. Conveniently, all
w
properties of T ,T−1 for scalar w(x) [1] carry over to matrix weights W(x) as well.
w w
Lemma 5. Let Π denote the component-wise L2 projection, and let W ∈(L∞)m. Then, T satisfies
N W
the following properties:
1. T−1T =Π
W W N
2. Π T−1 =T−1Π =T−1
N W W N W
3. (cid:13)(cid:13)TW−1(cid:13)(cid:13)L2(Dk) ≤w˜max.
4. (cid:0)T−1v,w(cid:1) forms an inner product on (cid:0)PN(cid:1)m×(cid:0)PN(cid:1)m, which is equivalent to the L2 inner
W L2(Dk)
product with equivalence constants C =w˜ ,C =w˜ .
1 min 2 max
Proof. The proofs of properties 1 and 2 are consequences of the definition of T , and are identical to
W
proofs for the scalar case. Property 3 is a straightforward extension from the scalar case. Let v ∈ (cid:0)PN(cid:1)m.
Then,
(cid:0)TW−1v,v(cid:1)L2(Dk) =(cid:0)W−1WTW−1v,v(cid:1)L2(Dk) ≤sup(cid:13)(cid:13)W−1(x)(cid:13)(cid:13)(cid:0)WTW−1v,v(cid:1)L2(Dk) =w˜max(cid:107)v(cid:107)2L2(Dk).
x
Property 4 then simply requires the lower bound
(cid:0)TW−1v,v(cid:1)L2(Dk) =(cid:0)W−1WTW−1v,v(cid:1)L2(Dk) ≥inxf(cid:13)(cid:13)W−1(x)(cid:13)(cid:13)(cid:107)v(cid:107)2L2(Dk) =w˜min(cid:107)v(cid:107)2L2(Dk).
Using Theorem 4, we may also show that the matrix-valued weight-adjusted inner product is also high
order accurate for sufficiently regular W.
10