The numerical continuation of steady states in thermohaline ocean circulation requires the solution of large systems. The (linearized) equations describing the thermohaline ocean circulation (for a detailed description of these equations see [5]) can be splitted in equations for velocities in longitudinal and latitudinal direction (x_uv), velocity in vertical direction (x_w), pressure (x_p) and salinity and temperature (x_ST). The equations are coupled in the following way,

where G_* represents the discrete gradient operator and D_* the discrete divergence operator. We assume that G_*^T=D_*, i.e. the discrete operators are each others transpose. Furthermore, A_ST is a convection-diffusion matrix, A_uv is also a convection-diffusion matrix, but includes a Coriolis force. Respectively the four equations describe conservation of momentum (in longitudinal and latitudinal direction), the hydrostatic pressure, conservation of mass and conservation of salt and heat.

Due to the size of the systems they are solved via a Krylov subspace method using a preconditioner to speed up the convergence. In general there are two ways to obtain a preconditioner for the matrix in the equation. One could try to exploit the structure of the system or apply a general black-box preconditioning method like an incomplete LU factorization. We will do both and compare the results.

In the current code for numerical simulation of thermohaline ocean circulation (THCM, see [5], the systems are solved with the MRILU [1] preconditioner. MRILU is a multilevel ILU method containing ideas similar to algebraic multigrid methods. However MRILU doesn't need any smoothing steps, because of sophisticated dropping- and lumping criteria. In THCM MRILU is applied to the clustered equations, that is the six unknowns (u,v,w,p,T,S), that belong to one and the same cell, are treated as one unknown. The method is able to construct a good factorization and the Krylov subspace method converges fast. Unfortunately the preconditioner requires a lot of construction time and memory, which becomes a problem if the size of the problem is increased.

In order to reduce the memory usage, we developed an alternative preconditioner that exploits the structure. The ingredients are: a splitting of the pressure in a depth-averaged pressure and a component perpendicular to that; a non-symmetric reordering of the matrix; dropping the block B_ST. This together gives a block-Gauss-Seidel preconditioner, that requires some matrix-vector products and the solution of three relatively easy systems: (i) a depth-averaged saddle point system based on A_uv,G_uv and D_uv, (ii) a convection-diffusion-Coriolis equation involving A_uv and (iii) a convection-diffusion equation involving A_ST. The last two problems can be handled very well by MRILU. For the saddle point system we can apply any of the preconditioners from literature [3,4] or the ones we developed ourselves based on artificial compressibility and grad-div-stabilization [2] respectively.

A comparison of the structured preconditioner with MRILU applied to the subsystems and MRILU applied to the clustered matrix at once shows that the required amount of memory is reduced with almost a factor 10. We have to pay for that in an increase of the number of iterations in the Krylov subspace method, but this factor is at most 2. Overall we see a serious speed up.

On the conference we will describe the block-Gauss-Seidel preconditioner and show the results of a comparison using both MRILU and an algebraic multigrid method as a subsolver in the preconditioner.