[petsc-users] KSPCG solution blow up

[Matthew Knepley]

On Mon, Apr 15, 2013 at 4:08 PM, Hugo Gagnon <
opensource.petsc at user.fastmail.fm> wrote:

> For the problem I'm describing my serial in-house solver does not work
> with ILU(0) but works with ILU(3).  I have no option to run Jacobi.  When I
> apply the same problem to PETSc's PC solver with ILU(3) in serial I
> get KSP_DIVERGED_INDEFINITE_PC on the first iteration (in MPI the solution
> somewhat converges but very slowly).
>

As Mark said, ILU(3) does not preserve either symmetry or definiteness.

   Matt


> call KSPGetPC(Pksp,Ppc,Pierr)
> call PCSetType(Ppc,PCILU,Pierr)
> call PCFactorSetLevels(Ppc,3,Pierr)
>
> This effectively changes the fill level from 0 to 3, right?
>
> --
>   Hugo Gagnon
>
> On 2013-04-15, at 4:35 PM, Mark F. Adams <mark.adams at columbia.edu> wrote:
>
> ILU is not guaranteed to stay positive and CG requires this.  PETSc's ILU
> is ILU(0).  If your in-house PCG solver works with ILU(0) then there is
> probably a bug in your construction of the PETSc matrix.  If your in-house
> code does not use ILU(0) then I would try PCJACOBI and verify that you get
> the same residaul history as your in-house solver (assuming that you can do
> Jacobi).
>
>
> On Apr 15, 2013, at 4:04 PM, Jed Brown <jedbrown at mcs.anl.gov> wrote:
>
> Hugo Gagnon <opensource.petsc at user.fastmail.fm> writes:
>
> Hi,
>
> I hope you don't mind me writing to you off list, I just don't think
> my problem would be helpful to others anyway.
>
>
> I'd rather keep discussions on the list.  We (PETSc developers) cannot
> scale to have private conversations with all users.
>
> I'm trying to converge a problem with KSPCG, which I know for sure
> that both A and b, as in Ax=b, are correct.  I can successfully
> converge this particular problem using our own serial in-house PCG
> solver, whereas with PETSc the solution blows up at the first
> iteration with error code -8.
>
>
> -ksp_converged_reason prints why, or use KSPConvergedReasons to get the
> string.
>
>              KSP_DIVERGED_INDEFINITE_PC       = -8,
>
>
> I've included parts of my code below.  Again, I've triple checked that
> I build A and b correctly.  Being a beginner with PETSc and parallel
> solvers in general, I'd appreciate if you could give me some simple
> pointers as to why my problem is blowing up and how can I improve my
> code.
>
>
>
> http://scicomp.stackexchange.com/questions/513/why-is-my-iterative-linear-solver-not-converging
>
> Note that this code does work for easier problems, although I've
> noticed that in general PETSc takes way more iterations to converge
> than our in-house solver.  Would you have an idea as to why this is
> so?
>
>
> You are either doing a different algorithm or evaluating convergence
> differently.  What preconditioner are you using in serial?  Why aren't
> you using PETSc's multigrid?  Send the diagnostics described in the link
> if you want more help.
>
> I know I'm asking quite a bit but I've been struggling with this
> problem for what I feel far too long!
>
> Thank you for your help,
> --
> Hugo Gagnon
>
>   !-- initialize petsc with the current communicator
>   PETSC_COMM_WORLD = COMM_CURRENT
>   call PetscInitialize(PETSC_NULL_CHARACTER,Pierr)
>
>   !-- instantiate the lhs and solution vectors
>   call VecCreate(PETSC_COMM_WORLD,Pb,Pierr)
>   call VecSetSizes(Pb,PETSC_DECIDE,numVar,Pierr))
>   call VecSetType(Pb,VECSTANDARD,Pierr)
>   call VecDuplicate(Pb,Psol,Pierr)
>
>   !-- instantiate the rhs
>   call MatCreate(PETSC_COMM_WORLD,Paa,Pierr)
>   call MatSetSizes(Paa,PETSC_DECIDE,PETSC_DECIDE,numVar,numVar,Pierr)
>   call MatSetType(Paa,MATAIJ,Pierr) !-- csr format
>
>   !-- preallocate memory
>   mnnz = 81 !-- knob
>   call MatSeqAIJSetPreallocation(Paa,mnnz,PETSC_NULL_INTEGER,Pierr)
>   call
> MatMPIAIJSetPreallocation(Paa,mnnz,PETSC_NULL_INTEGER,mnnz,PETSC_NULL_INTEGER,Pierr)
>
>   !-- instantiate the linear solver context
>   call KSPCreate(PETSC_COMM_WORLD,Pksp,Pierr)
>   call KSPSetType(Pksp,KSPCG,Pierr)
>
>   !-- use a custom monitor function
>   call
> KSPMonitorSet(Pksp,KSPMonitorPETSc,PETSC_NULL_OBJECT,PETSC_NULL_FUNCTION,Pierr)
>
>   !-- instantiate a scatterer so that root has access to full Psol
>   call VecScatterCreateToZero(Psol,Psols,Psol0,Pierr)
>
>   call MatZeroEntries(Paa,Pierr)
>
>   *** build Psol, Pb and Paa ****
>
>   !-- assemble (distribute) the initial solution, the lhs and the rhs
>   call VecAssemblyBegin(Psol,Pierr)
>   call VecAssemblyBegin(Pb,Pierr)
>   call MatAssemblyBegin(Paa,MAT_FINAL_ASSEMBLY,Pierr)
>   call VecAssemblyEnd(Psol,Pierr)
>   call VecAssemblyEnd(Pb,Pierr)
>   call MatAssemblyEnd(Paa,MAT_FINAL_ASSEMBLY,Pierr)
>   call MatSetOption(Paa,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE,Pierr)
>
>   !-- setup pcg, reusing the lhs as preconditioner
>   call KSPSetOperators(Pksp,Paa,Paa,SAME_NONZERO_PATTERN,Pierr)
>   call
> KSPSetTolerances(Pksp,tol,PETSC_DEFAULT_DOUBLE_PRECISION,PETSC_DEFAULT_DOUBLE_PRECISION,grid%pcg_its,Pierr)
>
>   !-- setup the preconditioner with the desired fill level
>   call KSPGetPC(Pksp,Ppc,Pierr)
>   if (nProc > 1) then
>     call PCSetType(Ppc,PCBJACOBI,Pierr)
>     call KSPSetup(Pksp,Pierr)
>     call PCBJacobiGetSubKSP(Ppc,j,j,Psubksp,Pierr)
>     call KSPGetPC(Psubksp(1),Ppc,Pierr)
>   end if
>   call PCSetType(Ppc,PCILU,Pierr)
>   call PCFactorSetLevels(Ppc,grid%pcg_fill,Pierr)
>
>   !-- do not zero out the solution vector
>   call KSPSetInitialGuessNonzero(Pksp,PETSC_TRUE,Pierr)
>   call KSPSetFromOptions(Pksp,Pierr)
>
>   !-- let petsc do its magic
>   call KSPSolve(Pksp,Pb,Psol,Pierr)
>   call KSPGetConvergedReason(Pksp,Pconv,Pierr)
>
>
>
>
>


-- 
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/20130415/df3a1d96/attachment.html>