Comparing between a PETSc matrix and a standard fortran array in compressed row format

[Barry Smith]

  At the "bad step" are the two matrices and right hand sides the same?
The the resulting solution vectors are different (a lot)? Is there any 
reason to think the matrix is (nearly) singular at that point?

   Barry



On Sun, 5 Aug 2007, Ben Tay wrote:

> Hi,
> 
> Now I'm comparing at every time step. The A matrix and rhs of the NSPCG and
> PETSc matrix are exactly the same, NORM is zero. The solutions given by both
> solvers are very similar for e.g. from 1-315 time step, and varying slowly
> since it is a oscillating body problem.  Then strangely, at time=316, PETSc's
> answer suddenly differs by a great difference. e.g. p=0.3 to 50. I must also
> add that the flowfield at that time step still "looks ok". However, the large
> pressure change of the pressure poisson eqn caused the computation of the
> lift/drag coefficient to be erroneous. Moreover, after that, the pressure
> slowly "returns back" to the same answer such that of the NSPCG solver after a
> few time steps ... then after a while the sudden change happens again. In the
> end, I got a cl/cd vs time graph with lots of spikes.
> 
> I hope that the thing can be solved because the NSPCG solver is much slower
> and it is not as stable.
> 
> Thanks
> 
> Barry Smith wrote:
> >    Ben,
> > 
> >     Are you comparing the two matrices at timestep 316? Also compare the
> > two right hand sides at 316. Are they similar, totally different??
> > 
> >    Barry
> > 
> > 
> > On Sat, 4 Aug 2007, Ben Tay wrote:
> > 
> >   
> > > Hi,
> > > 
> > > After using MatAXPY, I used MatGetRowMaxAbs, VecAb and VecMax to get the
> > > row
> > > with max value, get the absolute of that value and print out the
> > > location/value. I managed to find some minor mistakes and corrected them
> > > so
> > > that the 2 matrix are identical. Unfortunately, the same thing still
> > > happens!
> > > 
> > > I've checked for convergence and the iteration no. is 1 and 8 for direct
> > > solving and HYPRE's AMG respectively.  So there's no convergence problem.
> > > I've
> > > also tried to change to the same pc/ksp as that of NSPCG, which is jacobi
> > > presconditioning + KSPBCGS but the same thing happens.
> > > 
> > > Any idea what's happening?
> > > 
> > > Thanks
> > > 
> > > Barry Smith wrote:
> > >     
> > > >   Unless the matrix values are huge then this is a big difference. All
> > > > you
> > > > can do is print the difference matrix with ASCII MatView and
> > > > look for large entries. This will tell you the locations of the trouble
> > > > entries.
> > > > 
> > > >    Barry
> > > > 
> > > > 
> > > > On Sat, 4 Aug 2007, Ben Tay wrote:
> > > > 
> > > >         
> > > > > Hi,
> > > > > 
> > > > > I realised that the matrix storage format of NSPCG is actually ITPACK
> > > > > storage
> > > > > format. Anyway, I tried to insert the values into a PETSc matrix, use
> > > > > MatAXPY
> > > > > with the scalar a = -1 and then use MatNorm on the output matrix.
> > > > > 
> > > > > Using NORM_1, NORM_FROBENIUS and NORM_INFINITY
> > > > > <../Vec/NORM_FROBENIUS.html#NORM_FROBENIUS>, I got ~543, 3194 and 556
> > > > > respectively. Does this mean that the matrix are different? How "much"
> > > > > different does that means?
> > > > > 
> > > > > So what is the next step? Is there anyway to find out the location of
> > > > > the
> > > > > value which is different?
> > > > > 
> > > > > This is my comparison subroutine:
> > > > > 
> > > > > subroutine matrix_compare
> > > > > 
> > > > >    !compare big_A & PETSc matrix
> > > > > 
> > > > >    integer :: k,kk,II,JJ,ierr
> > > > > 
> > > > >    Vec    xx_test,b_rhs_test
> > > > >      Mat    A_mat_test
> > > > > 
> > > > >    PetscScalar aa
> > > > > 
> > > > >    PetscReal nrm
> > > > > 
> > > > >    aa=-1.
> > > > > 
> > > > >    call
> > > > > MatCreateSeqAIJ(PETSC_COMM_SELF,total_k,total_k,10,PETSC_NULL_INTEGER,A_mat_test,ierr)
> > > > > 
> > > > >    call VecCreateSeq(PETSC_COMM_SELF,total_k,b_rhs_test,ierr)
> > > > > 
> > > > >    call VecDuplicate(b_rhs_test,xx_test,ierr)
> > > > > 
> > > > >    call VecAssemblyBegin(b_rhs_test,ierr)
> > > > > 
> > > > >    call VecAssemblyEnd(b_rhs_test,ierr)
> > > > > 
> > > > >    call VecAssemblyBegin(xx_test,ierr)
> > > > > 
> > > > >    call VecAssemblyEnd(xx_test,ierr)
> > > > > 
> > > > >        do k=1,total_k
> > > > >          do kk=1,10
> > > > > 
> > > > >            II=k-1
> > > > > 
> > > > >            JJ=int_A(k,kk)-1
> > > > >                  call
> > > > > MatSetValues(A_mat_test,1,II,1,JJ,big_A(k,kk),INSERT_VALUES,ierr)
> > > > > call VecSetValue(b_rhs_test,II,q_p(k),INSERT_VALUES,ierr)
> > > > >              end do
> > > > > 
> > > > >    end do
> > > > > 
> > > > > !    call MatCopy(A_mat,A_mat_test,DIFFERENT_NONZERO_PATTERN ,ierr)
> > > > > 
> > > > >    call MatAssemblyBegin(A_mat_test,MAT_FINAL_ASSEMBLY,ierr)
> > > > >        call MatAssemblyEnd(A_mat_test,MAT_FINAL_ASSEMBLY,ierr)
> > > > > 
> > > > >    call MatAXPY(A_mat_test,aa,A_mat,SAME_NONZERO_PATTERN,ierr)
> > > > > 
> > > > >    call MatNorm(A_mat_test,NORM_INFINITY,nrm,ierr)
> > > > > 
> > > > >    print *, nrm
> > > > > 
> > > > >    end subroutine matrix_compare
> > > > > 
> > > > > Thanks!
> > > > > 
> > > > > 
> > > > > 
> > > > > Barry Smith wrote:
> > > > >             
> > > > > >   You can use MatCreateSeqAIJWithArrays() to "convert" the NSPCG
> > > > > > matrix
> > > > > > into
> > > > > > PETSc format and then MatAXPY() to difference them and then
> > > > > > MatNorm() to
> > > > > > see
> > > > > > how large the result is.   Make sure the PETSc or hypre solver is
> > > > > > always
> > > > > > converging. Run with
> > > > > > -ksp_converged_reason and or -ksp_monitor. My guess is that the
> > > > > > matrix
> > > > > > is
> > > > > > becoming very ill-conditioned so the solvers with the default
> > > > > > options
> > > > > > are
> > > > > > failing.
> > > > > > 
> > > > > >     Barry
> > > > > > 
> > > > > > 
> > > > > > On Fri, 3 Aug 2007, Ben Tay wrote:
> > > > > > 
> > > > > >                   
> > > > > > > Hi,
> > > > > > > 
> > > > > > > I used 2 packages to solve my poisson eqn, which is part of my NS
> > > > > > > unsteady
> > > > > > > solver. One is the NSPCG solver package which uses the compressed
> > > > > > > row
> > > > > > > format
> > > > > > > to store the A matrix. The other is PETSc. I found that using both
> > > > > > > solvers
> > > > > > > gave me similar answers for a number of time step. However, after
> > > > > > > that,
> > > > > > > the
> > > > > > > answer will suddenly change drastically for the PETSc case. This
> > > > > > > does
> > > > > > > not
> > > > > > > happen for the NSPCG solver.
> > > > > > > 
> > > > > > > For e.g. time step 1-315, oscillating airfoil case, pressure
> > > > > > > changes
> > > > > > > smoothly,
> > > > > > > similar answers in both cases
> > > > > > > 
> > > > > > > at time=316, pressure changes from -3.22 to -3.2 for NSPCG, but
> > > > > > > pressure
> > > > > > > changes from -3.21 to -60.2 for PETSc
> > > > > > > 
> > > > > > > This happens when I use HYPRE's AMG or PETSc's direct solver LU.
> > > > > > > 
> > > > > > > I have been trying to find out what's the cause and I can't find
> > > > > > > the
> > > > > > > answer in
> > > > > > > debugging. I would like to compare the values of the matrix of the
> > > > > > > 2
> > > > > > > different
> > > > > > > solvers and see if there's any difference. However, NSPCG's matrix
> > > > > > > is
> > > > > > > in
> > > > > > > compressed row format while PETSc's one is just an address and it
> > > > > > > can't be
> > > > > > > viewed easily. Moreover, it's a big matrix so it's not possible to
> > > > > > > check
> > > > > > > by
> > > > > > > inspection. I'm thinking of subtracting one matrix by the other
> > > > > > > and
> > > > > > > find
> > > > > > > if
> > > > > > > it's zero. What's the best way to solve this problem? Btw, I'm
> > > > > > > using
> > > > > > > fortran
> > > > > > > and there's no mpi
> > > > > > > 
> > > > > > > Thank you.
> > > > > > > 
> > > > > > > 
> > > > > > > 
> > > > > > >                         
> > > > > >                   
> > > > >             
> > > >         
> > >     
> > 
> > 
> >   
> 
>