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?
On Tue, May 19, 2020 at 1:17 PM Zhang, Hong <hongzh...@anl.gov> wrote: > Can you post your code (at least some essential snippets) for the residual > evaluation and the Jacobian evaluation? > > Hong (Mr.) > > On May 19, 2020, at 11:08 AM, zakaryah . <zakar...@gmail.com> wrote: > > 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. >> > > 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. > > On Tue, May 19, 2020 at 10:13 AM zakaryah . <zakar...@gmail.com> wrote: > >> 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. >> >> On Tue, May 19, 2020, 5:56 AM Matthew Knepley <knep...@gmail.com> wrote: >> >>> On Tue, May 19, 2020 at 1:57 AM zakaryah . <zakar...@gmail.com> wrote: >>> >>>> Hi all, >>>> >>>> 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. >>>> >>>> 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. >>>> >>>> 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. >>>> >>> >>> Since it only appears in parallel, I am guessing that your calculation >>> of global ordering does not take into account that we locally >>> reorder, rather than using lexicographic ordering, and you might have >>> periodic boundary conditions. >>> >>> Thanks, >>> >>> Matt >>> >>> -- >>> What most experimenters take for granted before they begin their >>> experiments is infinitely more interesting than any results to which their >>> experiments lead. >>> -- Norbert Wiener >>> >>> https://www.cse.buffalo.edu/~knepley/ >>> <http://www.cse.buffalo.edu/~knepley/> >>> >> >
