Some aspects of your questions come up frequenty, especially the unrealistic expectations from black-box iterative solvers. Hopefully this somewhat in-depth reply will be useful.
On Thu, 27 May 2010 09:20:55 +0200, Luca Mangani <luca.mangani at hslu.ch> wrote: > Good Morning, > > After many steps I'm not able to solve my coupled matrix with > iterative solver. The matrix is defined as the discretization of the > Navier Stokes equation in a pressure based form. This form is indefinite. While we often talk about there being no scalable black-box solvers, the statement is especially true for indefinite systems. As an example of the spectacular failure of black-box methods for elliptic systems when applied to other problems, I have observed BoomerAMG being massively singular when applied to a mixed finite element discretization of incompressible flow. This resulted in preconditioned residuals falling in essentially textbook manner, but the unpreconditioned residuals did not converge at all (after satisfying the boundary conditions). It is very important to check this. For *some* incompressible flow problems, it is possible to order unknowns such that incomplete factorization causes useful fill into the pressure-pressure block. Unfortunately, incomplete factorization has very unpredictable properties and is certainly not a robust solution. When playing with orderings, usually you would like all velocity unknowns in a subdomain to come before the pressure unknowns, a number of automatic orderings are also available with -pc_factor_mat_ordering_type. You may also try changing the shift behavior, see -pc_factor_shift_type, using an ILU with pivoting, or using sparse approximate inverse -pc_hypre_type parasails or -pc_type spai. While some of these black-box approaches will be better (worse) than others, and depending on your needs, perhaps even "good enough", none will scale to large problem size or be robust to parameters. They may all fail completely. There are essentially two scalable approaches: 1. Schur complement methods This class contains classical methods like Uzawa and SIMPLE, as well as newer approaches based on approximate commutators. See Elman's 2008 paper "A taxonomy and comparison of parallel block multi-level preconditioners for the incompressible Navier-Stokes equations" for a current overview. These methods break the problem into definite operators for which conventional preconditioners are suitable. It is relatively easy to implement using PCFieldSplit, but note that the approximation to the Schur complement is a crucial choice and involves problem-specific information. 2. Coupled multilevel methods This class contains multigrid methods based on Vanka smoothers as well as domain decomposition methods. It is crucial to define subdomains, interpolants, and smoothers/subdomain solvers that are compatible with the inf-sup condition. These stability requirements are highly dependent on your spatial discretization, therefore it really cannot be encapsulated in a black-box solver. For the multigrid end of the spectrum, see papers by Brandt in the finite difference context or Rannacher and Turek for certain finite element discretizations. At the DD end, you might like papers by Jing Li for the non-overlapping case or Xiao-Chuan Cai for overlapping cases. Note that 1-level Schwarz methods essentially only require a choice of stable subdomains which you can do with PCASMSetLocalSubdomains(). Many other methods can be posed as a multilevel Schwarz method in which case you can use PCMG and define your own problem-specific interpolants. There are currently no scalable algorithms with convergence rates that are independent of Reynolds number, the best cases with the latter class (coupled multigrid or DD) seem to be somewhat less sensitive than the Schur-complement methods. As usual, there are difficulties in the presence of large variation in coefficients or high aspect ratio, though some methods are indpendent or very weakly dependent on these parameters. The methods in class 1 almost always involve more frequent communication, therefore their scalability is likely to be somewhat worse than class 2, though some people have demonstrated quite good weak scalability using class 1. Benzi, Golub, and Liesen's 2005 Acta Numerica review paper may also be of interest. More things are likely to work with non-mixed discretizations, but those discretizations are generally less robust. A group at Sandia has had some luck applying ML to stabilized Q1-Q1 finite element methods, it is important to keep all fields coupled (PETSc will inform ML correctly if you order unknowns [u0,v0,w0,p0,u1,...] and call MatSetBlockSize or use BAIJ, note that this ordering is better for memory performance and smoother strength anyway). Their experience is that this approach is usually faster than approximate commutators, but sometimes fails to converge at very moderate parameter choices. Note that this Q1-Q1 discretization does not exactly enforce incompressibility and this divergence error can be quite large on under-resolved meshes or in the presence of non-smooth boundary conditions. Galerkin least squares methods offer a way to produce a positive definite system from a velocity-pressure formulation, therefore standard solvers can be applied. But note that it generally commits very large divergence errors (see 99% mass loss a few years ago, the newer methods can keep it to a few percent for simple problems), this can be circumvented by choosing inf-sup stable spaces, but then the solver issues become somewhat more complicated again. Once you have decided on a discretization and algorithm from the literature that you would like to implement (or your own variant), we can give suggestions. Please recognize that scalable solvers for incompressible flow is very much an active research topic, you cannot expect black-box solvers to produce good results. If you do not need very large problem sizes, you might want to stick with a direct solver like MUMPS; if do not need extreme parameter choices, you might get by with naive 1-level domain decomposition or coupled algebraic multigrid (especially with collocated discretization). Jed
