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 Afanasjew 1 Ralph-Uwe Börner 2 Oliver G. Ernst 1 Michael Eiermann 1 Stefan Güttel 1 Klaus Spitzer 2 1 Institut für Numerische Mathematik und Optimierung 2 Institut 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 1 23/08/2007, Harrachov Martin Afanasjew et al., TU Freiberg

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

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 3 23/08/2007, Harrachov Martin Afanasjew et al., TU Freiberg

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

Governing Equations Quasi-static Maxwell’s equations: � 1 � + ∂ t σ e = − ∂ t j e ∇× µ ∇× e (Maxwell) where e = e ( x , t ) is the electric field, µ = µ ( x ) is the magnetic permeability, σ = σ ( x ) is the electric conductivity, j e = j e ( x , t ) is the impressed source current density with x ∈ Ω ⊂ R 3 and t ∈ R . Problem Description Krylov Subspace Methods for TEM 5 23/08/2007, Harrachov Martin Afanasjew et al., TU Freiberg

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. of the form j e ( 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 6 23/08/2007, Harrachov Martin Afanasjew et al., TU Freiberg

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

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. x h y y h z e y h z h z e x e x e x z e y e z e z h y h y h x e z e z e z e x e x e y Spatial Discretization Krylov Subspace Methods for TEM 8 23/08/2007, Harrachov Martin Afanasjew et al., TU Freiberg

Example of a Graded Grid 2000 1000 0 z −1000 2000 1000 −2000 0 −2000 −1000 −1000 0 1000 −2000 2000 y x Spatial Discretization Krylov Subspace Methods for TEM 9 23/08/2007, Harrachov Martin Afanasjew et al., TU Freiberg

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

Computational Strategies ◮ We usually start with an initial solution to the electric field e 0 at a time t 0 > 0 . Our interest lies in computing e j = e ( t j ) at few times t j with t j − 1 < t j for 0 < j ≤ n . ◮ Depending on the used method we have different possibilities: ◮ Small steps calculating e j from e j − 1 . ◮ Steps of increasing size calculating e j from e 0 . ◮ One big step calculating e n from e 0 . t t 0 t 1 t 2 t 2 t n · · · Computational Strategies Krylov Subspace Methods for TEM 11 23/08/2007, Harrachov Martin Afanasjew et al., TU Freiberg

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

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 e j and magnetic fields h j in a leap-frog type iteration. e j − 1 h j − 1 e j h j e j +1 h j +1 e update: t h update: t j − 1 t j t j +1 Time-Stepping Methods Krylov Subspace Methods for TEM 13 23/08/2007, Harrachov Martin Afanasjew et al., TU Freiberg

Stability of the Du Fort-Frankel Method This method was shown to be stable if � µ min σ min t j ∆ t j = t j +1 − t j < ∆ x min 6 where ∆ x min 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 14 23/08/2007, Harrachov Martin Afanasjew et al., TU Freiberg

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

Rewriting the Problem (I) To apply Krylov subspace methods we have to rewrite our problem. Starting with (Maxwell) we get � 1 ∂ t e = − 1 � σ ∇× µ ∇× e . (PDE) This reduces to the solution of a linear first-order ordinary differential equation ∂ t e = A e , e ( t 0 ) = e 0 , (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 16 23/08/2007, Harrachov Martin Afanasjew et al., TU Freiberg

Rewriting the Problem (II) An explicit solution of (ODE) is given by e ( t ) = e ( t − t 0 ) A e 0 . 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 17 23/08/2007, Harrachov Martin Afanasjew et al., TU Freiberg

Arnoldi/Lanczos Procedure ◮ We define the m -th Krylov subspace as follows � � b , A b , A 2 b , . . . , A m − 1 b K m ( A , b ) := span . ◮ Generate an orthornormal basis V m = [ v 1 , v 2 , . . . , v m ] of K m ( A , b ) using a Gram-Schmidt procedure/three-term recurrence relation satisfying V ⊤ m AV m = H m with H m ∈ R m × m upper Hessenberg/tridiagonal. ◮ Calculate the Arnoldi approximation of order m f m := � b � V m f ( H m ) [1 , 0 , . . . , 0] ⊤ Krylov Subspace Methods Krylov Subspace Methods for TEM 18 23/08/2007, Harrachov Martin Afanasjew et al., TU Freiberg

Time-Stepping, Recycling, and Restarts (I) ◮ Remember: Given e 0 at time t 0 we are interested in evaluating the electric fields e 1 , e 2 , . . . , e n at times t 1 < t 2 < . . . < t n from an interval [ t 0 , t n ] . ◮ Time-stepped Arnoldi ◮ In each time step we compute f m j +1 ∈ K m ( A , f m j ) for f ( x ) = e ( t j +1 − t j ) x 0 = e 0 and m = m ( j ) ∼ � ( t j +1 − t j ) A � 1 / 2 . Such a with f m 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 19 23/08/2007, Harrachov Martin Afanasjew et al., TU Freiberg

Time-Stepping, Recycling, and Restarts (II) ◮ Arnoldi with Recycling ◮ Similarly to the time-stepped variant we compute in each step f m j ∈ K m ( A , e 0 ) for f ( x ) = e ( t j − t 0 ) x with m = m ( j ) ∼ � ( t j − t 0 ) A � 1 / 2 . This allows us to reuse basis vectors from the j − 1 -th step and only compute the additional basis vectors v m ( j )+1 , v m ( j )+2 , . . . , v m ( j +1) . ◮ Restarted Arnoldi Method ◮ If A cannot be symmetrized the unrestarted Arnoldi method might require a prohibitively high number of basis vectors. To overcome 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 20 23/08/2007, Harrachov Martin Afanasjew et al., TU Freiberg

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

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 � 10 − 6 , 10 − 3 � seconds. [ t 0 , t n ] = ◮ Everything implemented in pure MATLAB (Release 2007a). Numerical Examples Krylov Subspace Methods for TEM 22 23/08/2007, Harrachov Martin Afanasjew et al., TU Freiberg