[petsc-users] KSPCG solution blow up

[Hugo Gagnon]

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

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)
>> 
> 

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20130415/69e1655c/attachment-0001.html>