[petsc-users] discontinuous viscosity stokes equation 3D staggered grid

Matthew Knepley knepley at gmail.com
Sat Jul 27 07:14:28 CDT 2013 

On Sat, Jul 27, 2013 at 7:03 AM, Anton Popov <popov at uni-mainz.de> wrote:

>  Dave,
>
> coupled staggered grid possesses a non-trivial null space. That means it
> probably contains not just three translations, three rotations and a
> constant pressure. It would be nice if somebody finally figures out how it
> looks like. Any ideas? This information is valuable for building a proper
> coupled algebraic multigrid preconditioner.
>
> As for the fake dofs, I really don't like them. It's perfectly possible to
> design a numbering scheme that would only contain real dofs, and implement
> it with PETSc. I'm working on it.
>

My recommended way to do this is to use PetscSection() to define the dofs
over a DMDA. This currently works correctly in parallel
for vertex and cell dofs. The code for face dofs is in there, but I have
never tested it. You have to use the 'next' branch because this
is experimental stuff.

  Thanks,

      Matt


>
> Anton
>
>
> On 7/26/13 6:05 PM, Dave May wrote:
>
>  Yes, the nullspace is important.
>  There is definitely a pressure null space of constants which needs to be
> considered.
>  I don't believe ONLY using
>   -fieldsplit_1_ksp_constant_null_space
>  is not sufficient.
>
> You need to inform the KSP for the outer syetem (the coupled u-p system)
> that there is a null space of constants in the pressure system. This cannot
> (as far as I'm aware) be set via command line args. You need to write a
> null space removal function which accepts the complete (u,p) vector which
> only modifies the pressure dofs.
>
>  Take care though when you write this though. Because of the DMDA
> formulation you are using, you have introduced dummy/ficticious pressure
> dofs on the right/top/front faces of your mesh. Thus your null space
> removal function should ignore those dofs when your define the constant
> pressure constraint.
>
>
>
> Cheers,
>    Dave
>
>
> On 26 July 2013 17:42, Matthew Knepley <knepley at gmail.com> wrote:
>
>>  On Fri, Jul 26, 2013 at 10:13 AM, Bishesh Khanal <bisheshkh at gmail.com>wrote:
>>
>>>
>>>
>>>
>>>  On Fri, Jul 26, 2013 at 4:22 PM, Matthew Knepley <knepley at gmail.com>wrote:
>>>
>>>>  On Fri, Jul 26, 2013 at 9:11 AM, Bishesh Khanal <bisheshkh at gmail.com>wrote:
>>>>
>>>>>
>>>>>
>>>>>
>>>>>  On Fri, Jul 26, 2013 at 2:32 PM, Matthew Knepley <knepley at gmail.com>wrote:
>>>>>
>>>>>>  On Fri, Jul 26, 2013 at 7:28 AM, Bishesh Khanal <bisheshkh at gmail.com
>>>>>> > wrote:
>>>>>>
>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>>  On Wed, Jul 17, 2013 at 9:48 PM, Jed Brown <jedbrown at mcs.anl.gov>wrote:
>>>>>>>
>>>>>>>> Bishesh Khanal <bisheshkh at gmail.com> writes:
>>>>>>>>
>>>>>>>> > Now, I implemented two different approaches, each for both 2D and
>>>>>>>> 3D, in
>>>>>>>> > MATLAB. It works for the smaller sizes but I have problems
>>>>>>>> solving it for
>>>>>>>> > the problem size I need (250^3 grid size).
>>>>>>>> > I use staggered grid with p on cell centers, and components of v
>>>>>>>> on cell
>>>>>>>> > faces. Similar split up of K to cell center and faces to account
>>>>>>>> for the
>>>>>>>> > variable viscosity case)
>>>>>>>>
>>>>>>>>  Okay, you're using a staggered-grid finite difference
>>>>>>>> discretization of
>>>>>>>> variable-viscosity Stokes.  This is a common problem and I recommend
>>>>>>>> starting with PCFieldSplit with Schur complement reduction (make
>>>>>>>> that
>>>>>>>> work first, then switch to block preconditioner).
>>>>>>>
>>>>>>>
>>>>>>>  Ok, I made my 3D problem work with PCFieldSplit with Schur
>>>>>>> complement reduction using the options:
>>>>>>> -pc_fieldsplit_type schur -pc_fieldsplit_detect_saddle_point
>>>>>>> -fieldsplit_1_ksp_constant_null_space
>>>>>>>
>>>>>>>
>>>>>>>  You can use PCLSC or
>>>>>>>> (probably better for you), assemble a preconditioning matrix
>>>>>>>> containing
>>>>>>>> the inverse viscosity in the pressure-pressure block.  This diagonal
>>>>>>>> matrix is a spectrally equivalent (or nearly so, depending on
>>>>>>>> discretization) approximation of the Schur complement.  The velocity
>>>>>>>> block can be solved with algebraic multigrid.  Read the PCFieldSplit
>>>>>>>> docs (follow papers as appropriate) and let us know if you get
>>>>>>>> stuck.
>>>>>>>>
>>>>>>>
>>>>>>>  Now,  I got a little confused in how exactly to use command line
>>>>>>> options to use multigrid for the velocity bock and PCLS for the pressure
>>>>>>> block. After going through the manual I tried the following:
>>>>>>>
>>>>>>
>>>>>>  You want Algebraic Multigrid
>>>>>>
>>>>>>  -pc_type fieldsplit -pc_fieldsplit_detect_saddle_point
>>>>>> -pc_fieldsplit_type schur
>>>>>>   -fieldsplit_0_pc_type gamg
>>>>>>    -fieldsplit_1_ksp_type fgmres
>>>>>>    -fieldsplit_1_ksp_constant_null_space
>>>>>>    -fieldsplit_1_ksp_monitor_short
>>>>>>    -fieldsplit_1_pc_type lsc
>>>>>> -ksp_converged_reason
>>>>>>
>>>>>>    I tried the above set of options but the solution I get seem to
>>>>> be not correct. The velocity field I get are quite different than the one I
>>>>> got before without using gamg which were the expected one.
>>>>> Note: (Also, I had to add one extra option of
>>>>> -fieldsplit_1_ksp_gmres_restart 100 , because the fieldsplit_1_ksp residual
>>>>> norm did not converge within default 30 iterations before restarting).
>>>>>
>>>>
>>>>   These are all iterative solvers. You have to make sure everything
>>>> converges.
>>>>
>>>
>>>  When I set restart to 100, and do -ksp_monitor, it does converge (for
>>> the fieldsplit_1_ksp). Are you saying that in spite of having
>>> -ksp_converged_reason in the option and petsc completing the run with the
>>> message "Linear solve converged due to CONVERGED_RTOL .." not enough to
>>> make sure that everything is converging ? If that is the case what should I
>>> do for this particular problem ?
>>>
>>
>>   If your outer iteration converges, and you do not like the solution,
>> there are usually two possibilities:
>>
>>    1) Your tolerance is too high, start with it cranked down all the way
>> (1e-10) and slowly relax it
>>
>>    2) You have a null space that you are not accounting for
>>
>>    I have used the MAC scheme with indexing as shown in:  fig 7.5, page
>>> 96 of:
>>> http://books.google.co.uk/books?id=W83gxp165SkC&printsec=frontcover&dq=Introduction+to+Numerical+Geodynamic+Modelling&hl=en&sa=X&ei=v6TmUaP_L4PuOs3agJgE&ved=0CDIQ6AEwAA
>>>
>>> Thus I have a DM with 4 dof but there are several "ghost values" set to
>>> 0.
>>>  Would this cause any problem when using the multigrid ? (This has
>>> worked fine when not using the multigrid.)
>>>
>>
>>  I don't know exactly how you have implemented this. These should be
>> rows of the identity.
>>
>>
>>>     I do this problem in the tutorial with
>>>> a constant viscosity.
>>>>
>>>
>>>  Which tutorial are you referring to ? Could you please provide me the
>>> link please ?
>>>
>>
>>   There are a few on the PETSc Tutorials page, but you can look at this
>>
>>
>> http://www.geodynamics.org/cig/community/workinggroups/short/workshops/cdm2013/presentations/SessionIV_Solvers.pdf
>>
>>  for a step-by-step example of a saddle-point problem at the end.
>>
>>      Matt
>>
>>
>>>
>>>
>>>>      Matt
>>>>
>>>>
>>>>>          Matt
>>>>>>
>>>>>>
>>>>>>>  but I get the following errror:
>>>>>>>
>>>>>>> [0]PETSC ERROR: --------------------- Error Message
>>>>>>> ------------------------------------
>>>>>>> [0]PETSC ERROR: Null argument, when expecting valid pointer!
>>>>>>> [0]PETSC ERROR: Null Object: Parameter # 2!
>>>>>>> [0]PETSC ERROR:
>>>>>>> ------------------------------------------------------------------------
>>>>>>> [0]PETSC ERROR: Petsc Release Version 3.4.2, Jul, 02, 2013
>>>>>>> [0]PETSC ERROR: See docs/changes/index.html for recent updates.
>>>>>>> [0]PETSC ERROR: See docs/faq.html for hints about trouble shooting.
>>>>>>> [0]PETSC ERROR: See docs/index.html for manual pages.
>>>>>>> [0]PETSC ERROR:
>>>>>>> ------------------------------------------------------------------------
>>>>>>> [0]PETSC ERROR: src/AdLemMain on a arch-linux2-cxx-debug named
>>>>>>> edwards by bkhanal Fri Jul 26 14:23:40 2013
>>>>>>> [0]PETSC ERROR: Libraries linked from
>>>>>>> /home/bkhanal/Documents/softwares/petsc-3.4.2/arch-linux2-cxx-debug/lib
>>>>>>> [0]PETSC ERROR: Configure run at Fri Jul 19 14:25:01 2013
>>>>>>> [0]PETSC ERROR: Configure options --with-cc=gcc --with-fc=g77
>>>>>>> --with-cxx=g++ --download-f-blas-lapack=1 --download-mpich=1
>>>>>>> -with-clanguage=cxx --download-hypre=1
>>>>>>> [0]PETSC ERROR:
>>>>>>> ------------------------------------------------------------------------
>>>>>>> [0]PETSC ERROR: MatPtAP() line 8166 in
>>>>>>> /home/bkhanal/Documents/softwares/petsc-3.4.2/src/mat/interface/matrix.c
>>>>>>> [0]PETSC ERROR: PCSetUp_MG() line 628 in
>>>>>>> /home/bkhanal/Documents/softwares/petsc-3.4.2/src/ksp/pc/impls/mg/mg.c
>>>>>>> [0]PETSC ERROR: PCSetUp() line 890 in
>>>>>>> /home/bkhanal/Documents/softwares/petsc-3.4.2/src/ksp/pc/interface/precon.c
>>>>>>> [0]PETSC ERROR: KSPSetUp() line 278 in
>>>>>>> /home/bkhanal/Documents/softwares/petsc-3.4.2/src/ksp/ksp/interface/itfunc.c
>>>>>>> [0]PETSC ERROR: KSPSolve() line 399 in
>>>>>>> /home/bkhanal/Documents/softwares/petsc-3.4.2/src/ksp/ksp/interface/itfunc.c
>>>>>>> [0]PETSC ERROR: PCApply_FieldSplit_Schur() line 807 in
>>>>>>> /home/bkhanal/Documents/softwares/petsc-3.4.2/src/ksp/pc/impls/fieldsplit/fieldsplit.c
>>>>>>> [0]PETSC ERROR: PCApply() line 442 in
>>>>>>> /home/bkhanal/Documents/softwares/petsc-3.4.2/src/ksp/pc/interface/precon.c
>>>>>>> [0]PETSC ERROR: KSP_PCApply() line 227 in
>>>>>>> /home/bkhanal/Documents/softwares/petsc-3.4.2/include/petsc-private/kspimpl.h
>>>>>>> [0]PETSC ERROR: KSPInitialResidual() line 64 in
>>>>>>> /home/bkhanal/Documents/softwares/petsc-3.4.2/src/ksp/ksp/interface/itres.c
>>>>>>> [0]PETSC ERROR: KSPSolve_GMRES() line 239 in
>>>>>>> /home/bkhanal/Documents/softwares/petsc-3.4.2/src/ksp/ksp/impls/gmres/gmres.c
>>>>>>> [0]PETSC ERROR: KSPSolve() line 441 in
>>>>>>> /home/bkhanal/Documents/softwares/petsc-3.4.2/src/ksp/ksp/interface/itfunc.c
>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>
>>>>>>
>>>>>>  --
>>>>>> 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
>>>>>>
>>>>>
>>>>>
>>>>
>>>>
>>>>  --
>>>> 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
>>>>
>>>
>>>
>>
>>
>>  --
>> 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
>>
>
>
>


-- 
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/20130727/8cf94329/attachment-0001.html>

More information about the petsc-users mailing list