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

[Danyang Su]

Hi Matthew,

Thanks so much. It works out now.

Danyang

On 14/08/2013 11:14 AM, Matthew Knepley wrote:
> On Wed, Aug 14, 2013 at 12:28 PM, Danyang Su <danyang.su at gmail.com 
> <mailto:danyang.su at gmail.com>> wrote:
>
>     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)
>
>
>    Why are you doing this?
>
>         end if
>         !call KSPSetOperators(ksp,a,a2,SAME_PRECONDITIONER,ierr)
>
>
> Just use a, a for the matrices
>
>         call KSPSetOperators(ksp,a,a2,SAME_NONZERO_PATTERN,ierr)    
>         !These three patterns make no difference in current codes
>
>
> This DOES matter here if you are using the default PC which is ILU.
>
>    Matt
>
>         !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
>     !...
>     !***************************************************************************!
>
>
>
>
> -- 
> 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/20130814/17adfe3f/attachment.html>