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

Mon Jul 29 11:07:00 CDT 2013

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