Fast Sparse Spectral Methods for Higher Dimensional PDEs Jie Shen - PowerPoint PPT Presentation

Fast Sparse Spectral Methods for Higher Dimensional PDEs Jie Shen Purdue University Collaborators: Li-Lian Wang, Haijun Yu and Alexander Alekseenko Research supported by AFOSR and NSF ICERM workshop, June 3-7, 2013

Outlines Approximation by hyperbolic cross Sparse grids for bounded domains Sparse grids for unbounded domains Application to electronic Schr¨ odinger equation Adaptive sparse tensor products with multi-wavelets for the Boltzmann collision operator Concluding remarks

Motivations • Electronic Schr¨ odinger equation: find the eigenvalues and eigenfunctions of the Hamilton operator N N K N N � � � � � H = − 1 Z ν 1 ∆ i − | x i − a ν | + | x i − x j | , 2 i =1 i =1 ν =1 i =1 j > i where x 1 , · · · , x N ∈ R 3 are coordinates of N electrons. • Boltzmann Equation: ∂ f ∂ t + p m · ∇ x f = G ( f , f ) where x , p are position and momentum of particle in R 3 , m is the particle mass and f is the density distribution function.

• Fokker-Planck Navier-Stokes equations for FENE dumbbell model:  ∂ t u + ( u · ∇ x ) u + ∇ x p = ν Re∆ x u + 1 − ν   Re De ∇ x · τ p ,      ∇ x · u = 0 , ∂ t f + ∇ x · ( u f ) + ∇ q · ( ∇ x u · q f ) = 1 De ∇ q · ( F ( q ) f ) + 1 De∆ q f ,    �     τ p ( x , t ) = ( F ( q ) ⊗ q ) f ( x , q , t ) d q . W x ∈ Ω ⊂ R 3 , Ω is the fluid space; q ∈ W ⊂ R 3 , W is the configuration space; F ( q ) is the configuration force.

A model problem Consider the model elliptic equation: in Ω = ( − 1 , 1) d ; α u − ∆ u = f u | ∂ Ω = 0 . Given an approximation space V n and an interpolation operator I n , the (weighted) spectral-Galerkin method is to find u n ∈ V n such that ( ∇ u n , ∇ ( v n ω )) = ( I n f , v n ) ω , ∀ v n ∈ V n .

A model problem Consider the model elliptic equation: in Ω = ( − 1 , 1) d ; α u − ∆ u = f u | ∂ Ω = 0 . Given an approximation space V n and an interpolation operator I n , the (weighted) spectral-Galerkin method is to find u n ∈ V n such that ( ∇ u n , ∇ ( v n ω )) = ( I n f , v n ) ω , ∀ v n ∈ V n . Error estimate: �∇ ( u − u n ) � L 2 ω � inf v n ∈ V n �∇ ( u − v n ) � L 2 ω + � f − I n f � L 2 ω . approximation error interpolation error

Usual spectral method: Let V n = P n +1 ∩ H 1 0 (Ω) where P n +1 is the space of polynomials of degree ≤ n + 1 in EACH variable. N = dim( V n ) = n d ; ω � n 1 − s � u � H s ω � N (1 − s ) / d � u � H s inf v n ∈ V n �∇ ( u − v n ) � L 2 ω . The # of degree of freedom increases exponentially fast with d — the convergence rate deteriorates rapidly as d increases. This is the so called Curse of dimensionality!

Hyperbolic cross (Korobov, Babenko ’57) The hyperbolic cross is defined as � � d � k ∈ Z d K ( n ) = + : max( k i , 1) ≤ n . i =1 The cardinality of K ( n ) is of order n (log n ) d − 1 .

Hyperbolic cross (Korobov, Babenko ’57) The hyperbolic cross is defined as � � d � k ∈ Z d K ( n ) = + : max( k i , 1) ≤ n . i =1 The cardinality of K ( n ) is of order n (log n ) d − 1 . d d � � Set V n = span { φ k i ( x i ) : max( k i , 1) ≤ n } , i =1 i =1 in what functional spaces does V n provide good approximations?

Table : Comparison of grids and hyperbolic cross n d=2 d=3 d=4 |K ( n ) | |K ( n ) | Grids |K ( n ) | Grids |K ( n ) | Grids |K ( n ) | nlog ( n ) nlog ( n ) 2 9 100 42 2.1 1000 141 3.2 10000 424 13 196 64 1.9 2744 228 2.7 38416 720 17 324 87 1.8 5832 321 2.4 104976 1041 21 484 113 1.8 10648 435 2.2 234256 1457 25 676 138 1.7 17576 546 2.1 456976 1877 29 900 162 1.7 27000 646 2.0 810000 2229 33 1156 190 1.6 39304 778 1.9 1336336 2745 37 1444 217 1.6 54872 904 1.9 2085136 3239 41 1764 243 1.6 74088 1021 1.8 3111696 3683 45 2116 273 1.6 97336 1165 1.8 4477456 4243

Korobov spaces for periodic functions p (0 , 2 π ) d and its Fourier series Given f ∈ L 2 � � f ( k ) = 1 f ( k ) e i k · x ˆ with ˆ f ( x ) e − i k · x d x . f ( x ) = | Ω | Ω k ∈ Z d Let us define r ( k ) := � d j =1 max( | k j | , 1) , k = ( k 1 , . . . , k d ), then, the Korobov space of order α ∈ R is � � � f ( k ) | 2 r ( k ) 2 α < ∞ | ˆ K α f ∈ L 2 p := p (Ω) : k ∈ Z d � � ∂ q 1 + ... + q d f ∈ L 2 = p (Ω) : 0 ≤ q j ≤ α, 1 ≤ j ≤ d . ∂ x q 1 1 . . . ∂ x q d d

� f ( k ) e i k · x � f : f ( x ) = � r ( k ) ≤ n ˆ Setting X n := . Then, one can prove p � n r − s � v � K s v n ∈ X n � v − v n � K r inf p . (e.g. Griebel & Hamaekers 2007)

� f ( k ) e i k · x � f : f ( x ) = � r ( k ) ≤ n ˆ Setting X n := . Then, one can prove p � n r − s � v � K s v n ∈ X n � v − v n � K r inf p . (e.g. Griebel & Hamaekers 2007) Remarks: The estimate still depends on d since N = O ( n (log n ) d − 1 ). It is also possible to remove the dependence on d by using the optimized (energy based) hyperbolic cross space (Bungartz & Griebel ’04): � � d � k ∈ Z d max( k i , 1) | | k | − σ ∞ ≤ n 1 − σ K σ ( n ) = + : . i =1 Then, for σ ∈ (0 , 1), Card( K σ ( n )) = O ( n ). In this case, the estimate holds only for σ ≤ r s .

� f ( k ) e i k · x � f : f ( x ) = � r ( k ) ≤ n ˆ Setting X n := . Then, one can prove p � n r − s � v � K s v n ∈ X n � v − v n � K r inf p . (e.g. Griebel & Hamaekers 2007) Remarks: The estimate still depends on d since N = O ( n (log n ) d − 1 ). It is also possible to remove the dependence on d by using the optimized (energy based) hyperbolic cross space (Bungartz & Griebel ’04): � � d � k ∈ Z d max( k i , 1) | | k | − σ ∞ ≤ n 1 − σ K σ ( n ) = + : . i =1 Then, for σ ∈ (0 , 1), Card( K σ ( n )) = O ( n ). In this case, the estimate holds only for σ ≤ r s . These estimates were extended to the non-periodic case and unbounded domains in S. & Wang (SINUM ’10).

Korobov type spaces for non-periodic functions Let ω α,β = (1 − x ) α (1 + x ) β and J ( α,β ) be the Jacobi weight and n Jacobi polynomial with index α, β > − 1. We define the multi-dimensional Jacobi weights and Jacobi polynomials by α, ¯ α, ¯ i =1 ω α i ,β i ( x i ) , Φ ¯ β i =1 J α i ,β i ω ¯ β (¯ x ) = Π d x ) = Π d (¯ ( x i ) . ¯ k i k We then define (for any integer t ≥ 0) � β := { u : � u � 2 � ∂ ¯ r u � 2 K t β := r < ∞} . α, ¯ K t r , ¯ ¯ ω ¯ α +¯ β +¯ α, ¯ ¯ | ¯ r | ∞ ≤ t

Projection errors in the multi-D case α, ¯ α, ¯ Denote P ¯ β = span { Φ ¯ β x ) : Π d (¯ i =1 max( k i , 1) ≤ n } , and define n ¯ k α, ¯ α, ¯ Π ¯ β β → P ¯ β : L 2 be the orthogonal projector defined by n n α, ¯ ω ¯ α, ¯ α, ¯ ( u − Π ¯ β β = 0 , ∀ v n ∈ P ¯ β u , v K ) ω ¯ . α, ¯ n n Theorem. (S. & Wang ’10) For u ∈ K s β and 0 ≤ t < s , we have α, ¯ ¯ α, ¯ � u − Π ¯ β β � n t − s | u | K s u � K t β , n α, ¯ α, ¯ ¯ ¯ where � | u | 2 � ∂ ¯ r u � 2 β := r . K s r , ¯ ω ¯ α +¯ β +¯ α, ¯ ¯ | ¯ r | ∞ = s

Projection errors in the multi-D case α, ¯ α, ¯ Denote P ¯ β = span { Φ ¯ β x ) : Π d (¯ i =1 max( k i , 1) ≤ n } , and define n ¯ k α, ¯ α, ¯ Π ¯ β β → P ¯ β : L 2 be the orthogonal projector defined by n n α, ¯ ω ¯ α, ¯ α, ¯ ( u − Π ¯ β β = 0 , ∀ v n ∈ P ¯ β u , v K ) ω ¯ . α, ¯ n n Theorem. (S. & Wang ’10) For u ∈ K s β and 0 ≤ t < s , we have α, ¯ ¯ α, ¯ � u − Π ¯ β β � n t − s | u | K s u � K t β , n α, ¯ α, ¯ ¯ ¯ where � | u | 2 � ∂ ¯ r u � 2 β := r . K s r , ¯ ω ¯ α +¯ β +¯ α, ¯ ¯ | ¯ r | ∞ = s Approximation results in usual H 1 -norm and for elliptic equations are also established in S. & Wang ’10.

0 -5 -10 200 -15 400 -20 600 -25 -30 800 -35 1000 -40 0 200 400 600 800 1000 Figure : Spectrum of a highly anisotropic exact solution

Figure : Convergence history: comparison with the full grid

-5 10 -10 20 -15 30 -20 -25 40 -30 50 -35 60 -40 10 20 30 40 50 60 Figure : Spectrum of an isotropic exact solution

L2 error of solution d=2 full d=3 full 0.01 d=4 full d=5 full d=1 d=2 sparse 0.0001 d=3 sparse d=4 sparse d=5 sparse 1e-06 L2 error 1e-08 1e-10 1e-12 1e-14 1e-16 1 10 100 1000 10000 100000 1e+06 # of grid points Figure : Convergence history: comparison with the full grid

Sparse grids Sparse grid construction (Smolyak ’63): Start from a sequence of (preferably nested ) univariate quadrature/interpolation formulas: n i � f ( x i j ) ω i Q i f := j , ∆ i = Q i − Q i − 1 with Q 0 = 0 . j =1 Define a multidimensional quadrature rule by: � Q ( d ) f := (∆ i 1 ⊗ · · · ⊗ ∆ i d ) f . l d ≤| i |≤ l

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 Figure : Q (2) based on Chebyshev-Gauss-Lobatto quadrature. 5