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

[Jed Brown]

"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.
>>> >
>>> >
>>>
>>>
>>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 832 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20170905/becf663f/attachment.sig>