Improved numerical methods for elliptic problems in astrophysics

• Autores: Jose Enrique Adsuara Fuster
• Directores de la Tesis: Miguel-Ángel Aloy Toras (dir. tes.) , Pablo Cerdá Durán (codir. tes.)
• Lectura: En la Universitat de València ( España ) en 2017
• Idioma: español
• Tribunal Calificador de la Tesis: Michael Dumbser (presid.) , Isabel Cordero Carrión (secret.) , Alicia Cordero Barbero (voc.)
• Enlaces
  • Tesis en acceso abierto en: TESEO
• Resumen

  • Elliptic partial differential equations appear in a wide variety of areas of mathematics, physics and engineering. They are of particular interest in Astrophysics and appear, e.g., when computing the gravitational potential, in the solution of the Grad-Shafranov equation for force-free magnetospheres, to impose divergence free constraints in the numerical integration of MHD equations or when solving the constraint equations in General Relativity. Typically, elliptic equations must be solved numerically, which sets an ever-growing demand for efficient and highly parallel algorithms to tackle their computational solution.

    The Scheduled Relaxation Jacobi is a promising class of methods, atypical for combining simplicity and efficiency, that has been recently introduced for solving linear Poisson-like elliptic equations.

    It is an extension of the classical Jacobi iterative method to solve linear systems of equations (Au=b). It also inherits its robustness. Its methodology relies on computing the appropriate parameters of a multilevel approach with the goal of minimizing the number of iterations needed to cut down the residuals below specified tolerances.

    The efficiency in the reduction of the residual increases with the number of levels employed in the algorithm. Applying the original methodology to compute the algorithm parameters with more than 5 levels notably hinders obtaining optimal schemes, as the mixed (non-linear) algebraic-differential system of equations from which they result become notably stiff.

    On one hand, we present a new methodology for obtaining the parameters of Scheduled Relaxation Jacobi schemes that overcomes the limitations of the original algorithm and provides parameters for these schemes with up to 15 levels and resolutions of up to 215 points per dimension, allowing for acceleration factors larger than several hundreds with respect to the Jacobi method for typical resolutions and, in some high resolution cases, close to 1000. Most of the success in finding these optimal schemes with more than 10 levels is based on an analytic reduction of the complexity of the previously mentioned system of equations. Furthermore, we extend the original algorithm to apply it to certain systems of non-linear elliptic equations.

    On the other hand, in a typical Scheduled Relaxation Jacobi scheme, the former set of factors is employed in cycles of M consecutive iterations until a prescribed tolerance is reached. We present the analytic form for the optimal set of relaxation factors for the case in which all of them are strictly different, and find that the resulting algorithm is equivalent to a non-stationary generalized Richardson's method where the matrix of the system of equations is preconditioned multiplying it by D=diag(A). Our method to estimate the weights has the advantage that the explicit computation of the maximum and minimum eigenvalues of the matrix A (or the corresponding iteration matrix of the underlying weighted Jacobi scheme) is replaced by the (much easier) calculation of the maximum and minimum frequencies derived from a von Neumann analysis of the continuous elliptic operator. This set of weights is also the optimal one for the general problem, resulting in the fastest convergence of all possible Scheduled Relaxation Jacobi schemes for a given grid structure. We refer to it as the Chebyshev-Jacobi method. The amplification factor of the method can be found analytically and allows for the exact estimation of the number of iterations needed to achieve a desired tolerance.

    We also show that with the set of weights computed for the optimal SRJ scheme for a fixed cycle size it is possible to estimate numerically the optimal value of the relaxation factor in the successive overrelaxation method in some cases.

    We demonstrate with practical examples, with application in Astrophysics, that our method also works very well for Poisson-like problems in which a high-order discretization of the Laplacian operator is employed (e.g., a 9- or 17-points discretization). This is of interest since the former discretizations do not yield consistently ordered A matrices and, hence, the theory of Young cannot be used to predict the optimal value of the SOR parameter. Furthermore, the optimal SRJ schemes deduced here are advantageous over existing SOR implementations for high-order discretizations of the Laplacian operator in as much as they do not need to resort to multi-coloring schemes for their parallel implementation.

    We present the implementation of the Chebyshev-Jacobi method using a purely MPI implementation, an openMP / MPI hybrid implementation over both shared memory machines and distributed memory machines. They show ideal speed up. We also show how to reach a remarkable speed up when solving elliptic partial differential equations with finite differences thanks to the joint use of the Chebyshev-Jacobi method with high order discretizations and its parallel implementation over GPUs.

    Finally, we have tried to apply our methods beyond the realm of Astrophysics with limited success though. Specially, we have addressed the problem of finding normal modes of human eyeballs. This problem is ready for being solved with an improved variant of the methodology here presented. The improvement consists on extending the calculation of the optimal set of parameters to non positive-definite matrices. Our ideas on how to proceed in this field are sketched in the outlook of this thesis.