[petsc-users] Convergence problem in solving a symmetric positive definite matrix in a CFD code

[Xiao, Jianjun (IKET)]

Hello Jed and Matt,

Thanks for your reply.

I design a small 3D case with 10 cells in each direction. Please find the binary files generated with PETSc-dev in the attachment. The matrix should be symmetric and positive definite.

The calculation did not converge at the first time step with PCR+JACOBI. With PCG+JACOBI, it seems it converges in this case. I have very limited experience in using PETSc. In the past, I felt more flexible to change the combination of linear solvers and preconditioners. I could get divergence anyhow, and the main difference was the speed. But now, it seems only PCG+JACOBI works fine for my case. Many other combinations, such as PCR+JACOBI, PCR+BJACOBI,  PCG+BJACOBI, divergent at the first time step or take much longer time. I think there must be something wrong.

Below is how I set up the matrix and the solver.

       CALL DMSetMatType(da_ksp,MATMPISBAIJ,ierr)
       CALL DMCreateMatrix(da_ksp,gmat,ierr)
       CALL DMSetMatrixPreallocateOnly(da_ksp,PETSC_TRUE,ierr)
       CALL MatMPISBAIJSetPreallocation(gmat,1,4,PETSC_NULL_INTEGER,3,PETSC_NULL_INTEGER,ierr)
       CALL MatZeroEntries(gmat,ierr)

       CALL KSPCreate(PETSC_COMM_WORLD,gsolver,ierr)
!       CALL KSPSetType(gsolver,KSPCR,ierr)
       CALL KSPSetType(gsolver,KSPCG,ierr)

       CALL KSPGetPC(gsolver,gprec,ierr)
       CALL PCSetType(gprec,PCBJACOBI,ierr)
!       CALL PCSetType(gprec,PCBJACOBI,ierr)

********************************** matsetvalues 
       ...
           IF (k < kbar-1 .AND. j < jbar-1 .AND. i < ibar-1) THEN ! d-x-y-z

             valpm(1)     = value_ad_gfpp(m)
             colpm(MatStencil_i,1) = i
             colpm(MatStencil_j,1) = j
             colpm(MatStencil_k,1) = k

             valpm(2)     = value_ax_gfpp(m)
             colpm(MatStencil_i,2) = i+1
             colpm(MatStencil_j,2) = j
             colpm(MatStencil_k,2) = k

             valpm(3)     = value_ay_gfpp(m)
             colpm(MatStencil_i,3) = i
             colpm(MatStencil_j,3) = j+1
             colpm(MatStencil_k,3) = k

             valpm(4)     = value_az_gfpp(m)
             colpm(MatStencil_i,4) = i
             colpm(MatStencil_j,4) = j
             colpm(MatStencil_k,4) = k+1

             CALL MatSetValuesBlockedStencil(gmat,1,rowpm,4,colpm,valpm,INSERT_VALUES,ierr)
       ...

       ...
**********************************
       CALL MatAssemblyBegin(gmat,MAT_FINAL_ASSEMBLY,ierr)
       CALL MatAssemblyEnd(gmat,MAT_FINAL_ASSEMBLY,ierr)

       CALL KSPSetOperators(gsolver,gmat,gmat,SAME_NONZERO_PATTERN,ierr)

       epsi_rtol = 1.0e-50_dble
       epsi_atol = 1.0e-5_dble

       CALL KSPSetTolerances(gsolver,epsi_rtol, epsi_atol,&
           & PETSC_DEFAULT_DOUBLE_PRECISION,itmax,ierr)

       CALL KSPSetNormType(gsolver,KSP_NORM_PRECONDITIONED,ierr)

       CALL KSPSolve(gsolver,bvec0_gfpp,xsol0_gfpp,ierr)

Thank you for your help.

JJ


________________________________________
From: Jed Brown [jedbrown at mcs.anl.gov]
Sent: Friday, December 13, 2013 8:34 AM
To: Xiao, Jianjun (IKET); Matthew Knepley
Cc: petsc-users at mcs.anl.gov
Subject: Re: [petsc-users] Convergence problem in solving a symmetric positive definite matrix in a CFD code

"Xiao, Jianjun (IKET)" <jianjun.xiao at kit.edu> writes:

> Hi Matt,
>
> I am using PETSc-dev because only this version supports
> MatSetValuesBlockedStencil in Fortran. I tried many combinations of
> KSP and PC. None of them worked in my case. The same thing happened to
> PETSc-3.4.
>
> Then I changed back to 3.3-P7 and 3.2-P7 without using
> MatSetValuesBlockedStencil. I tried KSPCR+PCBJACOBI and many others,
> they worked fine.

We're going to need exact options and output for the cases you tried.
You should be able to use the same algorithm with both versions (with
suitable options).  And if you have a modest-sized problem, you can run
with -ksp_view_mat and send binaryoutput.  We'll use
src/ksp/ksp/examples/tutorials/ex10.c to reproduce the convergence you
see.
-------------- next part --------------
A non-text attachment was scrubbed...
Name: binaryoutput_PCR_JACOBI
Type: application/octet-stream
Size: 184320 bytes
Desc: binaryoutput_PCR_JACOBI
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20131213/70522c79/attachment-0002.obj>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: binaryoutput_PCG_JACOBI
Type: application/octet-stream
Size: 184320 bytes
Desc: binaryoutput_PCG_JACOBI
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20131213/70522c79/attachment-0003.obj>