<div dir="ltr"><div dir="ltr">On Tue, May 19, 2020 at 12:09 PM zakaryah . <<a href="mailto:zakaryah@gmail.com">zakaryah@gmail.com</a>> wrote:<br></div><div class="gmail_quote"><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div style="font-size:small">Thanks Matt. I should have said that the boundary is set to 
BOUNDARY_NONE. I am not sure what you mean about the ordering. I 
understand that the actual indexing, internal to PETSc, does not match 
the natural ordering, as explained in the manual. But, doesn't MatView 
always use the natural ordering? Also, my hand-coded Jacobian routine 
uses MatSetValuesStencil, and the corresponding output from 
snes_test_jacobian_display exactly matches what I expect to see - 
correct layout of nonzero terms according to natural ordering. It is 
only the finite difference Jacobian that contains the unexpected, 
off-stencil terms.</div></blockquote><div><br></div><div style="font-size:small">This is weird. I run the Jacobian test manually, by adding a small perturbation to the configuration vector at each index and calculating the function for the perturbed configuration, then the result looks like it should. This makes me think there is not a problem with my routine for calculating the function, but I can't explain why the Jacobian test by finite difference that PETSc calculates would be different, and have non-zero values outside the stencil.</div></div></blockquote><div><br></div><div>So it sounds like your residual function is the one generating the out of stencil changes. The FD just calls your residual with different inputs.</div><div><br></div><div>  Thanks,</div><div><br></div><div>     Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Tue, May 19, 2020 at 10:13 AM zakaryah . <<a href="mailto:zakaryah@gmail.com" target="_blank">zakaryah@gmail.com</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="auto"><div>Thanks Matt. I should have said that the boundary is set to BOUNDARY_NONE. I am not sure what you mean about the ordering. I understand that the actual indexing, internal to PETSc, does not match the natural ordering, as explained in the manual. But, doesn't MatView always use the natural ordering? Also, my hand-coded Jacobian routine uses MatSetValuesStencil, and the corresponding output from snes_test_jacobian_display exactly matches what I expect to see - correct layout of nonzero terms according to natural ordering. It is only the finite difference Jacobian that contains the unexpected, off-stencil terms.<br><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Tue, May 19, 2020, 5:56 AM Matthew Knepley <<a href="mailto:knepley@gmail.com" target="_blank">knepley@gmail.com</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div dir="ltr">On Tue, May 19, 2020 at 1:57 AM zakaryah . <<a href="mailto:zakaryah@gmail.com" rel="noreferrer" target="_blank">zakaryah@gmail.com</a>> wrote:<br></div><div class="gmail_quote"><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div style="font-size:small">Hi all,</div><div style="font-size:small"><br></div><div style="font-size:small">I'm debugging some convergence issues and I came across something strange. I am using a nonlinear solver on a structured grid. The DMDA is 3 dimensional, with 3 dof, and a box stencil of width 1. There are some small errors in the hand-coded Jacobian, which I am trying to sort out, but at least the fill pattern of the matrix is correct.</div><div style="font-size:small"><br></div><div style="font-size:small">However, when I run with -snes_test_jacobian -snes_test_jacobian_display -snes_compare_explicit, I see something very strange. The finite difference Jacobian has large terms outside the stencil. For example, for x,y,z,c = 0,0,0,0 (row 0), the columns 6, 7, 8, 12, 13, and 14 (column 6 => x=2,y=0,z=0,c=0, etc.) have large values, while columns 9 through 20 are calculated but equal to zero. The "correct" values, i.e., in the stencil, are calculated as well, and nearly agree with my hand-coded Jacobian. This issue does NOT occur in serial, but occurs for any number of processors greater than 1.<br></div><div style="font-size:small"><br></div><div style="font-size:small">I have checked the indexing by hand, carefully, and the memory access with valgrind, and the results were clean. Does anyone have an idea why the finite difference calculation of the Jacobian would produce large values outside the stencil? I am using PETSc 3.12.2 and openMPI 3.1.0. Thanks for your help.<br></div></div>
</blockquote></div><div><br></div>Since it only appears in parallel, I am guessing that your calculation of global ordering does not take into account that we locally<div>reorder, rather than using lexicographic ordering, and you might have periodic boundary conditions.</div><div><br></div><div>  Thanks,</div><div><br></div><div>     Matt<br clear="all"><div><br></div>-- <br><div dir="ltr"><div dir="ltr"><div><div dir="ltr"><div><div dir="ltr"><div>What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>-- Norbert Wiener</div><div><br></div><div><a href="http://www.cse.buffalo.edu/~knepley/" rel="noreferrer" target="_blank">https://www.cse.buffalo.edu/~knepley/</a><br></div></div></div></div></div></div></div></div></div>
</blockquote></div></div></div>
</blockquote></div>
</blockquote></div><br clear="all"><div><br></div>-- <br><div dir="ltr" class="gmail_signature"><div dir="ltr"><div><div dir="ltr"><div><div dir="ltr"><div>What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>-- Norbert Wiener</div><div><br></div><div><a href="http://www.cse.buffalo.edu/~knepley/" target="_blank">https://www.cse.buffalo.edu/~knepley/</a><br></div></div></div></div></div></div></div></div>