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

[Dave May]

On 4 February 2015 at 21:34, Fabian Gabel <gabel.fabian at gmail.com> wrote:

> Thank you for pointing me into the right direction. After some first
> tests on a test case with 2e6 cells (4dof) I could measure a slight
> improvement (25%) with respect to wall time, using a nested
> field split for the velocities:
>

Great. That's a good start. It can be made ever faster.



> -coupledsolve_pc_type fieldsplit
> -coupledsolve_pc_fieldsplit_0_fields 0,1,2
> -coupledsolve_pc_fieldsplit_1_fields 3
> -coupledsolve_pc_fieldsplit_type schur
> -coupledsolve_pc_fieldsplit_block_size 4
> -coupledsolve_fieldsplit_0_ksp_converged_reason
> -coupledsolve_fieldsplit_1_ksp_converged_reason
> -coupledsolve_fieldsplit_0_ksp_type gmres
> -coupledsolve_fieldsplit_0_pc_type fieldsplit
> -coupledsolve_fieldsplit_0_pc_fieldsplit_block_size 3
> -coupledsolve_fieldsplit_0_fieldsplit_0_pc_type ml
> -coupledsolve_fieldsplit_0_fieldsplit_1_pc_type ml
> -coupledsolve_fieldsplit_0_fieldsplit_2_pc_type ml
>
> Is it normal, that I have to explicitly specify the block size for each
> fieldsplit?
>

No. You should be able to just specify

-coupledsolve_fieldsplit_ksp_converged
-coupledsolve_fieldsplit_0_fieldsplit_pc_type ml

and same options will be applied to all splits (0,1,2).
Does this functionality not work?



> I attached the results (-converged_reason only for readability and another
> file
> solely for the output of -ksp_view). I am not sure if this result could
> be improved by modifying any of the solver options.
>

Yes, I believe they can.



>
> Are there any guidelines to follow that I could use to avoid taking wild
> guesses?
>

Sure. There are lots of papers published on how to construct robust block
preconditioners for saddle point problems arising from Navier Stokes.
I would start by looking at this book:

  Finite Elements and Fast Iterative Solvers
  Howard Elman, David Silvester and Andy Wathen
  Oxford University Press
  See chapters 6 and 8.

Your preconditioner is performing a full block LDU factorization with quite
accurate inner solves. This is probably overkill - but it will be robust.

The most relaxed approach would be :
  -coupledsolve_fieldsplit_0_ksp_type preonly
  -coupledsolve_fieldsplit_1_ksp_type preonly
  -coupledsolve_pc_fieldsplit_schur_fact_type DIAG

Something more aggressive (but less stringent that your original test)
would be:
  -coupledsolve_fieldsplit_0_ksp_type gmres
  -coupledsolve_fieldsplit_0_ksp_rtol 1.0e-2
  -coupledsolve_fieldsplit_1_ksp_type preonly
  -coupledsolve_pc_fieldsplit_schur_fact_type UPPER

When building the FS preconditioner, you can start with the absolute most
robust choices and then relax those choices to improve speed and hopefully
not destroy the convergence, or you can start with a light weight
preconditioner and make it stronger to improve convergence.
Where the balance lies is very much problem dependent.



>
> > Petsc has some support to generate approximate pressure schur
> > complements for you, but these will not be as good as the ones
> > specifically constructed for you particular discretization.
>
> I came across a tutorial (/snes/examples/tutorials/ex70.c), which shows
> 2 different approaches:
>
> 1- provide a Preconditioner \hat{S}p for the approximation of the true
> Schur complement
>
> 2- use another Matrix (in this case its the Matrix used for constructing
> the preconditioner in the former approach) as a new approximation of the
> Schur complement.
>
> Speaking in terms of the PETSc-manual p.87, looking at the factorization
> of the Schur field split preconditioner, approach 1 sets \hat{S}p while
> approach 2 furthermore sets \hat{S}. Is this correct?
>
>
No this is not correct.
\hat{S} is always constructed by PETSc as
  \hat{S} = A11 - A10 KSP(A00) A01
This is the definition of the pressure schur complement used by FieldSplit.
Note that it is inexact since the action
  y = inv(A00) x
is replaced by a Krylov solve, e.g. we solve A00 y = x for y

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



> >
> > [2] If you assembled a different operator for your preconditioner in
> > which the B_pp slot contained a pressure schur complement
> > approximation, you could use the simpler and likely more robust option
> > (assuming you know of a decent schur complement approximation for you
> > discretisation and physical problem)
>
> > -coupledsolve_pc_type fieldsplit
> > -coupledsolve_pc_fieldsplit_type MULTIPLICATIVE
> >
> > which include you U-p coupling, or just
> >
> > -coupledsolve_pc_fieldsplit_type ADDITIVE
> >
> >
> > which would define the following preconditioner
> >
> > inv(B) = diag( inv(B_uu,) , inv(B_vv) , inv(B_ww) , inv(B_pp) )
>
>
> What do you refer to with "B_pp slot"? I don't understand this approach
> completely. What would I need a Schur complement approximation for, if I
> don't use a Schur complement preconditioner?
>
>
I was referring to constructing an operator which approximates the schur
complement and inserting it into the pressure-pressure coupling block (pp
slot)



> > Option 2 would be better as your operator doesn't have an u_i-u_j, i !
> > = j coupling and you could use efficient AMG implementations for each
> > scalar terms associated with u-u, v-v, w-w coupled terms without
> > having to split again.
> >
> > Also, fieldsplit will not be aware of the fact that the Auu, Avv, Aww
> > blocks are all identical - thus it cannot do anything "smart" in order
> > to save memory. Accordingly, the KSP defined for each u,v,w split will
> > be a unique KSP object. If your A_ii are all identical and you want to
> > save memory, you could use MatNest but as Matt will always yell out,
> > "MatNest is ONLY a memory optimization and should be ONLY be used once
> > all solver exploration/testing is performed".
>
> Thanks, I will keep this in mind. Does this mean that I would only have
> to assemble one matrix for the velocities instead of three?
>

Yes, exactly.


>
> >
> >         - A_pp is defined as the matrix resulting from the
> >         discretization of the
> >         pressure equation that considers only the pressure related
> >         terms.
> >
> >
> > Hmm okay, i assumed for incompressible NS the pressure equation
> >
> > that the pressure equation would be just \div(u) = 0.
> >
>
> Indeed, many finite element(!) formulations I found while researching use
> this approach, which leads to the block A_pp being zero. I however use a
> collocated finite volume formulation and, to avoid checkerboarding of the
> pressure field, I deploy a pressure weighted interpolation method to
> approximate the velocities surging from the discretisation of \div{u}.
> This gives me an equation with the pressure as the dominant variable.
>

Ah okay, you stabilize the discretisation.
Now I understand why you have entries in the PP block.

Cheers
  Dave
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20150205/b659443d/attachment-0001.html>