[petsc-users] A number of questions about DMDA with SNES and Quasi-Newton methods

[Barry Smith]

> On Oct 3, 2017, at 5:54 AM, zakaryah . <zakaryah at gmail.com> 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 <jed at jedbrown.org> wrote:
> "zakaryah ." <zakaryah at gmail.com> 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 . <zakaryah at gmail.com> 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 <bsmith at mcs.anl.gov> wrote:
> >>
> >>>
> >>> > On Sep 4, 2017, at 4:09 PM, zakaryah . <zakaryah at gmail.com> 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 <bsmith at mcs.anl.gov> wrote:
> >>> >
> >>> > > On Sep 4, 2017, at 3:48 PM, zakaryah . <zakaryah at gmail.com> 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.
> >>> >
> >>> >
> >>>
> >>>
> >>
>