330 likes | 355 Views
Parallel solution of the Helmholtz equation with high frequency. Dan Gordon Computer Science University of Haifa. Rachel Gordon Aerospace Eng. Technion. OUTLINE. The wave equation The Kaczmarz algorithm ( KACZ ) KACZ CARP ( C omponent- A veraged R ow P rojections)
E N D
Parallel solution of the Helmholtz equation with high frequency Dan Gordon Computer Science University of Haifa Rachel Gordon Aerospace Eng. Technion ISCM-29, Technion, Haifa, Israel
OUTLINE The wave equation The Kaczmarz algorithm (KACZ) KACZ CARP (Component-Averaged Row Projections) CARP-CG: CG acceleration of CARP Sample results with the Helmholtz equation ISCM-29, Technion, Haifa, Israel
The Helmholtz Equation speed: c frequency: ν wave length: l = c/ν wave number: k = 2p/l = 2pν/c wave eqn: -Δu - k2u = f Discretization with uniform grid size h No. of grid pts per l: Ng = l/h = 2p/kh Considered desirable: Ng ≥ 8, but Ng = 6 also gave good results Linear system is complex and strongly indefinite ISCM-29, Technion, Haifa, Israel
The Helmholtz Equation Challenging problem when ν is high Shifted Laplacian approach: Bayliss, Goldstein & Turkel, 1983 introduced a shift into the Laplacian Erlangga, Vuik & Oosterlee, 2004/06 complex shift: -Δu - (1 -i b) k2 u = f Uses multigrid to solve the preconditioner Not simple for irregular grids ISCM-29, Technion, Haifa, Israel
The Helmholtz Equation Bollhöfer, Grote & Schenk, 2009: Introduced algebraic multilevel preconditioner Use symmetric max weight matchings and inverse-based pivoting Applied to heterogeneous 2D and 3D domains Can be parallelized in principle Apologies to many others! ISCM-29, Technion, Haifa, Israel
The Kaczmarz algorithm (KACZ) initial point eq. 1 eq. 2 eq. 3 ISCM-29, Technion, Haifa, Israel
KACZ with Relaxation Parameter KACZ can be used with a relaxation parameter w w=1: project exactly on the hyperplane w<1: project in front of hyperplane w>1: project beyond the hyperplane Cyclic relaxation: eq. i is assigned a relaxation parameter wi ISCM-29, Technion, Haifa, Israel
Algebraic formulation of KACZ Given the system Ax = b Consider the "normal equations" system AATy = b, x = ATy Well-known fact: KACZ is SOR applied to the normal equations The relaxation parameter of KACZ is the usual relax. par. of SOR ISCM-29, Technion, Haifa, Israel
CARP: Component-Averaged Row Projections A block-parallel version of KACZ The equations are divided into blocks (not necessarily disjoint) A variable shared by 2 or more blocks is "cloned" into its neighboring blocks. For each block (in parallel) do KACZ iterations Every shared variable becomes the average of its values in the different blocks Repeat until convergence ISCM-29, Technion, Haifa, Israel
CARP-CG: CG acceleration of CARP CARP is KACZ in some superspace (with cyclic relaxation parameters) Björck & Elfving (BIT '79): developed CGMN, which is a (sequential)CG-acceleration of KACZ (double sweep, fixed relax. parameter) We extended this result to allow cyclic relaxation parameters Result: CARP-CG ISCM-29, Technion, Haifa, Israel
CARP-CG: Properties On one processor, CARP-CG is identical to CGMN Particularly useful on systems with large off-diagonal elements example: convection-dominated PDEs Discontinuous coefficients are handled without requiring domain decomposition (DD) ISCM-29, Technion, Haifa, Israel
Robustness of CARP-CG KACZ inherently normalizes the eqns After normalization, the diagonal elements of AAT are larger than the off-diagonal ones (in each row) This is not diagonal dominance, but it makes the normal eqns manageable Normalization was also found to be useful for discontinuous coefficients ISCM-29, Technion, Haifa, Israel
Application of CARP-CG to Helmholtz equation A fixed relaxation parameter of 1.5 was used in all cases Domain: mostly unit square or unit cube 2nd order central difference scheme ISCM-29, Technion, Haifa, Israel
Homogeneous 2D problem Based on Erlangga et al. '04, §6.2 Eqn: Δu + k2u = 0 Domain: unit square [0,1] [0,1] Dirichlet bndry cond. on one side, with a discontinuity at midpoint 1st-order absorbing bndry cond. on other sides (Sommerfeld radiation condition) Grid points per l: Ng = 6, 8, 10 No. of processors: 1 – 32 k = (75), 150, 300, 450, 600 ISCM-29, Technion, Haifa, Israel
Heterogeneous 2D problem • 3-layer heterogeneous problem • Based on Erlangga et al. '04, §6.3 • Everything is identical to previous problem • EXCEPT: k=600 k=450 k=300 ISCM-29, Technion, Haifa, Israel
The Marmousi Model • Well-known benchmark for solvers of the Helmholtz equation • 6000m x 1600m vertical slice of earth surface, disturbance at top center • Highly heterogeneous and irregular • Speed of sound: 1500 m/s to 4000 m/s • Tested on 12 node infiniband machine • Each node: 2 quad CPUs ISCM-29, Technion, Haifa, Israel
Time (s) for rel-res<10-7, Ng ≥ 7.5 ISCM-29, Technion, Haifa, Israel
3D heterogeneous problem • Domain: [0,1] [0,1] [0,1] • Divided into 3 layers with k=60, 72, 90 • Point source in middle of one side • Sommerfeld radiation condition on bndry • Also tested with k=60, 90, 145 on the infiniband machine ISCM-29, Technion, Haifa, Israel
Time (s) for 3D hom. & het. problems,4 conv. goals, on 1 xeon E5520 proc. ISCM-29, Technion, Haifa, Israel
Summary – serial and parallel machines • When the frequency increases: • Faster convergence (on a fixed grid) • Improved scalability (on a fixed grid) • Improved speedup (with a fixed Ng) • For fixed Ng: no. of iterations is linear in k • Homogeneous & heterogeneous problems • Simple to implement • Generally useful for various problems with large off-diagonal elements and discontinuous coefficients ISCM-29, Technion, Haifa, Israel
Other Potential Applications Higher order schemes for the Helmholtz equation (good initial results) Maxwell equations? Saddle-point problems? Circuit problems? Linear solvers in some eigenvalue methods? Suggestions are welcome! ISCM-29, Technion, Haifa, Israel
Relevant Publications http://cs.haifa.ac.il/~gordon/pub.html CARP: SIAM J Sci Comp 2005 CGMN: ACM Trans Math Software 2008 Microscopy: J Parallel & Distr Comp 2008 Large convection + discontin coef: CMES 2009 CARP-CG: Parallel Comp 2010 Scaling for discont coef: J Comp & Appl Math 2010 Helmholtz equation: tech rept http://cs.haifa.ac.il/~gordon/helm.pdf CARP-CG SOFTWARE AVAILABLE ON REQUEST THANK YOU! ISCM-29, Technion, Haifa, Israel