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

[Matthew Knepley]

On Thu, Feb 5, 2015 at 9:32 AM, Dave May <dave.mayhem23 at gmail.com> wrote:

>
>
> 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
>

There are a few options here:


http://www.mcs.anl.gov/petsc/petsc-master/docs/manualpages/PC/PCFieldSplitSetSchurPre.html

  Thanks,

     Matt


>
>
>
>> >
>> > [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
>
>
>


-- 
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20150205/77cc87ad/attachment.html>