Time-parallel solution of the eddy current problem Iryna - - PowerPoint PPT Presentation

Time-parallel solution of the eddy current problem

Iryna Kulchytska-Ruchka1,2, Sebastian Schöps1,2, Herbert De Gersem1,2

1Graduate School of Computational Engineering, 2Institut für Theorie Elektromagnetischer Felder, Technische Universität Darmstadt

4th STEAM Collaboration Meeting, Darmstadt

www.graduate-school-ce.de September 21, 2017

Outline Introduction

Motivation The eddy current problem

Parallel-in-time solution

The Parareal method for IVPs Numerical example: coaxial cable model

Fourier basis for time-periodic systems

Coarse solution by spectral collocation Numerical results: coaxial cable model Systems with nonsmooth excitations

Conclusions and outlook

TU Darmstadt | GSC CE | Iryna Kulchytska | Time-parallel solution of the eddy current problem | 2/33

Overview Introduction

Motivation The eddy current problem

Parallel-in-time solution

The Parareal method for IVPs Numerical example: coaxial cable model

Fourier basis for time-periodic systems

Coarse solution by spectral collocation Numerical results: coaxial cable model Systems with nonsmooth excitations

Conclusions and outlook

TU Darmstadt | GSC CE | Iryna Kulchytska | Time-parallel solution of the eddy current problem | 3/33

Motivation → Transient FEM simulation

Fig.: Cross-section of an induction machine (J. Gyselinck).

TU Darmstadt | GSC CE | Iryna Kulchytska | Time-parallel solution of the eddy current problem | 4/33

Motivation → Transient FEM simulation

Fig.: Cross-section of an induction machine (J. Gyselinck).

→ System evolution in time

TU Darmstadt | GSC CE | Iryna Kulchytska | Time-parallel solution of the eddy current problem | 4/33

Motivation → Transient FEM simulation

Fig.: Cross-section of an induction machine (J. Gyselinck).

→ System evolution in time

TU Darmstadt | GSC CE | Iryna Kulchytska | Time-parallel solution of the eddy current problem | 4/33

Motivation → Transient FEM simulation

Fig.: Cross-section of an induction machine (J. Gyselinck).

→ System evolution in time

• Long settling time till the steady

state

TU Darmstadt | GSC CE | Iryna Kulchytska | Time-parallel solution of the eddy current problem | 4/33

Motivation → Transient FEM simulation

Fig.: Cross-section of an induction machine (J. Gyselinck).

→ System evolution in time

• Long settling time till the steady

state

• Many time steps =

⇒ time-consuming computation!

TU Darmstadt | GSC CE | Iryna Kulchytska | Time-parallel solution of the eddy current problem | 4/33

Overview Introduction

Motivation The eddy current problem

Parallel-in-time solution

The Parareal method for IVPs Numerical example: coaxial cable model

Fourier basis for time-periodic systems

Coarse solution by spectral collocation Numerical results: coaxial cable model Systems with nonsmooth excitations

Conclusions and outlook

TU Darmstadt | GSC CE | Iryna Kulchytska | Time-parallel solution of the eddy current problem | 5/33

The eddy current problem Fundamentals of electromagnetism: Maxwell’s equations. Assumptions:

Quasi-static regime:

|J| ≫

• ∂ D

∂t

• ;

Neglect hysteresis.

The eddy current equation: σ∂A ∂t + curl(ν(| curl A|) curl A) = Js, A − unknown magnetic vector potential; Js − impressed current density; σ, ν − electric conductivity and magnetic reluctivity.

TU Darmstadt | GSC CE | Iryna Kulchytska | Time-parallel solution of the eddy current problem | 6/33

Semi-discrete problem Solve the IVP in time: Mdtu(t) = f(t, u), t ∈ I := (0, T), u(0) = u0, where u : I → RNdof denotes the space-discretization of A .

TU Darmstadt | GSC CE | Iryna Kulchytska | Time-parallel solution of the eddy current problem | 7/33

Semi-discrete problem Solve the IVP in time: Mdtu(t) = f(t, u), t ∈ I := (0, T), u(0) = u0, where u : I → RNdof denotes the space-discretization of A .

• f(t, u) = −Ku(t) + g(t);
• M, K ∈ RNdof×Ndof − mass and stiffness matrices;
• g(t) − excitation (e.g., impressed current).

TU Darmstadt | GSC CE | Iryna Kulchytska | Time-parallel solution of the eddy current problem | 7/33

Overview Introduction

Motivation The eddy current problem

Parallel-in-time solution

The Parareal method for IVPs Numerical example: coaxial cable model

Fourier basis for time-periodic systems

Coarse solution by spectral collocation Numerical results: coaxial cable model Systems with nonsmooth excitations

Conclusions and outlook

TU Darmstadt | GSC CE | Iryna Kulchytska | Time-parallel solution of the eddy current problem | 8/33

The Parareal method: splitting of the time interval Partitioning the time interval into N windows (e.g., one per core) yields            Mdtu0 = f(t, u0), u0(T0) = U0, t ∈ (T0, T1], Mdtu1 = f(t, u1), u1(T1) = U1, t ∈ (T1, T2], . . . MdtuN−1 = f(t, uN−1), uN−1(TN−1) = UN−1, t ∈ (TN−1, TN],

TU Darmstadt | GSC CE | Iryna Kulchytska | Time-parallel solution of the eddy current problem | 9/33

The Parareal method: splitting of the time interval Partitioning the time interval into N windows (e.g., one per core) yields            Mdtu0 = f(t, u0), u0(T0) = U0, t ∈ (T0, T1], Mdtu1 = f(t, u1), u1(T1) = U1, t ∈ (T1, T2], . . . MdtuN−1 = f(t, uN−1), uN−1(TN−1) = UN−1, t ∈ (TN−1, TN], T0 T1 T2 T3 T4 T5 u(t) t

TU Darmstadt | GSC CE | Iryna Kulchytska | Time-parallel solution of the eddy current problem | 9/33

The Parareal method: splitting of the time interval Partitioning the time interval into N windows (e.g., one per core) yields            Mdtu0 = f(t, u0), u0(T0) = U0, t ∈ (T0, T1], Mdtu1 = f(t, u1), u1(T1) = U1, t ∈ (T1, T2], . . . MdtuN−1 = f(t, uN−1), uN−1(TN−1) = UN−1, t ∈ (TN−1, TN], T0 T1 T2 T3 T4 T5 u(t) t

TU Darmstadt | GSC CE | Iryna Kulchytska | Time-parallel solution of the eddy current problem | 9/33

Parareal as the Newton-Raphson method (I) Let F(t, Ti, U) be the solution operator of the IVP on (Ti, Ti+1], for i = 0, . . . , N − 1, which propagates the initial value U in time.

TU Darmstadt | GSC CE | Iryna Kulchytska | Time-parallel solution of the eddy current problem | 10/33

Parareal as the Newton-Raphson method (I) Let F(t, Ti, U) be the solution operator of the IVP on (Ti, Ti+1], for i = 0, . . . , N − 1, which propagates the initial value U in time. Matching conditions can be satisfied by root finding        U1 − F(T1, T0, U0) = 0, . . . UN−1 − F(TN−1, TN−2, UN−2) = 0

TU Darmstadt | GSC CE | Iryna Kulchytska | Time-parallel solution of the eddy current problem | 10/33

Parareal as the Newton-Raphson method (I) Let F(t, Ti, U) be the solution operator of the IVP on (Ti, Ti+1], for i = 0, . . . , N − 1, which propagates the initial value U in time. Matching conditions can be satisfied by root finding        U1 − F(T1, T0, U0) = 0, . . . UN−1 − F(TN−1, TN−2, UN−2) = 0

• r, equivalently,

F(U) = 0, with UT =

• UT

0 , UT 1 , ..., UT i , ..., UT N−1

• .
• M. J. Gander and E. Hairer, Nonlinear convergence analysis for the parareal algorithm, Domain Decomposition

Methods in Science and Engineering XVII, Springer Berlin Heidelberg, 2008. TU Darmstadt | GSC CE | Iryna Kulchytska | Time-parallel solution of the eddy current problem | 10/33

Parareal as the Newton-Raphson method (II) The Newton-Raphson iteration: for k = 0, 1, . . . U(k+1) = u0, U(k+1)

n

= F

• Tn, Tn−1, U(k)

n−1

• + ∂F(Tn, Tn−1, U)

∂U

• U(k+1)

n−1

− U(k)

n−1

• ,

where n = 1, . . . , N − 1.

TU Darmstadt | GSC CE | Iryna Kulchytska | Time-parallel solution of the eddy current problem | 11/33

Parareal as the Newton-Raphson method (II) The Newton-Raphson iteration: for k = 0, 1, . . . U(k+1) = u0, U(k+1)

n

= F

• Tn, Tn−1, U(k)

n−1

• + ∂F(Tn, Tn−1, U)

∂U

• U(k+1)

n−1

− U(k)

n−1

• ,

where n = 1, . . . , N − 1. How to calculate the derivative?

TU Darmstadt | GSC CE | Iryna Kulchytska | Time-parallel solution of the eddy current problem | 11/33

Parareal as the Newton-Raphson method (II) The Newton-Raphson iteration: for k = 0, 1, . . . U(k+1) = u0, U(k+1)

n

= F

• Tn, Tn−1, U(k)

n−1

• + ∂F(Tn, Tn−1, U)

∂U

• U(k+1)

n−1

− U(k)

n−1

• ,

where n = 1, . . . , N − 1. How to calculate the derivative? Cheap approximation by a coarse propagator G : ∂F(Tn, Tn−1, U) ∂U

• U(k+1)

n−1

− U(k)

n−1

• ≈

≈ G

• Tn, Tn−1, U(k+1)

n−1

• − G
• Tn, Tn−1, U(k)

n−1

• .

TU Darmstadt | GSC CE | Iryna Kulchytska | Time-parallel solution of the eddy current problem | 11/33

Parareal as the Newton-Raphson method (II) The Newton-Raphson iteration: for k = 0, 1, . . . U(k+1) = u0, U(k+1)

n

= F

• Tn, Tn−1, U(k)

n−1

• + ∂F(Tn, Tn−1, U)

∂U

• U(k+1)

n−1

− U(k)

n−1

• ,

where n = 1, . . . , N − 1. How to calculate the derivative? Cheap approximation by a coarse propagator G : ∂F(Tn, Tn−1, U) ∂U

• U(k+1)

n−1

− U(k)

n−1

• ≈

≈ G

• Tn, Tn−1, U(k+1)

n−1

• − G
• Tn, Tn−1, U(k)

n−1

• .

TU Darmstadt | GSC CE | Iryna Kulchytska | Time-parallel solution of the eddy current problem | 11/33

Parareal as the Newton-Raphson method (III) For k = 0, 1, . . . and n = 1, . . . , N solve U(k+1) = u0, U(k+1)

n

= F(Tn, Tn−1, U(k)

n−1) + G(Tn, Tn−1, U(k+1) n−1 ) − G(Tn, Tn−1, U(k) n−1).

TU Darmstadt | GSC CE | Iryna Kulchytska | Time-parallel solution of the eddy current problem | 12/33

Parareal as the Newton-Raphson method (III) For k = 0, 1, . . . and n = 1, . . . , N solve U(k+1) = u0, U(k+1)

n

= F(Tn, Tn−1, U(k)

n−1) + G(Tn, Tn−1, U(k+1) n−1 ) − G(Tn, Tn−1, U(k) n−1).

Propagators:

Fine ˜

U(k)

n

:= F

• Tn, Tn−1, U(k)

n−1

• : e.g., backward Euler’s method;

Coarse ¯

U(k)

n

:= G

• Tn, Tn−1, U(k)

n−1

• : lower-order scheme or the

same time-integrator as F but with coarser discretization.

TU Darmstadt | GSC CE | Iryna Kulchytska | Time-parallel solution of the eddy current problem | 12/33

Parareal as the Newton-Raphson method (III) For k = 0, 1, . . . and n = 1, . . . , N solve U(k+1) = u0, U(k+1)

n

= F(Tn, Tn−1, U(k)

n−1) + G(Tn, Tn−1, U(k+1) n−1 ) − G(Tn, Tn−1, U(k) n−1).

Propagators:

Fine ˜

U(k)

n

:= F

• Tn, Tn−1, U(k)

n−1

• : e.g., backward Euler’s method;

Coarse ¯

U(k)

n

:= G

• Tn, Tn−1, U(k)

n−1

• : lower-order scheme or the

same time-integrator as F but with coarser discretization. − → Parallel solution of fine problems;

TU Darmstadt | GSC CE | Iryna Kulchytska | Time-parallel solution of the eddy current problem | 12/33

Parareal as the Newton-Raphson method (III) For k = 0, 1, . . . and n = 1, . . . , N solve U(k+1) = u0, U(k+1)

n

= F(Tn, Tn−1, U(k)

n−1) + G(Tn, Tn−1, U(k+1) n−1 ) − G(Tn, Tn−1, U(k) n−1).

Propagators:

Fine ˜

U(k)

n

:= F

• Tn, Tn−1, U(k)

n−1

• : e.g., backward Euler’s method;

Coarse ¯

U(k)

n

:= G

• Tn, Tn−1, U(k)

n−1

• : lower-order scheme or the

same time-integrator as F but with coarser discretization. − → Parallel solution of fine problems; − → Sequential solution of coarse problems.

TU Darmstadt | GSC CE | Iryna Kulchytska | Time-parallel solution of the eddy current problem | 12/33

Considering the periodic constraint: PP-IC PP-IC: periodic parareal algorithm with initial value coarse problem: For k = 0, 1, . . . and n = 1, . . . , N solve U(k+1) = U(k)

N ,

U(k+1)

n

= F(Tn, Tn−1, U(k)

n−1) + G(Tn, Tn−1, U(k+1) n−1 ) − G(Tn, Tn−1, U(k) n−1).

TU Darmstadt | GSC CE | Iryna Kulchytska | Time-parallel solution of the eddy current problem | 13/33

Considering the periodic constraint: PP-IC PP-IC: periodic parareal algorithm with initial value coarse problem: For k = 0, 1, . . . and n = 1, . . . , N solve U(k+1) = U(k)

N ,

U(k+1)

n

= F(Tn, Tn−1, U(k)

n−1) + G(Tn, Tn−1, U(k+1) n−1 ) − G(Tn, Tn−1, U(k) n−1).

Propagators:

Fine ˜

U(k)

n

:= F

• Tn, Tn−1, U(k)

n−1

• : e.g., backward Euler’s method;

Coarse ¯

U(k)

n

:= G

• Tn, Tn−1, U(k)

n−1

• : lower-order scheme or the

same time-integrator as F but with coarser discretization.

TU Darmstadt | GSC CE | Iryna Kulchytska | Time-parallel solution of the eddy current problem | 13/33

Considering the periodic constraint: PP-IC PP-IC: periodic parareal algorithm with initial value coarse problem: For k = 0, 1, . . . and n = 1, . . . , N solve U(k+1) = U(k)

N ,

U(k+1)

n

= F(Tn, Tn−1, U(k)

n−1) + G(Tn, Tn−1, U(k+1) n−1 ) − G(Tn, Tn−1, U(k) n−1).

Propagators:

Fine ˜

U(k)

n

:= F

• Tn, Tn−1, U(k)

n−1

• : e.g., backward Euler’s method;

Coarse ¯

U(k)

n

:= G

• Tn, Tn−1, U(k)

n−1

• : lower-order scheme or the

same time-integrator as F but with coarser discretization. − → Parallel solution of fine problems;

TU Darmstadt | GSC CE | Iryna Kulchytska | Time-parallel solution of the eddy current problem | 13/33

Considering the periodic constraint: PP-IC PP-IC: periodic parareal algorithm with initial value coarse problem: For k = 0, 1, . . . and n = 1, . . . , N solve U(k+1) = U(k)

N ,

U(k+1)

n

= F(Tn, Tn−1, U(k)

n−1) + G(Tn, Tn−1, U(k+1) n−1 ) − G(Tn, Tn−1, U(k) n−1).

Propagators:

Fine ˜

U(k)

n

:= F

• Tn, Tn−1, U(k)

n−1

• : e.g., backward Euler’s method;

Coarse ¯

U(k)

n

:= G

• Tn, Tn−1, U(k)

n−1

• : lower-order scheme or the

same time-integrator as F but with coarser discretization. − → Parallel solution of fine problems; − → Sequential solution of coarse problems.

TU Darmstadt | GSC CE | Iryna Kulchytska | Time-parallel solution of the eddy current problem | 13/33

The PP-IC algorithm

1 initialize: U(0)

← U0 and k ← 0;

2 for n ← 1 to N do 3

U(0)

n

← ¯ U(0)

n

← G(Tn, Tn−1, U(0)

n−1);

4 end 5 while error (k) > tol do 6

parfor n ← 1 to N do

7

˜ U(k)

n

← F(Tn, Tn−1, U(k)

n−1);

8

end

9

U(k+1) ← U(k)

N ;

10

for n ← 1 to N do

11

¯ U(k+1)

n

← G(Tn, Tn−1, U(k+1)

n−1 );

12

U(k+1)

n

← ˜ U(k)

n

+ ¯ U(k+1)

n

− ¯ U(k)

n ;

13

end

14

increment: k ← k + 1;

15 end

TU Darmstadt | GSC CE | Iryna Kulchytska | Time-parallel solution of the eddy current problem | 14/33

The PP-IC algorithm

1 initialize: U(0)

← U0 and k ← 0;

2 for n ← 1 to N do 3

U(0)

n

← ¯ U(0)

n

← G(Tn, Tn−1, U(0)

n−1);

4 end 5 while error (k) > tol do 6

parfor n ← 1 to N do

7

˜ U(k)

n

← F(Tn, Tn−1, U(k)

n−1);

8

end

9

U(k+1) ← U(k)

N ;

10

for n ← 1 to N do

11

¯ U(k+1)

n

← G(Tn, Tn−1, U(k+1)

n−1 );

12

U(k+1)

n

← ˜ U(k)

n

+ ¯ U(k+1)

n

− ¯ U(k)

n ;

13

end

14

increment: k ← k + 1;

15 end

T0 T1 T2 T3 u(t) t

TU Darmstadt | GSC CE | Iryna Kulchytska | Time-parallel solution of the eddy current problem | 14/33

The PP-IC algorithm

1 initialize: U(0)

← U0 and k ← 0;

2 for n ← 1 to N do 3

U(0)

n

← ¯ U(0)

n

← G(Tn, Tn−1, U(0)

n−1);

4 end 5 while error (k) > tol do 6

parfor n ← 1 to N do

7

˜ U(k)

n

← F(Tn, Tn−1, U(k)

n−1);

8

end

9

U(k+1) ← U(k)

N ;

10

for n ← 1 to N do

11

¯ U(k+1)

n

← G(Tn, Tn−1, U(k+1)

n−1 );

12

U(k+1)

n

← ˜ U(k)

n

+ ¯ U(k+1)

n

− ¯ U(k)

n ;

13

end

14

increment: k ← k + 1;

15 end

T0 T1 T2 T3 u(t) t

TU Darmstadt | GSC CE | Iryna Kulchytska | Time-parallel solution of the eddy current problem | 14/33

The PP-IC algorithm

1 initialize: U(0)

← U0 and k ← 0;

2 for n ← 1 to N do 3

U(0)

n

← ¯ U(0)

n

← G(Tn, Tn−1, U(0)

n−1);

4 end 5 while error (k) > tol do 6

parfor n ← 1 to N do

7

˜ U(k)

n

← F(Tn, Tn−1, U(k)

n−1);

8

end

9

U(k+1) ← U(k)

N ;

10

for n ← 1 to N do

11

¯ U(k+1)

n

← G(Tn, Tn−1, U(k+1)

n−1 );

12

U(k+1)

n

← ˜ U(k)

n

+ ¯ U(k+1)

n

− ¯ U(k)

n ;

13

end

14

increment: k ← k + 1;

15 end

T0 T1 T2 T3 u(t) t

TU Darmstadt | GSC CE | Iryna Kulchytska | Time-parallel solution of the eddy current problem | 14/33

The PP-IC algorithm

1 initialize: U(0)

← U0 and k ← 0;

2 for n ← 1 to N do 3

U(0)

n

← ¯ U(0)

n

← G(Tn, Tn−1, U(0)

n−1);

4 end 5 while error (k) > tol do 6

parfor n ← 1 to N do

7

˜ U(k)

n

← F(Tn, Tn−1, U(k)

n−1);

8

end

9

U(k+1) ← U(k)

N ;

10

for n ← 1 to N do

11

¯ U(k+1)

n

← G(Tn, Tn−1, U(k+1)

n−1 );

12

U(k+1)

n

← ˜ U(k)

n

+ ¯ U(k+1)

n

− ¯ U(k)

n ;

13

end

14

increment: k ← k + 1;

15 end

T0 T1 T2 T3 u(t) t

TU Darmstadt | GSC CE | Iryna Kulchytska | Time-parallel solution of the eddy current problem | 14/33

The PP-IC algorithm

1 initialize: U(0)

← U0 and k ← 0;

2 for n ← 1 to N do 3

U(0)

n

← ¯ U(0)

n

← G(Tn, Tn−1, U(0)

n−1);

4 end 5 while error (k) > tol do 6

parfor n ← 1 to N do

7

˜ U(k)

n

← F(Tn, Tn−1, U(k)

n−1);

8

end

9

U(k+1) ← U(k)

N ;

10

for n ← 1 to N do

11

¯ U(k+1)

n

← G(Tn, Tn−1, U(k+1)

n−1 );

12

U(k+1)

n

← ˜ U(k)

n

+ ¯ U(k+1)

n

− ¯ U(k)

n ;

13

end

14

increment: k ← k + 1;

15 end

T0 T1 T2 T3 u(t) t

TU Darmstadt | GSC CE | Iryna Kulchytska | Time-parallel solution of the eddy current problem | 14/33

The PP-IC algorithm

1 initialize: U(0)

← U0 and k ← 0;

2 for n ← 1 to N do 3

U(0)

n

← ¯ U(0)

n

← G(Tn, Tn−1, U(0)

n−1);

4 end 5 while error (k) > tol do 6

parfor n ← 1 to N do

7

˜ U(k)

n

← F(Tn, Tn−1, U(k)

n−1);

8

end

9

U(k+1) ← U(k)

N ;

10

for n ← 1 to N do

11

¯ U(k+1)

n

← G(Tn, Tn−1, U(k+1)

n−1 );

12

U(k+1)

n

← ˜ U(k)

n

+ ¯ U(k+1)

n

− ¯ U(k)

n ;

13

end

14

increment: k ← k + 1;

15 end

T0 T1 T2 T3 u(t) t

TU Darmstadt | GSC CE | Iryna Kulchytska | Time-parallel solution of the eddy current problem | 14/33

The PP-IC algorithm

1 initialize: U(0)

← U0 and k ← 0;

2 for n ← 1 to N do 3

U(0)

n

← ¯ U(0)

n

← G(Tn, Tn−1, U(0)

n−1);

4 end 5 while error (k) > tol do 6

parfor n ← 1 to N do

7

˜ U(k)

n

← F(Tn, Tn−1, U(k)

n−1);

8

end

9

U(k+1) ← U(k)

N ;

10

for n ← 1 to N do

11

¯ U(k+1)

n

← G(Tn, Tn−1, U(k+1)

n−1 );

12

U(k+1)

n

← ˜ U(k)

n

+ ¯ U(k+1)

n

− ¯ U(k)

n ;

13

end

14

increment: k ← k + 1;

15 end

T0 T1 T2 T3 u(t) t

TU Darmstadt | GSC CE | Iryna Kulchytska | Time-parallel solution of the eddy current problem | 14/33

The PP-IC algorithm

1 initialize: U(0)

← U0 and k ← 0;

2 for n ← 1 to N do 3

U(0)

n

← ¯ U(0)

n

← G(Tn, Tn−1, U(0)

n−1);

4 end 5 while error (k) > tol do 6

parfor n ← 1 to N do

7

˜ U(k)

n

← F(Tn, Tn−1, U(k)

n−1);

8

end

9

U(k+1) ← U(k)

N ;

10

for n ← 1 to N do

11

¯ U(k+1)

n

← G(Tn, Tn−1, U(k+1)

n−1 );

12

U(k+1)

n

← ˜ U(k)

n

+ ¯ U(k+1)

n

− ¯ U(k)

n ;

13

end

14

increment: k ← k + 1;

15 end

T0 T1 T2 T3 u(t) t

TU Darmstadt | GSC CE | Iryna Kulchytska | Time-parallel solution of the eddy current problem | 14/33

The PP-IC algorithm

1 initialize: U(0)

← U0 and k ← 0;

2 for n ← 1 to N do 3

U(0)

n

← ¯ U(0)

n

← G(Tn, Tn−1, U(0)

n−1);

4 end 5 while error (k) > tol do 6

parfor n ← 1 to N do

7

˜ U(k)

n

← F(Tn, Tn−1, U(k)

n−1);

8

end

9

U(k+1) ← U(k)

N ;

10

for n ← 1 to N do

11

¯ U(k+1)

n

← G(Tn, Tn−1, U(k+1)

n−1 );

12

U(k+1)

n

← ˜ U(k)

n

+ ¯ U(k+1)

n

− ¯ U(k)

n ;

13

end

14

increment: k ← k + 1;

15 end

T0 T1 T2 T3 u(t) t

TU Darmstadt | GSC CE | Iryna Kulchytska | Time-parallel solution of the eddy current problem | 14/33

The PP-IC algorithm

1 initialize: U(0)

← U0 and k ← 0;

2 for n ← 1 to N do 3

U(0)

n

← ¯ U(0)

n

← G(Tn, Tn−1, U(0)

n−1);

4 end 5 while error (k) > tol do 6

parfor n ← 1 to N do

7

˜ U(k)

n

← F(Tn, Tn−1, U(k)

n−1);

8

end

9

U(k+1) ← U(k)

N ;

10

for n ← 1 to N do

11

¯ U(k+1)

n

← G(Tn, Tn−1, U(k+1)

n−1 );

12

U(k+1)

n

← ˜ U(k)

n

+ ¯ U(k+1)

n

− ¯ U(k)

n ;

13

end

14

increment: k ← k + 1;

15 end

T0 T1 T2 T3 u(t) t

TU Darmstadt | GSC CE | Iryna Kulchytska | Time-parallel solution of the eddy current problem | 14/33

The PP-IC algorithm

1 initialize: U(0)

← U0 and k ← 0;

2 for n ← 1 to N do 3

U(0)

n

← ¯ U(0)

n

← G(Tn, Tn−1, U(0)

n−1);

4 end 5 while error (k) > tol do 6

parfor n ← 1 to N do

7

˜ U(k)

n

← F(Tn, Tn−1, U(k)

n−1);

8

end

9

U(k+1) ← U(k)

N ;

10

for n ← 1 to N do

11

¯ U(k+1)

n

← G(Tn, Tn−1, U(k+1)

n−1 );

12

U(k+1)

n

← ˜ U(k)

n

+ ¯ U(k+1)

n

− ¯ U(k)

n ;

13

end

14

increment: k ← k + 1;

15 end

T0 T1 T2 T3 u(t) t

TU Darmstadt | GSC CE | Iryna Kulchytska | Time-parallel solution of the eddy current problem | 14/33

The PP-IC algorithm

1 initialize: U(0)

← U0 and k ← 0;

2 for n ← 1 to N do 3

U(0)

n

← ¯ U(0)

n

← G(Tn, Tn−1, U(0)

n−1);

4 end 5 while error (k) > tol do 6

parfor n ← 1 to N do

7

˜ U(k)

n

← F(Tn, Tn−1, U(k)

n−1);

8

end

9

U(k+1) ← U(k)

N ;

10

for n ← 1 to N do

11

¯ U(k+1)

n

← G(Tn, Tn−1, U(k+1)

n−1 );

12

U(k+1)

n

← ˜ U(k)

n

+ ¯ U(k+1)

n

− ¯ U(k)

n ;

13

end

14

increment: k ← k + 1;

15 end

T0 T1 T2 T3 u(t) t

TU Darmstadt | GSC CE | Iryna Kulchytska | Time-parallel solution of the eddy current problem | 14/33

The PP-IC algorithm

1 initialize: U(0)

← U0 and k ← 0;

2 for n ← 1 to N do 3

U(0)

n

← ¯ U(0)

n

← G(Tn, Tn−1, U(0)

n−1);

4 end 5 while error (k) > tol do 6

parfor n ← 1 to N do

7

˜ U(k)

n

← F(Tn, Tn−1, U(k)

n−1);

8

end

9

U(k+1) ← U(k)

N ;

10

for n ← 1 to N do

11

¯ U(k+1)

n

← G(Tn, Tn−1, U(k+1)

n−1 );

12

U(k+1)

n

← ˜ U(k)

n

+ ¯ U(k+1)

n

− ¯ U(k)

n ;

13

end

14

increment: k ← k + 1;

15 end

T0 T1 T2 T3 u(t) t

TU Darmstadt | GSC CE | Iryna Kulchytska | Time-parallel solution of the eddy current problem | 14/33

Overview Introduction

Motivation The eddy current problem

Parallel-in-time solution

The Parareal method for IVPs Numerical example: coaxial cable model

Fourier basis for time-periodic systems

Coarse solution by spectral collocation Numerical results: coaxial cable model Systems with nonsmooth excitations

Conclusions and outlook

TU Darmstadt | GSC CE | Iryna Kulchytska | Time-parallel solution of the eddy current problem | 15/33

Numerical example: coaxial cable model For t ∈ [0, T] find u(t) ∈ RNdof s.t. Mdtu(t) + Ku(t) = Xi(t), u(0) = u(T), space-discrete eddy current problem with Ndof = 2269. T = 0.02, ω = 2π/T, i(t) = 100 sin (ωt) Dimensions: rCu = 0.254cm, rAir = 1.27cm, rFe = 2.54cm.

Fig.: Wire inside of a steel tube.

TU Darmstadt | GSC CE | Iryna Kulchytska | Time-parallel solution of the eddy current problem | 16/33

Numerical example: coaxial cable model For t ∈ [0, T] find u(t) ∈ RNdof s.t. Mdtu(t) + Ku(t) = Xi(t), u(0) = u(T), space-discrete eddy current problem with Ndof = 2269. T = 0.02, ω = 2π/T, i(t) = 100 sin (ωt) Dimensions: rCu = 0.254cm, rAir = 1.27cm, rFe = 2.54cm.

Fig.: Wire inside of a steel tube.

TU Darmstadt | GSC CE | Iryna Kulchytska | Time-parallel solution of the eddy current problem | 16/33

Coaxial cable model: convergence PP-IC

Fig.: PP-IC: initial approximation. Fig.: Solution after the 1st iteration.

Iteration 1 50 125 200 260

• Rel. error

2.3 · 10−1 8.0 · 10−4 7.3 · 10−5 6.7 · 10−6 9.9 · 10−7

TU Darmstadt | GSC CE | Iryna Kulchytska | Time-parallel solution of the eddy current problem | 17/33

Coaxial cable model: speed-up PP-IC

Fig.: PP-IC: solution within tol ≈ 10−6. Fig.: Sequential steady-state solution.

Sequential time: 35.6 min; Number of periods: 261 PP-IC time: 1.5 min − → 23 times faster.

TU Darmstadt | GSC CE | Iryna Kulchytska | Time-parallel solution of the eddy current problem | 18/33

Overview Introduction

Motivation The eddy current problem

Parallel-in-time solution

The Parareal method for IVPs Numerical example: coaxial cable model

Fourier basis for time-periodic systems

Coarse solution by spectral collocation Numerical results: coaxial cable model Systems with nonsmooth excitations

Conclusions and outlook

TU Darmstadt | GSC CE | Iryna Kulchytska | Time-parallel solution of the eddy current problem | 19/33

Time-periodic problem to be solved: PP-PC PP-PC: periodic parareal algorithm with periodic coarse problem:      I . . . −G (TN, TN−1, ·) −G (T1, T0, ·) I . . . ... ... . . . . . . −G (TN−1, TN−2, ·) I            U(k+1) U(k+1)

1.

. . U(k+1)

N−1

      = =         F

• TN, TN−1, U(k)

N−1

• − G
• TN, TN−1, U(k)

N−1

• F
• T1, T0, U(k)
• − G
• T1, T0, U(k)
• .

. . F

• TN−1, TN−2, U(k)

N−2

• − G
• TN−1, TN−2, U(k)

N−2

• 

      

TU Darmstadt | GSC CE | Iryna Kulchytska | Time-parallel solution of the eddy current problem | 20/33

Time-periodic problem to be solved: PP-PC PP-PC: periodic parareal algorithm with periodic coarse problem:      I . . . −G (TN, TN−1, ·) −G (T1, T0, ·) I . . . ... ... . . . . . . −G (TN−1, TN−2, ·) I            U(k+1) U(k+1)

1.

. . U(k+1)

N−1

      = =         F

• TN, TN−1, U(k)

N−1

• − G
• TN, TN−1, U(k)

N−1

• F
• T1, T0, U(k)
• − G
• T1, T0, U(k)
• .

. . F

• TN−1, TN−2, U(k)

N−2

• − G
• TN−1, TN−2, U(k)

N−2

• 

       − → Large size and inconvenient structure!

TU Darmstadt | GSC CE | Iryna Kulchytska | Time-parallel solution of the eddy current problem | 20/33

PP-PC: backward Euler’s method as coarse solver Assume y(k+1)

n

:= G

• Tn, Tn−1, U(k+1)

n−1

• is defined by the backward

Euler’s method: M ∆t + K

• =:Q

y(k+1)

n

− M ∆t

• =:C

U(k+1)

n−1

= g(Tn) for n = 1, . . . , N with ∆t = T/N.

TU Darmstadt | GSC CE | Iryna Kulchytska | Time-parallel solution of the eddy current problem | 21/33

PP-PC: backward Euler’s method as coarse solver Assume y(k+1)

n

:= G

• Tn, Tn−1, U(k+1)

n−1

• is defined by the backward

Euler’s method: M ∆t + K

• =:Q

y(k+1)

n

− M ∆t

• =:C

U(k+1)

n−1

= g(Tn) for n = 1, . . . , N with ∆t = T/N.The system reads:      Q −C −C Q ... ... −C Q     

• =:G

      U(k+1) U(k+1)

1.

. . U(k+1)

N−1

      =         QF

• TN, TN−1, U(k)

N−1

• − CU(k)

N−1)

QF

• T1, T0, U(k)
• − CU(k)

. . . QF

• TN−1, TN−2, U(k)

N−2

• − CU(k)

N−2

       

• =:r(k)

TU Darmstadt | GSC CE | Iryna Kulchytska | Time-parallel solution of the eddy current problem | 21/33

PP-PC: backward Euler’s method as coarse solver Assume y(k+1)

n

:= G

• Tn, Tn−1, U(k+1)

n−1

• is defined by the backward

Euler’s method: M ∆t + K

• =:Q

y(k+1)

n

− M ∆t

• =:C

U(k+1)

n−1

= g(Tn) for n = 1, . . . , N with ∆t = T/N.The system reads:      Q −C −C Q ... ... −C Q     

• =:G

      U(k+1) U(k+1)

1.

. . U(k+1)

N−1

      =         QF

• TN, TN−1, U(k)

N−1

• − CU(k)

N−1)

QF

• T1, T0, U(k)
• − CU(k)

. . . QF

• TN−1, TN−2, U(k)

N−2

• − CU(k)

N−2

       

• =:r(k)

TU Darmstadt | GSC CE | Iryna Kulchytska | Time-parallel solution of the eddy current problem | 21/33

Coarse solution via the spectral basis U0, . . . UN−1 are values of a periodic function at t0, . . . , tN−1.

TU Darmstadt | GSC CE | Iryna Kulchytska | Time-parallel solution of the eddy current problem | 22/33

Coarse solution via the spectral basis U0, . . . UN−1 are values of a periodic function at t0, . . . , tN−1. − → can be expanded into the Fourier series over the period [0, T] :

TU Darmstadt | GSC CE | Iryna Kulchytska | Time-parallel solution of the eddy current problem | 22/33

Coarse solution via the spectral basis U0, . . . UN−1 are values of a periodic function at t0, . . . , tN−1. − → can be expanded into the Fourier series over the period [0, T] : U(t) =

• m∈M

ˆ Um exp (ıωmt), with the frequencies ωm = 2πm/T, and M := {−N/2 + 1, . . . , N/2}.

TU Darmstadt | GSC CE | Iryna Kulchytska | Time-parallel solution of the eddy current problem | 22/33

Coarse solution via the spectral basis U0, . . . UN−1 are values of a periodic function at t0, . . . , tN−1. − → can be expanded into the Fourier series over the period [0, T] : U(t) =

• m∈M

ˆ Um exp (ıωmt), with the frequencies ωm = 2πm/T, and M := {−N/2 + 1, . . . , N/2}. At the kth iteration U(k+1)

n

=

• m∈M

ˆ U(k+1)

m

exp (ıωmtn), n = 0, . . . , N − 1.

TU Darmstadt | GSC CE | Iryna Kulchytska | Time-parallel solution of the eddy current problem | 22/33

Coarse solution via the spectral basis U0, . . . UN−1 are values of a periodic function at t0, . . . , tN−1. − → can be expanded into the Fourier series over the period [0, T] : U(t) =

• m∈M

ˆ Um exp (ıωmt), with the frequencies ωm = 2πm/T, and M := {−N/2 + 1, . . . , N/2}. At the kth iteration U(k+1)

n

=

• m∈M

ˆ U(k+1)

m

exp (ıωmtn), n = 0, . . . , N − 1. − → into the PP-PC system.

TU Darmstadt | GSC CE | Iryna Kulchytska | Time-parallel solution of the eddy current problem | 22/33

Frequency domain solution of the coarse problem Let F denote the discrete Fourier transform matrix: Fpq = 1 √ N exp(−ıωptq) and ˜ F = F ⊗ I, with I ∈ RNdof ×Ndof .

TU Darmstadt | GSC CE | Iryna Kulchytska | Time-parallel solution of the eddy current problem | 23/33

Frequency domain solution of the coarse problem Let F denote the discrete Fourier transform matrix: Fpq = 1 √ N exp(−ıωptq) and ˜ F = F ⊗ I, with I ∈ RNdof ×Ndof . Restriction to the spectral basis in the discrete PP-PC system is equivalent to the solution of

• ˜

F G ˜ FH

• =:ˆ

G

ˆ U(k+1) = ˜ Fr(k)

• =:ˆ

r(k)

for ˆ U(k+1)T =

• ˆ

U(k+1)T

−N/2+1, . . . , ˆ

U(k+1)T

N/2

• in the frequency domain.

TU Darmstadt | GSC CE | Iryna Kulchytska | Time-parallel solution of the eddy current problem | 23/33

Frequency domain solution of the coarse problem Let F denote the discrete Fourier transform matrix: Fpq = 1 √ N exp(−ıωptq) and ˜ F = F ⊗ I, with I ∈ RNdof ×Ndof . Restriction to the spectral basis in the discrete PP-PC system is equivalent to the solution of

• ˜

F G ˜ FH

• =:ˆ

G

ˆ U(k+1) = ˜ Fr(k)

• =:ˆ

r(k)

for ˆ U(k+1)T =

• ˆ

U(k+1)T

−N/2+1, . . . , ˆ

U(k+1)T

N/2

• in the frequency domain.

G − block-circulant → ˆ G − block-diagonal:

TU Darmstadt | GSC CE | Iryna Kulchytska | Time-parallel solution of the eddy current problem | 23/33

Frequency domain solution of the coarse problem Let F denote the discrete Fourier transform matrix: Fpq = 1 √ N exp(−ıωptq) and ˜ F = F ⊗ I, with I ∈ RNdof ×Ndof . Restriction to the spectral basis in the discrete PP-PC system is equivalent to the solution of

• ˜

F G ˜ FH

• =:ˆ

G

ˆ U(k+1) = ˜ Fr(k)

• =:ˆ

r(k)

for ˆ U(k+1)T =

• ˆ

U(k+1)T

−N/2+1, . . . , ˆ

U(k+1)T

N/2

• in the frequency domain.

G − block-circulant → ˆ G − block-diagonal: ˆ Gpp = Q − C exp(−ı∆tωp).

TU Darmstadt | GSC CE | Iryna Kulchytska | Time-parallel solution of the eddy current problem | 23/33

Frequency domain solution of the coarse problem Solve for each harmonic component independently:

• Q − C exp(−ı∆tωp)
• ˆ

U(k+1)

p

= ˆ r(k)

p , p = −N/2 + 1, . . . , N/2.

TU Darmstadt | GSC CE | Iryna Kulchytska | Time-parallel solution of the eddy current problem | 24/33

Frequency domain solution of the coarse problem Solve for each harmonic component independently:

• Q − C exp(−ı∆tωp)
• ˆ

U(k+1)

p

= ˆ r(k)

p , p = −N/2 + 1, . . . , N/2.

− → Allows parallelization on the coarse level: one needs to solve N separate systems of Ndof equations.

TU Darmstadt | GSC CE | Iryna Kulchytska | Time-parallel solution of the eddy current problem | 24/33

Frequency domain solution of the coarse problem Solve for each harmonic component independently:

• Q − C exp(−ı∆tωp)
• ˆ

U(k+1)

p

= ˆ r(k)

p , p = −N/2 + 1, . . . , N/2.

− → Allows parallelization on the coarse level: one needs to solve N separate systems of Ndof equations. Solution in the time domain is obtained by the inverse Fourier transformation: U(k+1) = ˜ FH ˆ U(k+1).

TU Darmstadt | GSC CE | Iryna Kulchytska | Time-parallel solution of the eddy current problem | 24/33

Frequency domain solution of the coarse problem Solve for each harmonic component independently:

• Q − C exp(−ı∆tωp)
• ˆ

U(k+1)

p

= ˆ r(k)

p , p = −N/2 + 1, . . . , N/2.

− → Allows parallelization on the coarse level: one needs to solve N separate systems of Ndof equations. Solution in the time domain is obtained by the inverse Fourier transformation: U(k+1) = ˜ FH ˆ U(k+1). − → Can be calculated using the fast Fourier transform algorithm; − → Matrices ˜ F and ˜ FH do not have to be explicitly constructed.

TU Darmstadt | GSC CE | Iryna Kulchytska | Time-parallel solution of the eddy current problem | 24/33

Overview Introduction

Motivation The eddy current problem

Parallel-in-time solution

The Parareal method for IVPs Numerical example: coaxial cable model

Fourier basis for time-periodic systems

Coarse solution by spectral collocation Numerical results: coaxial cable model Systems with nonsmooth excitations

Conclusions and outlook

TU Darmstadt | GSC CE | Iryna Kulchytska | Time-parallel solution of the eddy current problem | 25/33

Coaxial cable model: convergence PP-PC

Fig.: Solution after the 1st iteration Fig.: Solution after the 10th iteration

Iteration 1 3 5 7 10

• Rel. error

2.0 · 10−1 6.3 · 10−3 3.2 · 10−4 2.0 · 10−5 3.7 · 10−7

TU Darmstadt | GSC CE | Iryna Kulchytska | Time-parallel solution of the eddy current problem | 26/33

Comparison of computational time: strong scaling

1 5 10 25 50 100 101 102 103

Number of CPUs Ncpu time / s

IVP PP-IC PP-PC LU PP-PC FT All the results are obtained for Ndof = 2269 degrees of freedom, rel. error ≈ 10−6.

TU Darmstadt | GSC CE | Iryna Kulchytska | Time-parallel solution of the eddy current problem | 27/33

Overview Introduction

Motivation The eddy current problem

Parallel-in-time solution

The Parareal method for IVPs Numerical example: coaxial cable model

Fourier basis for time-periodic systems

Coarse solution by spectral collocation Numerical results: coaxial cable model Systems with nonsmooth excitations

Conclusions and outlook

TU Darmstadt | GSC CE | Iryna Kulchytska | Time-parallel solution of the eddy current problem | 28/33

Systems with nonsmooth excitations PWM (pulse width modulation): high-order components in the frequency spectra of exciting currents/voltages. Mdtu(t) = ¯ f(t, u) +˜ f(t), t ∈ I u(0) = u(T),

¯

f smooth input;

˜

f piecewise continuous, square integrable on I. − → Very fine discretization required to capture the pulses.

Fig.: PWM signal with a switching frequency of 20 kHz and a sine wave

• f 50 Hz.

TU Darmstadt | GSC CE | Iryna Kulchytska | Time-parallel solution of the eddy current problem | 29/33

Comparison of computational time: PWM input

1 5 10 25 50 101 102 103 104 105

Number of CPUs Ncpu time / s

IVP PP-IC PP-PC LU PP-PC FT All the results are obtained for Ndof = 2269 degrees of freedom, rel. error ≈ 10−6.

TU Darmstadt | GSC CE | Iryna Kulchytska | Time-parallel solution of the eddy current problem | 30/33

Conclusions and outlook Conclusions

Parareal method significantly accelerates convergence to the

steady state;

TU Darmstadt | GSC CE | Iryna Kulchytska | Time-parallel solution of the eddy current problem | 31/33

Conclusions and outlook Conclusions

Parareal method significantly accelerates convergence to the

steady state; Outlook

Consider nonlinear problems; Application to an electrical machine; TU Darmstadt | GSC CE | Iryna Kulchytska | Time-parallel solution of the eddy current problem | 31/33

Thank you! ACKNOWLEDGEMENT This work has been supported by the Excellence Initiative of the German Federal and State Governments and the Graduate School of Computational Engineering at Technische Universität Darmstadt.

• S. Schöps, I. Niyonzima, and M. Clemens, Parallel-in-time Simulation of the

Eddy Current Problem using Parareal, 21st Conference on the Computation of Electromagnetic Fields (COMPUMAG 2017), Daejeon, Korea, June 2017.

• M. J. Gander, Y.-L. Jiang, B. Song, and H. Zhang, Analysis of two parareal

algorithms for time-periodic problems, SIAM J. Sci. Comput. 35 (5), 2013.

• J. Gyselinck et. al. A General Method for the Frequency Domain FE Modeling of

Rotating Electromagnetic Devices, IEEE Trans. Magn. 39 (3), 2003.

TU Darmstadt | GSC CE | Iryna Kulchytska | Time-parallel solution of the eddy current problem | 32/33