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

[Barry Smith]

On Aug 14, 2013, at 12:28 PM, Danyang Su <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:

  How are you saving the values and reloading them when you "restart the program"? Are you sure you are saving the exact values of everything in binary? Saving in ASCII and reading it in will introduce slight differences in the numerical values.

   Barry

> 
> 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
> !...
> !***************************************************************************!