Hey Max,
Without knowing anything about the specific application related to
your Stokes problem, or information about the mesh you are using, I
have a couple of questions and suggestions which might help.
The test case that I am working with is isoviscous convection, benchmark case 1a from Blankenbach 1989.
1) If A, is your stokes operator A = ( K,B ; B^T, 0 ), what is your
precondition operator?
Specifically, what is in the (2,2) slot in the precondioner? - i.e.
what matrix are you you applying -stokes_fieldsplit_1_pc_type jacobi
-stokes_fieldsplit_1_ksp_type preonly to?
Is it the identity as in the SpeedUp notes?
I think that this is the problem. The (2,2) slot in the LHS matrix is all zero (pressure does not appear in the continuity equation), so I think that the preconditioner is meaningless. I am still confused as to why this choice of preconditioner was suggested in the tutorial, and what is a better choice of preconditioner for this block? Should I be using one of the Schur complement methods instead of the additive or multiplicative field split?
2) This choice
-stokes_fieldsplit_0_pc_type ml -stokes_fieldsplit_0_ksp_type preonly
may simply not be a very effective and degrade the performance of the
outer solver.
I'd make the solver for the operator in the (1,1) slot much stronger,
for example
-stokes_fieldsplit_0_ksp_type gmres
-stokes_fieldsplit_0_ksp_rtol 1.0e-4
-stokes_fieldsplit_0_mg_levels_ksp_type gmres
-stokes_fieldsplit_0_mg_levels_pc_type bjacobi
-stokes_fieldsplit_0_mg_levels_ksp_max_it 4
Add a monitor on this solver (-stokes_fieldsplit_0_ksp_XXX) to see how
ML is doing.
3) Using -stokes_pc_fieldsplit_type MULTIPLICATIVE should reduce the
number of outer iterations by a factor of two, but it will use more
memory.
4) You should use a flexible Krylov method on the outer most solve
(-stokes_ksp_XXX) as the preconditioner is varying between each outer
iteration. Use -stokes_ksp_type fgmres or -stokes_ksp_type gcr
Thanks for pointing this out. I made that change.
5) Depending on how the physical problem is scaled
(non-dimensionalised), the size of the residuals associated with the
momentum and continuity equation make be quite different. You are
currently use the entire residual from (u,p) to determine when to stop
iterating. You might want to consider writing a monitor which examines
the these residuals independently.
I think that I have scaled the problem correctly. I (slowly) obtain a sufficiently accurate solution using as options only:
-stokes_ksp_atol 1e-5 -stokes_ksp_rtol 1e-5 -stokes_ksp_monitor_true_residual -stokes_ksp_norm_type UNPRECONDITIONED
Cheers,
Dave