Discrete Involutions, Resonance, and the Divergence Problem in MHD

Divergence Problem in MHD

Friedemann Kemm BTU Cottbus

kemm@math.tu-cottbus.de

0 0.05 0.1 0.15 0.2 0.25

0 20 40 60 80 100

central differences upwind differences transverse upwind differences

0 0.05 0.1 0.15 0.2 0.25

0 20 40 60 80 100

central differences upwind differences transverse upwind differences

The problem

Equations for the magnetic field:

B_t−∇ × (v × B) = 0

∇ · B = 0

⇒ No coupling of evolution and constraint!

The known solutions

Constrained Transport

• Apply special discretization

• Evans, Hawley, Balsara, Spicer, Torrilhon, Fey, Rossmanith, Kr¨oner, Besse . . .

• Discretized ∇ · B constant in time Divergence Cleaning

• Extend equations to control errors

• GLM-correction, Powell method, Transport correction

• With suitable chosen coefficients ∇ · B → 0 (t → ∞)

The surprise

Some schemes do not need any special technique:

• Zachary, Malagoli and Colella (1994), Upwind scheme

• Balbas and Tadmor (2006), Central scheme

Common to both schemes:

• No 1-d physics for inter-cell fluxes

Involutions

Conservation Law (F = (F₁, F₂, . . .))

q_t + ∇ · F(q) = 0

with Involution (M_i matrices)

X

i

M_i q_xi = 0

satisfied for every time whenever satisfied by initial condition Sufficient condition:

M_iF_j + M_jF_i = 0 i , j = 1,2, . . .

⇒ P

i M_iq_xi constant in time

Proof of Sufficiency

• Apply P

i M_{i ∂x}^∂

i to the conservation law

• Constant matrices commute with derivatives

• Partial derivatives commute with each other

• Due to condition all terms including fluxes vanish Result:

∂

∂t

X

i

M_i q_xi

= 0 Discrete Analogue?

Linear Schemes – Example

Discretized Conservation Law (2-d):

∂ˆ

∂tˆ q +

∂ˆ

∂xˆ F(q) +

∂ˆ

∂yˆ G(q) = 0 with constant coefficient vectors α, β, γ and

∂ˆ

∂tˆ h_{i ,j}ⁿ := X

m

α_m h^n+m_{i ,j}

∂ˆ

∂xˆ h_{i ,j}ⁿ := X

k ,l

β_{k ,l} h_iⁿ_{+k ,j+l}

∂ˆ

∂yˆ h_{i ,j}ⁿ := X

k ,l

γ_{k ,l} h_iⁿ_{+k ,j+l}

Linear Schemes – Simple Calculations

∂ˆ²

∂xˆ ∂tˆ h_{i ,j}ⁿ = X

k ,l

β_{k ,l} X

m

α_m h_i^n+m_{+k ,j+l} = X

k ,l ,m

β_{k ,l} α_m h_i^n+m_{+k ,j+l}

= X

m

α_m X

k ,l

β_{k ,l} h_i^n+m_{+k ,j+l} =

∂ˆ²

∂tˆ ∂xˆ h_{i ,j}ⁿ

M

∂ˆ

∂tˆ hⁿ_{i ,j} = X

m

α_m M h^n+m_{i ,j} = X

m

α_m (M h)^n+m_{i ,j} =

∂ˆ

∂tˆ (M h)ⁿ_{i ,j}

Other (pairs of) partial derivatives analogue

Discrete Proof of Sufficiency

• Apply P

i M_{i ˆ}^∂ˆ

∂x_i to the discrete conservation law

• Constant matrices commute with discrete derivatives

• Linear discrete partial derivatives commute with each other

• Due to condition all terms including fluxes vanish Result:

∂ˆ

∂tˆ

X

i

M_i

∂ˆ

∂xˆ _iq

= 0

Numerical example

Linearized induction equation (2d)

B_t − ∇ × (v × B) = 0 , v = (u, v)^T ≡ constant , u, v > 0

• Standard upwind (DCU)

• Transverse upwind (CTU):

∂ˆ

∂xˆ h = (1 − c_y)h_{i ,j} − h_i−1,j

∆x + c_yh_{i ,j−1} − h_i−1,j−1

∆x

∂ˆ

∂yˆ h = (1 − c_x)h_{i ,j} − h_{i ,j−1}

∆y + c_xh_i−1,j − h_i−1,j−1

∆y

• Forward Euler in time

Numerical Results

Initial condition

B₁ = cos(2πx + πy) B₂ = −2 cos(2πx + πy) L²-norm

0 0.05 0.1 0.15 0.2 0.25

0 20 40 60 80 100

central differences upwind differences transverse upwind differences

0 0.05 0.1 0.15 0.2 0.25

0 20 40 60 80 100

central differences upwind differences transverse upwind differences

Standard upwind Transverse upwind

Nonlinear Case

Finite differences for partial derivatives

∂ˆ

∂xˆ h_j = (h_x)_j + O(∆x^p)

∂ˆ

∂tˆ h_j = 1

∆t

X

i∈I

α_ih_j+i = (h_t)_j + O(∆t^q) with bounded, not necessarily constant, in general matrix valued α_j

∂ˆ

∂tˆ

∂ˆ

∂xˆ h

j = 1

∆t

X

i∈I

α_i

(h_x)_j+i + O(∆x^p)

=

∂ˆ

∂tˆ (h_x)_j + O ∆x^p

∆t

= (h_{x t})_j + O(∆t^q) + O ∆x^p

∆t

Commutators of finite differences

Space and time derivatives

∂ˆ

∂tˆ

∂ˆ

∂xˆ h

j − ∂ˆ

∂xˆ

∂ˆ

∂tˆ h

j = O ∆x^p

∆t

+ O ∆t^q

∆x

(1)

Space derivatives with each other

∂ˆ

∂yˆ

∂ˆ

∂xˆ h

j − ∂ˆ

∂xˆ

∂ˆ

∂yˆ h

j = O(∆x^p−1) . (2)

Partial derivatives with matrices M

∂ˆ

∂tˆ (h)_j − ∂ˆ

∂tˆ (Mh)_j = O(∆t^q) (3)

Discrete Nonlinear Proof of Sufficiency

• Use discrete differences of order O(∆x^p) in space and O(∆t^q) in time

• Apply P

i M_{i ˆ}^∂ˆ

∂x_i to the discrete conservation law

• Constant matrices commute with discrete derivatives with error of order O(∆t^q) (time) or O(∆x^p) (space)

• nonlinear discrete partial space derivatives commute with error according to (1) and (2)

• Due to condition all terms including fluxes vanish Result:

∂ˆ

∂tˆ

X

l

M_l

∂ˆ

∂xˆ_lq_j

= O ∆x^p

∆t

+ O ∆t^q

∆x

+ O(∆x^p−1)

Resonance

Weakly hyperbolic degeneracy for u = 0

B_1t + vB_1y = 0 (4)

B_2t − vB_1x = 0 (5)

Eq. (5) ODE for B₂ ⇒ grows (for B_1y = 0 linear) in time

∇ · B = B_1x + B_2y ⇒ B_1x = −B_2y + ∇ · B Inserted in (5):

B_2t + v B_2y = v(∇ · B)

• Advection in y-direction if divergence constrained fulfilled

• Otherwise, source, constant in time ⇒ Resonance

Resonance: Effects in Full MHD

• B grows by resonance

• Wave speeds depend on B

• Time step depends on leading wave speed

• Error estimate for the growth of ˆ∇ · B depends on ∆x^p/∆t

Result: Error estimate for divergence worthless

Resonance – Numerical Example

0 5 10 15 20 25

-0.1 -0.08 -0.06 -0.04 -0.02 0 0.02 0.04 0.06 0.08 0.1

B1 B2

0 5 10 15 20 25

-0.1 -0.08 -0.06 -0.04 -0.02 0 0.02 0.04 0.06 0.08 0.1

B1 B2

0 5 10 15 20 25

-0.1 -0.08 -0.06 -0.04 -0.02 0 0.02 0.04 0.06 0.08 0.1

B1 B2

Initial condition CIR-scheme Lax-Friedrichs

Consequences:

• Schemes resolving resonant waves fail

• Example System: all schemes taking care of wave-speeds coincide

Lax-Friedrichs Scheme

Starting point

qⁿ⁺¹_i − qⁿ_i

∆t + f (qⁿ_i+1) − f (qⁿ_i−1)

2∆x = 0

Unstable ⇒ Replace qⁿ_i by average of neighboring values LF-scheme

qⁿ⁺¹_i − ¹₂(qⁿ_i+1 + qⁿ_i−1)

∆t + f (qⁿ_i+1) − f (qⁿ_i−1)

2∆x = 0

Consequence for _ˆ^∂ˆ

∂th = 0

h_iⁿ⁺¹ − 1

2(h_iⁿ₊₁ + hⁿ_i−1) = 0

⇔ h_iⁿ⁺¹ = h_iⁿ + 1

2(h_iⁿ₊₁ − 2hⁿ_i + h_iⁿ₋₁)

= h_iⁿ + ∆t ∆x² 2∆t

h_iⁿ₊₁ − 2h_iⁿ + h_iⁿ₋₁

∆x²

HLL-scheme

Numerical flux function

g_HLL(q_r, q_l) = 1

2 f (q_r) + f (q_l)

− 1 2

S_R + S_L

S_R − S_L f (q_r) − f (q_l) + S_RS_L

S_R − S_L(q_r − q_l)

• Red term goes into time difference ⇒ central viscosity

• Vanishes if resonant wave exactly resolved (S_L = 0 or S_R = 0)

⇒ Impose lower bound on |S_L/R| ⇒ lower bound on central viscosity Problem with 1d physics for inter-cell fluxes: No control on resonant wave

Shallow Water MHD

Evolution equations

h_t + ∇ · [hv] = 0 (hv)_t + ∇ · [hvv^T−hBB^T + gh²

2 I] = 0 (hB)_t − ∇ × [v × (hB)] = 0

and the divergence constraint

∇ · (hB) = 0

Resonance e. g. for B k v

De Sterck Problem – Divergence for 1st Order

L^∞-norm of divergence

0 0.5 1 1.5 2 2.5 3

0 50 100 150 200 250 300

with Harten fix without entropy fix

Resonant Example

LLF (Rusanov) 2nd order, 1-d physics for inter-cell fluxes

-1 -0.5

0 0.5

1 -1 -0.5

0 0.5

1 0.5

1 1.5 2 2.5

-1 -0.5

0 0.5

1 -1 -0.5

0 0.5

1 0

0.5 1 1.5 2 2.5

-1 -0.5

0 0.5

1 -1 -0.5

0 0.5

1 0

20 40 60 80 100

B₁ after 4,5, and 6 time steps

De Sterck Problem (2nd order Roe with Harten entropy fix)

-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1

-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1

Initial state t=0.8

-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1

-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1

Distribution of Divergence

-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1

-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1

t=0.8 t=4.8

• No Divergence errors transported

• Computation survives even for long time

Outlook

0.001 0.01 0.1 1 10 100

100

Divergence error

Number of grid cells per direction Lax-Wendroff, dimensional splitting

CTU 2nd order

General criterion for ^∂_ˆ^ˆ

∂t

P

l M_l ^∂_ˆ^ˆ

∂x_lq_j

to be of the same order as the scheme itself?

Conclusions

• Divergence constraint prevents resonance

• Applying linear difference operators preserves discrete involution

• Approximate discrete involution for nonlinear systems

• Resonance can make the approximation worthless

• Stability by central viscosity on resonant wave

• 1-d physics for fluxes prevents control of viscosity on resonant wave

Problems with Initial Condition

B_l B^∗_l

B^∗_l

B_r B_r B_r

B^∗_r B^∗_r

B^∗_r

Divergence in Upwind Differences:

(Marked cell)

• Vanishes if u · v > 0

• Else

∇ · B = α

∆x(B_2r − B_2l)

∆x space step size

α fraction of marked cell right of discontinuity

Numerical Example

B₁ with different velocities

-2.5 -2 -1.5 -1 -0.5 0 0.5 1

0 50

100 150

200 250 300 0

10 20 30 40 50 60 70 80 90 100 -2.5

-2 -1.5 -1 -0.5 0 0.5 1

-2.5 -2 -1.5 -1 -0.5 0 0.5 1

0 50

100 150

200 250 300 0

10 20 30 40 50 60 70 80 90 100 -2.5

-2 -1.5 -1 -0.5 0 0.5 1

u, v > 0 u > 0, v < 0

Important Types of Involutions

Type I:

∂

∂t

X

i

M_i q_xi

= 0 Type II:

tlim→∞

X

i

M_i q_xi = 0 independent of initial condition

Examples:

MHD: ∇ · B Involution of Type I SMHD: ∇ · (hB) Involution of Type I GLM-MHD: ∇ · B Involution of Type II Powell-MHD: ∇ · B Involution of Type II

Resonance – Second Example (3d)

P_t + ∇(v · P) = 0 ∇ × P =: C Degeneracy for u = v = 0

P_t + w∇P₃ = 0 Insert Curl of P

P_t + wP_z = w





C₂

−C₁ 0





• Resonance if ∇ × P 6= 0 and not parallel to v

• Otherwise linear Advection in z-direction

Relation to GLM-approach

System including numerical viscosity introduced by LF

B_t − ∇ × (v × B) = ∆x²

2∆t∇ · (∇B) Parabolic GLM

B_t − ∇ × (v × B) + ∇ψ = 0 1

κψ + ∇ · B = 0 results in

B_t − ∇ × (v × B) = κ∇(∇ · B)