<div dir="ltr"><div class="gmail_default" style="font-size:small">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, just to make sure there are no mistakes there. I tried</div><div class="gmail_default" style="font-size:small"><br></div><div class="gmail_default" style="font-size:small">
<p class="gmail-p1"><span class="gmail-s1"><span class="gmail-Apple-converted-space"> </span>ierr = SNESSetType(mp->PETSc_snes, SNES</span><span class="gmail-s2">TEST</span><span class="gmail-s1">);CHKERRQ(ierr);<span class="gmail-Apple-converted-space"> </span></span></p>
<p class="gmail-p1"><span class="gmail-s1"><span class="gmail-Apple-converted-space"> </span>ierr = PetscOptionsSetValue(</span><span class="gmail-s3">NULL</span><span class="gmail-s1">, </span><span class="gmail-s4">"-snes_test_display"</span><span class="gmail-s1">, </span><span class="gmail-s4">""</span><span class="gmail-s1">); CHKERRQ(ierr);</span></p><p class="gmail-p1">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.</p></div></div><div class="gmail_extra"><br><div class="gmail_quote">On Tue, Sep 5, 2017 at 9:39 AM, Jed Brown <span dir="ltr"><<a href="mailto:jed@jedbrown.org" target="_blank">jed@jedbrown.org</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><span class="">"zakaryah ." <<a href="mailto:zakaryah@gmail.com">zakaryah@gmail.com</a>> writes:<br>
<br>
> OK - I've checked the Jacobian and function very thoroughly and I am<br>
> confident there are no errors.<br>
<br>
</span>Does Newton converge quadratically when you have a good enough initial guess?<br>
<br>
Globalization of large deformation elasticity is a persistent<br>
engineering challenge. The standard approach is to use a continuation,<br>
often in the form of load increments.<br>
<br>
Regarding trust region documentation, the man page says<br>
<br>
The basic algorithm is taken from "The Minpack Project", by More',<br>
Sorensen, Garbow, Hillstrom, pages 88-111 of "Sources and Development<br>
of Mathematical Software", Wayne Cowell, editor.<br>
<br>
You should be able to make sense of it reading from any other source on<br>
trust region methods.<br>
<div class="HOEnZb"><div class="h5"><br>
> I suspect that I am having problems with a bad starting point, and the SNES<br>
> cannot find the global minimum from there. I know that the global minimum<br>
> (with residual zero) exists in all cases but I would like the methods for<br>
> finding it to be more robust to the starting value.<br>
><br>
> The problem comes from the physics of finite deformations of elastic<br>
> materials. In short, I have a functional of smooth maps on a 3D domain to<br>
> itself. The functional contains two terms. The first term represents<br>
> forces which come from external data, and in the Lagrangian this term only<br>
> involves the values of the map at the point in question. The second term<br>
> penalizes fluctuations in the map, and can take various forms. The<br>
> simplest form is just the Dirichlet energy, but I'm also interested in the<br>
> infinitesimal strain energy and the finite strain energy. The first two<br>
> have terms in the Lagrangian which are first order in the second spatial<br>
> derivatives of the map, while the third (finite strain energy) has terms<br>
> which are up to third order in the first and second spatial derivatives of<br>
> the map. It is the finite strain energy term which has been problematic.<br>
><br>
> The Euler-Lagrange equations are discretized on a cubic grid, with equal<br>
> interval spacing in each dimension. The map is the dependent variable,<br>
> i.e. the x in F(x) = 0. I prefer Neumann boundary conditions. Because the<br>
> spatial derivatives of the map are usually small, the Jacobian typically<br>
> has large values in 3x3 blocks along the diagonal (which depend on the map<br>
> and the external data), and up to 54 values which are functions of the<br>
> spatial derivatives of the map and tend to be smaller.<br>
><br>
> Do you have any advice on diagnosing and improving situations in which<br>
> Newton's method finds a stationary point that is not the state with<br>
> globally minimal residual? One hint is that -snes_type newtonls does not<br>
> find as good a solution as -snes_type newtontr but I don't know much about<br>
> these trust region methods, or how to customize and assess them. I'm<br>
> grateful for any advice.<br>
><br>
> On Mon, Sep 4, 2017 at 5:44 PM, zakaryah . <<a href="mailto:zakaryah@gmail.com">zakaryah@gmail.com</a>> wrote:<br>
><br>
>> Yes, it looks like it IS the other way around, and I think the row is<br>
>><br>
>> 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<br>
>> [0,N-1], and r.k is in [0,P-1].<br>
>><br>
>> That matches the boundary conditions in the displayed Jacobian.<br>
>><br>
>> On Mon, Sep 4, 2017 at 5:33 PM, Barry Smith <<a href="mailto:bsmith@mcs.anl.gov">bsmith@mcs.anl.gov</a>> wrote:<br>
>><br>
>>><br>
>>> > On Sep 4, 2017, at 4:09 PM, zakaryah . <<a href="mailto:zakaryah@gmail.com">zakaryah@gmail.com</a>> wrote:<br>
>>> ><br>
>>> > OK that is super helpful. Just to be sure - for MxNxP, the row r in<br>
>>> the Jacobian is at r.i*P*N*3 + r.j*P*3 + r.k*3 + r.c?<br>
>>><br>
>>> It is that way, or the other way around r.k*M*N*3 + r.j*N*3 + r.k*3 +<br>
>>> r.c<br>
>>> ><br>
>>> ><br>
>>> > On Mon, Sep 4, 2017 at 4:58 PM, Barry Smith <<a href="mailto:bsmith@mcs.anl.gov">bsmith@mcs.anl.gov</a>> wrote:<br>
>>> ><br>
>>> > > On Sep 4, 2017, at 3:48 PM, zakaryah . <<a href="mailto:zakaryah@gmail.com">zakaryah@gmail.com</a>> wrote:<br>
>>> > ><br>
>>> > > One piece of information that would be useful is what ordering PETSc<br>
>>> uses for the Jacobian in the snes_test_display. Is it a natural ordering,<br>
>>> or the PETSc ordering? For debugging the Jacobian manually, the natural<br>
>>> ordering is much easier to work with.<br>
>>> ><br>
>>> > What is displayed is always the natural ordering (internally it is<br>
>>> 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<br>
>>> freedom, and the dimensions of my 3D DMDA are MxNxP, which row of the<br>
>>> 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<br>
>>> ordering and all the degrees of freedom for a single point are next to each<br>
>>> other in the vector/matrix.<br>
>>> ><br>
>>> ><br>
>>><br>
>>><br>
>><br>
</div></div></blockquote></div><br></div>