Beyond spatial scalability limitations with a massively parallel method for linear oscillatory problems

This is the author accepted manuscript. The final version is available from SAGE Publications via the DOI in this record.This paper presents, discusses and analyses a massively parallel-in-time solver for linear oscillatory PDEs, which is a key numerical component for evolving weather, ocean, climate and seismic models. The time parallelization in this solver allows us to significantly exceed the computing resources used by parallelization-in-space methods and results in a correspondingly significantly reduced wall-clock time. One of the major difficulties of achieving Exascale performance for weather prediction is that the strong scaling limit – the parallel performance for a fixed problem size with an increasing number of processors – saturates. A main avenue to circumvent this problem is to introduce new numerical techniques that take advantage of time parallelism. In this paper we use a time-parallel approximation that retains the frequency information of oscillatory problems. This approximation is based on (a) reformulating the original problem into a large set of independent terms and (b) solving each of these terms independently of each other which can now be accomplished on a large number of HPC resources. Our results are conducted on up to 3586 cores for problem sizes with the parallelization-in-space scalability limited already on a single node. We gain significant reductions in the time-to-solution of 118.3 for spectral methods and 1503.0 for finite-difference methods with the parallelizationin-time approach. A developed and calibrated performance model gives the scalability limitations a-priory for this new approach and allows us to extrapolate the performance method towards large-scale system. This work has the potential to contribute as a basic building block of parallelization-in-time approaches, with possible major implications in applied areas modelling oscillatory dominated problems.The authors gratefully acknowledge the Gauss Centre for Supercomputing e.V. (www.gauss-centre.eu) for funding this project by providing computing time on the GCS Supercomputer SuperMUC at Leibniz Supercomputing Centre (LRZ, www.lrz. de). We also acknowledge use of Hartree Centre resources in this work on which the early evaluation of the parallelization concepts were done

Numerical solution of 3-D electromagnetic problems in exploration geophysics and its implementation on massively parallel computers

The growing significance, technical development and employment of electromagnetic (EM) methods in exploration geophysics have led to the increasing need for reliable and fast techniques of interpretation of 3-D EM data sets acquired in complex geological environments. The first and most important step to creating an inversion method is the development of a solver for the forward problem. In order to create an efficient, reliable and practical 3-D EM inversion, it is necessary to have a 3-D EM modelling code that is highly accurate, robust and very fast. This thesis focuses precisely on this crucial and very demanding step to building a 3-D EM interpretation method. The thesis presents as its main contribution a highly accurate, robust, very fast and extremely scalable numerical method for 3-D EM modelling in geophysics that is based on finite elements (FE) and designed to run on massively parallel computing platforms. Thanks to the fact that the FE approach supports completely unstructured tetrahedral meshes as well as local mesh refinements, the presented solver is able to represent complex geometries of subsurface structures very precisely and thus improve the solution accuracy and avoid misleading artefacts in images. Consequently, it can be successfully used in geological environments of arbitrary geometrical complexities. The parallel implementation of the method, which is based on the domain decomposition and a hybrid MPI-OpenMP scheme, has proved to be highly scalable - the achieved speed-up is close to the linear for more than a thousand processors. Thanks to this, the code is able to deal with extremely large problems, which may have hundreds of millions of degrees of freedom, in a very efficient way. The importance of having this forward-problem solver lies in the fact that it is now possible to create a 3-D EM inversion that can deal with data obtained in extremely complex geological environments in a way that is realistic for practical use in industry. So far, such imaging tool has not been proposed due to a lack of efficient, parallel FE solutions as well as the limitations of efficient solvers based on finite differences. In addition, the thesis discusses physical, mathematical and numerical aspects and challenges of 3-D EM modelling, which have been studied during my research in order to properly design the presented software for EM field simulations on 3-D areas of the Earth. Through this work, a physical problem formulation based on the secondary Coulomb-gauged EM potentials has been validated, proving that it can be successfully used with the standard nodal FE method to give highly accurate numerical solutions. Also, this work has shown that Krylov subspace iterative methods are the best solution for solving linear systems that arise after FE discretisation of the problem under consideration. More precisely, it has been discovered empirically that the best iterative method for this kind of problems is biconjugate gradient stabilised with an elaborate preconditioner. Since most commonly used preconditioners proved to be either unable to improve the convergence of the implemented solvers to the desired extent, or impractical in the parallel context, I have proposed a preconditioning technique for Krylov methods that is based on algebraic multigrid. Tests for various problems with different conductivity structures and characteristics have shown that the new preconditioner greatly improves the convergence of different Krylov subspace methods, which significantly reduces the total execution time of the program and improves the solution quality. Furthermore, the preconditioner is very practical for parallel implementation. Finally, it has been concluded that there are not any restrictions in employing classical parallel programming models, MPI and OpenMP, for parallelisation of the presented FE solver. Moreover, they have proved to be enough to provide an excellent scalability for it

A parallel Heap-Cell Method for Eikonal equations

Numerous applications of Eikonal equations prompted the development of many efficient numerical algorithms. The Heap-Cell Method (HCM) is a recent serial two-scale technique that has been shown to have advantages over other serial state-of-the-art solvers for a wide range of problems. This paper presents a parallelization of HCM for a shared memory architecture. The numerical experiments in $R^3$ show that the parallel HCM exhibits good algorithmic behavior and scales well, resulting in a very fast and practical solver. We further explore the influence on performance and scaling of data precision, early termination criteria, and the hardware architecture. A shorter version of this manuscript (omitting these more detailed tests) has been submitted to SIAM Journal on Scientific Computing in 2012.Comment: (a minor update to address the reviewers' comments) 31 pages; 15 figures; this is an expanded version of a paper accepted by SIAM Journal on Scientific Computin

3D time-varying simulations of Ca2+ dynamics in arterial coupled cells: A massively parallel implementation

Preferential locations of atherosclerotic plaque are strongly associated with the areas of low wall shear stress and disturbed haemodynamic characteristics such as flow detachment, flow recirculation and oscillatory flow. The areas of low wall shear stress are also associated with the reduced production of adenosine triphosphate in the endothelial layer, as well as the resulting reduced production of inositol trisphosphate (IP3). The subsequent variation in Ca2+ signalling and nitric oxide synthesis could lead to the impairment of the atheroprotective function played by nitric oxide. In previous studies, it has been suggested that the reduced IP3 and Ca2+ signalling can explain the correlation of atherosclerosis with induced low WSS and disturbed flow characteristics. The massively parallel implementation described in this article provides insight into the dynamics of coupled smooth muscle cells and endothelial cells mapped onto the surface of an idealised arterial bifurcation. We show that variations in coupling parameters, which model normal and pathological conditions, provide vastly different smooth muscle cell Ca2+ dynamics and wave propagation profiles. The extensibility of the coupled cells model and scalability of the implementation provide a solid framework for in silico investigations of the interaction between complex cellular chemistry and the macro-scale processes determined by fluid dynamics

New numerical approaches for modeling thermochemical convection in a compositionally stratified fluid

Seismic imaging of the mantle has revealed large and small scale heterogeneities in the lower mantle; specifically structures known as large low shear velocity provinces (LLSVP) below Africa and the South Pacific. Most interpretations propose that the heterogeneities are compositional in nature, differing in composition from the overlying mantle, an interpretation that would be consistent with chemical geodynamic models. Numerical modeling of persistent compositional interfaces presents challenges, even to state-of-the-art numerical methodology. For example, some numerical algorithms for advecting the compositional interface cannot maintain a sharp compositional boundary as the fluid migrates and distorts with time dependent fingering due to the numerical diffusion that has been added in order to maintain the upper and lower bounds on the composition variable and the stability of the advection method. In this work we present two new algorithms for maintaining a sharper computational boundary than the advection methods that are currently openly available to the computational mantle convection community; namely, a Discontinuous Galerkin method with a Bound Preserving limiter and a Volume-of-Fluid interface tracking algorithm. We compare these two new methods with two approaches commonly used for modeling the advection of two distinct, thermally driven, compositional fields in mantle convection problems; namely, an approach based on a high-order accurate finite element method advection algorithm that employs an artificial viscosity technique to maintain the upper and lower bounds on the composition variable as well as the stability of the advection algorithm and the advection of particles that carry a scalar quantity representing the location of each compositional field. All four of these algorithms are implemented in the open source FEM code ASPECT