[petsc-users] Field Split PC for Fully-Coupled 3d stationary incompressible Navier-Stokes Solution Algorithm

[Fabian Gabel]

After some further testing I found, that the key to reduce the wall time
is to balance the number of iterations needed to apply the
preconditioner and the number of inner iterations of my solver for the
complete linear system (velocities and pressure). I found that the
following set of options led to a further improvement of wall time from
6083s to 2200s:

-coupledsolve_fieldsplit_0_fieldsplit_ksp_type preonly
-coupledsolve_fieldsplit_0_fieldsplit_pc_type ml
-coupledsolve_fieldsplit_0_ksp_rtol 1e-2 
-coupledsolve_fieldsplit_0_ksp_type gmres
-coupledsolve_fieldsplit_0_pc_fieldsplit_block_size 3
-coupledsolve_fieldsplit_0_pc_type fieldsplit
-coupledsolve_fieldsplit_1_ksp_rtol 1e-1 
-coupledsolve_fieldsplit_1_ksp_type gmres
-coupledsolve_fieldsplit_1_pc_type ilu
-coupledsolve_fieldsplit_ksp_converged_reason
-coupledsolve_fieldsplit_schur_precondition a11
-coupledsolve_ksp_type fgmres
-coupledsolve_pc_fieldsplit_0_fields 0,1,2
-coupledsolve_pc_fieldsplit_1_fields 3
-coupledsolve_pc_fieldsplit_block_size 4
-coupledsolve_pc_fieldsplit_type schur
-coupledsolve_pc_type fieldsplit
-coupledsolve_pc_fieldsplit_schur_fact_type lower

I will also attach the complete output at the end of my mail. 

On Fr, 2015-02-06 at 08:35 +0100, Dave May wrote:
>  
>         -coupledsolve_pc_fieldsplit_block_size 4
>         -coupledsolve_fieldsplit_0_pc_fieldsplit_block_size 3
>         
>         Without them, I get the error message
>         
>         [0]PETSC ERROR: PCFieldSplitSetDefaults() line 468
>         in /work/build/petsc/src/ksp/pc/impls/fieldsplit/fieldsplit.c
>         Unhandled
>         case, must have at least two fields, not 1
>         
>         I thought PETSc would already know, what I want to do, since I
>         initialized the fieldsplit with
>         
>         CALL PCFieldSplitSetIS(PRECON,PETSC_NULL_CHARACTER,ISU,IERR)
>         
> 
> 
> Huh. 
> If you called PCFieldSplitSetIS(), then certainly the FS knows you
> have four fields and you shouldn't need to set the block size 4
> option. I know this is true as I use it all the time. When you start
> grouping new splits (in your case u,v,w) I would have also thought
> that the block size 3 option would also be redundant - however I use
> this option less frequently.
> 
> 
> 
>  
>         As a matter of fact I spent the last days digging through
>         papers on the
>         regard of preconditioners or approximate Schur complements and
>         the names
>         Elman and Silvester have come up quite often.
>         
>         The problem I experience is, that, except for one publication,
>         all the
>         other ones I checked deal with finite element formulations.
>         Only
>         
>         Klaij, C. and Vuik, C. SIMPLE-type preconditioners for
>         cell-centered,
>         colocated finite volume discretization of incompressible
>         Reynolds-averaged Navier–Stokes equations
>         
>         presented an approach for finite volume methods. 
> 
> 
> The exact same analysis applies directly to any stable mixed u-p
> discretization (I agree with all Matt's comments as well). If your
> stablilization term is doing what is supposed to, then your
> discretization should be stable. 
> I don't know of any papers which show results using your exact
> discretization, but here are some preconditioning papers for Stokes
> which employ FV discretizations:
> 
> @article{olshanskii1999iterative,
>   title={An iterative solver for the Oseen problem and numerical solution of incompressible Navier--Stokes equations},
>   author={Olshanskii, Maxim A},
>   journal={Numerical linear algebra with applications},
>   volume={6},
>   number={5},
>   pages={353--378},
>   year={1999}
> }
> 
> @article{olshanskii2004grad,
>   title={Grad-div stablilization for Stokes equations},
>   author={Olshanskii, Maxim and Reusken, Arnold},
>   journal={Mathematics of Computation},
>   volume={73},
>   number={248},
>   pages={1699--1718},
>   year={2004}
> }
> @article{griffith2009accurate,
>   title={An accurate and efficient method for the incompressible Navier--Stokes equations using the projection method as a preconditioner},
>   author={Griffith, Boyce E},
>   journal={Journal of Computational Physics},
>   volume={228},
>   number={20},
>   pages={7565--7595},
>   year={2009},
>   publisher={Elsevier}
> }
> @article{furuichi2011development,
>   title={Development of a Stokes flow solver robust to large viscosity jumps using a Schur complement approach with mixed precision arithmetic},
>   author={Furuichi, Mikito and May, Dave A and Tackley, Paul J},
>   journal={Journal of Computational Physics},
>   volume={230},
>   number={24},
>   pages={8835--8851},
>   year={2011},
>   publisher={Elsevier}
> }
> 
> @article{cai2013efficient,
>   title={Efficient variable-coefficient finite-volume Stokes solvers},
>   author={Cai, Mingchao and Nonaka, AJ and Bell, John B and Griffith, Boyce E and Donev, Aleksandar},
>   journal={arXiv preprint arXiv:1308.4605},
>   year={2013}
> }
> 
> 
>         Furthermore, a lot of
>         literature is found on saddle point problems, since the linear
>         system
>         from stable finite element formulations comes with a 0 block
>         as pressure
>         matrix. I'm not sure how I can benefit from the work that has
>         already
>         been done for finite element methods, since I neither use
>         finite
>         elements nor I am trying to solve a saddle point problem (?).
> 
> 
> I would say you are trying to solve a saddle point system, only one
> which has been stabilized. I expect your stabilization term should
> vanish in the limit of h -> 0. The block preconditioners are directly
> applicable to what you are doing, as are all the issues associated
> with building preconditioners for schur complement approximations. I
> have used FS based preconditioners for stablized Q1-Q1 finite element
> discretizations for Stokes problems. Despite the stabilization term in
> the p-p coupling term, saddle point preconditioning techniques are
> appropriate. There are examples of this in
> src/ksp/ksp/examples/tutorials - see ex43.c ex42.c
> 
>  
> 
>  
>         But then what happens in this line from the
>         tutorial /snes/examples/tutorials/ex70.c
>         
>         ierr = KSPSetOperators(subksp[1], s->myS,
>         s->myS);CHKERRQ(ierr);
>         
>         It think the approximate Schur complement a (Matrix of type
>         Schur) gets
>         replaced by an explicitely formed Matrix (myS, of type
>         MPIAIJ).
> 
> 
> Oh yes, you are right, that is what is done in this example. But you
> should note that this is not the default way petsc's fieldsplit
> preconditioner will define the schur complement \hat{S}.  This
> particular piece of code lives in the users example. 
> 
> If you really wanted to do this, the same thing could be configured on
> the command line:
> 
>   -XXX_ksp_type preonly -XXX_pc_type ksp
> 
> 
>  
>         >
>         > You have two choices in how to define the preconditioned,
>         \hat{S_p}:
>         >
>         > [1] Assemble you own matrix (as is done in ex70)
>         >
>         > [2] Let PETSc build one. PETSc does this according to
>         >
>         >   \hat{S_p} = A11 - A10 inv(diag(A00)) A01
>         >
>         Regards,
>         Fabian
>         >
>         >
>         
>         
> 
>