<div dir="ltr"><div class="gmail_default" style="font-size:small">OK - I've checked the Jacobian and function very thoroughly and I am confident there are no errors.</div><div class="gmail_default" style="font-size:small"><br></div><div class="gmail_default" style="font-size:small">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.</div><div class="gmail_default" style="font-size:small"><br></div><div class="gmail_default" style="font-size:small"><div class="gmail_default">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.</div><div class="gmail_default"><br></div><div class="gmail_default">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.</div><div class="gmail_default"><br></div><div class="gmail_default">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.</div></div></div><div class="gmail_extra"><br><div class="gmail_quote">On Mon, Sep 4, 2017 at 5:44 PM, zakaryah . <span dir="ltr"><<a href="mailto:zakaryah@gmail.com" target="_blank">zakaryah@gmail.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><div class="gmail_default" style="font-size:small">Yes, it looks like it IS the other way around, and I think the row is</div><div class="gmail_default" style="font-size:small"><br></div><div class="gmail_default" style="font-size:small">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].</div><div class="gmail_default" style="font-size:small"><br></div><div class="gmail_default" style="font-size:small">That matches the boundary conditions in the displayed Jacobian.</div></div><div class="HOEnZb"><div class="h5"><div class="gmail_extra"><br><div class="gmail_quote">On Mon, Sep 4, 2017 at 5:33 PM, Barry Smith <span dir="ltr"><<a href="mailto:bsmith@mcs.anl.gov" target="_blank">bsmith@mcs.anl.gov</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><span><br>
> On Sep 4, 2017, at 4:09 PM, zakaryah . <<a href="mailto:zakaryah@gmail.com" target="_blank">zakaryah@gmail.com</a>> wrote:<br>
><br>
> 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?<br>
<br>
</span>  It is that way, or the other way around r.k*M*N*3 + r.j*N*3 + r.k*3 + r.c<br>
<div class="m_-787838728357351020HOEnZb"><div class="m_-787838728357351020h5">><br>
><br>
> On Mon, Sep 4, 2017 at 4:58 PM, Barry Smith <<a href="mailto:bsmith@mcs.anl.gov" target="_blank">bsmith@mcs.anl.gov</a>> wrote:<br>
><br>
> > On Sep 4, 2017, at 3:48 PM, zakaryah . <<a href="mailto:zakaryah@gmail.com" target="_blank">zakaryah@gmail.com</a>> wrote:<br>
> ><br>
> > 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.<br>
><br>
>    What is displayed is always the natural ordering (internally it is not the natural ordering).<br>
><br>
> >  For -n 1, are the orderings the same?<br>
><br>
>   yes<br>
><br>
><br>
><br>
> ><br>
> > 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?<br>
><br>
> 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.<br>
><br>
><br>
<br>
</div></div></blockquote></div><br></div>
</div></div></blockquote></div><br></div>