Model reduction of wave propagation via phase-preconditioned rational Krylov subspaces Delft University of Technology † and Schlumberger ∗ V. Druskin ∗ , R. Remis † , M. Zaslavsky ∗ , J¨ orn Zimmerling † November 8, 2017 J¨ orn Zimmerling (TU Delft) Model reduction of wave propagation November 8, 2017 1 / 37
Motivation • Second order wave equation with Receiver Source wave operator A PML 5500 Au − s 2 u = − δ ( x − x S ) 5000 500 4500 Wavespeed [m/s] 1000 • Assume N grid steps in every x-direction [m] 4000 spatial direction 1500 3500 • Scaling of surface seismic in 3D: 2000 3000 • # Grid points O ( N 3 ) 2500 2500 • # Sources O ( N 2 ) 2000 3000 1500 • # Frequencies O ( N ) 500 1000 1500 y-direction [m] • Overall O ( N 6 ) ψ ( N 3 ) (a) Section of wave speed profile of the Marmousi model. J¨ orn Zimmerling (TU Delft) Model reduction of wave propagation November 8, 2017 2 / 37
Goal of this work • Simulate and compress large scale wave fields in modern high performance computing environment ( parallel CPU and GPU environment ) • Use projection based model order reduction to • Approximate transfer function • reduce # of frequencies needed to solve • reduce # of sources to be considered • reduce # number of grid points needed J¨ orn Zimmerling (TU Delft) Model reduction of wave propagation November 8, 2017 3 / 37
Introduction • Simulating and compressing large scale wave fields Au [ l ] − s 2 u [ l ] = − δ ( x − x [ l ] S ) , (1) • With the wave operator given by A ≡ ν 2 ∆ , Laplace frequency s • We consider a Multiple-Input Multiple-Output problem • Define fields U = [ u [ 1 ] , u [ 2 ] , . . . , u [ N s ] ] and sources B = [ − δ ( x − x [ 1 ] S ) , − δ ( x − x [ 2 ] S ) , . . . , − δ ( x − x [ N s ] )] S • We are interested in the transfer function (Receivers and Sources coincide) � B H U ( s ) dx F( s ) = (2) • Open Domains J¨ orn Zimmerling (TU Delft) Model reduction of wave propagation November 8, 2017 4 / 37
Problem Formulation • After finite difference discretization with PML (A( s ) − s 2 I)U = ˆ B • Step sizes inside the PML h i = α i + β i s • Frequency dependent A( s ) caused by absorbing boundary Q( s )U = B with Q( s ) ∈ C N × N • Q( s ) propterties • sparse • complex symmetric (reciprocity) • Schwarz reflection principle Q( s ) = Q( s ) • passive (nonlinear NR 1 Re < 0) 1 NR: W { A ( s ) } = s ∈ C : x H A ( s ) x = 0 ∀ x ∈ C k \ 0 � � J¨ orn Zimmerling (TU Delft) Model reduction of wave propagation November 8, 2017 5 / 37
Problem Formulation • Transfer function from sources to receivers F(B , s ) = B H Q( s ) − 1 B • Reduced order modeling of transfer function over frequency range F m (B , s ) = V m B H (V H m Q( s )V m ) − 1 V H m B with V m ∈ C N × m • valid in a range of s , and easy to store Motivation • FD grid overdiscretized w.r.t. Nyquist • PML introduces losses • approximation F(B , s ) to noise level • limited I/O map J¨ orn Zimmerling (TU Delft) Model reduction of wave propagation November 8, 2017 6 / 37
Outline 1 Problem formulation 2 Model order reduction – Rational Krylov subspaces 3 Phase-Preconditioning 4 Numerical Experiments 5 Conclusions J¨ orn Zimmerling (TU Delft) Model reduction of wave propagation November 8, 2017 7 / 37
Structure preserving rational Krylov subspaces • Preserve: symmetry (w.r.t. matrix, frequency), passivity • Block rational Krylov subspace approach κ = s 1 , . . . , s m K m ( κ ) = span Q( s 1 ) − 1 B , . . . , Q( s m ) − 1 B � � K 2m Re = span {K m ( κ ) , K m ( κ ) } • Let V m be a (real) basis for K m Re then with reduced order model (via Galerkin condition) R m ( s ) = V H m Q( s )V m we approximate F m ( s ) = (V H m B) H R m ( s ) − 1 V H m B , J¨ orn Zimmerling (TU Delft) Model reduction of wave propagation November 8, 2017 8 / 37
Structure preserving rational Krylov subspaces • Numerical Range: W { R m ( s ) } ⊆ W { Q ( s ) } Proof: x H m R m ( s )x m = ( V m x m ) H Q ( s )( V m x m ) ⇒ R m ( s ) is passive • F m ( s ) is Hermite interpolant of F( s ) on κ ∪ κ Q( κ ) − 1 B ∈ K 2m Re + uniqueness of Galerkin • Schwarz reflection and symmetry hold aswell J¨ orn Zimmerling (TU Delft) Model reduction of wave propagation November 8, 2017 9 / 37
RKS for a Resonant cavity FDFD Response 6 RKS Response m=60 RKS Response m=20 5 1 4 Receiver Response [a.u] 20 Source PML 0.8 3 40 x-direction 0.6 60 2 80 0.4 1 100 0.2 120 0 0 0.1 0.2 0.3 0.4 0.5 0 20 40 60 80 100 120 Normalized Frequency y-direction • RKS has excellent convergence if singular Hankel values of system decay fast (few contributing eigenvectors) • F m ( s ) is a [2 m − 1 / 2 m ] rational function J¨ orn Zimmerling (TU Delft) Model reduction of wave propagation November 8, 2017 10 / 37
Problem of RKS in Geophysics - Nyquist Limit 0.2 Response [a.u.] 0.1 0 -0.1 -0.2 -0.3 0.05 0.1 0.15 Normalized Frequency • Long travel times ∗ δ ( t − T ) F − → exp ( − sT ) • Nyquist sampling of ∆ ω = π T � B H Q ( s ) − 1 Bdx is oscillatory • F( s ) = J¨ orn Zimmerling (TU Delft) Model reduction of wave propagation November 8, 2017 11 / 37
Filion Quadrature • Filion quadrature deals with oscillatory integral � F ( s ) = exp ( s t ) f ( t ) d t quadrature requires s ∆ t ≪ 1 � F ( s ) ≈ ∆ t a n exp ( s n ∆ t ) f ( n ∆ t ) n • Filion quadrature makes a n function of s ∆ t � F ( s ) ≈ ∆ t a n ( s ∆ t ) exp ( s n ∆ t ) f ( n ∆ t ) n • ⇒ Make projection basis s dependent • ⇒ Frequency dependence from asymptotic s → i ∞ (WKB) J¨ orn Zimmerling (TU Delft) Model reduction of wave propagation November 8, 2017 12 / 37
Phase-Preconditioning I - 1D • We can overcome the Nyquist demand by splitting the wavefield into oscillatory and smooth part u ( s j ) = exp ( − s j T eik ) c out ( s j ) + exp ( s j T eik ) c in ( s j ) . (3) • Oscillatory phase term obtained from high frequency asymptotics • Eikonal equation |∇ T [ l ] eik | 2 = 1 ν 2 • Amplitudes c out/in are smooth • Motivated by Filon quadrature • Handle oscillatory part analytically • Quadrature with smooth amplitudes • Note: Splitting not unique J¨ orn Zimmerling (TU Delft) Model reduction of wave propagation November 8, 2017 13 / 37
Phase-Preconditioning II • Projection on frequency dependent Reduced Order Basis K 2m EIK ( κ, s ) = span { exp ( − sT eik ) c out ( s 1 ) , . . . , exp ( − sT eik ) c out ( s m ) , exp ( sT eik ) c in ( s s ) , . . . , exp ( sT eik ) c in ( s m ) } • Preserve Schwartz reflection principle K 4m K 2m EIK ( κ, s ) , K 2m � � EIK;R ( κ, s ) = span EIK ( κ, s ) (4) • equivalent to changing Operator • Coefficients from Galerkin condition m m � � u m ( s ) = α i ( s ) exp ( − s T eik ) c out ( s i )+ β i ( s ) exp ( s T eik ) c in ( s i )+ . . . i =1 i =1 J¨ orn Zimmerling (TU Delft) Model reduction of wave propagation November 8, 2017 14 / 37
Phase-Preconditioning III • Non-uniqueness of splitting resolved by one-way WEQ c out ( s j ) = ν � s j ∂ � exp ( s j T eik ) ν u ( s j ) − ∂ | x − x S | u ( s j ) , (5) 2s j � s j � c in ( s j ) = ν ∂ exp ( − s j T eik ) ν u ( s j ) + ∂ | x − x S | u ( s j ) . (6) 2s j Effects of Phase preconditioning on • Number of Interpolation points • Size of the computational Grid • MIMO problems • Computational Scheme J¨ orn Zimmerling (TU Delft) Model reduction of wave propagation November 8, 2017 15 / 37
Projection on frequency dependent space • Let V m ;EIK ( s ) be a real basis of K 4 m EIK;R ( κ, s ) • The reduced order model is given by R m ;EIK ( s ) = V H m ;EIK ( s ) Q ( s ) V m ;EIK ( s ) . • large inner products on GPU • This preserves • symmetry • Schwarz reflection principle • passivity • Interpolation J¨ orn Zimmerling (TU Delft) Model reduction of wave propagation November 8, 2017 16 / 37
Number of interpolation points needed • Double interpolation of transfer function still holds F m ( s ) = F m ( s ) and d ds F m ( s ) = d ds F( s ) with s ∈ κ ∪ κ. (7) • Number of interpolation point needed dependent on complexity medium , not latest arrival • Proposition: Let a 1D problem have ℓ homogenous layers . Then there exist m ≤ ℓ + 1 non-coinciding interpolation points, such that the solution u m ;EIK ( s ) = u . J¨ orn Zimmerling (TU Delft) Model reduction of wave propagation November 8, 2017 17 / 37
Illustration of Proposition Wavespeed 1 0.5 0 Source L1 L2 L3 L4 L5 L6 L7 Real part field 10 0 -10 Source L1 L2 L3 L4 L5 L6 L7 Imag(C) outgooing 10 5 0 Source L1 L2 L3 L4 L5 L6 L7 Imag(C) incoming 4 2 0 Source L1 L2 L3 L4 L5 L6 L7 • Amplitudes are constants in layers + left and right of source • Basis is complete after ℓ + 1 iterations J¨ orn Zimmerling (TU Delft) Model reduction of wave propagation November 8, 2017 18 / 37
Recommend
More recommend
