[petsc-users] KSPCG solution blow up

[Jed Brown]

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)