Exactly divergence-free methods for induction and related equations - - PowerPoint PPT Presentation
Exactly divergence-free methods for induction and related equations Praveen Chandrashekar praveen@math.tifrbng.res.in Center for Applicable Mathematics Tata Institute of Fundamental Research Bangalore-560065, India http://cpraveen.github.io
Praveen Chandrashekar praveen@math.tifrbng.res.in
Center for Applicable Mathematics Tata Institute of Fundamental Research Bangalore-560065, India http://cpraveen.github.io
7 May 2018
Supported by Airbus Foundation Chair at TIFR-CAM, Bangalore http://math.tifrbng.res.in/airbus-chair
1 / 41
1 DG scheme using Raviart-Thomas polynomials 2 Divergence-free reconstruction approach
Joint work with Dinshaw S. Balsara, Univ. of Notre Dame Rakesh Kumar, TIFR-CAM Arijit Hazra, TIFR-CAM
2 / 41
Linear hyperbolic system
∂B ∂t + ∇ × E = 0, ∂D ∂t − ∇ × H = −J B = magnetic flux density D = electric flux density E = electric field H = magnetic field J = electric current density B = µH, D = εE, J = σE µ, ε ∈ R3×3 symmetric ε = permittivity tensor µ = magnetic permeability tensor σ = conductivity ∇ · B = 0, ∇ · D = ρ (electric charge density), ∂ρ ∂t + ∇ · J = 0
3 / 41
Nonlinear hyperbolic system
Compressible Euler equations with Lorentz force ∂ρ ∂t + ∇ · (ρv) = ∂ρv ∂t + ∇ · (pT I + ρv ⊗ v − B ⊗ B) = ∂E ∂t + ∇ · ((E + pT )v − (v · B)B) = ∂B ∂t − ∇ × (v × B) = Magnetic monopoles do not exist: = ⇒ ∇ · B = 0
4 / 41
5 / 41
Find B ∈ H(div, Ω) H(div, Ω) = {B ∈ L2(Ω) : div(B) ∈ L2(Ω)} To approximate functions in H(div, Ω) on a mesh Th, the normal trace B · n must be continuous. Pk(x), Pk(y): 1-D polynomials of degree at most k Qr,s(x, y): tensor product polynomials Qr,s(x, y) = span{xiyj, 0 ≤ i ≤ r, 0 ≤ j ≤ s} For k ≥ 0, the Raviart-Thomas space of vector functions is defined as RTk = Qk+1,k × Qk,k+1, dim(RTk) = 2(k + 1)(k + 2)
6 / 41
Two consequences:
div(Bh) ∈ Qk,k(x, y) =: Qk(x, y)
Bh
x(±∆x/2, y) ∈ Pk(y),
Bh
y (x, ±∆y/2) ∈ Pk(x)
For doing the numerical computations, it is useful to map each cell to a reference cell. {ξi, 0 ≤ i ≤ k + 1} = Gauss-Lobatto-Legendre (GLL) nodes {ˆ ξi, 0 ≤ i ≤ k} = Gauss-Legendre (GL) nodes
7 / 41
Let φi and ˆ φi be the corresponding 1-D Lagrange polynomials. Then the magnetic field is given by Bh
x(ξ, η) = k+1
k
(Bx)ijφi(ξ)ˆ φj(η), Bh
y (ξ, η) = k
k+1
(By)ij ˆ φi(ξ)φj(η)
Bx ∈ Q1,0 By ∈ Q0,1
Location of dofs of Raviart-Thomas polynomial for k = 0
8 / 41
Bx ∈ Q2,1 By ∈ Q1,2
Location of dofs of Raviart-Thomas polynomial for k = 1
Bx ∈ Q3,2 By ∈ Q2,3
Location of dofs of Raviart-Thomas polynomial for k = 2
9 / 41
Our choice of nodes ensures that the normal component of the magnetic field is continuous on the cell faces. We have the error estimates on Cartesian meshes [?], [?] B − BhL2(Ω) ≤ Chk+1|B|Hk+1(Ω) div(B) − div(Bh)L2(Ω) ≤ Chk+1|div(B)|Hk+1(Ω)
10 / 41
e+
x
e−
x
e−
y
e+
y
C 1 2 3
The edge moments are given by
x
Bh
xφdy
∀φ ∈ Pk(y) and
y
Bh
y φdx
∀φ ∈ Pk(x) The cell moments are given by
Bh
xψdxdy
∀ψ ∈ ∂xQk(x, y) := Qk−1,k(x, y)
11 / 41
and
Bh
y ψdxdy
∀ψ ∈ ∂yQk(x, y) := Qk,k−1(x, y) Note that dim Pk(x) = dim Pk(y) = k + 1 and dim ∂xQk(x, y) = dim ∂yQk(x, y) = k(k + 1) so that we have in total 4(k + 1) + 2k(k + 1) = 2(k + 1)(k + 2) = dim RTk The moments on the edges e∓
x uniquely determine the restriction of
Bh
x on those edges, and similarly the moments on e∓ y uniquely
determine the restriction of Bh
y on the corresponding edges. This
ensures continuity of the normal component of Bh on all the edges.
12 / 41
Theorem
If all the moments are zero for any cell C, then Bh ≡ 0 inside that cell.
Theorem
Let Bh ∈ RTk satisfy the moments
x
Bh
xφdy
=
x
Bxφdy ∀φ ∈ Pk(y) (1)
y
Bh
y φdx
=
y
Byφdx ∀φ ∈ Pk(x) (2)
Bh
xψdxdy
=
Bxψdxdy ∀ψ ∈ ∂xQk(x, y) (3)
Bh
y ψdxdy
=
Byψdxdy ∀ψ ∈ ∂yQk(x, y) (4) for a given vector field B ∈ H(div, Ω). If div(B) ≡ 0 then div(Bh) ≡ 0.
13 / 41
∂B ∂t + ∇ × E = −M Divergence evolves according to ∂ ∂t(divB) + ∇ · M = 0 (5) If M ≡ 0, then div(B) = constant. In 2-D ∂Bx ∂t + ∂Ez ∂y = −Mx, ∂By ∂t − ∂Ez ∂x = −My and for the induction equation Ez = vyBx − vxBy
14 / 41
Constructing Bh from the edge and cell moments allowed us to get divergence-free approximation. We will construct a scheme to evolve the same moments in time.
e+
x
e−
x
e−
y
e+
y
C 1 2 3
Consider the right edge e+
x . Multiply by test
function φ ∈ Pk(y)
x
∂Bx ∂t + ∂E ∂y
x
Mxφdy and integrate by parts in second term
x
∂Bx ∂t φdy −
x
E ∂φ ∂y dy + (Eφ)3 − (Eφ)1 = −
x
Mxφdy
15 / 41
Edge moments are evolved by
x
∂Bh
x
∂t φdy −
x
ˆ E ∂φ ∂y dy + [ ˜ Eφ]e∓
x = −
x
ˆ Mxφdy, ∀φ ∈ Pk(y)
y
∂Bh
y
∂t φdx +
y
ˆ E ∂φ ∂xdx − [ ˜ Eφ]e∓
y = −
y
ˆ Myφdy, ∀φ ∈ Pk(x) where ˆ E = numerical flux from a 1-D Riemann solver required on the faces ˜ E = numerical flux from a multi-D Riemann solver needed at vertices [ ˜ Eφ]e−
x = ( ˜
Eφ)2 − ( ˜ Eφ)0, [ ˜ Eφ]e+
x = ( ˜
Eφ)3 − ( ˜ Eφ)1 [ ˜ Eφ]e−
y = ( ˜
Eφ)1 − ( ˜ Eφ)0, [ ˜ Eφ]e+
y = ( ˜
Eφ)3 − ( ˜ Eφ)2
16 / 41
The cells moments are evolved by the following standard DG scheme
∂Bh
x
∂t ψdxdy −
E ∂ψ ∂y dxdy +
ˆ Eψnyds = −
Mxψdxdy, ∀ψ ∈ ∂xQk(x, y)
∂Bh
y
∂t ψdxdy +
E ∂ψ ∂x dxdy −
ˆ Eψnxds = −
Myψdxdy, ∀ψ ∈ ∂yQk(x, y)
Note that the same 1-D numerical flux ˆ E is used in both the edge and cell moment equations whereas the vertex numerical flux ˜ E is needed only in the edge moment equations.
17 / 41
Theorem
If M = 0, the DG scheme preserves the divergence of the magnetic field. Proof: For any φ ∈ Qk take ψ = ∂xφ and ψ = ∂yφ in the two cell moment equations respectively and add them together to obtain
x
∂t ∂xφ + ∂Bh
y
∂t ∂yφ
ˆ E(nx∂yφ − ny∂xφ)ds = 0 Two of the cell integrals cancel since ∂x∂yφ = ∂y∂xφ. Performing an integration by parts in the first term, we obtain −
φ ∂ ∂tdiv(Bh)dxdy+
φ ∂ ∂t(Bh·n)ds−
ˆ E(nx∂yφ−ny∂xφ)ds = 0 The last two terms can be re-arranged as
x
φ∂Bh
x
∂t dy −
x
φ∂Bh
x
∂t dy +
y
φ∂Bh
y
∂t dx −
y
φ∂Bh
y
∂t dx −
x
ˆ E∂yφdy +
x
ˆ E∂yφdy +
y
ˆ E∂xφdx −
y
ˆ E∂xφdx
18 / 41
The restriction of φ on each edge is a one dimensional polynomial of degree k. We can use the edge moment equations to simplify as
= − [ ˜ Eφ]e+
x + [ ˜
Eφ]e−
x + [ ˜
Eφ]e+
y − [ ˜
Eφ]e−
y
= − ( ˜ Eφ)3 + ( ˜ Eφ)1 + ( ˜ Eφ)2 − ( ˜ Eφ)0 + ( ˜ Eφ)3 − ( ˜ Eφ)2 − ( ˜ Eφ)1 + ( ˜ Eφ)0 = 0
Hence we have
φ ∂ ∂tdiv(Bh)dxdy = 0 ∀φ ∈ Qk(x, y) Since div(Bh) ∈ Qk(x, y) we conclude that the divergence is preserved by the numerical scheme.
19 / 41
Theorem (M = 0)
The divergence evolves consistently with equation (5) in the sense that
φ ∂ ∂tdiv(Bh)dxdy −
M ·∇φdxdy +
φ ˆ M ·nds = 0, ∀φ ∈ Qk Since div(Bh) ∈ Qk, we can expect div(Bh) to be accurate to O(hk+1). In case of Maxwell equations, if we want to compute charge density ρ = ∇ · D we can get ρ to same accuracy as the vector fields.
20 / 41
Using the zero divergence condition, rewrite the induction equation ∂Bx ∂t +v·∇Bx+Bx ∂vy ∂y −By ∂vx ∂y = 0, ∂By ∂t +v·∇By+By ∂vx ∂x −Bx ∂vy ∂x = 0 The characteristics are the integral curves of v. The 1-D numerical flux is given by ˆ E =
if v · n > 0 ER
For example, across the face e∓
x , the flux is given by
ˆ E =
y
if vx > 0 vyBx − vxBR
y
21 / 41
BDL = (BD
x , BL y )
BDR = (BD
x , BR y )
BUL = (BU
x , BL y )
BUR = (BU
x , BR y )
The corner flux is given by ˜ E = EDL if vx > 0, vy > 0 EUL if vx > 0, vy < 0 EDR if vx < 0, vy > 0 EUR if vx < 0, vy < 0 which can be written in compact form as ˜ E = vy 2 (BU
x + BD x ) − vx
2 (BL
y + BR y ) − |vy|
2
x − BD x
2
y − BL y
An equivalent expression is given by [?] ˜ E =vy 4 (BUL
x
+ BUR
x
+ BDL
x
+ BDR
x
) − vx 4 (BUL
y
+ BUR
y
+ BDL
y
+ BDR
y
) − |vy| 2 BUL
x
+ BUR
x
2 − BDL
x
+ BDR
x
2
2
y
+ BDR
y
2 − BUL
y
+ BDL
y
2
x
= BDR
x
, etc.
23 / 41
ˆ E: Given boundary condition B∗ ˆ E =
v · n > 0 E(B∗, v)
˜ E: Corner flux for boundary point: use numerical flux
BDL = (B∗
x, B∗ y)
BDR = (BD
x , BR y )
BUL = (B∗
x, B∗ y)
BUR = (BU
x , BR y )
(a) BDL = (BD
x , BL y )
BDR = (BD
x , BL y )
BUL = (BU
x , BL y )
BUR = (BU
x , BL y )
(b)
Vertex states at boundary: (a) inflow vertex on left side of domain, (b)
24 / 41
4 1 5 8 10 9 11 2 6 3 7 Bx By
Mx Mx My My Nx
l
Nx
r
Qx Ny
b
Ny
t
Qy d dt[dofs] = rhs
25 / 41
Mx d dt[Bx]e = Re = ⇒ d dt[Bx]e = ˜ Re := (Mx)−1Re
Nx
l
d dt[Bx]e−
x + Nx
r
d dt[Bx]e+
x + Qx d
dt[Bx]int = Rint ⇓ d dt[Bx]int = (Qx)−1[Rx
int − Nx l ˜
Re−
x − Nx
r ˜
Re+
x ]
26 / 41
27 / 41
The initial condition is given by B0 = (∂yΦ, −∂xΦ) where Φ(x, y) = 1 10 exp[−20((x − 1/2)2 + y2)] and the velocity field is v = (y, −x). The exact solution is a pure rotation
B(r, t) = R(t)B0(R(−t)r), R(t) = cos t − sin t sin t cos t
28 / 41
We compute the numerical solution on the computational domain [−1, +1] × [−1, +1] upto a final time of T = 2π at which time the solution comes back to the initial condition. (a) (b) (c) Contour of |Bh| for Test 2a, 10 contours between 0 and 0.3867: (a) initial, (b) final, k = 1, (c) final, k = 2
29 / 41
h Bh − BL2(Ω) div(Bh)L2(Ω) 0.0312 2.1427e-03
0.0156 3.2571e-04 2.71 1.8566e-13 0.0078 5.9640e-05 2.45 5.8486e-13 0.0039 1.3209e-05 2.17 1.8853e-12 k = 1 h Bh − BL2(Ω) div(Bh)L2(Ω) 0.0625 2.4003e-04
0.0312 2.5212e-05 3.25 1.4299e-13 0.0156 3.0946e-06 3.02 4.5663e-13 0.0078 3.8448e-07 3.00 1.5058e-12 k = 2
30 / 41
We compute the numerical solution on the computational domain [0, 1] × [0, 1] upto a final time of T = π/4. (a) (b) (c) Contour of |Bh| for Test 2b, 10 contours between 0 and 0.3867: (a) initial, (b) final, k = 1, (c) final, k = 2
31 / 41
h Bh − BL2(Ω) div(Bh)L2(Ω) 0.0312 6.5882e-04
0.0156 1.4979e-04 2.13 9.8666e-14 0.0078 3.6394e-05 2.04 3.2902e-13 0.0039 9.0308e-06 2.01 1.1356e-12 k = 1 h Bh − BL2(Ω) div(Bh)L2(Ω) 0.0625 1.4110e-04
0.0312 1.7238e-05 3.03 7.9129e-14 0.0156 2.1442e-06 3.00 2.5910e-13 0.0078 2.6749e-07 3.00 9.2720e-13 k = 2
32 / 41
The exact solution is taken to be B(x, y, t) = cos t − sin t sin t cos t
where B0 = ∇φ, φ = 1 10 exp(−20(x2 + y2)) so that ∇ · B = 0, and the velocity field is taken as v = ∇⊤ψ, ψ = 1 π sin(πx) sin(πy) The right hand side source term M is computed from the above solution using the formula M = − ∂B
∂t + ∇ × (v × B). The problem is solved on
the domain [−1, +1] × [−1, +1] until a final time of T = 2π.
33 / 41
(a) (b) (c) Contour plot of Bx for Test 3 showing 16 contours between -0.3838 and +0.3838: (a) initial condition, (b) final, k = 1 and (c) final, k = 2
34 / 41
h Bh − BL2(Ω) div(B) − div(Bh)L2(Ω) 0.0312 8.5550e-04
1.8915e-04 2.17 1.7299e-03 1.99 0.0078 3.8730e-05 2.29 4.3267e-04 1.99 0.0039 7.8346e-06 2.30 1.0818e-04 1.99 k = 1 h Bh − BL2(Ω) div(B) − div(Bh)L2(Ω) 0.0625 3.4775e-04
3.3408e-05 3.38 2.3550e-04 2.99 0.0156 3.0287e-06 3.46 2.9491e-05 2.99 0.0078 2.7345e-07 3.47 3.6881e-06 2.99 k = 2
35 / 41
We take the potential Φ(x, y) =
x > y
and the velocity field is v = (1, 2). This leads to a discontinuous magnetic field with the discontinuity along the line x = y B0 =
x > y (0, 0) x < y The exact solution is given by B(x, y, t) = B0(x − t, y − 2t) Animation
36 / 41
(a) IC (b) k = 0 (c) k = 1 (d) k = 2
37 / 41
38 / 41
39 / 41
∂B ∂t + ∇ × E = −∇ × (µh∇ × B)
◮ Maxwell equations ◮ Magnetohydrodynamics 40 / 41
41 / 41
Define U = [ρ, ρv, ρe, Bz]⊤, B = (Bx, By) the ideal MHD equations are ∂U ∂t + ∇ · F (U, B) = 0, ∂Bx ∂t + ∂Ez ∂y = 0, ∂By ∂t − ∂Ez ∂x = 0 where Ez = vyBx − vxBy Magnetic monopoles do not exist: = ⇒ ∇ · B = 0
1 / 18
B · n on the faces must be continuous: directly approximate B · n by polynomial functions Pk on the faces.
B+
x (η)
B−
x (η)
B−
y (ξ)
B+
y (ξ)
e+
x
e−
x
e−
y
e+
y
C 1 2 3
∆x, ∆y → [− 1
2, + 1 2] × [− 1 2, + 1 2]
ξ = x − xc ∆x , η = y − yc ∆y B±
x (η) = k
a±
j φj(η)
B±
y (ξ) = k
b±
j φj(ξ)
2 / 18
φj are mutually orthogonal polynomials on the interval [− 1
2, 1 2] given by
φ0(ξ) = 1, φ1(ξ) = ξ, φ2(ξ) = ξ2 − 1
12,
φ3(ξ) = ξ3 − 3
20ξ
φ4(ξ) = ξ4 − 3
14ξ2 + 3 560,
etc. Note
2
− 1
2
B±
x (η)dη = a± 0 ,
2
− 1
2
B±
y (ξ)dξ = b±
Given the data B±
x , B± y on the faces, we want to recontruct a
divergence free vector field (Bx(ξ, η), By(ξ, η)) inside the cell. Necessary condition
B · n = (a+
0 − a− 0 )∆y + (b+ 0 − b− 0 )∆x = 0
(1)
3 / 18
B+
x (η)
B−
x (η)
B−
y (ξ)
B+
y (ξ)
Bx(ξ, η) =? By(ξ, η) =?
We can try to satisfy the conditions Bx(± 1
2, η) = B± x (η),
∀η ∈ [− 1
2, + 1 2]
By(ξ, ± 1
2) = B± y (ξ),
∀ξ ∈ [− 1
2, + 1 2]
In addition, we have to make the vector field divergence-free ∇ · B = 1 ∆x ∂Bx ∂ξ + 1 ∆y ∂By ∂η = 0, ∀ ξ, η ∈ [− 1
2, + 1 2]
4 / 18
h where V k h = RT(k)/BDM(k)/BDFM(k)
h
◮ We can remove basis functions in V k
h \ Pk
◮ face solution determines reconstruction ◮ all three spaces yield same reconstruction
5 / 18
Constant approximation on the faces B±
x (η) = a± 0 ,
B±
y (ξ) = b±
Vector field inside the cell (dim = 2 + 2 = 4) Bx(ξ, η) = a00 + a10ξ, By(ξ, η) = b00 + b01η Solution is given by a00 = 1 2(a−
0 + a+ 0 ),
b00 = 1 2(b−
0 + b+ 0 ),
a10 = a+
0 − a− 0 ,
b01 = b+
0 − b−
6 / 18
The face solution is a polynomial of degree one B±
x (η) = a± 0 + a± 1 η,
B±
x (η) = b± 0 + b± 1 ξ
Vector field inside the cell (dim = 5 + 5 = 10) Bx(ξ, η) = a00 + a10ξ + a01η + a20(ξ2 − 1
12) + a11ξη
By(ξ, η) = b00 + b10ξ + b01η + b11ξη + b02(η2 − 1
12)
Reconstruction is given by
a00 = 1 2(a−
0 + a+ 0 ) + 1 12(b+ 1 − b− 1 )∆x
∆y a10 = a+
0 − a−
a01 = 1 2(a−
1 + a+ 1 )
a20 = 1 2(b−
1 − b+ 1 )∆x
∆y a11 = a+
1 − a− 1
b00 = 1 2(b−
0 + b+ 0 ) + 1 12(a+ 1 − a− 1 ) ∆y
∆x b10 = 1 2(b−
1 + b+ 1 )
b01 = b+
0 − b−
b02 = 1 2(a−
1 − a+ 1 ) ∆y
∆x b11 = b+
1 − b− 1 7 / 18
The face solution is a polynomial of degree two B±
x (η) = a± 0 + a± 1 η + a± 2 (η2 − 1 12),
B±
y (ξ) = b± 0 + b± 1 ξ + b± 2 (ξ2 − 1 12)
The field inside the cell is approximated by (dim = 8 + 8 = 16) Bx(ξ, η) = a00 + a10ξ + a01η + a20(ξ2 − 1
12) + a11ξη + a02(η2 − 1 12)+
a30(ξ3 − 3
20ξ) + a12ξ(η2 − 1 12)
By(ξ, η) = b00 + b10ξ + b01η + b20(ξ2 − 1
12) + b11ξη + b02(η2 − 1 12)+
b21(ξ2 − 1
12)η + b03(η3 − 3 20η)
8 / 18
The reconstruction is given by
a00 = 1 2 (a−
0 + a+ 0 ) + 1 12 (b+ 1 − b− 1 )
∆x ∆y a10 = a+
0 − a− 0 +
1 30 (b+
2 − b− 2 )
∆x ∆y a11 = a+
1 − a− 1
a02 = 1 2 (a−
2 + a+ 2 )
a20 = − 1 2 (b+
1 − b− 1 )
∆x ∆y a12 = a+
2 − a− 2
a30 = − 1 3 (b+
2 − b− 2 )
∆x ∆y a01 = 1 2 (a−
1 + a+ 1 )
b00 = 1 2 (b−
0 + b+ 0 ) + 1 12 (a+ 1 − a− 1 )
∆y ∆x b01 = b+
0 − b− 0 +
1 30 (a+
2 − a− 2 )
∆y ∆x b11 = b+
1 − b− 1
b20 = 1 2 (b−
2 + b+ 2 )
b02 = − 1 2 (a+
1 − a− 1 )
∆y ∆x b03 = − 1 3 (a+
2 − a− 2 )
∆y ∆x b21 = b+
2 − b− 2
b10 = 1 2 (b−
1 + b+ 1 )
9 / 18
The face solution is a polynomial of degree three B±
x (η) = a± 0 φ0(η) + a± 1 φ1(η) + a± 2 φ2(η) + a± 3 φ3(η)
B±
x (η) = b± 0 φ0(ξ) + b± 1 φ1(ξ) + b± 2 φ2(ξ) + b± 3 φ3(ξ)
The vector field has the form (dim = 12 + 12 = 24)
Bx(ξ, η) = a00 + a10ξ + a01η + a20(ξ2 −
1 12) + a11ξη + a02(η2 − 1 12)+
a30(ξ3 −
3 20ξ) + a21(ξ2 − 1 12)η + a12ξ(η2 − 1 12) + a03(η3 − 3 20η)+
a40(ξ4 −
3 14ξ2 + 3 560) + a13ξ(η3 − 3 20η)
By(ξ, η) = b00 + b10ξ + b01η + b20(ξ2 −
1 12) + b11ξη + b02(η2 − 1 12)+
b30(ξ3 −
3 20ξ) + b21(ξ2 − 1 12)η + b12ξ(η2 − 1 12) + b03(η3 − 3 20η)+
b31(ξ3 −
3 20ξ)η + b04(η4 − 3 14η2 + 3 560)
In addition, we assume that ω = b10 − a01 is given.
10 / 18
The reconstruction is given by
a00 = 1 2 (a−
0 + a+ 0 ) +
1 12 (b+
1 − b− 1 )
∆x ∆y a10 = a+
0 − a− 0 +
1 30 (b+
2 − b− 2 )
∆x ∆y a20 = − 1 2 (b+
1 − b− 1 )
∆x ∆y + 3 140 (b+
3 − b− 3 )
∆x ∆y a30 = − 1 3 (b+
2 − b− 2 )
∆x ∆y a03 = 1 2 (a−
3 + a+ 3 )
a12 = a+
2 − a− 2
a40 = − 1 4 (b+
3 − b− 3 )
∆x ∆y a13 = a+
3 − a− 3
a02 = 1 2 (a−
2 + a+ 2 )
a11 = a+
1 − a− 1
a01 = 1 1 + ∆y
∆x
∆y ∆x + r2 − ω
= 6(r1 − a01) b00 = 1 2 (b−
0 + b+ 0 ) +
1 12 (a+
1 − a− 1 )
∆y ∆x b01 = b+
0 − b− 0 +
1 30 (a+
2 − a− 2 )
∆y ∆x b02 = − 1 2 (a+
1 − a− 1 )
∆y ∆x + 3 140 (a+
3 − a− 3 )
∆y ∆x b30 = 1 2 (b−
3 + b+ 3 )
b03 = − 1 3 (a+
2 − a− 2 )
∆y ∆x b21 = b+
2 − b− 2
b31 = b+
3 − b− 3
b04 = − 1 4 (a+
3 − a− 3 )
∆y ∆x b20 = 1 2 (b−
2 + b+ 2 )
b11 = b+
1 − b− 1
b10 = ω + a01 b12 = 6(r2 − b10) 11 / 18
∂t + ∂Ez ∂y = 0
Choose k + 1 Gauss-Legendre points {ηi, 0 ≤ i ≤ k} on the face, and interpolate Ez using Lagrange polynomials Eδ
z(η) = k
( ˆ Ez)jℓj(η) ˆ Ez is obtained by solving 2-state Riemann problem for full MHD system with data (U L, Bx, BL
y ) and (U R, Bx, BR y ).
˜ Et
z
˜ Eb
z
( ˆ Ez)0 ( ˆ Ez)1 ( ˆ Ez)2
BL
y , BR y obtained from div-free
reconstruction. ˜ Ez is obtained from solving 4-state Riemann problem for full MHD.
12 / 18
∂t + ∂Ez ∂y = 0
We correct the above interpolant as Ec
z(η) = Eδ z(η) + [ ˜
Eb
z − Eδ z(−1/2)]gl(η) + [ ˜
Et
z − Eδ z(+1/2)]gr(η)
where the correction functions gl, gr are polynomials of degree k + 1 and have the property gl(−1/2) = gr(+1/2) = 1, gl(+1/2) = gl(−1/2) = 0 This implies that Ec
z(−1/2) = ˜
Eb
z,
Ec
z(+1/2) = ˜
Et
z,
Ec
z ∈ Pk+1
Following Huynh, the correction functions will be taken to be Radau polynomials of degree k + 1 given by gl(η) = (−1)k 2 [Lk(2η) − Lk+1(2η)], gr(η) = 1 2[Lk(2η) + Lk+1(2η)]
13 / 18
∂t + ∂Ez ∂y = 0
where Lk : [−1, +1] → R is the Legendre polynomial of degree k. Normal component Bx updated using a collocation approach ∂Bx ∂t (ηi) = − 1 ∆y ∂Ec
z
∂η (ηi), 0 ≤ i ≤ k (2) We can use DG on the faces, but FR is useful to deal with quadtree and
hanging nodes.
14 / 18
ω = b10 − a01 = 12 (Byξ − Bxη)dξdη Derive an evolution equation for ω using the induction equation
1 12 dω dt = ∂By ∂t ξ − ∂Bx ∂t η
1 ∆x ∂Ez ∂ξ ξ + 1 ∆y ∂Ez ∂η η
Performing an integration by parts and using a numerical flux on the faces which is based on a 1-D Riemann solver, we obtain a semi-discrete DG scheme 1 12 dω dt = 1 ∆x
2
2
− 1
2
ˆ Ex−
z dη + 1
2
2
− 1
2
ˆ Ex+
z dη −
2
− 1
2
2
− 1
2
Ezdξdη
∆y
2
2
− 1
2
ˆ Ey−
z dξ + 1
2
2
− 1
2
ˆ Ey+
z dξ −
2
− 1
2
2
− 1
2
Ezdξdη
B+
x (η, t)
B−
x (η, t)
B−
y (ξ, t)
B+
y (ξ, t)
U(ξ, η, t) ω(t)
On the faces B±
x ∈ Pk,
B±
y ∈ Pk
At fourth order ω ∈ R Approximate U inside cells by
U(ξ, η, t) =
Uj(t)Φj(ξ, η) ∈ Pk
16 / 18
C
∂U ∂t Φidξdη −
C
Fx(U, B) ∆x ∂Φi ∂ξ + Fy(U, B) ∆y ∂Φi ∂η
+ 1 ∆x 1
2 − 1 2
ˆ F x+
x
Φi( 1
2, η)dη − 1
∆x 1
2 − 1 2
ˆ F x−
x
Φi(− 1
2, η)dη
+ 1 ∆y 1
2 − 1 2
ˆ F y+
y
Φi(ξ, 1
2)dξ − 1
∆y 1
2 − 1 2
ˆ F y−
y
Φi(ξ, − 1
2)dξ = 0
B is obtained from div-free reconstruction. ˆ Fx, ˆ Fy are obtained from 2-state Riemann problem on the faces. Choose (k + 1)-point GL quadrature on faces, re-use numerical flux from FR scheme on the faces. WENO-type limiter for solution on faces, and solution inside the cell.
17 / 18
Algorithm 1: Divergence-free reconstruction scheme
Allocate memory for all variables; Set initial condition on the faces and cells; t = 0; while t < T do Copy current solution into old solution; Compute time step ∆t; for each RK stage do Loop over cells and do divergence-free reconstruction; Loop over vertices and compute vertex flux; Loop over faces and compute all face integrals (FR and DG); Loop over cells and compute all cell integrals; Update solution to next stage; Apply WENO limiter; end t = t + ∆t; end
18 / 18