> On Oct 3, 2017, at 5:54 AM, zakaryah . <[email protected]> wrote: > > I'm still working on this. I've made some progress, and it looks like the > issue is with the KSP, at least for now. The Jacobian may be > ill-conditioned. Is it possible to use -snes_test_display during an > intermediate step of the analysis? I would like to inspect the Jacobian > after several solves have already completed,
No, our currently code for testing Jacobians is poor quality and poorly organized. Needs a major refactoring to do things properly. Sorry Barry > just to make sure there are no mistakes there. I tried > > ierr = SNESSetType(mp->PETSc_snes, SNESTEST);CHKERRQ(ierr); > > > ierr = PetscOptionsSetValue(NULL, "-snes_test_display", ""); > CHKERRQ(ierr); > > and the first line works, of course, but the second line doesn't seem to > activate the printing of the Jacobian. I also tried it with "true" in the > last argument and that didn't work either. > > > On Tue, Sep 5, 2017 at 9:39 AM, Jed Brown <[email protected]> wrote: > "zakaryah ." <[email protected]> writes: > > > OK - I've checked the Jacobian and function very thoroughly and I am > > confident there are no errors. > > Does Newton converge quadratically when you have a good enough initial guess? > > Globalization of large deformation elasticity is a persistent > engineering challenge. The standard approach is to use a continuation, > often in the form of load increments. > > Regarding trust region documentation, the man page says > > The basic algorithm is taken from "The Minpack Project", by More', > Sorensen, Garbow, Hillstrom, pages 88-111 of "Sources and Development > of Mathematical Software", Wayne Cowell, editor. > > You should be able to make sense of it reading from any other source on > trust region methods. > > > I suspect that I am having problems with a bad starting point, and the SNES > > cannot find the global minimum from there. I know that the global minimum > > (with residual zero) exists in all cases but I would like the methods for > > finding it to be more robust to the starting value. > > > > The problem comes from the physics of finite deformations of elastic > > materials. In short, I have a functional of smooth maps on a 3D domain to > > itself. The functional contains two terms. The first term represents > > forces which come from external data, and in the Lagrangian this term only > > involves the values of the map at the point in question. The second term > > penalizes fluctuations in the map, and can take various forms. The > > simplest form is just the Dirichlet energy, but I'm also interested in the > > infinitesimal strain energy and the finite strain energy. The first two > > have terms in the Lagrangian which are first order in the second spatial > > derivatives of the map, while the third (finite strain energy) has terms > > which are up to third order in the first and second spatial derivatives of > > the map. It is the finite strain energy term which has been problematic. > > > > The Euler-Lagrange equations are discretized on a cubic grid, with equal > > interval spacing in each dimension. The map is the dependent variable, > > i.e. the x in F(x) = 0. I prefer Neumann boundary conditions. Because the > > spatial derivatives of the map are usually small, the Jacobian typically > > has large values in 3x3 blocks along the diagonal (which depend on the map > > and the external data), and up to 54 values which are functions of the > > spatial derivatives of the map and tend to be smaller. > > > > Do you have any advice on diagnosing and improving situations in which > > Newton's method finds a stationary point that is not the state with > > globally minimal residual? One hint is that -snes_type newtonls does not > > find as good a solution as -snes_type newtontr but I don't know much about > > these trust region methods, or how to customize and assess them. I'm > > grateful for any advice. > > > > On Mon, Sep 4, 2017 at 5:44 PM, zakaryah . <[email protected]> wrote: > > > >> Yes, it looks like it IS the other way around, and I think the row is > >> > >> r.c + r.i*3 + r.j*3*M + r.k*3*M*N, where r.i is in [0,M-1], r.j is in > >> [0,N-1], and r.k is in [0,P-1]. > >> > >> That matches the boundary conditions in the displayed Jacobian. > >> > >> On Mon, Sep 4, 2017 at 5:33 PM, Barry Smith <[email protected]> wrote: > >> > >>> > >>> > On Sep 4, 2017, at 4:09 PM, zakaryah . <[email protected]> wrote: > >>> > > >>> > OK that is super helpful. Just to be sure - for MxNxP, the row r in > >>> the Jacobian is at r.i*P*N*3 + r.j*P*3 + r.k*3 + r.c? > >>> > >>> It is that way, or the other way around r.k*M*N*3 + r.j*N*3 + r.k*3 + > >>> r.c > >>> > > >>> > > >>> > On Mon, Sep 4, 2017 at 4:58 PM, Barry Smith <[email protected]> wrote: > >>> > > >>> > > On Sep 4, 2017, at 3:48 PM, zakaryah . <[email protected]> wrote: > >>> > > > >>> > > One piece of information that would be useful is what ordering PETSc > >>> uses for the Jacobian in the snes_test_display. Is it a natural ordering, > >>> or the PETSc ordering? For debugging the Jacobian manually, the natural > >>> ordering is much easier to work with. > >>> > > >>> > What is displayed is always the natural ordering (internally it is > >>> not the natural ordering). > >>> > > >>> > > For -n 1, are the orderings the same? > >>> > > >>> > yes > >>> > > >>> > > >>> > > >>> > > > >>> > > If I use a MatStencil r to represent a field with 3 degrees of > >>> freedom, and the dimensions of my 3D DMDA are MxNxP, which row of the > >>> Jacobian corresponds to r.i=x, r.j=y, r.k=z, r.c=f? > >>> > > >>> > Internally it is complicated but for any viewing it is just the natural > >>> ordering and all the degrees of freedom for a single point are next to > >>> each > >>> other in the vector/matrix. > >>> > > >>> > > >>> > >>> > >> >
