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

[zakaryah .]

OK - I've checked the Jacobian and function very thoroughly and I am
confident there are no errors.

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.
>> >
>> >
>>
>>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20170904/040444f7/attachment-0001.html>