Krylov Subspace Methods for an Initial Value Problem Arising in TEM - - PowerPoint PPT Presentation

Krylov Subspace Methods for an Initial Value Problem Arising in TEM

Martin Afanasjew1 Ralph-Uwe Börner2 Oliver G. Ernst1 Michael Eiermann1 Stefan Güttel1 Klaus Spitzer2

1Institut für Numerische Mathematik und Optimierung 2Institut für Geophysik

TU Bergakademie Freiberg, Germany

August 23, 2007 Computational Methods with Applications, Harrachov

Outline

Problem Description Spatial Discretization Computational Strategies Time-Stepping Methods Krylov Subspace Methods Numerical Examples Summary and Future Work

Outline Krylov Subspace Methods for TEM 23/08/2007, Harrachov Martin Afanasjew et al., TU Freiberg 1

Problem Description

Problem Description Krylov Subspace Methods for TEM 23/08/2007, Harrachov Martin Afanasjew et al., TU Freiberg 2

Our Aim

Our aim:

◮ Simulate propagation of transient electromagnetic fields

(TEM) in the subsurface.

◮ Fields are a response to controlled electromagnetic sources.

◮ Here: Vertical magnetic dipole.

◮ Solve the forward problem.

Practical aspects:

◮ TEM is an important method in geophysical exploration to

infer properties of the subsurface.

Problem Description Krylov Subspace Methods for TEM 23/08/2007, Harrachov Martin Afanasjew et al., TU Freiberg 3

Typical Setup

Problem Description Krylov Subspace Methods for TEM 23/08/2007, Harrachov Martin Afanasjew et al., TU Freiberg 4

Governing Equations

Quasi-static Maxwell’s equations: ∇×

1

µ ∇× e

• + ∂tσe = −∂tje

(Maxwell) where e = e(x, t) is the electric field, µ = µ(x) is the magnetic permeability, σ = σ(x) is the electric conductivity, je = je(x, t) is the impressed source current density with x ∈ Ω ⊂ R3 and t ∈ R.

Problem Description Krylov Subspace Methods for TEM 23/08/2007, Harrachov Martin Afanasjew et al., TU Freiberg 5

Further Assumptions

◮ Typically, the spatial domain Ω is a parallelepiped with its

upper boundary at ground level or above it.

◮ We assume the perfect conductor boundary condition

n × e = 0 on ∂Ω.

◮ The impressed source current is typically of shut-off type, i.e.

• f the form

je(x, t) = q(x) H(−t) with H denoting the Heaviside unit step function and the vector field q being the spatial current pattern.

◮ In our case the right-hand side of (Maxwell) vanishes since we

are looking at times t > 0.

Problem Description Krylov Subspace Methods for TEM 23/08/2007, Harrachov Martin Afanasjew et al., TU Freiberg 6

Spatial Discretization

Spatial Discretization Krylov Subspace Methods for TEM 23/08/2007, Harrachov Martin Afanasjew et al., TU Freiberg 7

Spatial Discretization

Subdivision of spatial domain Ω:

◮ Graded grid. ◮ Increasing cell size as we move away from the source. ◮ Staggered grid (Yee grid).

◮ Electric components e at the center of the edges. ◮ Magnetic components h at the center of the faces. ◮ System of elementary electric and magnetic loops.

y z x ex ex ex ey ey ey ez ez ez hx hy hz ex ex ez ez hy hy hz hz

Spatial Discretization Krylov Subspace Methods for TEM 23/08/2007, Harrachov Martin Afanasjew et al., TU Freiberg 8

Example of a Graded Grid

−2000 −1000 1000 2000 −2000 −1000 1000 2000 −2000 −1000 1000 2000 y x z

Spatial Discretization Krylov Subspace Methods for TEM 23/08/2007, Harrachov Martin Afanasjew et al., TU Freiberg 9

Computational Strategies

Computational Strategies Krylov Subspace Methods for TEM 23/08/2007, Harrachov Martin Afanasjew et al., TU Freiberg 10

Computational Strategies

◮ We usually start with an initial solution to the electric field e0

at a time t0 > 0. Our interest lies in computing ej = e(tj) at few times tj with tj−1 < tj for 0 < j ≤ n.

◮ Depending on the used method we have different possibilities:

◮ Small steps calculating ej from ej−1. ◮ Steps of increasing size calculating ej from e0. ◮ One big step calculating en from e0.

t

t0 t1 t2 t2 · · · tn

Computational Strategies Krylov Subspace Methods for TEM 23/08/2007, Harrachov Martin Afanasjew et al., TU Freiberg 11

Time-Stepping Methods

Time-Stepping Methods Krylov Subspace Methods for TEM 23/08/2007, Harrachov Martin Afanasjew et al., TU Freiberg 12

Du Fort-Frankel Method

◮ Proposed by [Wang & Hohmann, 1993]. ◮ Explicit. ◮ Solves coupled first-order Maxwell’s equations. ◮ Uses Yee grid for spatial dicretization. ◮ Computes time-interleaved electric fields ej and magnetic

fields hj in a leap-frog type iteration. t

e update: h update: ej−1 tj−1 hj−1 ej tj hj ej+1 tj+1 hj+1

Time-Stepping Methods Krylov Subspace Methods for TEM 23/08/2007, Harrachov Martin Afanasjew et al., TU Freiberg 13

Stability of the Du Fort-Frankel Method

This method was shown to be stable if ∆tj = tj+1 − tj < ∆xmin

µminσmintj

6 where ∆xmin is the smallest grid size, µmin is the minimal magnetic permeability, σmin is the minimal electric conductivity.

Time-Stepping Methods Krylov Subspace Methods for TEM 23/08/2007, Harrachov Martin Afanasjew et al., TU Freiberg 14

Krylov Subspace Methods

Krylov Subspace Methods Krylov Subspace Methods for TEM 23/08/2007, Harrachov Martin Afanasjew et al., TU Freiberg 15

Rewriting the Problem (I)

To apply Krylov subspace methods we have to rewrite our problem. Starting with (Maxwell) we get ∂te = − 1 σ ∇×

1

µ ∇× e

• .

(PDE) This reduces to the solution of a linear first-order ordinary differential equation ∂te = Ae, e(t0) = e0, (ODE) where A represents the discrete action of −1/σ ∇× (1/µ ∇× ·) on the spatial discretization of the electric field e.

Krylov Subspace Methods Krylov Subspace Methods for TEM 23/08/2007, Harrachov Martin Afanasjew et al., TU Freiberg 16

Rewriting the Problem (II)

An explicit solution of (ODE) is given by e(t) = e(t−t0)A e0. Thus, we have to evaluate the exponential function for a sparse matrix times a vector, which is what Krylov subspace methods are well suited for.

Krylov Subspace Methods Krylov Subspace Methods for TEM 23/08/2007, Harrachov Martin Afanasjew et al., TU Freiberg 17

Arnoldi/Lanczos Procedure

◮ We define the m-th Krylov subspace as follows

Km(A, b) := span

• b, Ab, A2b, . . . , Am−1b
• .

◮ Generate an orthornormal basis Vm = [v1, v2, . . . , vm] of

Km(A, b) using a Gram-Schmidt procedure/three-term recurrence relation satisfying V ⊤

m AVm = Hm

with Hm ∈ Rm×m upper Hessenberg/tridiagonal.

◮ Calculate the Arnoldi approximation of order m

f m := b Vmf (Hm) [1, 0, . . . , 0]⊤

Krylov Subspace Methods Krylov Subspace Methods for TEM 23/08/2007, Harrachov Martin Afanasjew et al., TU Freiberg 18

Time-Stepping, Recycling, and Restarts (I)

◮ Remember: Given e0 at time t0 we are interested in evaluating

the electric fields e1, e2, . . . , en at times t1 < t2 < . . . < tn from an interval [t0, tn].

◮ Time-stepped Arnoldi

◮ In each time step we compute

f m

j+1 ∈ Km(A, f m j ) for f (x) = e(tj+1−tj)x

with f m

0 = e0 and m = m(j) ∼ (tj+1 − tj) A1/2. Such a

choice of m(j) guarantees a certain relative error for our approximation.

◮ This approach requires us to build a new Krylov subspace in

every time step.

Krylov Subspace Methods Krylov Subspace Methods for TEM 23/08/2007, Harrachov Martin Afanasjew et al., TU Freiberg 19

Time-Stepping, Recycling, and Restarts (II)

◮ Arnoldi with Recycling

◮ Similarly to the time-stepped variant we compute in each step

f m

j ∈ Km(A, e0) for f (x) = e(tj−t0)x

with m = m(j) ∼ (tj − t0) A1/2. This allows us to reuse basis vectors from the j − 1-th step and only compute the additional basis vectors vm(j)+1, vm(j)+2, . . . , vm(j+1).

◮ Restarted Arnoldi Method

◮ If A cannot be symmetrized the unrestarted Arnoldi method

might require a prohibitively high number of basis vectors. To

• vercome this problem a restared Arnoldi method was

proposed in [Eiermann & Ernst, 2006] where we discard all but the last basis vector after every m steps and start to build a new Krylov subspace starting with this last vector.

◮ See Stefan Güttel’s talk today at 18:20 (same session).

Krylov Subspace Methods Krylov Subspace Methods for TEM 23/08/2007, Harrachov Martin Afanasjew et al., TU Freiberg 20

Numerical Examples

Numerical Examples Krylov Subspace Methods for TEM 23/08/2007, Harrachov Martin Afanasjew et al., TU Freiberg 21

Setup

◮ Vertical magnetic dipole of unit strength located at the origin. ◮ Yee discretization of (PDE). ◮ Grid with 58 × 58 × 58 cells, i.e. 565326 unknowns. ◮ Constant coefficients µ = 1.26 · 10−6, σ = 1.00 · 10−1. ◮ 24 logarithmically equidistant time steps from the interval

[t0, tn] =

10−6, 10−3 seconds.

◮ Everything implemented in pure MATLAB (Release 2007a). Numerical Examples Krylov Subspace Methods for TEM 23/08/2007, Harrachov Martin Afanasjew et al., TU Freiberg 22

Du Fort-Frankel Method vs. ROCK4

5 10 15 20 25 10 20 30 40 Computation time in seconds Time step number

• verall time in s: 266.2893

Du Fort−Frankel 5 10 15 20 25 10 20 30 40 Computation time in seconds Time step number

• verall time in s: 300.873

ROCK4

Numerical Examples Krylov Subspace Methods for TEM 23/08/2007, Harrachov Martin Afanasjew et al., TU Freiberg 23

Lanczos Method: Time-Stepping vs. Recycling (I)

10

−6

10

−5

10

−4

10

−3

10

−3

10

−2

10

−1

10 10

1

Time step Relative error of Krylov approx. dim(K) = 25 dim(K) = 50 dim(K) = 100 dim(K) = 150

Numerical Examples Krylov Subspace Methods for TEM 23/08/2007, Harrachov Martin Afanasjew et al., TU Freiberg 24

Lanczos Method: Time-Stepping vs. Recycling (II)

5 10 15 20 25 5 10 15 20 25 30 Computation time in seconds Time step number

• verall time in s: 159.4429
• verall time in s: 108.8866

Lanczos with time−stepping Lanczos with recycling 5 10 15 20 25 100 200 300 400 500 600 Dimension of Krylov space Time step number Lanczos with recycling Lanczos with time−stepping

Numerical Examples Krylov Subspace Methods for TEM 23/08/2007, Harrachov Martin Afanasjew et al., TU Freiberg 25

Restarted Arnoldi (I)

◮ Algorithm 1: Proposed in [Eiermann & Ernst, 2006]. ◮ Algorithm 2: Improved variant of Algorithm 1. ◮ m = ∞ (I): Standard Lanczos (unrestarted). ◮ m = ∞ (II): Two-pass Lanczos (unrestarted).

Algorithm 1 Algorithm 2 m time [s] mvp acc. time [s] mvp acc. ∞ (I) 118 1072 9.93e-13 86 1072 9.93e-13 ∞ (II) 176 2144 9.93e-13 144 2144 9.93e-13 90 273 1350 1.92e-13 118 1350 2.01e-13 70 339 1400 3.28e-13 112 1400 9.13e-13 50 613 1600 2.10e-13 very slow convergence 30 2014 2040 5.64e-13 divergence

Numerical Examples Krylov Subspace Methods for TEM 23/08/2007, Harrachov Martin Afanasjew et al., TU Freiberg 26

Restarted Arnoldi (II)

500 1000 1500 2000 10

−16

10

−14

10

−12

10

−10

10

−8

10

−6

10

−4

matrix−vector products ||error||

• Alg. 1, m=∞
• Alg. 2, m=∞
• Alg. 1, m=90
• Alg. 2, m=90
• Alg. 1, m=70
• Alg. 2, m=70
• Alg. 1, m=50
• Alg. 2, m=50
• Alg. 1, m=30
• Alg. 2, m=30

stopping error

Numerical Examples Krylov Subspace Methods for TEM 23/08/2007, Harrachov Martin Afanasjew et al., TU Freiberg 27

Summary and Future Work

Summary and Future Work Krylov Subspace Methods for TEM 23/08/2007, Harrachov Martin Afanasjew et al., TU Freiberg 28

Summary and Future Work

◮ Summary

◮ We successfully applied various methods to our model problem. ◮ Restarted Krylov subspace methods might become a viable

alternative to the well known Du Fort-Frankel and Spectral Lanczos Decomposition methods.

◮ Future work

◮ Implicit time stepping. ◮ FE discretization. ◮ Heterogeneous materials. ◮ Error estimates.

◮ References

◮ Implementation of a Restarted Krylov Subspace Method for

the Evaluation of Matrix Functions [A, Eiermann, Ernst, & Güttel, 2007].

Summary and Future Work Krylov Subspace Methods for TEM 23/08/2007, Harrachov Martin Afanasjew et al., TU Freiberg 29