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

Fri Jul 26 10:42:35 CDT 2013

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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20130726/7eee271a/attachment.html>