[petsc-users] Debugging KSP output error

[Dave May]

On 20 August 2015 at 11:01, TAY wee-beng <zonexo at gmail.com> wrote:

>
> On 20/8/2015 3:29 PM, Dave May wrote:
>
>
>
> On 20 August 2015 at 05:28, TAY wee-beng <zonexo at gmail.com> wrote:
>
>> Hi,
>>
>> I run my code on 1, 2 and 3 procs. KSP is used to solve the Poisson eqn.
>>
>> Using MatView and VecView, I found that my LHS matrix and RHS vec are the
>> same for 1,2 and 3 procs.
>>
>> However, my pressure (ans) output is the almost the same (due to
>> truncation err) for 1,2 procs.
>>
>> But for 3 procs, the output is the same as for the 1,2 procs for all
>> values except:
>>
>> 1. the last few values for procs 0
>>
>> 2. the first and last few values for procs 1 and 2.
>>
>> Shouldn't the output be the same when the LHS matrix and RHS vec are the
>> same? How can I debug to find the err?
>>
>>
> It's a bit hard to say much without knowing exactly what solver
> configuration you actually ran and without seeing the difference in the
> solution you are referring too.
>
> Some preconditioners have different behaviour in serial and parallel.
> Thus, the convergence of the solver and the residual history (and thus the
> answer) can look slightly different. This difference will become smaller as
> you solve the system more accurately.
> Do you solve the system accurately? e.g. something like -ksp_rtol 1.0e-10
>
> To avoid the problem mentioned above, try using -pc_type jacobi. This PC
> is the same in serial and parallel. Thus, if your A and b are identical  on
> 1,2,3 procs, then the residuals and solution will also be identical on
> 1,2,3 procs (upto machine precision).
>
> Hi Dave,
>
> I tried using jacobi and it's the same result. I found out that the error
> is actually due to mismatched size between DMDACreate3d and
> MatGetOwnershipRange.
>
> Using
>
> *call
> DMDACreate3d(MPI_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,size_x,size_y,&*
>
>
> *size_z,1,PETSC_DECIDE,PETSC_DECIDE,1,stencil_width,lx,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,da_w,ierr)*
>
> *call
> DMDAGetCorners(da_u,start_ijk(1),start_ijk(2),start_ijk(3),width_ijk(1),width_ijk(2),width_ijk(3),ierr)*
>
> and
>
> *call
> MatCreateAIJ(MPI_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,size_x*size_y*size_z,size_x*size_y*size_z,7,PETSC_NULL_INTEGER,7,PETSC_NULL_INTEGER,A_mat,ierr)*
>
> *call MatGetOwnershipRange(A_mat,ijk_sta_p,ijk_end_p,ierr)*
>
> Is this possible? Or is there an error somewhere? It happens when using 3
> procs, instead of 1 or 2.
>
>
Sure it is possible you get a mismatch in the local sizes if you create the
matrix this way as the matrix created knows nothing about the DMDA, and
specifically, it does not know how it has been spatially decomposed.

If you want to ensure consistency between the DMDA and the matrix,
you should always use
  DMCreateMatrix()
to create the matrix.
Any subsequent calls to MatGetOwnershipRange() will then be consistent with
the DMDA parallel layout.



> For my size_x,size_y,size_z = 4,8,10, it was partitioned along z direction
> with 1->4, 5->7, 8->10 using 3 procs with DMDACreate3d which should give
> ownership (with Fortran index + 1)  of:
>
> myid,ijk_sta_p,ijk_end_p           1         129         192
>  myid,ijk_sta_p,ijk_end_p           0           1         128
>  myid,ijk_sta_p,ijk_end_p           2         193         320
>
> But with MatGetOwnershipRange, I got
>
> myid,ijk_sta_p,ijk_end_p           1         108         214
>  myid,ijk_sta_p,ijk_end_p           0           1         107
>  myid,ijk_sta_p,ijk_end_p           2         215         320
>
> Thanks,
>   Dave
>
>
>
>> --
>> Thank you
>>
>> Yours sincerely,
>>
>> TAY wee-beng
>>
>>
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20150820/5c0316a1/attachment.html>