On Mon, Jul 29, 2013 at 3:42 PM, Bishesh Khanal <[email protected]> wrote:
> > > > On Fri, Jul 26, 2013 at 6:05 PM, Dave May <[email protected]> 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 <[email protected]> wrote: >> >>> On Fri, Jul 26, 2013 at 10:13 AM, Bishesh Khanal <[email protected]>wrote: >>> >>>> >>>> >>>> >>>> On Fri, Jul 26, 2013 at 4:22 PM, Matthew Knepley <[email protected]>wrote: >>>> >>>>> On Fri, Jul 26, 2013 at 9:11 AM, Bishesh Khanal >>>>> <[email protected]>wrote: >>>>> >>>>>> >>>>>> >>>>>> >>>>>> On Fri, Jul 26, 2013 at 2:32 PM, Matthew Knepley >>>>>> <[email protected]>wrote: >>>>>> >>>>>>> On Fri, Jul 26, 2013 at 7:28 AM, Bishesh Khanal <[email protected] >>>>>>> > wrote: >>>>>>> >>>>>>>> >>>>>>>> >>>>>>>> >>>>>>>> On Wed, Jul 17, 2013 at 9:48 PM, Jed Brown <[email protected]>wrote: >>>>>>>> >>>>>>>>> Bishesh Khanal <[email protected]> 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 >>> >> >> >
