헬름홀츠 방정식을 풀기 위한 개선된 shifted Laplace preconditioner

Abstract

헬름홀츠 방정식은 주파수 성분의 파동장을 지배하는 편미분 방정식(PDE)이다. 이는 타원형 PDE이기 때문에, 소스 정보를 포함한 주어진 우변에 대해 선형 시스템, 즉 임피던스 매트릭스를 계산함으로써 해을 얻을 수 있다. 헬름홀츠 방정식의 상호상관성(reciorocity) 때문에, 이산화된 임피던스 행렬은 복소 대칭 행렬로 자연스럽게 표현될 수 있으며, 이를 통해 Cholesky 분해를 적용하여 메모리 요구량과 산술 연산의 수를 절반으로 줄일 수 있다. 하지만 PML 경계조건을 사용하면서 행렬의 대칭 특성을 보존하는 방법은 널리 사용되지 않았다. 따라서 본 연구에서는 PML 경계조건 층을 포함한 전체 공간 노드에 대한 대칭 행렬로서 헬름홀츠 방정식을 이산화하는 방법을 제안한다. 또한 파장당 격자 수가 4보다 큰 범위 내 에서 0.5% 미만의 분산 오차를 만족하는 최적화된 이산화 계수를 제안한다. 선형 시스템은 직접법을 사용하여 2차원 혹은 작은 3차원 문제를 풀 수 있지만, 복잡도를 고려하면 일반적인 3차원 문제를 풀기에는 시간이 소요량, 행렬의 분해시 저장공간 요구량 그리고 임피던스 행렬을 계산하는데 극히 높은 비용이 필요하게 된다. 이 경우, 반복법, 특히, conjugate gradient(CG), GMRES 등과 같은 크릴로브 부분공간법이 효과적으로 사용될 수 있다. 본 연구에서는, 여러 크릴로 브 부분공간법들을 대칭 행렬에 적용해 실험적으로 어떤 방법이 임피던스 행렬을 계산하는데 적합한지 확인하였다. 이를 통해, 다양한 kh 조건에 대해 conjugate residual(CR)이 가장 빠르게 해을 제공하는 것을 관찰하였다. 수렴 속도의를 개선 하기 위해, preconditioner 행렬의 복소 시프트(ϵ)를 조절하여 최적의 고유값 군집화를 보여주는 Shifted Laplace preconditioner가 사용되었다. 일반적으로 이러한 preconditioner를 적용하기 위해 불완전 분해가 주로 사용된다. 본 연구에서는, 분해된 행렬의 p개의 가장 큰 값만 남기는 ICT(p), 불완전 Cholesky 분해를 사용하는 것을 제안한다. 주어진 kh 조건에 대해 최적의 ϵ과 p를 제안하고, 이를 인공합성 속도모델을 사용한 수치 실험을 통해 검증한다. 여기서 직접법과 비교하여 파동장의 해를 구하는데 적은 시간과 무시할 수 있는 오차를 보여주었다. 특히, 169X169X50 크기의 3차 문제의 경우, 고도로 최적화된 프론탈 솔버인 UMFPACK의 연산 시 간의 약 5%인 것으로 나타났다.

The Helmholtz equation is a partial differential equation (PDE) that governs the wavefield of a frequency component. Since it is an elliptic type PDE, it is expressed in the form that the solution can be obtained by solving a linear system, which is composed of the impedance matrix and a given right hand side including source information. Due to the reciprocity property of the Helmholtz equation, the discretized impedance matrix can be naturally expressed as a complex symmetric matrix that can reduce a memory requirement and the number of arithmetic operations by half by applying the Cholesky factorization. However, it is not widely disseminated how to reserve such symmetric property of the matrix using the PML boundary condition. Therefore, firstly, we propose a method that discretizes the Helmholtz equation as a symmetric matrix for the entire spatial nodes including the PML layer. We also suggest the optimized discrete coefficient to generate impedance matrix which satisfies the dispersion error of less than 0.5\% within the range that the number of grids per wavelength is greater than 4. The system can be solved by using the direct matrix solver for 2D or small 3D problems, however, time taken and memory requirement to factorize and solve the impedance matrix is extremely high to solve general 3D problems considering the complexities. In this case, iterative methods, in particular, the Krylov subspace methods, such as conjugate gradient method, general minimal residual method, etc., can be effectively applied. In this study, many Krylov iterative solvers that can be applied to the symmetric matrix are implemented to experimentally investigate which method is advantageous for solving the impedance matrix. As a result, it was observed that for various kh conditions, the conjugate residual method solves the matrix in the fastest time consistently. To shorten the convergence speed, we applied the shifted Laplacian preconditioner that can optimally modify the eigenvalue clustering by adjusting the complex shift(ϵ) of the preconditioner matrix. To apply this preconditioner, incomplete factorization method is generally used. In this study, we suggest incomplete Cholesky decomposition, especially ICT(p) method, that leaves only p largest values of the factorized matrix. Optimal pairs of ϵ and p for given kh condition are also suggested, which is validated by numerical examples using the synthetic velocity models. It is shown that only few minutes are taken to solve the problem to yield the wavefield that has negligibly small error compared to that obtained by the direct method, the error is negligibly small. Especially, in the case of a 3D problem with a size of 169X169X50, it is shown that it is about 5\% of the computation time of UMFPACK, a highly optimized directive multifrontal solver.