# Research

My research focusses on the following main themes:

- Solution of linear systems with large sparse matrices (MRILU)
- Eigenvalue computations of large sparse matrices
- Discretization of PDEs
- Parallel and vector computing
- Fluid flow computations
- Bifurcation analysis
- Stability study of Ocean Flow Circulations

Description:

**Solution of large sparse linear systems** In CFD problems often sparse linear systems have to be solved, e.g. in explicit methods for incompressible flow a pressure Poisson equation and in implicit methods the fully coupled Navier-Stokes equations. The run time of codes for CFD problems is usually dominated by the time needed for solving such a system. We were able to construct two efficient methods for convection-diffusion equations: NGILU (Van der Ploeg et al. 1996) and MRILU (Botta &W. 1999). Both are incomplete LU factorizations of the system matrix with special orderings for the unknowns and sophisticated dropping criteria, so the product of the matrix factors L and U is not precisely the system matrix. NGILU was constructed for matrices that arise from discretizations on regular cartesian grids and MRILU can handle matrices arising from discretizations of unstructured grids since it chooses self the ordering of the unknowns based on the magnitude of the entries of the matrix. Like the almost classic multi-grid methods, both show convergence independent of the number of unknowns. So the amount of work to solve the system is linear in the number of unknowns. In benchmark problems posed for a workshop at the university of Groningen, we have shown that NGILU and MRILU combined with a Krylov subspace solver are competitive to other solvers, including multi-grid methods (Botta et al. 1997).

An interesting application of these methods are the ocean climate models as developed by IMAU. Since a decade we apply MRILU successfully to these problems (Wijer et al. 2003). Recently we worked out a special purpose method for the ocean climate model which speeded up the computations significantly; MRILU is used as a subsolver in this process.

**Eigenvalue computations** To study the stability of steady flows often a (generalized) eigenvalue problem has to be solved. One approach for this problem is the Jacobi-Davidson method developed at Utrecht University. Together with MRILU the convergence of this method can be drastically improved. Moreover, it is possible to make a good initial guess for the iteration process by constructing a small replacement of the eigenvalue problem using the MRILU factorization (Sleijpen & W. 2003)

**Well-posed discretizations for PDEs** A discretization of a PDE is called well-posed if the solution of the discretization is unique and varies continuously with the coefficients in the PDE and associated boundary and initial conditions. Moreover the discretization should be consistent to the PDE. In many cases discretization of a well-posed partial differential equation including boundary conditions is not a trivial task. The discretization may be ill-posed leading to instabilities during the numerical solution due to approximation errors and round-off errors. This already starts with second-order elliptic PDEs with suitable boundary conditions. Though by itself these are well-posed and all eigenvalues of the associated operator together with the homogeneous form of the boundary conditions are in the right half plane, the discretization may lead to a matrix which has eigenvalues left and right from zero or even may be singular. The start of a stable discretization lies in a proper choice of the formulation of the PDE. Usually one should start from a divergence form. From this form it is easy to see whether the PDE, or the part of it where we are looking at, is well-posed. Then the discretization should be performed in such a way that this property is not lost, which leads to a stable discretization. Such situations occurred in the ocean model THCM (developed at Utrecht University) for the discretization of a complex mixing term. As a numerical analyst one should also be aware of the possibility that the given PDEs are not well-posed. We encountered these during the numerical treatment of the Boussinesq equations for shallow-water flow. The original equations had to be adapted to become well-posed (Van der Veen & W. 1995). Recently, we also encountered a similar problem in an ocean model with a free surface: the free surface condition conflicted with mass conservation. Accuracy of approximation and stability of the discretization can also be conflicting requirements. In this respect we showed that for second-order accuracy in oceanography applications it is better to use a so-called Arakawa B-grid instead of a C-grid (Wubs et al. 2005) for the horizontal directions.