<div dir="ltr"><div class="gmail_default" style="font-size:small">I think I figured it out. It seems that only the hand-coded Jacobian is displayed in natural ordering, while the finite difference (and hand-coded minus finite difference) are displayed in global ordering. I guess the finite difference Jacobian doesn't inherit the natural ordering from the DM. Thanks for your help!<br></div></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Tue, May 19, 2020 at 5:06 PM Zhang, Hong <<a href="mailto:hongzhang@anl.gov">hongzhang@anl.gov</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 style="overflow-wrap: break-word;">
<br>
<div><br>
<blockquote type="cite">
<div>On May 19, 2020, at 2:51 PM, zakaryah . <<a href="mailto:zakaryah@gmail.com" target="_blank">zakaryah@gmail.com</a>> wrote:</div>
<br>
<div>
<div dir="ltr">
<div class="gmail_default" style="font-size:small">I'd like to have one more pass at debugging it myself first. I think Matt is probably right about the ordering. I need to make sure I understand the ordering used in VecView. If I use DMCreateGlobalVector to
create the vector for the residual, calculate the residual, then call VecView on the residual vector, how do I ensure that the output is in natural ordering?<br>
</div>
</div>
</div>
</blockquote>
<div><br>
</div>
<div>The output is in natural ordering by default. You can use PETSC_VIEWER_NATIVE to make output printed in PETSc ordering and have a comparison.</div>
<div><br>
</div>
<span>PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_NATIVE);<br>
VecView(your_global_vector,PETSC_VIEWER_STDOUT_WORLD);</span></div>
<div><span>PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD);<br>
</span><span><br>
Hong (Mr.)<br>
<br>
</span>
<blockquote type="cite">
<div><br>
<div class="gmail_quote">
<div dir="ltr" class="gmail_attr">On Tue, May 19, 2020 at 1:17 PM Zhang, Hong <<a href="mailto:hongzhang@anl.gov" target="_blank">hongzhang@anl.gov</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>Can you post your code (at least some essential snippets) for the residual evaluation and the Jacobian evaluation?
<div><br>
</div>
<div>Hong (Mr.) <br>
<div><br>
<blockquote type="cite">
<div>On May 19, 2020, at 11:08 AM, zakaryah . <<a href="mailto:zakaryah@gmail.com" target="_blank">zakaryah@gmail.com</a>> wrote:</div>
<br>
<div>
<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 class="gmail_default" 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" class="gmail_default">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.<br>
</div>
</div>
<br>
<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>
</div>
</blockquote>
</div>
<br>
</div>
</div>
</blockquote>
</div>
</div>
</blockquote>
</div>
<br>
</div>
</blockquote></div>