<div dir="ltr"><br><div class="gmail_extra"><br><br><div class="gmail_quote">On Mon, Jul 29, 2013 at 3:42 PM, Bishesh Khanal <span dir="ltr"><<a href="mailto:bisheshkh@gmail.com" target="_blank">bisheshkh@gmail.com</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><br><div class="gmail_extra"><br><br><div class="gmail_quote"><div class="im">On Fri, Jul 26, 2013 at 6:05 PM, Dave May <span dir="ltr"><<a href="mailto:dave.mayhem23@gmail.com" target="_blank">dave.mayhem23@gmail.com</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div>Yes, the nullspace is important.<br></div><div>There is definitely a pressure null space of constants which needs to be considered.<br>
</div><div>I don't believe ONLY using<br> -fieldsplit_1_ksp_constant_null_space<br>
</div><div>is not sufficient. <br><br>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.<br>
<br></div><div>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.<br>
<br></div></div></blockquote><div><br></div></div><div>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!<br>
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.<br>
</div><div><br>The relevant lines in my solver function are (mKsp and mDa are my relevant KSP and DM objects):<br><br> ierr = DMKSPSetComputeRHS(mDa,computeRHSTaras3D,this);CHKERRQ(ierr); <br> ierr = DMKSPSetComputeOperators(mDa,computeMatrixTaras3D,this);CHKERRQ(ierr);<br>
ierr = KSPSetDM(mKsp,mDa);CHKERRQ(ierr); <br> ierr = KSPSetNullSpace(mKsp,mNullSpace);CHKERRQ(ierr);<br> ierr = KSPSetFromOptions(mKsp);CHKERRQ(ierr);<br>
ierr = KSPSetUp(mKsp);CHKERRQ(ierr);<br> ierr = KSPSolve(mKsp,NULL,NULL);CHKERRQ(ierr);<br><br></div><div>And a corresponding addition in the function computeMatrixTaras3D is something equivalent to:<br> ierr = MatNullSpaceRemove(mNullSpace,b,NULL);CHKERRQ(ierr); //where b is the vector that got updated with the rhs values before this line. <br>
<br></div><div>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.<br>
</div></div></div></div></blockquote><div><br></div><div>Never mind, I found out that I had forgotten to normalize the mNullBasis to form the basis of the null space.<br></div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<div dir="ltr"><div class="gmail_extra"><div class="gmail_quote"><div>
</div><div>Am I doing something wrong or completely missing some point here with nullspace removal idea ? <br></div><div>(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)<br>
</div></div></div></div></blockquote><div><br></div><div>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.<br>
</div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><div class="gmail_extra"><div class="gmail_quote"><div>
<br></div><div><div class="h5"><div><br><br></div><div><br><br></div><div><br></div><div> <br> <br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">
<div dir="ltr"><div>
<br></div><div><br>Cheers,<br></div> Dave<br></div><div><div><div class="gmail_extra"><br><br><div class="gmail_quote">On 26 July 2013 17:42, Matthew Knepley <span dir="ltr"><<a href="mailto:knepley@gmail.com" target="_blank">knepley@gmail.com</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div><div>On Fri, Jul 26, 2013 at 10:13 AM, Bishesh Khanal <span dir="ltr"><<a href="mailto:bisheshkh@gmail.com" target="_blank">bisheshkh@gmail.com</a>></span> wrote:<br>
</div></div><div class="gmail_extra"><div class="gmail_quote"><div><div>
<blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><br><div class="gmail_extra"><br><br><div class="gmail_quote">
On Fri, Jul 26, 2013 at 4:22 PM, Matthew Knepley <span dir="ltr"><<a href="mailto:knepley@gmail.com" target="_blank">knepley@gmail.com</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div><div>On Fri, Jul 26, 2013 at 9:11 AM, Bishesh Khanal <span dir="ltr"><<a href="mailto:bisheshkh@gmail.com" target="_blank">bisheshkh@gmail.com</a>></span> wrote:<br>
</div></div><div class="gmail_extra"><div class="gmail_quote"><div><div>
<blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><br><div class="gmail_extra"><br><br><div class="gmail_quote">
On Fri, Jul 26, 2013 at 2:32 PM, Matthew Knepley <span dir="ltr"><<a href="mailto:knepley@gmail.com" target="_blank">knepley@gmail.com</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div>On Fri, Jul 26, 2013 at 7:28 AM, Bishesh Khanal <span dir="ltr"><<a href="mailto:bisheshkh@gmail.com" target="_blank">bisheshkh@gmail.com</a>></span> wrote:<br>
</div><div class="gmail_extra"><div class="gmail_quote"><div>
<blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><br><div class="gmail_extra"><br><br><div class="gmail_quote">
On Wed, Jul 17, 2013 at 9:48 PM, Jed Brown <span dir="ltr"><<a href="mailto:jedbrown@mcs.anl.gov" target="_blank">jedbrown@mcs.anl.gov</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div>Bishesh Khanal <<a href="mailto:bisheshkh@gmail.com" target="_blank">bisheshkh@gmail.com</a>> writes:<br>
<br>
> Now, I implemented two different approaches, each for both 2D and 3D, in<br>
> MATLAB. It works for the smaller sizes but I have problems solving it for<br>
> the problem size I need (250^3 grid size).<br>
> I use staggered grid with p on cell centers, and components of v on cell<br>
> faces. Similar split up of K to cell center and faces to account for the<br>
> variable viscosity case)<br>
<br>
</div>Okay, you're using a staggered-grid finite difference discretization of<br>
variable-viscosity Stokes. This is a common problem and I recommend<br>
starting with PCFieldSplit with Schur complement reduction (make that<br>
work first, then switch to block preconditioner). </blockquote><div><br></div><div>Ok, I made my 3D problem work with PCFieldSplit with Schur complement reduction using the options:<br>-pc_fieldsplit_type schur -pc_fieldsplit_detect_saddle_point -fieldsplit_1_ksp_constant_null_space<br>
</div><div> <br><br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">You can use PCLSC or<br>
(probably better for you), assemble a preconditioning matrix containing<br>
the inverse viscosity in the pressure-pressure block. This diagonal<br>
matrix is a spectrally equivalent (or nearly so, depending on<br>
discretization) approximation of the Schur complement. The velocity<br>
block can be solved with algebraic multigrid. Read the PCFieldSplit<br>
docs (follow papers as appropriate) and let us know if you get stuck.<br>
</blockquote></div><br></div><div class="gmail_extra">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:<br>
</div></div></blockquote><div> </div></div><div>You want Algebraic Multigrid</div><div><br></div><div>-pc_type fieldsplit -pc_fieldsplit_detect_saddle_point</div><div>-pc_fieldsplit_type schur</div><div> -fieldsplit_0_pc_type gamg</div>
<div>
<div> -fieldsplit_1_ksp_type fgmres</div><div> -fieldsplit_1_ksp_constant_null_space</div><div> -fieldsplit_1_ksp_monitor_short</div><div> -fieldsplit_1_pc_type lsc</div><div>-ksp_converged_reason<br></div><div><br>
</div></div></div></div></div></blockquote><div>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.<br>
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).</div></div></div></div></blockquote>
<div><br></div></div></div><div>These are all iterative solvers. You have to make sure everything converges. </div></div></div></div></blockquote><div><br></div><div>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 ?<br>
</div></div></div></div></blockquote><div><br></div></div></div><div>If your outer iteration converges, and you do not like the solution, there are usually two possibilities:</div><div><br></div><div> 1) Your tolerance is too high, start with it cranked down all the way (1e-10) and slowly relax it</div>
<div><br></div><div> 2) You have a null space that you are not accounting for</div><div><div><br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">
<div dir="ltr"><div class="gmail_extra"><div class="gmail_quote"><div></div><div>I have used the MAC scheme with indexing as shown in: fig 7.5, page 96 of: <a href="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" target="_blank">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</a><br>
<br>Thus I have a DM with 4 dof but there are several "ghost values" set to 0.<br></div><div>Would this cause any problem when using the multigrid ? (This has worked fine when not using the multigrid.)</div></div>
</div></div></blockquote><div><br></div></div><div>I don't know exactly how you have implemented this. These should be rows of the identity.</div><div><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">
<div dir="ltr"><div class="gmail_extra"><div class="gmail_quote"><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">
<div dir="ltr"><div class="gmail_extra"><div class="gmail_quote"><div>I do this problem in the tutorial with</div>
<div>a constant viscosity.</div></div></div></div></blockquote><div><br></div><div>Which tutorial are you referring to ? Could you please provide me the link please ? </div></div></div></div></blockquote><div><br></div></div>
<div>
There are a few on the PETSc Tutorials page, but you can look at this</div><div><br></div><div> <a href="http://www.geodynamics.org/cig/community/workinggroups/short/workshops/cdm2013/presentations/SessionIV_Solvers.pdf" target="_blank">http://www.geodynamics.org/cig/community/workinggroups/short/workshops/cdm2013/presentations/SessionIV_Solvers.pdf</a></div>
<div><br></div><div>for a step-by-step example of a saddle-point problem at the end.</div><div><br></div><div> Matt</div><div><div><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">
<div dir="ltr"><div class="gmail_extra"><div class="gmail_quote"><br><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">
<div dir="ltr"><div class="gmail_extra"><div class="gmail_quote"><div><br></div><div> Matt</div><div><div><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">
<div dir="ltr"><div class="gmail_extra"><div class="gmail_quote"><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">
<div dir="ltr"><div class="gmail_extra"><div class="gmail_quote"><div><div>
</div></div><div> Matt</div><div><div><div><br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">
<div dir="ltr"><div class="gmail_extra">
<br></div><div class="gmail_extra">but I get the following errror:<br><br>[0]PETSC ERROR: --------------------- Error Message ------------------------------------<br>[0]PETSC ERROR: Null argument, when expecting valid pointer!<br>
[0]PETSC ERROR: Null Object: Parameter # 2!<br>[0]PETSC ERROR: ------------------------------------------------------------------------<br>[0]PETSC ERROR: Petsc Release Version 3.4.2, Jul, 02, 2013 <br>[0]PETSC ERROR: See docs/changes/index.html for recent updates.<br>
[0]PETSC ERROR: See docs/faq.html for hints about trouble shooting.<br>[0]PETSC ERROR: See docs/index.html for manual pages.<br>[0]PETSC ERROR: ------------------------------------------------------------------------<br>
[0]PETSC ERROR: src/AdLemMain on a arch-linux2-cxx-debug named edwards by bkhanal Fri Jul 26 14:23:40 2013<br>
[0]PETSC ERROR: Libraries linked from /home/bkhanal/Documents/softwares/petsc-3.4.2/arch-linux2-cxx-debug/lib<br>[0]PETSC ERROR: Configure run at Fri Jul 19 14:25:01 2013<br>[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<br>
[0]PETSC ERROR: ------------------------------------------------------------------------<br>[0]PETSC ERROR: MatPtAP() line 8166 in /home/bkhanal/Documents/softwares/petsc-3.4.2/src/mat/interface/matrix.c<br>[0]PETSC ERROR: PCSetUp_MG() line 628 in /home/bkhanal/Documents/softwares/petsc-3.4.2/src/ksp/pc/impls/mg/mg.c<br>
[0]PETSC ERROR: PCSetUp() line 890 in /home/bkhanal/Documents/softwares/petsc-3.4.2/src/ksp/pc/interface/precon.c<br>[0]PETSC ERROR: KSPSetUp() line 278 in /home/bkhanal/Documents/softwares/petsc-3.4.2/src/ksp/ksp/interface/itfunc.c<br>
[0]PETSC ERROR: KSPSolve() line 399 in /home/bkhanal/Documents/softwares/petsc-3.4.2/src/ksp/ksp/interface/itfunc.c<br>[0]PETSC ERROR: PCApply_FieldSplit_Schur() line 807 in /home/bkhanal/Documents/softwares/petsc-3.4.2/src/ksp/pc/impls/fieldsplit/fieldsplit.c<br>
[0]PETSC ERROR: PCApply() line 442 in /home/bkhanal/Documents/softwares/petsc-3.4.2/src/ksp/pc/interface/precon.c<br>[0]PETSC ERROR: KSP_PCApply() line 227 in /home/bkhanal/Documents/softwares/petsc-3.4.2/include/petsc-private/kspimpl.h<br>
[0]PETSC ERROR: KSPInitialResidual() line 64 in /home/bkhanal/Documents/softwares/petsc-3.4.2/src/ksp/ksp/interface/itres.c<br>[0]PETSC ERROR: KSPSolve_GMRES() line 239 in /home/bkhanal/Documents/softwares/petsc-3.4.2/src/ksp/ksp/impls/gmres/gmres.c<br>
[0]PETSC ERROR: KSPSolve() line 441 in /home/bkhanal/Documents/softwares/petsc-3.4.2/src/ksp/ksp/interface/itfunc.c<br><br></div><div class="gmail_extra"><br></div></div>
</blockquote></div></div></div><br><br clear="all"><span><font color="#888888"><span><font color="#888888"><div><div><br></div>-- <br>What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>
-- Norbert Wiener
</div></font></span></font></span></div></div><span><font color="#888888">
</font></span></blockquote></div><span><font color="#888888"><br></font></span></div></div><span><font color="#888888">
</font></span></blockquote></div></div></div><span><font color="#888888"><div><div><br><br clear="all"><div><br></div>-- <br>What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>
-- Norbert Wiener
</div></div></font></span></div></div>
</blockquote></div><br></div></div>
</blockquote></div></div></div><div><div><br><br clear="all"><div><br></div>-- <br>What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>
-- Norbert Wiener
</div></div></div></div>
</blockquote></div><br></div>
</div></div></blockquote></div></div></div><br></div></div>
</blockquote></div><br></div></div>