Solving the Helmholtz equation via row- projections Tristan ¡van ¡Leeuwen ¡(Univ. ¡of ¡BC, ¡Canada) Dan ¡Gordon ¡(Univ. ¡of ¡Haifa, ¡Israel) Rachel ¡Gordon ¡(Technion, ¡Israel) Felix ¡J. ¡Herrmann ¡(Univ. ¡of ¡BC, ¡Canada)
Modelling ¡engine ¡for ¡ 3D ¡Frequency-‑domain ¡FWI: • work ¡with ¡few ¡sources/ frequencies ¡at ¡each ¡iteration • flexibility ¡in ¡type ¡of ¡equation • robust ¡ • parallel
3D ¡Helmholtz ¡equation: • large, ¡sparse, ¡indefinite ¡system • direct ¡factorization ¡not ¡feasible • `standard’ ¡preconditioners ¡often ¡ fail • successful ¡preconditioners ¡often ¡ tailored ¡to ¡specific ¡wave ¡equation [Ernst ¡& ¡MarLn, ¡2011]
VS. simple, ¡robust, ¡... fast, ¡complicated,..
Overview • Kaczmarz ¡preconditioning • Examples • Parallelization • 3D ¡Benchmark • Inversion • Conclusions
Kaczmarz The ¡Kaczmarz ¡method ¡solves ¡a ¡ system ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡by ¡successive ¡row ¡ A x = b projections λ i b i − a T � � a i , x := x + i x || a i || 2 2 with ¡relaxation ¡parameter 0 < λ i < 2 [Kaczmarz, ¡’37]
Kaczmarz λ < 1 λ = 1 λ > 1
Kaczmarz rewrite: ✓ ◆ λ i λ i a i a T I − b i a i x := x + i || a i || 2 || a i || 2 2 2 | {z } Q i a ¡double ¡sweep ¡yields x := ( Q 1 Q 2 . . . Q n Q n . . . Q 1 ) x + ( . . . ) b | {z } |{z} Q R
Kaczmarz Find ¡a ¡fixed ¡point ¡by ¡solving ( I − Q ) x = R b where ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡is ¡symmetric ¡and ¡ I − Q positive ¡semidefinite, ¡so ¡we ¡can ¡ use ¡CG ¡(CGMN). [Bjork ¡& ¡Elfving, ¡’79]
Kaczmarz We ¡never ¡form ¡the ¡matrix ¡ explicitly, ¡but ¡compute ¡its ¡action: Algorithm 1 DKSWP ( A, x , b , λ ) = Q x + R b forward sweep for i = 1 to n do x := x + λ ( b i − a T i x ) a i / || a i || 2 2 end for backward sweep for i = n to 1 do x := x + λ ( b i − a T i x ) a i / || a i || 2 2 end for return x
Kaczmarz • low ¡complexity ¡ • low ¡memory ¡(same ¡as ¡original ¡matrix) • no ¡setup ¡time
1D results 1D ¡profile, ¡varying ¡k, ¡10 ¡p/wavelength 1.5 0.01 1.4 0.005 1.3 1.2 p(x) 0 u 1.1 1 − 0.005 0.9 0.8 − 0.01 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 x x
1D results eigenvalues # ¡of ¡CG ¡iteraLons 1 2000 A T A AA* WAA*W* A T W T WA I − Q w=0.5 0.8 I − Q w=0.5 1500 I − Q w=1 I − Q w=1 # of CG iterations I − Q w=1.5 I − Q w=1.5 0.6 1000 � i 0.4 500 0.2 0 0 50 100 150 200 500 1000 1500 2000 2500 k 0 i
2D results Marmousi, ¡304 ¡x ¡1100, ¡f=20, ¡h=10 • CG ¡+ ¡Kaczmarz ¡(CGMN) • BiCGstab ¡+ ¡ILU(0) • SQMR ¡+ ¡ML [Bollhofer ¡et ¡al, ¡’08]
2D results ¡solve ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡starting ¡from ¡random ¡vector A r = 0 r 1 = ( I − M − 1 A ) r 0 r 1 = Q r 0 2 − 0.05 − 0.05 − 0.05 1.5 k z [1/m] 0 0 0 1 0.5 0.05 0.05 0.05 0 − 0.05 0 0.05 − 0.05 0 0.05 − 0.05 0 0.05 k x [1/m] k x [1/m] k x [1/m] ILU(0) ML Kaczmarz
2D results iteraLons Lme ¡[s]* CG ¡+ ¡Kaczmarz 5542 603 BiCGstab ¡+ ¡ILU(0) div. div. SQMR ¡+ ¡ML 514 379 0 1000 z [m] 2000 3000 0 2000 4000 6000 8000 10000 x [m]
Parallelization • divide ¡domain ¡in ¡blocks • Kaczmarz ¡sweeps ¡on ¡blocks ¡are ¡ done ¡in ¡parallel ¡(CARP) • average ¡boundary ¡elements ¡ between ¡each ¡sweep • convergence ¡guaranteed [Gordon ¡& ¡Gordon, ¡’10]
SEG/EAGE salt 7-‑point ¡stencil, ¡ABC
SEG/EAGE salt 0 10 -1 10 h=160; f= 1.25 f h iterations -2 h= 80; f= 2.5 10 h= 40; f= 5. h= 20; f=10. 1.25 160 310 -3 10 Res/Res(0) 2.5 80 510 -4 10 5 40 760 -5 10 10 20 1780 -6 10 on ¡1 ¡processor, ¡ ✏ = 10 − 4 -7 10 -8 10 0 1000 2000 3000 4000 5000 6000 7000 8000 iter on ¡12 ¡processors
SEG/EAGE salt grid: ¡105 ¡x ¡338 ¡x ¡338, ¡h=40, ¡f=5, ✏ = 10 − 4 np iter time (s) e ffi ciency 1 621 4444.90 1.00 2 619 3091.10 0.72 4 593 1335.00 0.83 8 599 737.90 0.75
Overthrust 27 ¡point ¡stencil ¡(2 nd ¡order), ¡PML
Overthrust 0 10 f=1.5, ¡h=200 f=1.5 Hz. f=3 Hz. f=3, ¡ ¡ ¡ ¡h=100 f=6 Hz. f=6, ¡ ¡ ¡ ¡h=50 − 1 10 rel. residual − 2 10 − 3 10 − 4 10 100 200 300 400 500 600 iteration on ¡1 ¡processor
Overthrust grid: ¡47x201x201, ¡h=100, ¡f=3 ¡Hz, ¡ ✏ = 10 − 4 e ffi ciency np iter time 1 659 20785.40 1.00 2 657 11306.90 0.92 4 596 4882.50 0.96 8 603 3960.10 0.60
Inversion Camembert ¡model ¡ ¡in ¡2D ¡... 3500 velocity [m/s] 3000 2500 1000 1000 500 500 0 0 z [m] x [m]
Inversion EDAM ¡model ¡ ¡in ¡3D ¡! 0 z [km] 500 1000 1000 1000 500 500 0 0 y [km] x [km]
Inversion transmission ¡setup, ¡9 ¡sources, ¡ 3 ¡frequencies ¡1.0 0 10 � = 10 − 2 0 � = 10 − 3 ¡0.9 rel. model error z [km] 500 ¡0.8 1000 1000 1000 ¡0.7 500 500 0 2 4 6 8 10 iteration 0 0 y [km] x [km]
Conclusions • simple, ¡robust ¡and ¡generic ¡ preconditioner • no ¡overhead, ¡cheap ¡to ¡apply • easy ¡to ¡parallelize
Future plans • efficient ¡ ¡implementation ¡of ¡ sweeps ¡using ¡multi-‑threading • investigate ¡BlockCG • incorporate ¡in ¡inversion • high-‑order ¡schemes
http://cs.haifa.ac.il/~gordon/soft.html
Acknowledgements This ¡work ¡was ¡in ¡part ¡financially ¡supported ¡by ¡the ¡Natural ¡Sciences ¡and ¡Engineering ¡ Research ¡Council ¡of ¡Canada ¡Discovery ¡Grant ¡(22R81254) ¡and ¡the ¡Collaborative ¡ Research ¡and ¡Development ¡Grant ¡DNOISE ¡II ¡(375142-‑08). ¡This ¡research ¡was ¡carried ¡ out ¡as ¡part ¡of ¡the ¡SINBAD ¡II ¡project ¡with ¡support ¡from ¡the ¡following ¡organizations: ¡ BG ¡Group, ¡BP, ¡Chevron, ¡ConocoPhillips, ¡Petrobras, ¡Total ¡SA, ¡and ¡WesternGeco.
References Bollhöfer, ¡M., ¡Grote, ¡M. ¡J., ¡& ¡Schenk, ¡O. ¡(2009). ¡Algebraic ¡MulLlevel ¡PrecondiLoner ¡for ¡ the ¡Helmholtz ¡EquaLon ¡in ¡Heterogeneous ¡Media. ¡ SIAM ¡J. ¡on ¡Sc. ¡Comp. , ¡ 31 (5), ¡3781. ¡ Gordon, ¡D. ¡and ¡Gordon, ¡R. ¡(2010). ¡CARP-‑CG: ¡A ¡robust ¡and ¡efficient ¡parallel ¡solver ¡for ¡ linear ¡systems, ¡applied ¡to ¡strongly ¡convecLon ¡dominated ¡PDEs. ¡ Parallel ¡CompuTng, ¡ 36 , ¡495–515. Björck, ¡Å. ¡and ¡Elfving, ¡T. ¡(1979) ¡Accelerated ¡projecLon ¡methods ¡for ¡compuLng ¡ pseudoinverse ¡soluLons ¡of ¡systems ¡of ¡linear ¡equaLons. ¡ BIT, ¡19 , ¡145–163. Kaczmarz, ¡S. ¡[1937] ¡Angenäherte ¡auflösung ¡von ¡systemen ¡linearer ¡gleichungen. ¡ BulleTn ¡InternaTonal ¡de ¡l’Académie ¡Polonaise ¡des ¡Sciences ¡et ¡des ¡LeYres, ¡35 , ¡355–357. Ernst, ¡O.G. ¡and ¡MarLn ¡J. ¡(2011) ¡Why ¡it ¡is ¡Difficult ¡to ¡Solve ¡Helmholtz ¡Problems ¡with ¡ Classical ¡IteraLve ¡Methods. ¡ Unpublished.
Recommend
More recommend
