Tridiagonal Matrix Solver Research Articles

PittPack: An open-source Poisson’s equation solver for extreme-scale computing with accelerators

We present a parallel implementation of a direct solver for Poisson’s equation on extreme-scale supercomputers with accelerators. We introduce a chunked-pencil decomposition as the domain-decomposition strategy to distribute work among processing elements to achieve improved scalability at high counts of accelerators. Chunked-pencil decomposition enables overlapping MPI communication and data transfer between the central processing units (CPUs) and the graphics processing units (GPUs). It enables contiguous message transfer among the nodes and improves data locality by keeping neighboring elements in adjacent memory locations while permitting the use of shared memory for certain segments of the algorithm when possible. We study two different communication patterns within the chunked-pencil decomposition. The first pattern fully overlaps the communication with data transfer and aims to speedup the overall turnaround time. The second pattern concentrates on low memory usage and is more network friendly than the first pattern for computations at extreme scale. In our parallel implementation, we interleave OpenACC with MPI to support computations on the GPU or the CPU. The numerical solution and its formal second order of accuracy is verified using the method of manufactured solutions for various combinations of boundary conditions. Additionally, we used PittPack within an incompressible flow solver to further validate its accuracy and as well as demonstrate its versatility as a software package. We performed weak scaling analysis with up to 1.1 trillion Cartesian mesh points distributed over 16384 GPUs on a petascale leadership class supercomputer. Program SummaryProgram Title: PittPackProgram Files doi:http://dx.doi.org/10.17632/59ktwdby4r.1Licensing provisions: GNU General Public License 3Programming language: C++Nature of problem: An elliptic Poisson’s equation for pressure arises in the numerical solution of incompressible fluid flow. This equation needs to be solved accurately to enforce conservation of mass principle.Solution method: The Poisson’s equation is discretized with a second order accurate finite difference scheme. The resulting system of linear equations is solved with a direct method that employs fast Fourier transforms in the horizontal plane and a tridiagonal matrix solver in the third direction.Additional comments including restrictions and unusual features: The mesh has to be directionally uniform. Due to adoption of FFT in the horizontal plane, the boundary conditions in each direction must be of the same type. Each direction in the horizontal plane has to be partitioned by an equal integer number such that n2 parallel processes are deployed at run-time.

A Fast Poisson Solver of Second-order Accuracy for Isolated Systems in Three-dimensional Cartesian and Cylindrical Coordinates

Abstract We present an accurate and efficient method to calculate the gravitational potential of an isolated system in 3D Cartesian and cylindrical coordinates subject to vacuum (open) boundary conditions. Our method consists of two parts: an interior solver and a boundary solver. The interior solver adopts an eigenfunction expansion method together with a tridiagonal matrix solver to solve the Poisson equation subject to the zero boundary condition. The boundary solver employs James’s method to calculate the boundary potential due to the screening charges required to keep the zero boundary condition for the interior solver. A full computation of gravitational potential requires running the interior solver twice and the boundary solver once. We develop a method to compute the discrete Green’s function in cylindrical coordinates, which is an integral part of the James algorithm to maintain second-order accuracy. We implement our method in the Athena++ magnetohydrodynamics code and perform various tests to check that our solver is second-order accurate and exhibits good parallel performance.

A locally modified single-phase model for analyzing magnetohydrodynamic boundary layer flow and heat transfer of nanofluids over nonlinearly stretching sheet with chemical reaction

The problem of boundary layer flow and heat transfer of nanofluids over nonlinear stretching of a flat sheet in the presence of a magnetic field and chemical reaction is investigated numerically. In this paper, a new locally modified single-phase model for the analysis is introduced. In this model, the effective viscosity, density and thermal conductivity of the solid-liquid mixtures (nanofluids) which are commonly utilized in the homogenous single-phase model, are locally combined with the prevalent single-phase model. Similarity transformation is used to convert the governing equations into three coupled nonlinear ordinary differential equations. These equations depend on five local functions of the nanoparticle volume fraction viz., local viscosity ratio, magnetic, Prandtl, Brownian motion and thermophoresis functions. The equations are solved using Newton’s method and a block tridiagonal matrix solver. The results are compared to the prevalent single-phase model. In addition, the effect of important governing parameters on the velocity, temperature, volume fraction distribution and the heat and mass transfer rates are examined.

Research on tridiagonal matrix solver design based on a combination of processors

Large-scale tridiagonal matrix solvers based on heterogeneous systems currently cannot balance computational efficiency and numerical stability when solving a non-diagonally dominant matrix. A tridiagonal solver combined central processing unit with graphics processing unit is proposed, based on SPIKE2 as a solver framework, a simplified SPIKE algorithm as a central processing unit component, and a diagonal pivot algorithm as a graphics processing unit component. The solver performance is further improved by using a data-layout-transformation mechanism to obtain continuous addresses, reducing memory communication using constant memory to store unchanged data in the calculation process, and employing a kernel-fusion mechanism to reduce power consumption of graphics processing unit. For a diagonally dominant matrix, extended Thomas algorithms and cycle reduction to replace the graphics processing unit component are proposed in the solver. Experimental results show that the tridiagonal matrix solver in this paper can effectively consider both numerical stability and computational efficiency, and reduce total power consumption while improving memory efficiency.

The analysis of MHD blood flows through porous arteries using a locally modified homogenous nanofluids model

In this paper, magneto-hydrodynamic blood flows through porous arteries are numerically simulated using a locally modified homogenous nanofluids model. Blood is taken into account as the third-grade non-Newtonian fluid containing nanoparticles. In the modified nanofluids model, the viscosity, density, and thermal conductivity of the solid-liquid mixture (nanofluids) which are commonly utilized as an effective value, are locally combined with the prevalent single-phase model. The modified governing equations are solved numerically using Newton's method and a block tridiagonal matrix solver. The results are compared to the prevalent nanofluids single-phase model. In addition, the efficacies of important physical parameters such as pressure gradient, Brownian motion parameter, thermophoresis parameter, magnetic-field parameter, porosity parameter, and etc. on temperature, velocity and nanoparticles concentration profiles are examined.

Parallelisation of Implicit Time Domain Methods: Progress with ADI-FDTD

The parallelisation of implicit time domain solvers is essential for applications requiring large, highly over-sampled meshes that exceed the memory limitations of standalone workstations (1). Prime application examples include large arrays of small antennas (2), and on- chip interconnects for state of the art integrated circuits (IC) (3). The need for highly over sampled meshes in these applications can be illustrated by example. Current on-chip IC interconnect widths can be as low as 22nm or less, yet cover an area > 100(mm) 2 and transmit signals with fundamental frequencies in the range DC-10GHz. Typically, transient efiects involving high frequencies are of the most interest, but even this can still require meshes with ¢x < (‚0=10 6 ), i.e., four to flve orders of magnitude smaller than for typical FDTD meshes (neglecting considerations relating to material properties). It is not practical to use explicit FDTD solvers for such meshes, even in parallel, because the Courant-Friedrichs-Lewy (CFL) stability criteria enforces a commensurate reduction in the time step size for numerical reasons. While implicit FDTD solvers, such as ADI-FDTD, ofier freedom from the CFL stability criteria, it has been presumed that parallel implementations on all but the most specialised architectures (4) would be of little beneflt due to the high communication overhead. However, we have been able to show that this is not the case (5). In this Paper we will present an overview of our work to date on the parallelisation of implicit time domain methods, including parallel ADI-FDTD (1,5) and parallel ADI-BOR-FDTD (6) on both symmetric multiprocessor (SMP) and distributed memory computer clusters (DMCC). We do not parallelise the tri-diagonal matrix solver itself, because each 2D matrix must have more than 4,000 (40,000) elements per direction before this becomes e-cient for SMP (DMCC) ma- chines. This requires 3D domains at are at or beyond current memory limits for state of the art machines. Instead, we employ domain decomposition and solve multiple, smaller, tri-diagonal ma- trix systems in parallel. Carefully organised data exchanges avoid unnecessary double-handling of data during communication, improving the parallel algorithm performance. We will describe our domain decomposition scheme, and show results for small and large domains, comparing par- allel FDTD and parallel ADI-FDTD and parallel BOR-FDTD with parallel ADI-BOR-FDTD. We demonstrate e-cient solutions of large full 3D meshes with 8 billion mesh cells for paral- lel ADI-FDTD. Since machines with SMP and DMCC architectures are widely available, our demonstration of parallel speed up represents an important step forward for the application of implicit time domain solvers for large, highly oversampled meshes with ¢x < 10 i2 ‚. We expect that our parallelisation approach can be adopted for related implicit FDTD methods.

An ADI-FDTD method for periodic structures

In this communication, the alternating-direction implicit finite-difference time-domain (ADI-FDTD) method is extended to analyze periodic structures. In the implicit updates of the ADI-FDTD method, the periodic boundary condition leads to a cyclic matrix. Instead of inverting the cyclic matrix directly, the problem is converted into two auxiliary linear systems that can be solved using the tridiagonal matrix solver. Consequently, only 7n arithmetic operations are required for each implicit update and the efficiency of the ADI-FDTD method is retained. Numerical examples further demonstrate the effectiveness of this periodic ADI-FDTD method.

Cascade Flow Calculations Using the k - w; Turbulence Model with Explicit-Implicit Solver

A four-stage Runge-Kutta scheme, which is a typical explicit method, is adopted to solve the compressible Navier-Stokes equations. An implicit approximate factorization scheme, which can be solved by a tridiagonal matrix solver, is developed.

On an efficient numerical method for modeling sea ice dynamics

A computationally efficient numerical method for the solution of nonlinear sea ice dynamics models employing viscous‐plastic rheologies is presented. The method is based on a semi‐implicit decoupling of the x and y ice momentum equations into a form having better convergence properties than the coupled equations. While this decoupled form also speeds up solutions employing point relaxation methods, a line successive overrelaxation technique combined with a tridiagonal matrix solver procedure was found to converge particularly rapidly. The procedure is also applicable to the ice dynamics equations in orthogonal curvilinear coordinates which are given in explicit form for the special case of spherical coordinates.

Thermal characteristics of packed bed storage system

The thermal behaviour of a packed bed storage system charged with hot air is modelled using two partial differential equations representing the energy conservation in the air and solid phases constituting the bed. These two equations are coupled through the heat exchange process between the two phases. A fully implicit numerical scheme based on forward, upwind and central differencing for the time, first and second space derivatives, respectively, is used to solve the modelling equations. Marching technique is used for the air equation and a tri-diagonal matrix solver is employed to solve the solid equation. The solution yields the thermal structure of the bed, namely the air and solid temperature distribution inside the bed at any particular time, and the variation of total energy stored in the bed with time. The effect of bed length, solid diameter and void fraction on the thermal characteristics of the packed bed is studied. Further, the performance of the bed under variable inlet air temperature and mass flow rate is investigated.

Plasma turbulence calculations on the intel iPSC/860 (RX) hypercube

One approach to improving the real-time efficiency of plasma turbulence calculations is to use a parallel algorithm. A serial algorithm used for plasma turbulence calculations is modified to allocate a radial region in each node. In this way, convolutions at a fixed radius are performed in parallel, and communication is limited to boundary values for each radial region. For a semi-implicit numerical scheme (tridiagonal matrix solver), there is a factor of 2 improvement in efficiency with the Intel iPSC/860 machine using 64 processors over a single-processor CRAY-II. For block-tridiagonal matrix cases (fully implicit code), a second parallelization takes place. The Fourier components are distributed in nodes. In each node, the block-tridiagonal matrix is inverted for each of the allocated Fourier components. The algorithm for this second case has similar efficiency in both machines.