Fast Poisson Solvers Research Articles

Difference potentials method for models with dynamic boundary conditions and bulk-surface problems

In this work, we consider parabolic models with dynamic boundary conditions and parabolic bulk-surface problems in 3D. Such partial differential equations–based models describe phenomena that happen both on the surface and in the bulk/domain. These problems may appear in many applications, ranging from cell dynamics in biology, to grain growth models in polycrystalline materials. Using difference potentials framework, we develop novel numerical algorithms for the approximation of the problems. The constructed algorithms efficiently and accurately handle the coupling of the models in the bulk and on the surface, approximate 3D irregular geometry in the bulk by the use of only Cartesian meshes, employ fast Poisson solvers, and utilize spectral approximation on the surface. Several numerical tests are given to illustrate the robustness of the developed numerical algorithms.

FFT-based high order central difference schemes for three-dimensional Poisson's equation with various types of boundary conditions

In this paper, a unified approach is introduced to implement high order central difference schemes for solving Poisson's equation via the fast Fourier transform (FFT). Popular high order fast Poisson solvers in the literature include compact finite differences and spectral methods. However, FFT-based high order central difference schemes have never been developed for Poisson problems, because with long stencils, central differences require fictitious nodes outside the boundary, which poses a challenge to integrate boundary conditions in FFT computations. To overcome this difficulty, several layers of exterior grid lines are introduced to convert the problem to an immersed boundary problem with zero-padding solutions beyond the original cubic domain. Over the boundary of the enlarged cubic domain, the anti-symmetric property is naturally satisfied so that the FFT fast inversion is feasible, while the immersed boundary problem can be efficiently solved by the proposed augmented matched interface and boundary (AMIB) method. As the first fast Poisson solver based on high order central differences, the AMIB method can be easily implemented in any dimension, due to its tensor product nature of the discretization. As a systematical approach, the AMIB method can be made to arbitrarily high order in principle, and can handle the Dirichlet, Neumann, Robin or any combination of boundary conditions. The accuracy, efficiency, and robustness of the proposed AMIB method are numerically validated by considering various Poisson problems in two and three dimensions.

Fast electrostatic solvers for kinetic Monte Carlo simulations

Kinetic Monte Carlo (KMC) is an important computational tool in theoretical physics and chemistry. In contrast to standard Monte Carlo, KMC permits the description of time dependent dynamical processes and is not restricted to systems in equilibrium. Compared to Molecular Dynamics, it allows simulations over significantly longer timescales. Recently KMC has been applied successfully in modelling of novel energy materials such as Lithium-ion batteries and organic/perovskite solar cells. Motivated by this, we consider general solid state systems which contain free, interacting particles which can hop between localised sites in the material. The KMC transition rates for those hops depend on the change in total potential energy of the system. For charged particles this requires the frequent calculation of electrostatic interactions, which is usually the bottleneck of the simulation. To avoid this issue and obtain results in reasonable times, many studies replace the long-range potential by a phenomenological short range approximation. This, however, leads to systematic errors and unphysical results. On the other hand standard electrostatic solvers such as Ewald summation or fast Poisson solvers are highly inefficient in the KMC setup or introduce uncontrollable systematic errors at high resolution.In this paper we describe how the Fast Multipole Method by Greengard and Rokhlin can be adapted to overcome this issue by dramatically reducing computational costs. We exploit the fact that each update in the transition rate calculation corresponds to a single particle move and changes the configuration only by a small amount. This allows us to construct an algorithm which scales linearly in the number of charges for each KMC step, something which had not been deemed to be possible before.We demonstrate the performance and parallel scalability of the method by implementing it in a performance portable software library, which was recently developed in our group. We describe the high-level Python interface of the code which makes it easy to adapt to specific use cases.

Analysis of a constant-coefficient pressure equation method for fast computations of two-phase flows at high density ratios

An analysis of a modified pressure-correction formulation for fast simulations of fully resolved incompressible two-phase flows has been carried out. By splitting of the density weighted pressure gradient, the pressure equation is reduced to a constant-coefficient Poisson equation, for which efficient linear solvers can be used. While the gain in speed-up is well documented, the error introduced by the temporal extrapolation of the pressure gradient requires further investigations. In this paper it is shown that the modified pressure equation can lead to unphysical pressure oscillations and large errors. By appropriately combining the extrapolated pressure gradient with a matching volume fraction gradient grid convergence at high density ratios could be recovered. The cases of a one-dimensional front and a sphere translating at uniform velocity were first considered, allowing to decouple the pressure equation from the momentum equation. Subsequently, the case of a rising bubble in an upflow is analysed for which the full set of governing equations is solved. The pressure jump extrapolation error has been found dependent on the density ratio and the CFL number. Ultimately, the gain in the computational time, made possible by the use of fast Poisson solvers, should be weighted by the additional computational time the reduction of the aforementioned error may require.

Efficient Numerical Algorithms Based on Difference Potentials for Chemotaxis Systems in 3D

In this work, we propose efficient and accurate numerical algorithms based on Difference Potentials Method for numerical solution of chemotaxis systems and related models in 3D. The developed algorithms handle 3D irregular geometry with the use of only Cartesian meshes and employ Fast Poisson Solvers. In addition, to further enhance computational efficiency of the methods, we design a Difference-Potentials-based domain decomposition approach which allows mesh adaptivity and easy parallelization of the algorithm in space. Extensive numerical experiments are presented to illustrate the accuracy, efficiency and robustness of the developed numerical algorithms.

FROM IIM TO AUGMENTED IIM: A POWERFUL TOOL FOR COMPLEX PROBLEMS USING CARTESIAN MESHES

The immersed interface method (IIM) ?rst proposed in is an accurate numerical method for solving elliptic interface problems on Cartesian meshes. It is a sharp interface method that was intended to improve accuracy of the immersed boundary (IB) method. The IIM is second order accurate in the maximum norm (pointwise, strongest) while the IB method is ?rst order accurate. The ?rst IIM paper is one of the most downloaded one from the SIAM website and is one of the most cited papers. While IIM provided a way of accurate discretization of the partial differential equations (PDEs) with discontinuous coefficients, the augmented IIM ?rst proposed in made the IIM much more efficient and faster by utilizing existing fast Poisson solvers. More important is that the augmented IIM provides an efficient way for multi-physics models with different governing equations, problems on irregular domains, multi-scales and multi-connected domains. A brie?y introduction of the augmented strategy including some recently progress is presented in this article.

POSTPROCESSING MIXED FINITE ELEMENT METHODS FOR SOLVING CAHN-HILLIARD EQUATION: METHODS AND ERROR ANALYSIS.

A postprocessing technique for mixed finite element methods for the Cahn-Hilliard equation is developed and analyzed. Once the mixed finite element approximations have been computed at a fixed time on the coarser mesh, the approximations are postprocessed by solving two decoupled Poisson equations in an enriched finite element space (either on a finer grid or a higher-order space) for which many fast Poisson solvers can be applied. The nonlinear iteration is only applied to a much smaller size problem and the computational cost using Newton and direct solvers is negligible compared with the cost of the linear problem. The analysis presented here shows that this technique remains the optimal rate of convergence for both the concentration and the chemical potential approximations. The corresponding error estimate obtained in our paper, especially the negative norm error estimates, are non-trivial and different with the existing results in the literatures.

Large Eddy Simulation of Wind Turbine Wakes

We present the coupling of a vortex particle-mesh method with immersed lifting lines for the Large Eddy Simulation of wind turbine wakes. The method relies on the Lagrangian discretization of the Navier–Stokes equations in vorticity-velocity formulation. Advection is handled by the particles while the mesh allows the evaluation of the differential operators and the use of fast Poisson solvers. We use a Fourier-based fast Poisson solver which simultaneously allows unbounded directions and inlet/outlet boundaries. The method also allows the feeding of a turbulent incoming flow. We apply this methodology to the study of large scale aerodynamics and wake behavior of tandem wind turbines. We analyze the generators performance, unsteady power, loads and aerodynamics they are subjected to. The average flow field of the wakes is also computed and turbulence statistics are extracted. In particular, we investigate the influence of the type of turbulent inflow used—precomputed or synthetic—, and study wake meandering.

Asymptotic-preserving schemes for kinetic-fluid modeling of disperse two-phase flows

We consider a system coupling the incompressible Navier–Stokes equations to the Vlasov–Fokker–Planck equation. Such a problem arises in the description of particulate flows. We design a numerical scheme to simulate the behavior of the system. This scheme is asymptotic-preserving, thus efficient in both the kinetic and hydrodynamic regimes. It has a numerical stability condition controlled by the non-stiff convection operator, with an implicit treatment of the stiff drag term and the Fokker–Planck operator. Yet, consistent to a standard asymptotic-preserving Fokker–Planck solver or an incompressible Navier–Stokes solver, only the conjugate–gradient method and fast Poisson and Helmholtz solvers are needed. Numerical experiments are presented to demonstrate the accuracy and asymptotic behavior of the scheme, with several interesting applications.

A NOVEL FAST SOLVER FOR POISSON'S EQUATION WITH NEUMANN BOUNDARY CONDITION

In this paper, we present a novel fast method to solve Poisson's equation in an arbitrary two dimensional region with Neumann boundary condition, which are frequently encountered in solving electrostatic boundary problems. The basic idea is to solve the original Poisson's equation by a two-step procedure. In the flrst stage, we expand the electric fleld of interest by a set of tree basis functions and solve it with a fast tree solver in O(N) operations. The fleld such obtained, however, fails to expand the exact fleld because the tree basis is not curl-free. Despite of this, we can retrieve the correct electric fleld by purging the divergence-free fleld. Next, for the second stage, we flnd the potential distribution rapidly with a same fast solution of O(N) complexity. As a result, the proposed method dramatically reduces solution time compared with traditional FEM methods. In addition, it is the flrst time that the loop-tree decomposition technique has been introduced to develop fast Poisson solvers. Numerical examples including electrostatic simulations are presented to demonstrate the e-ciency of the proposed method.

Fast poisson solvers for thermal analysis

Accurate and efficient thermal analysis for a VLSI chip is crucial, both for sign-off reliability verification and for design-time circuit optimization. To determine an accurate temperature profile, it is important to simulate a die together with its thermal mounts: this requires solving Poisson's equation on a nonrectangular 3D domain. This article presents a class of eigendecomposition-based Fast Poisson Solvers (FPS) for chip-level thermal analysis. We start with a solver that solves a rectangular 3D domain with mixed boundary conditions in O( N ⋅ log N ) time, where N is the dimension of the finite difference matrix. Then we reveal, for the first time in the literature, a strong relation between fast Poisson solvers and Green-function-based methods. Finally, we propose an FPS method that leverages the preconditioned conjugate gradient method to solve nonrectangular 3D domains efficiently. We demonstrate this approach on thermal analysis of an industrial microprocessor, showing accurate results verified by a commercial tool, and that it solves a system of dimension 4.54e6 in only 13 conjugate gradient iterations, with a runtime of 65 seconds, a 15X speedup over the popular ICCG solver.

Upwind-Difference Potentials Method for Patlak-Keller-Segel Chemotaxis Model

We develop a novel upwind-difference potentials method for the Patlak-Keller-Segel chemotaxis model that can be used to approximate problems in complex geometries. The chemotaxis model under consideration is described by a system of two nonlinear PDEs: a convection-diffusion equation for the cell density coupled with a reaction-diffusion equation for the chemoattractant concentration. Chemotaxis is an important process in many medical and biological applications, including bacteria/cell aggregation and pattern formation mechanisms, as well as tumor growth. Furthermore modeling of real biomedical problems often has to deal with the complex structure of computational domains. There is consequently a need for accurate, fast, and computationally efficient numerical methods for different chemotaxis models that can handle arbitrary geometries. The upwind-difference potentials method proposed here handles complex domains with the use of only Cartesian meshes, and can be easily combined with fast Poisson solvers. In the numerical tests presented below we demonstrate the robustness of the proposed scheme.

A formulation to compute mass-consistent models of hydrodynamic flows

Standard interpolation methods of measured data of an incompressible fluid yield a non-solenoidal field v 0. A formulation to estimate a solenoidal field from v 0, is proposed. Variational calculus reduces the problem to the solution of an elliptic equation for a Lagrange multiplier. Examples illustrate how boundary conditions improve the mass-balance of velocity fields obtained in meteorology with similar approaches. The elliptic equation is separable in meteorological problems over a complex orography. This allows the use of fast-Poisson solvers. It is shown how the flow-rate can be used to define a low-pass filter which improves the results given by the Fast Fourier Algorithm.

A Fractional Step Immersed Boundary Method for Stokes Flow with an Inextensible Interface Enclosing a Solid Particle

In this paper, we develop a fractional step method based on the immersed boundary (IB) formulation for Stokes flow with an inextensible (incompressible) interface enclosing a solid particle. In addition to solving for the fluid variables such as the velocity and pressure, the present problem involves finding an extra unknown elastic tension such that the surface divergence of the velocity is zero along the interface, and an extra unknown particle surface force such that the velocity satisfies the no-slip boundary condition along the particle surface. While the interface moves with local fluid velocity, the enclosed particle hereby undergoes a rigid body motion, and the system is closed by the force-free and torque-free conditions along the particle surface. The equations are then discretized by standard centered difference schemes on a staggered grid, and the interactions between the interface and particle with the fluid are discretized using a discrete delta function as in the IB method. The resultant linear system of equations is symmetric and can be solved by fractional steps so that only fast Poisson solvers are involved. The present method can be extended to Navier--Stokes flow with moderate Reynolds number by treating the nonlinear advection terms explicitly for the time integration. The convergent tests for a Stokes solver with or without an inextensible interface are performed and confirm the desired accuracy. The tank-treading to tumbling motion for an inextensible interface enclosing a solid particle with different filling fractions under a simple shear flow has been studied extensively, and the results here are in good agreement with those obtained in literature. (A corrected version is attached.)

Gene H. Golub 1932–2007

Gene Howard Golub, professor in the Computer Science Department of Stanford University (USA) passed away suddenly from an acute myeloid leukemia in the Stanford Hospital on November 16, 2007. Gene was one of the prominent figures of numerical linear algebra or, as he preferred to call it, matrix computations. Gene was born in Chicago (Illinois) on February 29, 1932, from parents who emigrated from Ukraine and Latvia to the United States in 1923. He did his graduate studies first in the University of Chicago and then, from 1953 to 1959, in the University of Illinois at Urbana-Champaign. There, he studied matrix theory and statistics. He had a partial time job in the computer center and learn how to program one of the first parallel computers, the ILLIAC IV. For a while he was considering doing a Ph.D. in statistics but he became a student of Professor Abraham Taub, an applied mathematician, who oriented him toward the study of Chebyshev polynomials for solving linear systems, starting from John von Neumann’s works. In 1959 Taub invited Richard Varga to the University of Illinois and Gene discovered that Varga was working on the same topic. This lead to the writing of a seminal joint paper that was published in 1961. After his Ph.D. thesis Gene received an NSF fellowship that allowed him to spend 15 months in Cambridge (UK) where he met for the first time Jim Wilkinson who was working at the National Physical Laboratory. Jim became a prominent expert in rounding error analysis in numerical linear algebra algorithms and a frequent visitor to Stanford. Gene met also Cornelius Lanczos in the UK. When he was back to the US, he worked for several industrial companies and then decided to back to academics. In 1962, George Forsythe offered him an assistant professorship in Stanford and soon a permanent position. It was in these times that Gene met some of Forsythe’s students that were about to become famous like Cleve Moler (founder of Matlab and present SIAM president) and Beresford Parlett. The Stanford University Computer Science Department was founded by Forsythe in 1966 and Gene naturally get a position in the department together with John McCarthy, Donald Knuth and other well–known researchers. He was chairman of the department in 1982. During his brilliant career Gene authored or co–authored more than 180 papers published in the best journals and many contributed or invited papers to an uncountable number of international conferences. It is almost impossible to summarize all of his many important contributions to applied mathematics. Perhaps the most well– known ones are the use of QR factorizations to solve least squares problems and the invention of direct and iterative algorithms for the computation of the singular value decomposition (SVD) together with W. Kahan. SVD was a great passion of Gene and he used it in many applications. He was so in love with SVD that the license plate of his last car was “Pr. SVD”. Noticeable also is the discovery of the cyclic reduction algorithm for solving certain structured linear systems leading to the construction of fast Poisson solvers. Gene was one of the most active evangelists for the preconditioned conjugate gradient algorithm through his joint papers with P. Concus and D.P. O’Leary. More recently he worked on algorithms for computing bounds or estimates of bilinear forms u f(A)v where u and v are vectors, A is a square

Performance of Multigrid in the Context of Beam Dynamics Simulations

Precise and fast 3D space charge calculations for bunches or clouds of charged particles are of growing importance in design studies for future linear accelerators and light sources. One of the possible approaches is the computation of the potential of the bunch in the rest frame by means of Poisson's equation. The software package MOEVE has been developed for space charge calculations on non-equidistant grids. It consists of several iterative Poisson solvers (MOEVE Poisson solvers), among them the state-of-the-art multigrid Poisson solver. Furthermore, the MOEVE Poisson solvers have been implemented in the tracking code Astra and GPT. In this paper the algorithms of the software package MOEVE will be described and the performance will be tested for very large linear systems of equations. The numerical results will show that the conditions for the Poisson solvers have to be chosen carefully in order to achieve optimal results for real life applications. The design of future light sources and colliders requires increasingly precise 3D beam dynamics sim- ulations. In so-called tracking simulations the particle trajectory is determined which is described by the relativistic equation of motion. The equation of motion is solved by means of an appropriate time integration scheme. In regimes of rather low energy, this implies that the space charge flelds have to be taken into account in each time step of the numerical integration. Recently, the e-cient calculation of 3D space charge flelds gained particular importance in the context of electron cloud studies for the ILC (International Linear Collider) damping rings. Based on the geometric multigrid technique fast Poisson solvers have been developed and suc- cessfully applied for 3D space charge simulations. In theory, these multigrid Poisson solvers have optimal performance, i.e., the numerical efiort depends linearly on the number of mesh points. Unfortunately, this optimal convergence rate can sometimes not be achieved in simulations of real life problems. In this paper the performance of the geometric multigrid technique is investigated in the context of space charge calculations, in particular, for huge numbers of mesh points, i.e., up to 4 million. Another problem is the behavior of the multigrid Poisson solver within the particle tracking procedure. Since space charge flelds have to be computed in each time step of the numer- ical integration, the calculated flelds of the previous time step can be used as initial guess for the iterative Poisson solver. With this approach the efiort for the new space charge calculation can be reduced. The question under investigation is how multigrid can compete with other iterative algorithms for real life tracking simulations. The iterative Poisson solvers are available as software package MOEVE 2.0 (MOEVE: Multigrid for non-equidistant grids to solve Poisson's equation) (8). Furthermore, these Poisson solvers are implemented in the tracking code Astra (DESY, Hamburg, Germany) (3) and the tracking code GPT (Pulsar Physics, Eindhoven, The Netherlands) (13). 2. MATHEMATICAL MODEL Space charge calculations for beam dynamics studies are performed within a tracking procedure. The tracking is a method to determine the trajectories of the particles which are described by the relativistic equations of motion (1). The equations of motion are solved by means of an appropriate time integration scheme. The space charge flelds have to be taken into account in each time step of the numerical integration. The space charge calculations are performed in the rest frame of the bunch by means of Poisson's equation given by i¢' = % 0 in › ‰R 3 ;

Sobolev gradient preconditioning for the electrostatic potential equation

Sobolev gradient type preconditioning is proposed for the numerical solution of the electrostatic potential equation. A constructive representation of the gradients leads to efficient Laplacian preconditioners in the iteration thanks to the available fast Poisson solvers. Convergence is then verified for the corresponding sequence in Sobolev space, implying mesh independent convergence results for the discretized problems. A particular study is devoted to the case of a ball: due to the radial symmetry of this domain, a direct realization without discretization is feasible. The simplicity of the algorithm and the fast linear convergence are finally illustrated in a numerical test example.

A parallel two-level hybrid method for tridiagonal systems and its application to fast poisson solvers

A new method, namely, the parallel two-level hybrid (PTH) method, is developed to solve tridiagonal systems on parallel computers. PTH has two levels of parallelism. The first level is based on algorithms developed from the Sherman-Morrison modification formula, and the second level can choose different parallel tridiagonal solvers for different applications. By choosing different outer and inner solvers and by controlling its two-level partition, PTH can deliver better performance for different applications on different machine ensembles and problem sizes. In an extreme case, the two levels of parallelism can be merged into one, and PTH can be the best algorithm otherwise available. Theoretical analyses and numerical experiments indicate that PTH is significantly better than existing methods on massively parallel computers. For instance, using PTH in a fast Poisson solver results in a 2-folds speedup compared to a conventional parallel Poisson solver on a 512 nodes IBM machine. When only the tridiagonal solver is considered, PTH is over 10 times faster than the currently used implementation.

A Fourth Order Scheme for Incompressible Boussinesq Equations

A fourth order finite difference method is presented for the 2D unsteady viscous incompressible Boussinesq equations in vorticity-stream function formulation. The method is especially suitable for moderate to large Reynolds number flows. The momentum equation is discretized by a compact fourth order scheme with the no-slip boundary condition enforced using a local vorticity boundary condition. Fourth order long-stencil discretizations are used for the temperature transport equation with one-sided extrapolation applied near the boundary. The time stepping scheme for both equations is classical fourth order Runge–Kutta. The method is highly efficient. The main computation consists of the solution of two Poisson-like equations at each Runge–Kutta time stage for which standard FFT based fast Poisson solvers are used. An example of Lorenz flow is presented, in which the full fourth order accuracy is checked. The numerical simulation of a strong shear flow induced by a temperature jump, is resolved by two perfectly matching resolutions. Additionally, we present benchmark quality simulations of a differentially-heated cavity problem. This flow was the focus of a special session at the first MIT conference on Computational Fluid and Solid Mechanics in June 2001.

Parallel Distributed FFT-Based Solvers for 3-D Poisson Problems in Meso-Scale Atmospheric Simulations

Fast Poisson solvers based on fast Fourier transform (FFT) computations are among the fastest techniques to solve Poisson equation on uniform grids. In this paper, the authors present two parallel distributed implementations of a 3-D fast Poisson solver in the context of the atmospheric simulation code. The first approach consists of implementing data movement between each computational step so that no elementary computational routine involves communication. The second approach aims at reducing the global data movement and requires the parallelization of one step of the fast Poisson solver. Experimental results on a 128-node Cray T3E are presented to illustrate the advantages and drawbacks of each approach.