[petsc-users] petsc-users Digest, Vol 55, Issue 70
Bishesh Khanal
bisheshkh at gmail.com
Mon Jul 29 12:05:12 CDT 2013
On Mon, Jul 29, 2013 at 6:07 PM, <petsc-users-request at mcs.anl.gov> wrote:
> Send petsc-users mailing list submissions to
> petsc-users at mcs.anl.gov
>
> To subscribe or unsubscribe via the World Wide Web, visit
> https://lists.mcs.anl.gov/mailman/listinfo/petsc-users
> or, via email, send a message with subject or body 'help' to
> petsc-users-request at mcs.anl.gov
>
> You can reach the person managing the list at
> petsc-users-owner at mcs.anl.gov
>
> When replying, please edit your Subject line so it is more specific
> than "Re: Contents of petsc-users digest..."
>
>
> Today's Topics:
>
> 1. Re: discontinuous viscosity stokes equation 3D staggered
> grid (Bishesh Khanal)
> 2. Re: discontinuous viscosity stokes equation 3D staggered
> grid (Anton Popov)
>
>
> ----------------------------------------------------------------------
>
> Message: 1
> Date: Mon, 29 Jul 2013 18:07:00 +0200
> From: Bishesh Khanal <bisheshkh at gmail.com>
> To: Dave May <dave.mayhem23 at gmail.com>, Jed Brown
> <jedbrown at mcs.anl.gov>
> Cc: PETSc users list <petsc-users at mcs.anl.gov>
> Subject: Re: [petsc-users] discontinuous viscosity stokes equation 3D
> staggered grid
> Message-ID:
> <
> CAEhex8icJ1wqR+Eef74R+_b0vDUGP3jMDay72rpUFhOOvC6NbQ at mail.gmail.com>
> Content-Type: text/plain; charset="iso-8859-1"
>
> On Mon, Jul 29, 2013 at 3:42 PM, Bishesh Khanal <bisheshkh at gmail.com>
> wrote:
>
> >
> >
> >
> > On Fri, Jul 26, 2013 at 6:05 PM, Dave May <dave.mayhem23 at gmail.com>
> 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.
> >>
> >>
> > I'm aware of the constant pressure null space but did not realize that
> > just having -fieldsplit_1_ksp_constant_null_space would work, thanks.
> Now I
> > tried doing what you said but I get some problems!
> > I created a null basis vector say, mNullBasis (creating a global vector
> > using the DMDA formulation). I set all it's values to zero except the
> > non-dummy pressure dofs whose values are set to one. Now I use
> > MatNullSpaceCreate to create a MatNullSpace say mNullSpace using the
> > mNullBasis vector.
> >
> > The relevant lines in my solver function are (mKsp and mDa are my
> relevant
> > KSP and DM objects):
> >
> > ierr =
> > DMKSPSetComputeRHS(mDa,computeRHSTaras3D,this);CHKERRQ(ierr);
> > ierr =
> > DMKSPSetComputeOperators(mDa,computeMatrixTaras3D,this);CHKERRQ(ierr);
> > ierr =
> > KSPSetDM(mKsp,mDa);CHKERRQ(ierr);
> >
> > ierr = KSPSetNullSpace(mKsp,mNullSpace);CHKERRQ(ierr);
> > ierr = KSPSetFromOptions(mKsp);CHKERRQ(ierr);
> > ierr = KSPSetUp(mKsp);CHKERRQ(ierr);
> > ierr = KSPSolve(mKsp,NULL,NULL);CHKERRQ(ierr);
> >
> > And a corresponding addition in the function computeMatrixTaras3D is
> > something equivalent to:
> > ierr = MatNullSpaceRemove(mNullSpace,b,NULL);CHKERRQ(ierr); //where
> b
> > is the vector that got updated with the rhs values before this line.
> >
> > When using the -ksp_view in the option I see that a null space now gets
> > attached to the outer solver, while I still need to keep
> > -fieldsplit_1_ksp_constant_null_space option as it is to have a null
> space
> > attached to the inner pressure block solver. However the problem is I do
> > not get the expected result with gamg and resulting solution blows up.
> >
>
> Never mind, I found out that I had forgotten to normalize the mNullBasis to
> form the basis of the null space.
>
>
> > Am I doing something wrong or completely missing some point here with
> > nullspace removal idea ?
> > (Note that I am imposing a Dirichlet boundary condition, so in my case I
> > think the only null space is indeed the constant pressure and not the
> > translation and rigid body rotation which would be present if using for
> > e.g. Neumann boundary condition)
> >
>
> So now I think the PCFieldsplit with the Algebraic multigrid on the
> velocity block works fine for the constant viscosity case. Now I'll look
> for the variable viscosity case. I'll have a look at Jed's suggestions for
> this in the beginning of the thread and come back with few more questions,
> thanks.
>
> >
> >
> >
> >
> >
> >
> >
> >
> >
> >>
> >>
> >> 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
> >>>
> >>
> >>
> >
> -------------- next part --------------
> An HTML attachment was scrubbed...
> URL: <
> http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20130729/0935faa2/attachment-0001.html
> >
>
> ------------------------------
>
> Message: 2
> Date: Mon, 29 Jul 2013 18:06:53 +0200
> From: Anton Popov <popov at uni-mainz.de>
> To: <petsc-users at mcs.anl.gov>
> Subject: Re: [petsc-users] discontinuous viscosity stokes equation 3D
> staggered grid
> Message-ID: <51F6931D.4070700 at uni-mainz.de>
> Content-Type: text/plain; charset="iso-8859-1"; Format="flowed"
>
> On 7/29/13 3:42 PM, Bishesh Khanal wrote:
> >
> >
> >
> > On Fri, Jul 26, 2013 at 6:05 PM, Dave May <dave.mayhem23 at gmail.com
> > <mailto:dave.mayhem23 at gmail.com>> 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.
> >
> >
> > I'm aware of the constant pressure null space but did not realize that
> > just having -fieldsplit_1_ksp_constant_null_space would work, thanks.
> > Now I tried doing what you said but I get some problems!
> > I created a null basis vector say, mNullBasis (creating a global
> > vector using the DMDA formulation). I set all it's values to zero
> > except the non-dummy pressure dofs whose values are set to one. Now I
> > use MatNullSpaceCreate to create a MatNullSpace say mNullSpace using
> > the mNullBasis vector.
> >
> > The relevant lines in my solver function are (mKsp and mDa are my
> > relevant KSP and DM objects):
> >
> > ierr = DMKSPSetComputeRHS(mDa,computeRHSTaras3D,this);CHKERRQ(ierr);
> > ierr =
> > DMKSPSetComputeOperators(mDa,computeMatrixTaras3D,this);CHKERRQ(ierr);
> > ierr = KSPSetDM(mKsp,mDa);CHKERRQ(ierr);
> > ierr = KSPSetNullSpace(mKsp,mNullSpace);CHKERRQ(ierr);
> > ierr = KSPSetFromOptions(mKsp);CHKERRQ(ierr);
> > ierr = KSPSetUp(mKsp);CHKERRQ(ierr);
> > ierr = KSPSolve(mKsp,NULL,NULL);CHKERRQ(ierr);
> >
> > And a corresponding addition in the function computeMatrixTaras3D is
> > something equivalent to:
> > ierr = MatNullSpaceRemove(mNullSpace,b,NULL);CHKERRQ(ierr);
> > //where b is the vector that got updated with the rhs values before
> > this line.
> >
> > When using the -ksp_view in the option I see that a null space now
> > gets attached to the outer solver, while I still need to keep
> > -fieldsplit_1_ksp_constant_null_space option as it is to have a null
> > space attached to the inner pressure block solver. However the problem
> > is I do not get the expected result with gamg and resulting solution
> > blows up.
> > Am I doing something wrong or completely missing some point here with
> > nullspace removal idea ?
> > (Note that I am imposing a Dirichlet boundary condition, so in my case
> > I think the only null space is indeed the constant pressure and not
> > the translation and rigid body rotation which would be present if
> > using for e.g. Neumann boundary condition)
> >
> Bishesh,
>
> yes, pressure is the only filed that is defined up to a constant in the
> incompressible case, if you properly set the velocity Dirichlet BC.
>
> What I was saying is that you need to know the full null space of your
> linear system matrix if you want to setup an efficient algebraic
> multigrid preconditioner. (Near) null space is a starting point for the
> tentative prolongator. I don't know how it looks like for staggered
> grid, and don't know how to make PETSc use this information (the latter
> should be easy I suspect). Geometric multigrid for staggered grid is
> also not a straightforward option, because you need to define and pass
> to PETSc a custom interpolation-restriction operator with
> PCMGSetRestriction or PCMGSetInterpolation . Without efficient multigrid
> you code will be severely limited for large-scale 3D problems
> (especially with variable viscosity, if this is your goal).
>
I do need to solve the system for a variable viscosity case. Yes, the null
space does not seem to be trivial when we have the staggered grid case.
>
> Which problem with variable viscosity are you solving?
I am working in modelling some biological phenomenon and with certain
assumptions I ended up deriving a set of equations which happened to be
this stokes equation except for a non-zero divergence. (This led me to this
whole "new" world of {fluid dynamics, mechanics, solving pde's numerically,
large-scale linear systems}!! :) )
Actually I do not assume the material to be incompressible, hence the
non-zero divergence field. But it must satisfy this "compressibility
constraint" of:
div(v) = f2, everywhere in the domain strictly.
It basically results from the constrained minimization of some sort of
linear elastic energy, and for my problem it is important that it satisfies
the "compressibility constraint" everywhere. And this deforming object has
different parts with very different viscosity, (or actually they are Lame's
parameters if you look at it from elastic material point of view).
To my opinion the
> easiest way to remove the pressure null space in your case is to add
> minor compressibility. In many cases you can easily get rid of strict
> incompressibility assumption, because everything is compressible (at
> least slightly).
>
> If you still want to project out constant pressure, then Dave is right,
> -fieldsplit_1_ksp_constant_null_space will most likely not work, because
> of dummy DOFs. Another complication is that you should define the null
> space only for the pressure. You probably should retrieve KSP for the
> pressure Schur complement using PCFieldSplitGetSubKSP after setup of
> fieldsplit preconditioner, and try to set your custom null space
> directly for it. I'm not sure however if this will work out.
It seems to work for me (Only for the constant viscosity) when I create a
global vector using the DMDA I have and set to non-zero constant the
components corresponding to only those non-dummy pressure dofs. Then I
normalize it and use as the nullspace basis for the MatNullSpace context.
This I use in KSPSetNullSpace() in the outermost KSP solver, then I use
MatNullSpaceRemove() in the rhs computing function. Then for the pressure
Schur complement, I could simply use the run-time option
fieldsplit_1_ksp_constant_null_space.
> We have a
> version of staggered grid solver without dummy pressure DOFs. For us it
> works both with and without ..._ksp_constant_null_space for pressure.
>
Does your version work with strict incompressibilty constraint ? Do you use
your own specialized multigrid preconditioner or the petsc provided one for
this case ?
And do you use the PetscSection as suggested in the thread before (I'm not
familiar with PetscSection yet) for creating staggered grid without dummy
pressure ?
Thanks,
Bishesh
>
> Anton
> >
> >
> >
> >
> >
> >
> >
> >
> >
> > Cheers,
> > Dave
> >
> >
> > On 26 July 2013 17:42, Matthew Knepley <knepley at gmail.com
> > <mailto:knepley at gmail.com>> wrote:
> >
> > On Fri, Jul 26, 2013 at 10:13 AM, Bishesh Khanal
> > <bisheshkh at gmail.com <mailto:bisheshkh at gmail.com>> wrote:
> >
> >
> >
> >
> > On Fri, Jul 26, 2013 at 4:22 PM, Matthew Knepley
> > <knepley at gmail.com <mailto:knepley at gmail.com>> wrote:
> >
> > On Fri, Jul 26, 2013 at 9:11 AM, Bishesh Khanal
> > <bisheshkh at gmail.com <mailto:bisheshkh at gmail.com>>
> wrote:
> >
> >
> >
> >
> > On Fri, Jul 26, 2013 at 2:32 PM, Matthew Knepley
> > <knepley at gmail.com <mailto:knepley at gmail.com>>
> wrote:
> >
> > On Fri, Jul 26, 2013 at 7:28 AM, Bishesh
> > Khanal <bisheshkh at gmail.com
> > <mailto:bisheshkh at gmail.com>> wrote:
> >
> >
> >
> >
> > On Wed, Jul 17, 2013 at 9:48 PM, Jed Brown
> > <jedbrown at mcs.anl.gov
> > <mailto:jedbrown at mcs.anl.gov>> wrote:
> >
> > Bishesh Khanal <bisheshkh at gmail.com
> > <mailto: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/20130729/f0583dd1/attachment.html
> >
>
> ------------------------------
>
> _______________________________________________
> petsc-users mailing list
> petsc-users at mcs.anl.gov
> https://lists.mcs.anl.gov/mailman/listinfo/petsc-users
>
>
> End of petsc-users Digest, Vol 55, Issue 70
> *******************************************
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20130729/0a65e0a9/attachment-0001.html>
More information about the petsc-users
mailing list