[petsc-users] discontinuous viscosity stokes equation 3D staggered grid
Anton Popov
popov at uni-mainz.de
Mon Jul 29 11:06:53 CDT 2013
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).
Which problem with variable viscosity are you solving? 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. 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.
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-0001.html>
More information about the petsc-users
mailing list