Discontinuous Galerkin Methods for anisotropic diffusion Praveen - - PowerPoint PPT Presentation
Discontinuous Galerkin Methods for anisotropic diffusion Praveen Chandrashekar praveen@tifrbng.res.in Center for Applicable Mathematics Tata Institute of Fundamental Research Bangalore-560065, India http://cpraveen.github.io Oberseminar,
Praveen Chandrashekar praveen@tifrbng.res.in
Center for Applicable Mathematics Tata Institute of Fundamental Research Bangalore-560065, India http://cpraveen.github.io
Oberseminar, Dept. of Mathematics, Univ. of W¨ urzburg 21 July, 2016 Supported by Airbus Foundation Chair at TIFR-CAM, Bangalore http://math.tifrbng.res.in/airbus-chair
1 / 57
2.1. Equations of MHD with Anisotropic Heat Conduction The equations of MHD with the addition of the heat flux, Q, and a vertical gravitational acceleration, g = −g0ˆ z are ∂ρ ∂t + ∇ · (ρv) = 0, (1) ∂(ρv) ∂t + ∇ ·
8π
4π
(2) ∂B ∂t + ∇ × (v × B) = 0, (3) ∂E ∂t + ∇ ·
8π
4π
(4) where the symbols have their usual meaning. The total energy E is given as E = + ρv · v 2 + B · B 8π , (5) with the internal energy, = p/(γ − 1). For this paper, we assume γ = 5/3 throughout.
Parrish & Stone, http://arxiv.org/abs/astro-ph/0507212
2 / 57
The heat flux contains contributions both from electron motions (which are constrained to move primarily along field lines) and from isotropic transport which may arise due to photons or particle collisions which drive cross field diffusion. Thus, Q = QC + QR, where QC = −χCˆ bˆ b · ∇T, (6) QR = −χR∇T, (7) where χC is the Spitzer Coulombic conductivity (Spitzer 1962), ˆ b is a unit vector in the direction of the magnetic field, and χR is the coefficient of isotropic conductivity, ostensibly due to radiation. We consider both χC and χR as free parameters in the problem and will vary both independently.
Parrish & Stone, http://arxiv.org/abs/astro-ph/0507212
3 / 57
2 THERMAL CONDUCTION We start with the energy conservation equation which can be written as ρ@u @t + r · j = 0, (4) where u is the gas internal energy per unit mass, j is the di- rectional heat flux and ρ is the gas density. If the conduction is isotropic, then the heat flux j is opposite to the direction
j = −κsp(T) rT. (5) In the presence of the magnetic fields, electrons prefer- entially move along the magnetic field lines. Therefore the heat flux is modified as j = −κ[b(b · rT)], (10) where b = B/|B|. The coefficient of conduction perpendic- ular to the magnetic fields is set to zero (see Arth et al. 2014 for a discussion on non vanishing perpendicular conduction coefficients). We can also rewrite Equation (10) as j = − κ cv [b(b · ru)], (11) by replacing temperature for an ideal gas with the internal energy per unit mass u = kBT (γ − 1)µ = cvT, (12) where γ is the adiabatic index, and µ is the mean molec- ular weight. The value of µ depends on the temperature, ionization state and the composition of the gas. Since we are dealing mainly with a high temperature, low metallic- ity plasma of primordial composition we assume that the gas is fully ionised and consequently set the mean molecular weight to the constant value µ ∼ 0.6 mp. The final equation for anisotropic thermal conduction then becomes @u @t = 1 cvρr · [κb(b · r)u], (13) and below we focus on studying numerical solutions schemes for this partial differential equation.
Kannan et al., Accurately simulating anisotropic thermal conduction on a moving mesh (2015) http://arxiv.org/abs/1512.03053
4 / 57
∂θ ∂t = χ∂2θ ∂x2 Usual finite difference θn+1
j
= θn
j + χ∆t
h2 (θn
j−1 − 2θn j + θn j+1)
Stable in L2 and L∞ norms if σ = χ∆t h2 ≤ 1 2
5 / 57
∂θ ∂t + ∇ · q = 0, in Ω (1) where the heat flux vector is given by [7], [6] q(θ) = −(χ − χ⊥)b(b · ∇θ) − χ⊥∇θ The vector field b has unit magnitude and is assumed to be given, and χ⊥ ≥ 0, χ − χ⊥ > 0 are taken to be constant In many applications χ ≫ χ⊥ and the diffusion is mainly along the vector field b which is refered to as anisotropic diffusion, while the term containing χ⊥ leads to isotropic diffusion. For simplicity of notation, let us define κ = χ − χ⊥ > 0, η = χ⊥ ≥ 0
6 / 57
The heat flux vector can be written as q = −D∇θ where D = ηI + κbb⊤ The diffusion matrix is positive definite and satisfies η|p|2 ≤ (Dp, p) ≤ (η + κ)|p|2 ∀ p ∈ R2 If η = 0, then the diffusion matrix is only positive semi-definite. Boundary conditions θ = Θ
Γd (2) q · n = g
Γn (3) Total energy
d dt
1 2θ2dx+
κ|b·∇θ|2dx+
η|∇θ|2dx+
(q(θ)·n)Θds+
gθds = 0
If Θ = 0, g = 0, we observe that the energy decreases with time.
7 / 57
Consider b = (1, 0) and initial condition θ(x, y, 0) = θ0(x, y) =
y < 0 θl y > 0 On the surface y = 0, n = (0, 1), heat flux q · n = −(b · n)(b · ∇θ) = 0 since b · n = 0. Hence θ(x, y, t) = θ0(x, y)
8 / 57
The initial temperature is given by θ(r, φ, 0) =
if 0.5 < r < 0.7 and
11 12π < φ < 13 12π
10,
The “magnetic field” is b = (− sin φ, cos φ).
9 / 57
7046 cells, 3624 vertices CG(1) using BDF2 and σ = 1 at t = 200
10 / 57
Heat flux q = (qx, qy) qx = −κBx(Bx∂xθ + By∂yθ), qy = −κBy(Bx∂xθ + By∂yθ) Finite volume method θn+1
i,j
− θn
i,j
∆t + (qx)i+ 1
2 ,j − (qx)i− 1 2 ,j
∆x + (qy)i,j+ 1
2 − (qy)i,j− 1 2
∆y = 0 Approximation of fluxes: centered asymmetric scheme [6] (qx)i+ 1
2 ,j = −κBx
Bx θi+1,j − θi,j ∆x + By ∂θ ∂y
2 ,j
∂θ ∂y
2 ,j
= 1 4∆y(θi,j+1 + θi+1,j+1 − θi,j−1 − θi+1,j−1)
11 / 57
i,j+1 i,j i+1,j
T =0.1
i+1,j+1
T =10 T =0.1 T =0.1 2 4 6 8 10 −0.5 0.5 1 1.5 2 2.5 3 time
temperature at (i,j)
with bx ¼ by ¼ 1= ffiffiffi 2 p . With the asymmetric method heat flows out of the third quadrant which is already a temperature minimum, resulting in a negative temperature Ti,j. However due to numerical perpendicular diffusion, at late times the temperature becomes positive
with slope limiters (dashed line; both give the same result), and symmetric with entropy limiting (dot dashed line).
Sharma & Hammett [6]
12 / 57
Th = triangulation of the domain by disjoint triangles, Ω = ∪K∈ThK. Ei = interior edges, Ed = Dirichlet edges, En = Neumann edges he = length of edge e hK = diameter of K ∈ Th = length of largest side of K Diameter of the mesh h = max
K∈Th hK
ρK = diameter of the largest circle that can be inscribed in triangle K The mesh is shape regular, i.e., there is a constant α > 0 such that hK ρK ≤ α hmin = length of the smallest edge in the triangulation
13 / 57
Space of broken polynomials V k
h = {ϕ ∈ L2(Ω) : ϕ|K ∈ Pk(K), ∀K ∈ Th}
and their vector version W k
h = V k h × V k h
For each interior edge { {q} } = 1 2(q− + q+) ϕ = ϕ+ − ϕ− ϕn = ϕ−n− + ϕ+n+ On a boundary edge { {q} } = q ϕ = ϕ ϕn = ϕn
n− n+ ϕ− ϕ+
14 / 57
Inverse inequality
h
1 2
e ϕe ≤ CtϕK,
∀ϕ ∈ Pk(K) where Ct =
Inverse Poincare-type inequality
hK∇ϕK ≤ CpϕK, ∀ϕ ∈ Pk(K) with Cp = √ 2α, where α is the shape regularity parameter.
15 / 57
n− K ϕ−
Multiply the heat equation (1) by a test function ϕh ∈ V k
h and integrate
ϕh∂tθhdx +
D∇θh · ∇ϕhdx −
ˆ D{ {∇θh} } · ϕ−
h n− = 0
where ˆ D = ηI + κ{ {b} }{ {b} }⊤
16 / 57
Sum over all cells K ∈ Th
ϕh∂tθhdx +
D∇θh · ∇ϕhdx −
ˆ D{ {∇θh} } · ϕhnds −
(D∇θh · n)ϕhds +
gϕhds = 0 However such a scheme is unstable !!! IPDG scheme: find θh ∈ V k
h such that for all ϕh ∈ V k h
(∂tθh, ϕh) + A(κ,η)
h
(θh, ϕh) = ℓΘ(ϕh) + ℓg(ϕh) (4)
17 / 57
where the bilinear form is given by A(κ,η)
h
(θh, ϕh) =
D∇θh · ∇ϕhdx −
D{ {∇θh} } · ϕhn ± ˆ D{ {∇ϕh} } · θhn
+
Cip he
ˆ Dθhn · ϕhnds (5) and the linear functionals due to the boundary conditions are given by ℓΘ(ϕh) = ∓
(D∇ϕh · n)Θhds +
Cip he
(Dn · n)Θϕhds (6) ℓg(ϕh) = −
gϕhds (7)
18 / 57
± signs in the bilinear form lead to the symmetric and non-symmetric interior penalty schemes, respectively. Cip is a dimension-less positive penalty parameter Penalty terms make the bilinear form associated with the scheme to be coercive, and hence stable. Interior penalty term on an edge
ˆ Dθhn · ϕhnds =
(η + κ|{ {b} } · n|2)θhϕhds The anisotropic diffusion is active across the edge only if { {b} } · n = 0 which is the physically correct behaviour. The penalty terms are similar to those used in Pestiaux et al. [4].
19 / 57
Define the mesh dependent norm ϕ2
h =
|∇ϕ|2dx +
Cip he ϕ2ds +
Cip he ϕ2ds (8) This is a norm only if the Dirichlet part of the boundary is non-empty,
20 / 57
Coercivity of A(κ,η)
h
(·, ·)
(1) Assume that η ≥ 0. The non-symmetric form is positive semi-definite A(κ,η)
h
(ϕ, ϕ) ≥ 0 ∀φ ∈ V k
h
(2) Assume that η > 0. The non-symmetric form is coercive A(κ,η)
h
(ϕ, ϕ) ≥ ηϕ2
h
∀φ ∈ V k
h
for any Cip > 0. (3) Assume that η > 0. The symmetric form is coercive A(κ,η)
h
(ϕ, ϕ) ≥ 1 2ηϕ2
h
∀φ ∈ V k
h
provided the penalty parameter is sufficiently large Cip ≥ 4C2
t n0
η + κ η 2 where n0 = 3 for triangles and n0 = 4 for quadrilaterals.
21 / 57
The above results can be proved by adapting the proof of the coercivity found in [5]. For the symmetric scheme to be coercive, the penalty constant Cip has to be chosen sufficiently large Cip ≥ 4C2
t n0
χ χ⊥ 2 to compensate for the lack of positive-definiteness of the anisotropic diffusion. Dirichlet boundary conditions are imposed in a weak manner through the fluxes and penalty terms, and are not built into the approximation space.
Conservation property
The DG methods are locally conservative. If we take the test function ϕh = 1 for (x, y) ∈ K and zero outside K, then we obtain
∂tθhdx + flux across ∂K = 0
22 / 57
Stability of semi-discrete IPDG scheme
Consider the non-symmetric scheme with η ≥ 0, Cip > 0 or the symmetric scheme with η > 0 and Cip sufficiently large. If Θ = 0, g = 0, the semi-discrete scheme (4) is stable in L2 norm. Proof: If we take ϕh = θh in (4), we obtain the energy equation d dt
1 2θ2
hdx + A(κ,η) h
(θh, θh) = 0 and under the conditions stated in the theorem A(κ,η)
h
(θh, θh) ≥ 0.
23 / 57
The backward Euler scheme is given by 1 ∆t
h
− θn
h, ϕh
h
(θn+1
h
, ϕh) = 0, ∀ϕh ∈ V k
h
(9) and the Crank-Nicholson scheme is given by 1 ∆t
h
− θn
h, ϕh
h
(θ
n+ 1
2
h
, ϕh) = 0, ∀ϕh ∈ V k
h
(10) where θ
n+ 1
2
h
= 1
2(θn h + θn+1 h
). We will say that the scheme is stable if θn+1
h
≤ θn
h
for the case of homogeneous boundary conditions.
24 / 57
Stability of BE/CN - IPDG
Consider the non-symmetric scheme with η ≥ 0, Cip > 0 or the symmetric scheme with η > 0 and Cip sufficiently large. If Θ = 0, g = 0, the backward Euler and Crank-Nicholson schemes are stable for any time step. Proof: Taking ϕh = θn+1
h
in (9), we get 1 ∆tθn+1
h
2 = 1 ∆t
h, θn+1 h
h
(θn+1
h
, θn+1
h
) ≤ 1 ∆t
h, θn+1 h
1 2∆tθn
h2 +
1 2∆tθn+1
h
2 which implies that θn+1
h
≤ θn
h and hence the backward Euler scheme
is stable. Taking ϕh = θ
n+ 1
2
h
in (10), we get 1 2∆tθn+1
h
2 − 1 2∆tθn
h2 = −A(κ,η) h
(θ
n+ 1
2
h
, θ
n+ 1
2
h
) ≤ 0 which shows that θn+1
h
≤ θn
h and hence the Crank-Nicholson scheme
is stable.
25 / 57
The simplest time integration scheme is the forward Euler scheme 1 ∆t
h
, ϕh
∆t (θn
h, ϕh) − A(κ,η) h
(θn
h, ϕh)
(11) which is explicit in nature. We show that this scheme is stable only under very restrictive time step conditions. To do this, we need the continuity of the bilinear form.
26 / 57
For both the symmetric and non-symmetric forms, we have A(κ,η)
h
(θ, ϕ) ≤ Cc(η + κ)θhϕh (12) for some constant Cc > 0 that depends only on the constant in the inverse
θh ≤ Cb hmin θ (13) for some constant Cb > 0.
27 / 57
Stability of IPDG with forward Euler
Consider the non-symmetric scheme with η > 0, Cip > 0 or the symmetric scheme with η > 0 and Cip sufficiently large. The forward Euler scheme (11) is stable under the time step restriction χ∆t h2
min
≤ α χ⊥ χ
28 / 57
This case corresponds to only anisotropic diffusion being present. In this case, the symmetric scheme is not even positive-semidefinite and hence we consider only the non-symmetric scheme. The forward Euler scheme is given by 1 ∆t
h
, ϕh
∆t (θn
h, ϕh) − A(κ,0) h
(θn
h, ϕh)
(14) Taking ϕh = θn+1
h
we get 1 ∆tθn+1
h
2 = 1 ∆t
h, θn+1 h
h
(θn
h, θn+1 h
) We can rewrite this as 1 2∆tθn+1
h
2+ 1 2∆tθn+1
h
− θn
h2 + A(κ,0) h
(θn+1
h
, θn+1
h
) = 1 2∆tθn
h2 + A(κ,0) h
(θn+1
h
− θn
h, θn+1 h
)
29 / 57
The last term on the right has to be controlled by the second and third terms on the left hand side. However this is not possible since A(κ,0)
h
(·, ·) is not coercive and we are unable to establish the stability of the forward Euler scheme in this case. We can study the stability of the scheme through numerical approach. For homogeneous boundary conditions, the semi-discrete scheme can be written as M dθh dt = Aθh All the eigenvalues of (A, M) must be non-positive. Define λM = max
i
|λi| Stability of the forward Euler scheme requires ∆t ≤ 2 λM
30 / 57
Compute the parameter σ = χ∆t h2
min
≤ 2χ λMh2
min
Mesh 1 (structured) Mesh 2 (unstructured)
31 / 57
We use two types of meshes consisting of 50 edges on each side of the square domain. For each φ ∈ [0, 2π], take a constant field b = (cos φ, sin φ) Compute largest eigenvalue λM numerically using SLEPc and the upper bound on σ
1 2 3 4 5 6 Angle φ 0.01 0.02 0.03 0.04 0.05 0.06 0.07 σ
Mesh 1 Mesh 2
Dependence of σ on the angle of b for NIPDG scheme and forward Euler time stepping minimum at φ = π/4, 5π/4 maximum at φ = 3π/4, 7π/4
32 / 57
This indicates the time step should approximately satisfy χ∆t h2
min
1 50 This limit may be specific to the grid and b field that we have assumed and we cannot claim this to be a universal upper bound.
33 / 57
Main idea: Write as system of first order equations, apply a DG scheme Matrix D is positive (semi)-definite and has eigenvalues λ1 = η, λ2 = η + κ|b|2 Its square root is well defined and given by D
1 2 =
1 |b|2 √λ1b2
2 + √λ2b2 1
(√λ2 − √λ1)b1b2 (√λ2 − √λ1)b1b2 √λ1b2
1 + √λ2b2 2
p = −D
1 2 ∇θ
(15) so that the heat equation can be written as ∂θ ∂t + ∇ · (D
1 2 p) = 0
(16)
34 / 57
Mutliplying by test functions ψh ∈ W k
h and ϕh ∈ V k h , respectively, and
performing an integration by parts on a single cell K, we obtain
ϕh∂tθhdx −
D
1 2 ph · ∇ϕhdx +
ϕ−
h ˆ
D
1 2 ˆ
ph · n−ds = 0
ph · ψhdx +
D
1 2 ∇θh · ψhdx +
( ˆ D
1 2 n− · ψ−
h )(ˆ
θh − θ−
h )ds = 0
where ˆ θ, ˆ p are numerical fluxes that have to specified. LDG scheme: Find (θh, ph) ∈ V k
h × W k h such that
(∂tθh, ϕh) + Ah(ph, ϕh) = ℓg(ϕh), ∀ ϕh ∈ V k
h
(17) (ph, ψh) + Bh(θh, ψh) = ℓΘ(ψh), ∀ ψh ∈ W k
h
(18)
35 / 57
where Ah(ph, ϕh) = −
D
1 2 ph · ∇ϕhdx +
ˆ D
1 2 ˆ
ph · ϕhnds Bh(θh, ψh) =
D
1 2 ∇θh · ψhdx +
θh ˆ D
1 2 ψhn − θh ˆ
D
1 2 ψhn
−
(D
1 2 n · ψh)θhds
ℓΘ(ψh) = −
(D
1 2 n · ψh)Θds 36 / 57
On the interior edges, the numerical fluxes are taken in an alternate manner (Cockburn & Shu [3]), e.g., ˆ p = p−, ˆ θ = θ+ (19)
ˆ p = { {p} }, ˆ θ = { {θ} } (20) while on the Dirichlet edges ˆ p = p, ˆ θ = Θ On the Neumann boundaries, we do not require the numerical flux. For other numerical fluxes, see [1], but the above two numerical fluxes were the only ones which allowed us to prove stability of the forward Euler scheme in the case η = 0 and leads to the usual time-step conditions.
37 / 57
Stability of semi-discrete LDG
If Θ = 0, g = 0, the semi-discrete LDG scheme is stable. Proof: Take ϕh = θh in (17) and ψh = ph in (18), and adding the two equations, we obtain the energy equation d dt
1 2θ2
hdx + (ph, ph) + Ah(ph, θh) + Bh(θh, ph) = ℓg(θh) + ℓΘ(ph)
and
Ah(ph, θh)+Bh(θh, ph) =
D
1 2 ˆ
ph · θhn + ˆ θh ˆ D
1 2 phn − θh ˆ
D
1 2 phn
Using the numerical fluxes (19) or (20) ˆ D
1 2 ˆ
ph · ϕhn + ˆ θh ˆ D
1 2 phn − θh ˆ
D
1 2 phn = 0 38 / 57
The energy equation becomes d dt
1 2θ2
hdx +
|ph|2dx +
(D
1 2 ph · n)Θds +
gθhds = 0 If Θ = 0, g = 0, then d dt
1 2θ2
hdx = −
|ph|2dx ≤ 0 which shows the stability of the semi-discrete scheme.
39 / 57
The backward Euler scheme is given by 1 ∆t
h
− θn
h, ϕh
h
, ϕh) = ℓg(ϕh), ∀ ϕh ∈ V k
h
(21)
h
, ψh
h
, ψh) = ℓΘ(ψh), ∀ ψh ∈ W k
h (22)
while the Crank-Nicholson scheme is given by 1 ∆t
h
− θn
h, ϕh
n+ 1
2
h
, ϕh) = ℓg(ϕh), ∀ ϕh ∈ V k
h
(23)
n+ 1
2
h
, ψh
n+ 1
2
h
, ψh) = ℓΘ(ψh), ∀ ψh ∈ W k
h (24)
Stability of implicit LDG
The backward Euler and Crank-Nicholson schemes are stable in L2 norm without any restriction on the time step.
40 / 57
1 ∆t
h
− θn
h, ϕh
h, ϕh)
= ℓg(ϕh), ∀ ϕh ∈ V k
h
(pn
h, ψh) + Bh(θn h, ψh)
= ℓΘ(ψh), ∀ ψh ∈ W k
h
The second equation can solved on each cell to obtain pn
h which is used in
the first equation to update θh.
Stability
If Θ = 0, g = 0, the forward Euler scheme is stable provided the time step satisfies χ∆t h2
min
≤ α for some α > 0.
41 / 57
Proof: Taking ϕh = θn+1
h
and ψh = pn
h and adding the two equations
1 ∆t
h
− θn
h, θn+1 h
h, θn+1 h
) + (pn
h, pn h) + Bh(θn h, pn h) = 0
This can be rewritten as 1 2∆tθn+1
h
2+ 1 2∆tθn+1
h
− θn
h2 + pn h2 + Ah(pn h, θn h) + Bh(θn h, pn h)
= 1 2∆tθn
h2 − Ah(pn h, θn+1 h
− θn
h)
(25)
42 / 57
We can establish a bound for Ah(p, ϕ) using the inverse inequality and inverse Poincare-type inequalities as follows. Ah(p, ϕ) ≤
D
1 2 pK∇ϕK +
ˆ D
1 2 peϕne
≤
K∈Th
D
1 2 p2
K
1 2
K∈Th
∇ϕ2
K
1 2
+
he ˆ D
1 2 p2
e
1
2
≤ (κ + η)
1 2 p · Cp
hmin ϕ +
2Ct(κ + η)
1 2 p ·
2 Ct hmin ϕ ≤ Cl(κ + η)
1 2
hmin pϕ, Cl = Cp + 3 2C2
t
43 / 57
Using this in (25), we get the inequality 1 2∆tθn+1
h
2 + 1 2∆tθn+1
h
− θn
h2 + pn h2
≤ 1 2∆tθn
h2 + Cl(κ + η)
hmin phθn+1
h
− θn
h
≤ 1 2∆tθn
h2 + ε
2ph2 + C2
l (κ + η)
2εh2
min
θn+1
h
− θn
h2
for some ε > 0, which can be rearranged as 1 2∆tθn+1
h
2+ 1 2∆t − C2
l (κ + η)
2εh2
min
h
− θn
h2 +
2
≤ 1 2∆tθn
h2
44 / 57
The largest time step is obtained by choosing ε = 2, and the scheme is stable in L2 norm if the time step satisfies the condition in the theorem with α = 2/C2
l .
Remark
The LDG scheme together with forward Euler time integration is thus stable under a time step restriction even in the case of χ⊥ = 0. Moreover the time step restriction does not depend on the ratio of χ to χ⊥ which was the case in the IPDG scheme, thus making LDG scheme preferable for explicit time integration.
45 / 57
46 / 57
47 / 57
48 / 57
49 / 57
50 / 57
Domain Ω = [− 1
2, 1 2] × [− 1 2, 1 2]
∂θ ∂t + ∇ · q = s(x, y), s(x, y) = 2π2 cos(πx) cos(πy) Both χ and χ⊥ are taken to be non-zero and b = (+ cos(πx) sin(πy), − sin(πx) cos(πy)) Steady solution θs(x, y) = χ−1
⊥ cos(πx) cos(πy)
is independent of χ since b · ∇θs = 0. Initial condition θ = 0
Ω × {t = 0} Boundary condition θ = 0
∂Ω × (0, ∞)
51 / 57
θiso
h
= solution with χ/χ⊥ = 1 (isotropic diffusion) θh = solution with χ/χ⊥ > 1 Numerical cross-diffusion [6] χ⊥,num =
θh(0, 0) − 1 θiso
h (0, 0)
10−1 h 10−6 10−5 10−4 10−3 10−2 χ⊥,num/χ⊥
χ||/χ⊥ = 10 χ||/χ⊥ = 100 52 / 57
◮ Explicit scheme has usual CFL condition ◮ Average flux can give exact resolution of discontinuity if mesh aligned
with b
◮ Better schemes to reduce χ⊥,num for χ ≫ 1 ? ◮ Temperature dependent coefficients χ, χ⊥ ? ◮ Positivity preserving scheme ? ◮ Coupling with full MHD simulation ? ◮ Need good linear solver in parallel ? 53 / 57
◮ Explicit scheme has usual CFL condition ◮ Average flux can give exact resolution of discontinuity if mesh aligned
with b
◮ Better schemes to reduce χ⊥,num for χ ≫ 1 ? ◮ Temperature dependent coefficients χ, χ⊥ ? ◮ Positivity preserving scheme ? ◮ Coupling with full MHD simulation ? ◮ Need good linear solver in parallel ?
54 / 57
[1] D. N. Arnold, F. Brezzi, B. Cockburn, and L. D. Marini, Unified analysis of discontinuous Galerkin methods for elliptic problems, SIAM Journal on Numerical Analysis, 39 (2002), pp. pp. 1749–1779. [2] F. Bassi and S. Rebay, A high-order accurate discontinuous finite element method for the numerical solution of the compressible Navier-Stokes equations, Journal of Computational Physics, 131 (1997), pp. 267 – 279. [3] B. Cockburn and C.-W. Shu, The local discontinuous Galerkin method for time-dependent convection-diffusion systems, SIAM J.
55 / 57
[4] A. Pestiaux, S. Melchior, J. Remacle, T. K¨ arn¨ a,
element discretization of a strongly anisotropic diffusion operator, International Journal for Numerical Methods in Fluids, 75 (2014),
[5] B. Rivi` ere, Discontinuous Galerkin Methods for Solving Elliptic and Parabolic Equations, Society for Industrial and Applied Mathematics, 2008. [6] P. Sharma and G. W. Hammett, Preserving monotonicity in anisotropic diffusion, J. Comput. Phys., 227 (2007), pp. 123–142.
56 / 57
[7] C. Sovinec, A. Glasser, T. Gianakon, D. Barnes,
and M. Chu, Nonlinear magnetohydrodynamics simulation using high-order finite elements, Journal of Computational Physics, 195 (2004), pp. 355 – 386.
57 / 57