[petsc-users] Result error in repeatedly solving linear equations

[Danyang Su]

Hi All,

I have many linear equations with the same matrix structure (same 
non-zero entries) that are derived from a flow problem at different time 
steps. I feel puzzled that the results are a little different when the 
solver run repeatedly and one by one. Say, I have three equations, I can 
get the following results if running three equations together

Equation 1: Iterations     1 norm  0.9457E-02    Result error PETSc vs 
Solver2, max  0.4362E-02 min -0.2277E-04 norm  0.9458E-02
Equation 2: Iterations     2 norm 0.2994E-05    Result error PETSc vs 
Solver2, max  0.1381E-05 min -0.7209E-08 norm  0.2994E-05
Equation 3: Iterations     2 norm  0.3919E-04    Result error PETSc vs 
Solver2, max  0.9435E-07 min -0.1808E-04 norm  0.3919E-04

But if I solve only one equation every time, then restart the program to 
run another one, the results are like this:

Equation 1: Iterations     1 norm  0.9457E-02    Result error PETSc vs 
Solver2, max  0.4362E-02 min -0.2277E-04 norm  0.9458E-02
Equation 2: Iterations     1 norm 0.7949E-05    Result error PETSc vs 
Solver2, max  0.3501E-05 min -0.8377E-06 norm  0.7949E-05
Equation 3: Iterations     1 norm  0.1980E-04    Result error PETSc vs 
Solver2, max  0.4168E-08 min -0.9085E-05 norm  0.1980E-04

Note: Solver2 is the original sequential solver used in this flow model.

Though there are no big difference in the solution for the above 
equations, I want to know why?

For another large linear equations with more than 400,000 unknowns and 
10,000,000 non-zero entries, if the equations are solved repeatedly, 
they need a lot of iterations or fail, but if the equations are solved 
one by one, it only needs 1 to 2 iterations.

How does this difference come from?

The sample codes are attached bellow.

Thanks and regards,

Danyang

!***************************************************************************!
!Create matrix, rhs and solver
call MatCreateAIJ(Petsc_Comm_World, Petsc_Decide, Petsc_Decide, nb, nb, 
nd_nzrow, &
                   Petsc_Null_Integer, nd_nzrow, Petsc_Null_Integer, a, 
ierr)
call MatSetOption(a,Mat_New_Nonzero_Allocation_Err,Petsc_False,ierr)
call VecCreateMPI(Petsc_Comm_World, Petsc_Decide, nb, b, ierr)
call VecDuplicate(b, x, ierr)
call VecDuplicate(x, u, ierr)
call KSPCreate(Petsc_Comm_World,ksp,ierr)
call KSPSetTolerances(ksp,tol,                                 &
                       PETSC_DEFAULT_DOUBLE_PRECISION,          &
                       PETSC_DEFAULT_DOUBLE_PRECISION,          &
                       100,ierr)
call KSPSetFromOptions(ksp,ierr)

!Do time loop
do i = 1, nTimeStep
     call MatGetOwnershipRange(a,istart,iend,ierr)
     do i = istart, iend - 1
        ii = ia_in(i+1)
        jj = ia_in(i+2)
        call MatSetValues(a, ione, i, jj-ii, ja_in(ii:jj-1)-1, 
a_in(ii:jj-1), Insert_Values, ierr)
     end do
     call MatAssemblyBegin(a, Mat_Final_Assembly, ierr)
     call MatAssemblyEnd(a, Mat_Final_Assembly, ierr)

     call VecGetOwnershipRange(b,istart,iend,ierr)
     call VecSetValues(b, iend-istart, ix(istart+1:iend), 
b_in(istart+1:iend), Insert_Values, ierr)
     call VecAssemblyBegin(b,ierr)
     call VecAssemblyEnd(b,ierr)

     if(i == 1) then
         call MatConvert(a,MATSAME,MAT_INITIAL_MATRIX,a2,ierr)
     end if
     !call KSPSetOperators(ksp,a,a2,SAME_PRECONDITIONER,ierr)
     call KSPSetOperators(ksp,a,a2,SAME_NONZERO_PATTERN,ierr)     !These 
three patterns make no difference in current codes
     !call KSPSetOperators(ksp,a,a2,DIFFERENT_NONZERO_PATTERN,ierr)

     call KSPSolve(ksp,b,x,ierr)

     call KSPGetResidualNorm(ksp,norm,ierr)
     call KSPGetIterationNumber(ksp,its,ierr)
end do

!Destroy objects
!...
!***************************************************************************!
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20130814/eb5ee97d/attachment.html>