The solution of partial differential equations (PDEs) plays a fundamental role in numerous scientific and engineering applications. Helmholtz equations hold particular significance among the various types of PDEs due to their wide-ranging applicability in fields such as acoustics, electromagnetics, quantum mechanics, and seismology, which is the scientific study of earthquakes and elastic waves \cite{juraev2023applications}. However, obtaining accurate and efficient solutions for these equations can be challenging, especially when dealing with large-scale problems or complex geometries. Therefore, iterative Krylov subspace methods \cite{liesen2013krylov} have emerged as powerful techniques for solving linear systems arising from discretized PDEs. One of the simplest methods within this class is the \cg{} (CG) method \cite{shewchuk1994introduction}, renowned for its robustness and efficiency in the case of symmetric positive definite matrices. In the case of unsymmetric matrices, one of the fastest Krylov subspace methods is the GMRES Method (Generalized Minimal Residual Method) \cite{saad1993flexible}, which can also be used in the case of more complex preconditioning techniques, like the \ddm{} (DDM), that we will describe later. 

The \cg{} method is an iterative algorithm to solve symmetric positive definite systems of linear equations. It is beneficial for solving large, sparse systems, which arise in many scientific and engineering applications. Some properties of the CG method are given below.
\begin{itemize}
    \item \textbf{Efficiency} -  The \cg{} method is efficient, and it usually requires fewer iterations than other iterative methods to achieve a given level of accuracy. 
    \item \textbf{Symmetry} - The \cg{} method is suitable for solving symmetric positive definite systems. The symmetric positive definite systems appear in many scientific and engineering applications. In such cases, the \cg{} method is often faster and more accurate compared to other methods not specifically designed for symmetric systems.
    \item \textbf{Parallelizability} - The \cg{} method is easily parallelizable, so it can leverage modern computer architectures and solve large systems efficiently. This parallelizability makes it a good candidate for GPGPU use.
    \item \textbf{Low memory footprint} - The \cg{} method requires less memory than other methods. That results in CG being ideal for solving large systems on computational devices with limited memory capacity that cannot store the whole system in memory simultaneously.
\end{itemize}

When discussing iterative methods for solving linear systems of equations, it is critical to consider using preconditioners. The preconditioners are used in iterative methods for solving linear systems of equations to improve their convergence properties. They are generally used when the coefficient matrix is large, sparse, and poorly conditioned. A preconditioner is a matrix or an operator that approximates the inverse of the coefficient matrix and is applied to the system of equations before or during the usage of the iterative method. The goal of the preconditioner is to transform the system of equations into an equivalent system that is easier to solve iteratively by obtaining a new matrix with a lower condition number. 

For instance, the Additive Schwarz Method \cite{dryja1989additive} can be a preconditioner for the Krylov subspace methods like GMRES (Generalized Minimal Residual) to improve their convergence properties. In this approach, the domain of the problem is divided into a set of overlapping subdomains. Each subdomain is solved independently with a direct or an inexact solver like ILU (Incomplete Cholesky Factorisation) or CG. At the same time, the outer iterative method is the Flexible GMRES (FGMRES) method, allowing the robust convergence of outer iterations with different preconditions on each iteration, see \cite{graham2017domain}. During the preconditioning step, the solutions obtained from each subdomain are combined using the Additive Schwarz method, which updates the global solution by adding up the local (possibly approximate) solutions from each subdomain. The outer iterative process is repeated until the global solution converges to the desired level of accuracy. The key advantage of this approach is that the computational cost of the local preconditioner is significantly reduced compared to solving the entire system using GMRES with direct subdomain solvers with sparse matrix factorization.

In this thesis, the background section (see Section \ref{background}) explores GPU programming and its benefits for computational tasks. It also provides an overview of Krylov subspace and domain decomposition methods. The data and methods section (see Section \ref{methods}) details the implementation of GPU-accelerated solvers for the Helmholtz equations, while the results section (see Section \ref{results}) analyzes the solver's performance across various metrics. Finally, the discussion section (see Section \ref{discussion}) examines the results' interpretation, discusses the methods' limitations, and suggests potential improvements.

While writing this thesis, the author used ChatGPT \cite{chatgpt} to get feedback on the content and correct language errors. The feedback helped improve the text and fix the mistakes regarding the language.




\subsection{Problem}
The solution of Helmholtz equations presents numerous challenges, making it a complex and demanding task. These equations describe the behavior of wave phenomena, such as acoustic or electromagnetic waves, as they propagate through a medium. Wave propagation brings unique obstacles originating from diffraction, interference, and resonance. Accurate wave behavior representation necessitates meticulous treatment of boundary conditions, computational domain size, and discretization techniques.

The Helmholtz operator can also exhibit severe ill-conditioning, mainly when dealing with large wave numbers or complex geometries. Ill-conditioning refers to situations where slight variations in the input or data can lead to significant changes in the solution. Furthermore, Helmholtz equations arise in various applications involving large-scale systems, such as wave propagation in realistic environments or the modeling of complex structures. Solving such large-scale problems calls for efficient computational techniques capable of managing vast amounts of data and computational operations while maintaining accuracy and reasonable time complexity.

Addressing these challenges requires developing and utilizing advanced numerical techniques, including iterative Krylov subspace methods like the Flexible Generalized Minimal Residual method (FGMRES) and the \cg{} (CG) method.  These methods aim to overcome the limitations presented by direct solvers by efficiently approximating the solution through iterative means, making them well-suited for handling the complexities associated with Helmholtz equations.

\subsection{Contribution}
In this thesis, we study the solution to Helmholtz's problem using the FGMRES method as an outer Krylov subspace method preconditioned with Restricted Averaging Additive Schwarz (RAAS) to improve the convergence of FGMRES. The \cg{} (CG) is used as an inner method for equal-sized submatrices obtained by applying an RAAS preconditioner to the original domain. This makes it possible to solve for multiple subdomains on a single GPU simultaneously – using multiple right-hand side vectors since the matrices are precisely the same on each subdomain. The aim of the thesis is to shed light on the properties, strengths, and limitations of those methods.

In order to enhance the computational speed of these operations, advanced techniques are employed, such as utilizing the power of Graphical Processing Units (GPUs) to accelerate matrix operations using the OpenCL programming framework.

The choice of an appropriate preconditioner can significantly impact the performance of the CG method when solving Helmholtz equations. Preconditioners aim to transform the original system into a more well-conditioned form, thereby reducing the number of iterations required for convergence. By exploring and comparing different preconditioning strategies, this thesis aims to identify the most effective and efficient preconditioners tailored specifically for Helmholtz equations.

Furthermore, numerical experiments will be conducted to evaluate the performance of different solutions in various scenarios using the resources from UT HPC. This will involve systematically varying parameters such as the total number of domains, domain size, and number of inner iterations (and possibly the presence of heterogeneities - as future work) to assess the preconditioners' robustness and scalability. Comparisons will be made based on convergence behavior, computational efficiency, and overall solution accuracy. The outcomes of this investigation will provide the performance gains achieved by utilizing the parallel processing capabilities of GPUs and explore the integration of \ddm{s} to enhance the solution process's efficiency further. The objective is to achieve significant speedup and reduce the time-to-solution for large-scale Helmholtz equation problems.