Preconditioning flow- and wave propagation problems using fundamental solutions

Motivation

Recent advances in computer technology and numerical algorithms have opened new possibilities for realistic computer simulations based on advanced mathematical models. An integral part of the effort of making applied mathematics and numerical simulation techniques available and widely used is further development of efficient and robust numerical algorithms. For many problems occurring in applications, efficient techniques exist. However, in some important situations, there is still a need for improvements when aiming at solving application problems.

Goal

The overall goal of the proposed project is to provide an efficient numerical technique, suitable for high-performance parallel computers, for solving flow- and wave propagation problems. Two important applications are studied; the flow of a viscous, incompressible fluid, and the propagation of electromagnetic waves. The new solution methods will be applied to problems studied in the GEMS (General ElectroMagnetic Solver) project within the center of excellence PSCI (NUTEK/VINNOVA), and to problems arising in direct numerical simulation of turbulent flows. The project builds on previous experience of solving similar problems, but the new methods will specifically be aimed at solving problems where the governing partial differential equation (PDE) is discretized on an unstructured or curvilinear grid.

It is important to note that the numerical methods that will be developed are non-domain specific; apart from the problem areas mentioned above, the same type of solution methodology could be applied to equations describing inviscid and/or compressible flow and wave propagation in other media, e.g., in magnetohydrodynamics and different regimes of acoustics.

Organization and environment

The research at the Department of Scientific Computing (TDB) is organized in informal research groups. The project described below is part of the work planned in the group Solving systems of equations for PDE problems, or in short-hand notation, the EQSYS group.

The new project will primarily be performed by two faculty members (Holmgren and Brandén) and one PhD student (Sundqvist). Holmgren is the coordinator of the work in the EQSYS group, and was one of the original members when the group was formed. Brandén received his PhD degree within the group in February 2001, and has been active in the project described here since it started with the MSc thesis work by Sundqvist during winter 1999/2000, where Brandén was the primary advisor. Sundqvist was admitted as a PhD student in April 2000, and has performed initial studies within the project since then. So far, no funding has been available for senior researches within the project. Still, significant progress has been made, resulting in one submitted paper [?] (Brandén, Sundqvist), one manuscript [?] (Sundqvist, Brandén and Holmgren), and three oral presentations at international conferences by Sundqvist. For the continuation of the project according to the plan below, additional recourses according to the budget are required. Holmgren has applied for funding for Sundqvist from the Graduate School in Mathematics and Computing at Uppsala University. The descision if this grant will be approved or not will be taken 010514.

The EQSYS group has a well developed collaboration with the numerical analysis group at Oxford University (Dr Andy Wathen, Prof. Nick Trefethen). The Oxford group has word-leading expertise in using unstructured grids for industrial problems, which will be essential for the project described here. During autumn 1999, Brandén visited Oxford University for three months. During the same time period, also Sundqvist made a shorter visit in Oxford. The intention is that Sundqvist will visit Oxford for another short period closer to the end of the project. The EQSYS group also has collaboration with Dr Alison Ramage, University of Strathclyde, Scotland. Ramage is also studying iterative methods for problems discretized on unstructured grids, and her expertise will be valuable for the project.

The project will benefit from the results derived in the proposed continuation of the “ramprogram” project High-performance computing methods for accurate simulation of wave propagation and turbulence, which is coordinated by Gustafsson and Lötstedt. However, note that the research considered in the project described below is not covered within the “ramprogram”.

The EQSYS group has established contacts with industry within the NUTEK/VINNOVA center of excellence Parallel and Scientific Computation Institute (PSCI). Of special interest here is the involvement in the General Electromagnetic Solver (GEMS) project, see

http://www.psci.kth.se/Programs/GEMS/.

Holmgren presently has a 30% grant as course coordinator for the National Graduate School in Scientific Computing (NGSSC), and previously had a grant as HPDR lecturer. Thus, Holmgren is well informed concerning the current development in the large-scale computation arena in Sweden and internationally, and has many useful contacts.

Project plan

Background

Since the EQSYS group was formed in the end of the 1980:s, considerable competence in the field of numerical solution of large-scale flow and wave propagation problems, including parallel implementations, has been acquired. Problems of immediate relevance for applications in industry and technological/natural sciences have been approached, and in the project described here this line of research will be pursued. We will combine the experience and competence gained from previous projects with inside knowledge on applications and exciting new ideas, studying new challenging issues.

For many PDE problem settings, the major part of the computational work consist of solving large, sparse systems of linear equations. The state of the art for solving these systems is to employ some Krylov subspace method [?]. In order to achieve an acceptable rate of convergence and, more importantly, a short total execution time, it is crucial to construct effective, parallelizable preconditioners. For the type of PDEs that we are considering, standard preconditioning techniques often perform poorly. A viable alternative is to use preconditioners based on fast transforms, which has been the main research area for the EQSYS group. Originally, the preconditioners were defined using a linear algebra approach, applying circulant approximations as preconditioners for Toeplitz matrices [?], [?]. The preconditioning technique has been generalized to non-Toeplitz systems [?, and the references therein] and other transforms [?, ?, ?]. For PDE problems, members of the EQSYS group have contributed to the development [?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?]. In the work within the group, the fast-transform based preconditioners have also been generalized to PDE problems arising from applications in fluid mechanics, electromagnetics and acoustics [?, ?, ?, ?, ?, ?]. This has involved treating nonlinear problems, variable coefficients, and relatively complex geometries. Also, the technique has been applied in domain decomposition settings [?, ?, ?, ?], and implemented on parallel computers [?, ?, ?, ?].

The Green’s function approach

The fast transform-based preconditioners described above are intimately connected to using structured grids for the spatial discretization of the PDE. The preconditioners are defined using the spatial difference operator, i.e., after the discretization is performed. In both the preconditioner set-up and the preconditioner solve, the algorithms employ that the grid is structured. In many cases, complex problem settings can be handled using multi-block structured grids. However, in some situations, it is advantageous to use an unstructured or highly stretched grid. The goal for the project described here is to develop effective preconditioners also for such problems, using a technique where the preconditioner is defined using the differential operator in space, which is of course independent of the discretization. Both first and second order PDE will be considered.

First order PDE

The application we will consider here is electromagnetic wave propagation. The problem is governed by a time-dependent first order system of PDEs, the Maxwell equations. Usually, problems of this type are solved using an explicit time marching method, i.e., for problems in elecromagnetics, the Yee scheme [?] is often employed. However, for some problems of interest in applications, the cell sizes in the mesh has to vary by orders of magnitude to resolve small geometrical features present in the problem [?]. The smallest cell size is much smaller than required for resolving the wavelength of the electromagnetic waves. If an explicit scheme is employed, the time step has to be unrealistically small for stability reasons. Another problem is that there is often a problem with long time instability with the explicit schemes [?]. A solution to these dilemmas is to use an implicit method such as a backward differentiation formula (BDF) for advancing the solution in time [?]. Then, a system of linear equations has to be solved in every time step. Formulating the system of PDEs as a first order system implies that these systems will have positive eigenvalues, i.e., they are definite. For problems in 2D and 3D, only iterative methods result in reasonable memory consumption, and such an iteration needs an effective preconditioner. In the project, we will apply the approach described below for the implicit solution of the Maxwell equations in 2D in the transverse electric or transverse magnetic case. Here, we will use problem settings of interest for the GEMS project in PSCI, see above.

In [?], preliminary studies of preconditioners for first order PDE is performed, and this approach will be pursued in this part of the project. Here, a general PDE

P u = f, x ∈ Ω ⊂ ℜ^{d},

with boundary conditions

B_{i}u = b_{i}, x ∈ Γ_{i} ⊂ ∂Ω, i = 1, …, q,

is considered. By introducing

B = ∑_{i=1}^{q} δ_{Γ_{i}}B_{i}, b = ∑_{i=1}^{q} δ_{Γ_{i}}b_{i},

the problem may be written as a single equation

(P + B)u = f + b, x ∈ \overline{Ω},  (1)

where δ_{Γ_{i}} denotes the Dirac measure with support on Γ_{i}.

Equation (1) will be discretized by a a finite volume or finite element method. Let P̃ be a linear differential operator with constant coefficients approximating P. Then, P̃ has

a fundamental solution, i.e. there is a function g satisfying P g = δ, where δ is the Dirac measure. We will employ

Ku(x) = ∫\overline{Ω} g(x − y)u(y) dy, x ∈ \overline{Ω},

as a preconditioner for Equation (1). Note that, when constructing the preconditioner, the properties of the of the discretized problem, i.e., the difference operator, are not used.

On an infinite domain, without boundary conditions, K would be the exact inverse of P̃, suggesting that K approximates the inverse of P + B and yields a good preconditioner.

The preconditioned problem

K(P + B)u = K(f + b),

should be easier to solve. The preliminary studies show that the operators K and K(P + B) have some interesting properties:

1. The integral operator K is a convolution operator, and it is possible to compute the action of K using the Fourier transform. In the discrete case, FFT makes a fast application possible on a Cartesian grid. Also, the memory requirements are low, since K can be represented by the fundamental solution only.

2. Theoretical results show that K(P +B) is a bounded operator when P is a first order operator. This indicates that grid independent convergence is possible, and numerical experiments for model problems show [?] that this can be achieved. For linear problems with constant coefficients it may even be possible to construct an exact inverse.

For higher order differential operators, K(P +B) is unbounded. This observation is consistent with results in [?], [?] and [?].

In [?], results for a few scalar model problems are given. In the project, we will generalize the technique to implicit solution of systems of PDE like the Maxwell equations, that includes first order derivatives only. We will consider the following topics:

1. Although existence of fundamental solutions for linear partial differential operators with constant coefficients has been proven, they may be difficult to derive in physical space. However, note that in Fourier space, the derivation of a fundamental solution is trivial. Since we intend to apply the integral operator K using Fourier technique, a representation in Fourier space is sufficient. It is even possible to automatize the computation. Furthermore, a fundamental solution of a discretized operator may also easily be computed in Fourier space.

2. A fast application of K, using the fast Fourier transform, is straight forward on structured grids. For unstructured grids, an efficient interpolation procedure, mapping values from a structured “Fourier grid” to an unstructured computational grid will be developed. Fast Fourier transforms for nonequispaced data have been investigated in [?].

3. For problems with rapidly varying or discontinuous coefficients, a domain decomposition approach will be used to decompose the domain into several domains with mildly varying or constant coefficients.

The preconditioning technique we suggest is, to the best of our knowledge, a new approach. However, work based on the same basic observation that a preconditioner should be related to a fundamental solution or a Green’s function has been published in [?] and [?]. The authors observe that the Green’s tensor of the Oseen equations can be expressed in terms of Green’s functions of the Poisson equation and the advection diffusion equation. By using efficient solvers for these problems, they derive a good preconditioner for the Oseen equations, which is a simplification of the Navier–Stokes equations.

Second order equations

As an application problem, we will here consider the incompressible Navier-Stokes equations, which governs laminar and turbulent flow of an incompressible fluid. The space derivatives are approximated by a finite difference, finite volume or finite element method, and the time marching method is a normally a mixed implicit-explicit scheme. To allow for large Reynolds numbers and a small grid size, the viscous part of the operator is treated implicitly and the nonlinear convective part is advanced explicitly. The velocity at the new time level should be divergence-free. The result is a system of linear equations to solve in each time step for the velocity and the pressure. The system, which is almost symmetric (depending on the boundary conditions) but indefinite, is the same as we obtain when we solve Stokes’ equations.

In [?], the Stokes equations are solved using the Uzawa algorithm and different preconditioners. The linearized Oseen equations are solved with a preconditioner in [?, ?]. In [?], it is shown how the system of equations for the Navier-Stokes equations can be solved approximately in three steps with second order accuracy. The first equation is a discretized heat equation for an intermediate velocity. The second equation is almost a Poisson equation for the pressure. Then the new velocity can be updated at the new time level. The PDE has constant coefficients. We intend to solve the equations in the first two steps iteratively using preconditioners based on approximate Green’s functions. This idea is related to the work presented in [?, ?, ?, ?].

As remarked in the previous section, the operator K(P + B) is unbounded for PDE problems of higher order than one. This implies that, for a second order PDE, the preconditioning technique should take the boundary conditions of the PDE into account. It is not sufficient to use only the fundamental solution. The solution u(x) to the Poisson equation with constant coefficients in the domain Ω with Dirichlet boundary conditions on the boundary ∂Ω and the right hand side f can be written [?]

u(x) = ∫_{Ω} f(ξ)g(x − y)dy + ∫_{∂Ω} σ(y)g_{n}(x − y)dy,  (2)

where g is the fundamental solution of the Laplace operator and a subscript n denotes the normal derivative at ∂Ω. The doublet density σ is determined by an integral equation such that the boundary conditions are satisfied by u. This integral equation will be solved using well established techniques, see below. We have chosen to use the fundamental solution in the expression in (2) since it is explicitly known. Using the Green’s function, the second integral is not needed in (2) but for arbitrary boundaries ∂Ω, the Green’s function is general not known.

After discretization in space of the PDE we have a system of linear equations to solve for the unknowns. An accurate approximation of the inverse of the system matrix is the best preconditioner of the system. The integral relation in (2) provides us with such an approximation.

It is important that the evaluation of the first integral in (2) is efficient. Assume that the number of unknowns in 2D are of O(n2). The fast multipole method [?] and the panel clustering technique [?] are methods that compute the integral to high precision in O(n2) or O(n2 log n) operations. Also the integral equation for O(n) components in σ in (2) can be solved in at most O(n2) operations using the same type of techniques. In our case we need only an approximation of the right hand side and we can reduce the computational cost and programming efforts substantially by simplifying the calculations. How precise this evaluation has to be for an efficient iterative method will be investigated in the project. In [?] the Poisson equation is solved with a fast solver with special boundary conditions (thus replacing the need to evaluate the first integral in (2)) and then a boundary correction is added with a boundary integral as in [?].

The fundamental solution g is known explicitly only for constant coefficient problems. If our coefficients were constant, then we could compute the solution using a boundary element method derived from (2). In our application, the coefficients may vary slowly over the grid. Since we need only approximate inverses, we will use the fundamental solution for the the Poisson equation with locally or globally frozen coefficients. The extreme choices here are the point values and the averaged values over Ω.

The same technique will be applied also to the heat equation, which defines the first step in the procedure to solve the Navier-Stokes equation.

Parallel implementations

As an integral part of the work, we will implement and study the performance of the new solution methods on high-performance computers. We intend to implement parallel versions of the preconditioner solves using shared memory parallelization techniques. Within the EQSYS group, there is significant experience of parallel computing, resent results which are relevant for the project described here are presented in [?] and [?]. The department has recently acquired a new parallel computer system with 32 processors, which is currently beeing upgraded to a completely new 48 processor system with a large amount of primary memory. This will be a significant computer resource for this project, but the group also has access to, and experience from, several other larger parallel computers, e.g., at PDC, NSC, and Parallab in Bergen.