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

[Matthew Knepley]

On Wed, 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:
>
> 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/c8d2aa84/attachment-0001.html>