[petsc-users] a strange bug on using DMDA_STENCIL_BOX

[Barry Smith]

On Jul 10, 2014, at 10:48 PM, Sun, Hui <hus003 at ucsd.edu> wrote:

> I see. Should I use the command MatMPIAIJSetPreallocation to preallocate the memory? 
> 
> By the way, it seems that using DMDA_STENTIL_BOX with stencil width 3 will have 4 times as much memory allocation than using DMDA_STENTIL_STAR, is that supposed to be the problem? 

  It is not really a “problem”. Since the box stencil includes the corners it naturally will use many more locations in the vector. You should just use the stencil that matches the stencil you are using so if you are using a stencil that only sticks out directly on each side and doesn’t use the corners you should use the star stencil. It uses much less memory. If your discretization requires using the corners then you need the box stencil.

> 
> However, here is another code, 2D poisson, with 3 degrees of freedom on each node. I tried the same thing, both DMDA_STENTIL_BOX and DMDA_STENTIL_STAR with stencil width 3. But this one doesn't seem to produce the same kind of problem. 
  
  Yes ILU(0) is very problem dependent, it is mostly used for Poisson type problems (though it is much slower for large problems than using good algorithms like multigrid).

  Barry

> 
> Best,
> Hui
> 
> 
> ________________________________________
> From: Barry Smith [bsmith at mcs.anl.gov]
> Sent: Thursday, July 10, 2014 8:34 PM
> To: Sun, Hui
> Cc: petsc-users at mcs.anl.gov
> Subject: Re: [petsc-users] a strange bug on using DMDA_STENCIL_BOX
> 
> On Jul 10, 2014, at 10:07 PM, Sun, Hui <hus003 at ucsd.edu> wrote:
> 
>> I have also noticed those unnecessary zeros, in both cases DMDA_STENTIL_BOX and DMDA_STENTIL_STAR there are many unnecessary zeros. Do you have any idea how to get rid of those entries?
> 
>   These zeros are for all the locations allowed by the stencil (assuming a full stencil for every variable). The only way to get rid of them is to use a custom matrix allocation.
> 
>   The box stencil will have many extra of these zeros since it includes a lot of “coupling” in the matrix that doesn’t exist.
> 
>  Barry
> 
>> 
>> 
>> ________________________________________
>> From: Barry Smith [bsmith at mcs.anl.gov]
>> Sent: Thursday, July 10, 2014 7:46 PM
>> To: Sun, Hui
>> Cc: petsc-users at mcs.anl.gov
>> Subject: Re: [petsc-users] a strange bug on using DMDA_STENCIL_BOX
>> 
>>  Hmm. Because the box stencil puts extra zeros in the original matrix (due to the “extra” stencil locations) ilu(0) will converge differently with the box stencil than the star stencil with the linear solver.  Thus one will actually get a different convergence history. With Jacobi (which doesn’t depend on the non-zero structure of the matrix) the convergence won’t change. I checked with Jacobi and -da_refine 1 and they converge the same. Are you sure there is an issue with box stencil? You wrote originally
>> 
>>>>> However, if I change it to DMDA_STENTIL_BOX, it does not converge, the residue \|Ax-b\|_2 remains large.
>> 
>> can you explain to me how to reproduce this “does not converge, the residue \|Ax-b\|_2 remains large.”  business.
>> 
>>  Barry
>> 
>> On Jul 10, 2014, at 8:46 PM, Sun, Hui <hus003 at ucsd.edu> wrote:
>> 
>>> Thanks for the hint. That is a possible bug. However, I still haven't caught the bug after rethinking about it. I've got three degrees of freedom on each node, MAC grid, and std 5 points laplacian, central differencing for derivatives. Here I've also attached my code.
>>> 
>>> 
>>> ________________________________________
>>> From: Barry Smith [bsmith at mcs.anl.gov]
>>> Sent: Thursday, July 10, 2014 6:14 PM
>>> To: Sun, Hui
>>> Cc: petsc-users at mcs.anl.gov
>>> Subject: Re: [petsc-users] a strange bug on using DMDA_STENCIL_BOX
>>> 
>>> Perhaps the code is not correct for the star stencil? That is perhaps you are actually setting values in the “corner” locations that are not being communicated with stencil but are communicated with star?  If your finite difference stencil involves the “corner” locations then you need to always use the box stencil type.
>>> 
>>> What is exactly your finite difference stencil? How many degrees of freedom per point?
>>> 
>>> 
>>> Barry
>>> 
>>> On Jul 10, 2014, at 8:07 PM, Sun, Hui <hus003 at ucsd.edu> wrote:
>>> 
>>>> If I run:
>>>> ./test_29 -da_grid_x 6 -da_grid_y 6 -ksp_type bcgs -ksp_monitor_short -ksp_converged_reason
>>>> 
>>>> For DMDA_STENCIL_BOX, I get:
>>>> 
>>>> 0 KSP Residual norm 17.5431
>>>> 1 KSP Residual norm 3.23193
>>>> 2 KSP Residual norm 0.399307
>>>> 3 KSP Residual norm 0.0103842
>>>> 4 KSP Residual norm 0.000124391
>>>> Linear solve converged due to CONVERGED_RTOL iterations 4
>>>> 
>>>> For DMDA_STENCIL_STAR, I get:
>>>> 
>>>> 0 KSP Residual norm 60.7996
>>>> 1 KSP Residual norm 13.7279
>>>> 2 KSP Residual norm 4.15583
>>>> 3 KSP Residual norm 0.430568
>>>> 4 KSP Residual norm 0.0999957
>>>> 5 KSP Residual norm 0.00260724
>>>> 6 KSP Residual norm 0.00191392
>>>> 7 KSP Residual norm 5.4681e-05
>>>> Linear solve converged due to CONVERGED_RTOL iterations 7
>>>> 
>>>> The history seems different.
>>>> 
>>>> 
>>>> ________________________________________
>>>> From: Barry Smith [bsmith at mcs.anl.gov]
>>>> Sent: Thursday, July 10, 2014 5:54 PM
>>>> To: Sun, Hui
>>>> Cc: petsc-users at mcs.anl.gov
>>>> Subject: Re: [petsc-users] a strange bug on using DMDA_STENCIL_BOX
>>>> 
>>>> Your interpretation is correct. If DMDA_STENTIL_STAR works then one absolutely expects that DMDA_STENTIL_BOX would produce the same answer.  If you are running with linux could you run with valgrind to see if there is some strange memory corruption or lack of initialization http://www.mcs.anl.gov/petsc/documentation/faq.html#valgrind
>>>> 
>>>> What happens with -ksp_monitor on both cases? Same convergence history or different?
>>>> 
>>>> Barry
>>>> 
>>>> 
>>>> 
>>>> On Jul 10, 2014, at 7:49 PM, Sun, Hui <hus003 at ucsd.edu> wrote:
>>>> 
>>>>> I'm using KSP for 2D Navier Stokes, finite difference. I have DM object created with DMDA_STENTIL_STAR. I'm only using one single process to run the code. It works fine, the linear system converges to the correct solution with correct order, as expected. However, if I change it to DMDA_STENTIL_BOX, it does not converge, the residue \|Ax-b\|_2 remains large.
>>>>> 
>>>>> I have checked that the rhs for both cases are exactly the same, and the matrices also seem to be the same because I have checked the output of A*y, where y is the analytic solution, and both cases give exactly the same answer. If the linear system is exactly the same, then it seems very strange to me that DMDA_STENTIL_STAR works while DMDA_STENTIL_BOX does not work, for I use the same solver, there shouldn't be any difference.
>>>>> 
>>>>> If you happen to encounter this kind of bug before and have some idea what's going on, please let me know. Thank you very much! If you want to look at the code I've written, feel free to ask.
>>>> 
>>> 
>>> <test_29_cp.c>
>> 
> 
> <2D_poisson_3dof.c>