An asymptotic preserving scheme for strongly anisotropic elliptic problems

In this article we introduce an asymptotic preserving scheme designed to compute the solution of a two dimensional elliptic equation presenting large anisotropies. We focus on an anisotropy aligned with one direction, the dominant part of the elliptic operator being supplemented with Neumann boundary conditions. A new scheme is introduced which allows an accurate resolution of this elliptic equation for an arbitrary anisotropy ratio.Comment: 21 page

Efficient Multigrid Preconditioners for Atmospheric Flow Simulations at High Aspect Ratio

Many problems in fluid modelling require the efficient solution of highly anisotropic elliptic partial differential equations (PDEs) in "flat" domains. For example, in numerical weather- and climate-prediction an elliptic PDE for the pressure correction has to be solved at every time step in a thin spherical shell representing the global atmosphere. This elliptic solve can be one of the computationally most demanding components in semi-implicit semi-Lagrangian time stepping methods which are very popular as they allow for larger model time steps and better overall performance. With increasing model resolution, algorithmically efficient and scalable algorithms are essential to run the code under tight operational time constraints. We discuss the theory and practical application of bespoke geometric multigrid preconditioners for equations of this type. The algorithms deal with the strong anisotropy in the vertical direction by using the tensor-product approach originally analysed by B\"{o}rm and Hiptmair [Numer. Algorithms, 26/3 (2001), pp. 219-234]. We extend the analysis to three dimensions under slightly weakened assumptions, and numerically demonstrate its efficiency for the solution of the elliptic PDE for the global pressure correction in atmospheric forecast models. For this we compare the performance of different multigrid preconditioners on a tensor-product grid with a semi-structured and quasi-uniform horizontal mesh and a one dimensional vertical grid. The code is implemented in the Distributed and Unified Numerics Environment (DUNE), which provides an easy-to-use and scalable environment for algorithms operating on tensor-product grids. Parallel scalability of our solvers on up to 20,480 cores is demonstrated on the HECToR supercomputer.Comment: 22 pages, 6 Figures, 2 Table

Duality-based Asymptotic-Preserving method for highly anisotropic diffusion equations

The present paper introduces an efficient and accurate numerical scheme for the solution of a highly anisotropic elliptic equation, the anisotropy direction being given by a variable vector field. This scheme is based on an asymptotic preserving reformulation of the original system, permitting an accurate resolution independently of the anisotropy strength and without the need of a mesh adapted to this anisotropy. The counterpart of this original procedure is the larger system size, enlarged by adding auxiliary variables and Lagrange multipliers. This Asymptotic-Preserving method generalizes the method investigated in a previous paper [arXiv:0903.4984v2] to the case of an arbitrary anisotropy direction field

Degenerate anisotropic elliptic problems and magnetized plasma simulations

This paper is devoted to the numerical approximation of a degenerate anisotropic elliptic problem. The numerical method is designed for arbitrary space-dependent anisotropy directions and does not require any specially adapted coordinate system. It is also designed to be equally accurate in the strongly and the mildly anisotropic cases. The method is applied to the Euler-Lorentz system, in the drift-fluid limit. This system provides a model for magnetized plasmas

Matrix-free GPU implementation of a preconditioned conjugate gradient solver for anisotropic elliptic PDEs

Many problems in geophysical and atmospheric modelling require the fast solution of elliptic partial differential equations (PDEs) in "flat" three dimensional geometries. In particular, an anisotropic elliptic PDE for the pressure correction has to be solved at every time step in the dynamical core of many numerical weather prediction models, and equations of a very similar structure arise in global ocean models, subsurface flow simulations and gas and oil reservoir modelling. The elliptic solve is often the bottleneck of the forecast, and an algorithmically optimal method has to be used and implemented efficiently. Graphics Processing Units have been shown to be highly efficient for a wide range of applications in scientific computing, and recently iterative solvers have been parallelised on these architectures. We describe the GPU implementation and optimisation of a Preconditioned Conjugate Gradient (PCG) algorithm for the solution of a three dimensional anisotropic elliptic PDE for the pressure correction in NWP. Our implementation exploits the strong vertical anisotropy of the elliptic operator in the construction of a suitable preconditioner. As the algorithm is memory bound, performance can be improved significantly by reducing the amount of global memory access. We achieve this by using a matrix-free implementation which does not require explicit storage of the matrix and instead recalculates the local stencil. Global memory access can also be reduced by rewriting the algorithm using loop fusion and we show that this further reduces the runtime on the GPU. We demonstrate the performance of our matrix-free GPU code by comparing it to a sequential CPU implementation and to a matrix-explicit GPU code which uses existing libraries. The absolute performance of the algorithm for different problem sizes is quantified in terms of floating point throughput and global memory bandwidth.Comment: 18 pages, 7 figure

A multiresolution space-time adaptive scheme for the bidomain model in electrocardiology

This work deals with the numerical solution of the monodomain and bidomain models of electrical activity of myocardial tissue. The bidomain model is a system consisting of a possibly degenerate parabolic PDE coupled with an elliptic PDE for the transmembrane and extracellular potentials, respectively. This system of two scalar PDEs is supplemented by a time-dependent ODE modeling the evolution of the so-called gating variable. In the simpler sub-case of the monodomain model, the elliptic PDE reduces to an algebraic equation. Two simple models for the membrane and ionic currents are considered, the Mitchell-Schaeffer model and the simpler FitzHugh-Nagumo model. Since typical solutions of the bidomain and monodomain models exhibit wavefronts with steep gradients, we propose a finite volume scheme enriched by a fully adaptive multiresolution method, whose basic purpose is to concentrate computational effort on zones of strong variation of the solution. Time adaptivity is achieved by two alternative devices, namely locally varying time stepping and a Runge-Kutta-Fehlberg-type adaptive time integration. A series of numerical examples demonstrates thatthese methods are efficient and sufficiently accurate to simulate the electrical activity in myocardial tissue with affordable effort. In addition, an optimalthreshold for discarding non-significant information in the multiresolution representation of the solution is derived, and the numerical efficiency and accuracy of the method is measured in terms of CPU time speed-up, memory compression, and errors in different norms.Comment: 25 pages, 41 figure