Introduction Geometric properties Volume-preserving integrators Linear stability of the volume-preserving methods Numerical Exp
Volume-preserving algorithms for charged particle
dynamics
Yang He
University of Science and Technology of China
March. 2015, Hefei
Yang He
Volume-preserving Lorentz force equation
Introduction Geometric properties Volume-preserving integrators Linear stability of the volume-preserving methods Numerical Exp
Single particle model
Newton–lorentz equation
x˙ = v,
p˙ = q (E(x, t) + v × B(x, t))
p = mv in non-relativistic case,
p= p
m0 v
1 − v 2 /c2
in relativistic case.
Numerical algorithms
Due to the multi-scale nature of plasma dynamics, numerical
algorithms need to have long term accuracy.
Geometric integrators: to preserve
conservative volume in phase space: det
∂(x(t),p(t))
∂(x0 ,p0 )
constants of motion (energy, momentum, etc.),
etc.
Yang He
Volume-preserving Lorentz force equation
= 1,
Introduction Geometric properties Volume-preserving integrators Linear stability of the volume-preserving methods Numerical Exp
Splitting method
Yang He
Volume-preserving Lorentz force equation
Introduction Geometric properties Volume-preserving integrators Linear stability of the volume-preserving methods Numerical Exp
Two second order symmetric methods
F2
F1
1
2
3
RG2h := φF
◦ φF
◦ φF
h ◦ φh/2 ◦ φh/2 ;
h/2
h/2
xk+ 1 = xk +
2
h
vk ,
2
F2
F1
1
2
3
CG2h := φF
◦φF
◦ΦF
h ◦φh/2 ◦φh/2 ;
h/2
h/2
xk+ 1 = xk +
2
h
vk ,
2
hq
p− = pk +
E(xk+ 1 , tk ),
2
2
h
+
ˆ
,
t
)
p− ,
p = exp
B(x
1
1
k+ 2 k+ 2
m0 γ(p− )
hq
pk+1 = p+ +
E(xk+ 1 , tk+1 ),
2
2
hq
E(xk+ 1 , tk ),
2
2
h p+ + p−
+
ˆ
p =
× B(x
, tk+ 1 ),
k+ 1
2
2
2m0 γ(p− )
hq
pk+1 = p+ +
E(xk+ 1 , tk+1 ),
2
2
vk+1 = pk+1 /(m0 γ(pk+1 )),
vk+1 = pk+1 /(m0 γ(pk+1 )),
xk+1 = xk+ 1 +
2
h
vk+1 .
2
r
where γ(p) =
1+
p2
2
m2
0c
p− = pk +
xk+1 = xk+ 1 +
2
h
vk+1 .
2
in relativistic case, and γ(p) is replaced by 1 in
non-relativistic case.
Yang He
Volume-preserving Lorentz force equation
Introduction Geometric properties Volume-preserving integrators Linear stability of the volume-preserving methods Numerical Exp
Contents of this talk
Analyze the performances of the two volume-preserving
methods, in
linear stability,
accuracy,
cost,
preservation of the invariants.
Optimizing the methods by designing higher order
volume-preserving methods.
Yang He
Volume-preserving Lorentz force equation
Introduction Geometric properties Volume-preserving integrators Linear stability of the volume-preserving methods Numerical Exp
Linear stability: how to choose the proper step size
Normalized relativistic Lorentz force equation
x˙ =
p
,
γ(p)
p˙ = E(x) +
where v = p/γ(p),
p
× B(x),
γ(p)
γ(p) =
p
1 + p2
The test linearized equation
x˙ =
γ0 =
√
1
p,
γ0
p˙ = Λx +
ω
p×b
γ0
1 + s2 > 1, where 0 ≤ s ≤ p(t).
Λx = −(λ2x xex + λ2y yey + λ2z zez ), λ > 0. b = [0, 0, 1].
λ2x corresponds to kx E(x0 ), ω corresponds to B(x0 ).
Yang He
Volume-preserving Lorentz force equation
Introduction Geometric properties Volume-preserving integrators Linear stability of the volume-preserving methods Numerical Exp
Linear stability
The test equation
x˙ =
1
p,
γ0
p˙ = −λ2 x +
where = ±1, x = [x, y], v = [vx , vy ], J =
0 1
−1 0
ω
Jp
γ0
, λ2 ∼ kE(x0 ), ω ∼ B(x0 ).
The numerical solution
γ0 x
hp
!n
2
= [Mm ((hλ) , hω)]
n
γ0 x(0)
!
hp(0)
Proposition
For the n-dimensional source-free system, a volume-preserving
method is linear stable iff
ρ(Mm ) = 1 and Mm is diagonalizable.
Yang He
Volume-preserving Lorentz force equation
Introduction Geometric properties Volume-preserving integrators Linear stability of the volume-preserving methods Numerical Exp
Linear stability domain
(a) RG2h
(b) CG2h
0)
In relativistic case: λ2 ∼ k E(x
,ω∼
B0 c
B(x0 )
,
B0
γ0 =
√
1 + s2 , 0 ≤ s ≤
p(t)
.
m0 c
0)
0)
In non-relativistic case: λ2 ∼ k E(x
, ω ∼ B(x
, γ0 = 1.
B0 v0
B0
√
The two schemes are stable for hλ/ γ0 < 2.
To simulate accurately the Larmor cyclotron it is needed that hω/γ0 ≤ π.
Yang He
Volume-preserving Lorentz force equation
Introduction Geometric properties Volume-preserving integrators Linear stability of the volume-preserving methods Numerical Exp
Other properties of the two volume-preserving methods
Accuracy
They are of second order, i.e. |x1 − x(0)| ≤ Ch2 .
Cost
They are explicit, only needs 1 evaluation of E(x) and B(x)
per step.
Preservation of invariants
They can preserve the kinetic energy of the particle exactly
with static magnetic field;
They can bound the error of invariants up to O(h2 ).
Yang He
Volume-preserving Lorentz force equation
Introduction Geometric properties Volume-preserving integrators Linear stability of the volume-preserving methods Numerical Exp
Higher order methods
Two fourth order method designed by composition
3-stage (Yoshida, 90): G4h = Φa1 h ◦ Φa2 h ◦ Φa1 h
where a1 = (2 − 21/3 )−1 , a2 = 1 − 2a1 ;
5-stage (Suzuki, 92): G4h = Φb1 h ◦ Φb2 h ◦ Φb3 h ◦ Φb2 h ◦ Φb1 h
where b1 = b2 = (4 − 41/3 )−1 , b3 = 1 − 2(b1 + b2 );
A fourth order method designed by processing(Blanes, 99)
G4h = χh ◦ Φh ◦ χ−1
h
1
23
1
23
1
The kernel: Φh = φ1a1 h ◦ φ23
b1 h ◦ φa2 h ◦ φb2 h ◦ φa2 h ◦ φb1 h ◦ φa1 h
The processor: χh =
5
Y
φ1zi h ◦ φ23
yi h
i=1
a1 = 0.6762,
Yang He
b1 = −0.175,
a2 = −0.1762,
Volume-preserving Lorentz force equation
b2 = 1.35
Introduction Geometric properties Volume-preserving integrators Linear stability of the volume-preserving methods Numerical Exp
Linear stability domain
S
U
S
S
(c) BG2h
S
U
U
U
(d) CG2h
S
U
S
U
(e) by Yoshida composi- (f) by Suzuki composition (g) by processing(only for
tion of BG2h
of BG2h
non-relativistic case)
Yang He
Volume-preserving Lorentz force equation
Introduction Geometric properties Volume-preserving integrators Linear stability of the volume-preserving methods Numerical Exp
Accuracy and cost
Truncation error
Error coefficients
Evaluations of E(x)
after N steps
Effective error constants
G4h Y
C1 h4
3.92e-2
3N
G4h S
C 2 h4
6.86e-4
5N
G4h P
C3 h4
1.3e-3
3N + 10
1.3349
0.8092
0.5696
Table: G4h Y denotes the Yoshida’s composition method and G4h S
denotes the Suzuki’s composition method. G4h P denotes the processing
method for the non-relativistic motion equation.
Yang He
Volume-preserving Lorentz force equation
Introduction Geometric properties Volume-preserving integrators Linear stability of the volume-preserving methods Numerical Exp
Numerical experiment 1: non-relativistic case
Static electromagnetic field
B = B0 ez ,
E = El (
x
y
ex +
ey ),
R0
R0
where B0 = 3T , E0 = 3V /m, R0 = 1m.
−4
1.5
10
1
10
−5
−6
10
∆ H/H0
y*ωce
0.5
0
−0.5
−7
10
−8
10
−9
10
−10
−1
10
−1.5
−1.5
10
2nd order
4th order
−11
−1
−0.5
0
x*ωce
0.5
1
1.5
0
1000 2000 3000 4000 5000 6000 7000 8000
Steps
Figure: Take the initial momentum as vk (0) = 1v0 , v⊥ (0) = 5v0 ,
v0 = 0.1m/s. Choose the step size as h = 0.1π, run for 104 steps.
Yang He
Volume-preserving Lorentz force equation
Introduction Geometric properties Volume-preserving integrators Linear stability of the volume-preserving methods Numerical Exp
Ex1: relative errors of the solution
0
−1
10
10
0.067h2
−2
log(phase error) at the time 39.17
relative error of (x,v)
−2
10
−4
10
G2h I
−6
10
Suzuki
processing
Yoshida
−8
10
4
6.8e−4h
10
−3
10
−4
10
G2hI
−5
10
Suzuki of G2hI
−6
10
processing of G2hI
−7
Yoshida of G2hI
10
−10
10
1
10
0
−1
10
10
−2
10
time step h
−8
10
1
10
2
3
10
10
4
10
log (effort)
(a) as functions of h
(b) as functions of efforts (evaluations of E(x))
Yang He
Volume-preserving Lorentz force equation
Introduction Geometric properties Volume-preserving integrators Linear stability of the volume-preserving methods Numerical Exp
Numerical experiment 1: relativistic case
Static electromagnetic field
B = B0 ez ,
E = El (
x
y
ex +
ey ),
R0
R0
where B0 = 3T , E0 = 3V /m, R0 = 1m.
−3
2
x 10
−15
10
t=0s
1.5
−16
10
1
−17
∆ H/H
y
0
0.5
0
−0.5
−1
t=6e−8 s
10
−18
10
2nd order
4th order Yoshida
4th order Suzuki
−19
10
−1.5
−20
−2
1.05 1.0505 1.051 1.0515 1.052 1.0525 1.053 1.0535 1.054
x
10
0
1
2
3
Steps
4
5
6
Figure: Take the initial momentum as pk0 = 1m0 c, p⊥0 = 5m0 c.
Choose the step size as h = 0.1π, run for 105 steps.
Yang He
4
x 10
Volume-preserving Lorentz force equation
Introduction Geometric properties Volume-preserving integrators Linear stability of the volume-preserving methods Numerical Exp
Ex1: relative errors of the solution
2
2
10
10
CG2h
0
10
0
10
4th order Yoshida
4th order Suzuki
−2
relative error of x/m
relative error of x/m
10
−4
10
−6
10
−8
10
CG2h
−6
10
10
−10
−12
10
−4
10
−8
4th order Yoshida
4th order Suzuki
−10
10
−2
10
1
10
0
10
−1
10
time step h
−2
10
−3
10
(a) as functions of h
10
1
10
2
10
3
10
log(effort)
4
10
5
10
(b) as functions of efforts (evaluations of φ3h )
Yang He
Volume-preserving Lorentz force equation
Introduction Geometric properties Volume-preserving integrators Linear stability of the volume-preserving methods Numerical Exp
Numerical experiment 2: non-relativistic case
Axi-symmetric tokamak geometry
B0 R0
B0 r
B=∇×A=
eξ +
eθ ,
R
qR
where q = 2 and R0 = 1m.
ϕ(x) = 0,
−14
0.08
0.06
relative error
0.04
z
0.02
0
−0.02
−0.04
−0.06
−0.08
1
1.02
1.04
1.06
1.08
1.1
R
(c) Banana Orbit
x 10
8
7
6
5
4
3
2
1
0
0
2
volume−preserving methods
4
6
8
10
time
12
14
16
4
x 10
(d) Relative error of energy
Figure: Take the initial momentum as vk0 = 1m0 v0 , v⊥0 = 5m0 v0 .
Choose the step size as h = 0.1π, run for 2 × 106 steps.
Yang He
Volume-preserving Lorentz force equation
Introduction Geometric properties Volume-preserving integrators Linear stability of the volume-preserving methods Numerical Exp
Numerical experiment 2: relativistic case
Axi-symmetric tokamak geometry
B=∇×A=
B0 R0
B0 r
eξ +
eθ ,
R
qR
ϕ(x) = ϕ0
R0
,
R
where q = 2 and R0 = 1m.
−16
0.08
20
x 10
2nd order
4th order
0.06
15
0.04
∆ H/H0
z
0.02
0
−0.02
−0.04
10
5
0
−0.06
−0.08
0.92 0.94 0.96 0.98
1
1.02 1.04 1.06
−5
0
R
2
4
6
Steps
8
10
12
Figure: Take the initial momentum as pk0 = 5m0 c, p⊥0 = 1m0 c.
Choose the step size as h = 0.1π, run for 105 steps.
Yang He
4
x 10
Volume-preserving Lorentz force equation
Introduction Geometric properties Volume-preserving integrators Linear stability of the volume-preserving methods Numerical Exp
Summary
We analyze the performances of two volume-preserving
algorithms when simulating charged particle trajectories in
both relativistic and non relativistic cases.
Higher order methods are constructed by composition and by
processing.
Numerical analysis and experiments show that the two higher
order methods are more efficient than the second order
methods for more strict error restriction.
Ongoing work
To construct more efficient higher order methods for
simulating relativistic dynamics of charged particles.
To analyze the linear stability of the higher order methods for
relativistic case.
Yang He
Volume-preserving Lorentz force equation