<div dir="ltr"><div class="gmail_default" style="font-size:small">Thanks Lukasz, this is extremely helpful.</div><div class="gmail_default" style="font-size:small"><br></div><div class="gmail_default" style="font-size:small">By dynamic relaxation, do you mean converting the elliptic PDEs (Euler-Lagrange eqs.) to parabolic PDEs by introducing time, i.e. instead of first variation = 0, solve dx/dt = -nu first variation, where nu is like a viscosity? It seems to me that this should approach an extremum of the Lagrangian. What do you mean by not consistent?</div><div class="gmail_default" style="font-size:small"><br></div><div class="gmail_default" style="font-size:small">I will look into this and the other ideas you mentioned - thanks again!</div><div class="gmail_default" style="font-size:small"><br></div><div class="gmail_default" style="font-size:small">Cheers, Zak</div></div><div class="gmail_extra"><br><div class="gmail_quote">On Tue, Oct 10, 2017 at 1:06 PM, Lukasz Kaczmarczyk <span dir="ltr"><<a href="mailto:Lukasz.Kaczmarczyk@glasgow.ac.uk" target="_blank">Lukasz.Kaczmarczyk@glasgow.ac.uk</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<div style="word-wrap:break-word">
<br>
<div><br>
</div>
<div><br>
<div><span class="">
<blockquote type="cite">
<div>On 10 Oct 2017, at 17:08, zakaryah . <<a href="mailto:zakaryah@gmail.com" target="_blank">zakaryah@gmail.com</a>> wrote:</div>
<br class="m_-4930927281515451023Apple-interchange-newline">
<div>
<div dir="ltr">
<div class="gmail_default" style="font-size:small">Thanks for clearing that up.</div>
<div class="gmail_default" style="font-size:small"><br>
</div>
<div class="gmail_default" style="font-size:small">I'd appreciate any further help. Here's a summary:</div>
<div class="gmail_default" style="font-size:small"><br>
</div>
<div class="gmail_default" style="font-size:small">My ultimate goal is to find a vector field which minimizes an action. The action is a (nonlinear) function of the field and its first spatial derivatives.</div>
<div class="gmail_default" style="font-size:small"><br>
</div>
<div class="gmail_default" style="font-size:small">My current approach is to derive the (continuous) Euler-Lagrange equations, which results in a nonlinear PDE that the minimizing field must satisfy. These Euler-Lagrange equations are then discretized, and
I'm trying to use an SNES to solve them.</div>
<div class="gmail_default" style="font-size:small"><br>
</div>
<div class="gmail_default" style="font-size:small">The problem is that the solver seems to reach a point at which the Jacobian (this corresponds to the second variation of the action, which is like a Hessian of the energy) becomes nearly singular, but where
the residual (RHS of PDE) is not close to zero. The residual does not decrease over additional SNES iterations, and the line search results in tiny step sizes. My interpretation is that this point of stagnation is a critical point.</div>
<div class="gmail_default" style="font-size:small"><br>
</div>
<div class="gmail_default" style="font-size:small">I have checked the hand-coded Jacobian very carefully and I am confident that it is correct.</div>
<div class="gmail_default" style="font-size:small"><br>
</div>
<div class="gmail_default" style="font-size:small">I am guessing that such a situation is well-known in the field, but I don't know the lingo or literature. If anyone has suggestions I'd be thrilled. Are there documentation/methodologies within PETSc for
this type of situation?</div>
<div class="gmail_default" style="font-size:small"><br>
</div>
<div class="gmail_default" style="font-size:small">Is there any advantage to discretizing the action itself and using the optimization routines? With minor modifications I'll have the gradient and Hessian calculations coded. Are the optimization routines
likely to stagnate in the same way as the nonlinear solver, or can they take advantage of the structure of the problem to overcome this?</div>
<div class="gmail_default" style="font-size:small"><br>
</div>
<div class="gmail_default" style="font-size:small">Thanks a lot in advance for any help.</div>
</div>
</div>
</blockquote>
<div><br>
</div>
<div><br>
</div>
</span><div>Hello,
<div><br>
</div>
<div>
<div>The problem similar to yours, i.e. singular tangent matrix is well known in structural and solid mechanics when the structure becomes unstable and buckling. There are good and bad methods to solve this. </div>
<div><br>
</div>
<div>One is dynamic relaxation, which is not consistent. A good method is to use control equation, which controls external forces/boundary conditions. Search for spherical arc-length control and many derivatives of this method.</div>
<div><br>
</div>
<div>Control equation is a scalar equation, and add "moustaches" to the matrix, you can make shell preconditioner which solves such system efficiently with iterative solvers. I using this for some time with petsc with great success.</div>
</div>
<div><br>
</div>
<div>Regards,</div>
<div>Lukasz</div>
</div><div><div class="h5">
<div><br>
</div>
<br>
<blockquote type="cite">
<div>
<div class="gmail_extra"><br>
<div class="gmail_quote">On Sun, Oct 8, 2017 at 5:57 AM, 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">
<br>
There is apparently confusing in understanding the ordering. Is this all on one process that you get funny results? Are you using MatSetValuesStencil() to provide the matrix (it is generally easier than providing it yourself). In parallel MatView() always
maps the rows and columns to the natural ordering before printing, if you use a matrix created from the DMDA. If you create the matrix yourself it has a different MatView in parallel that is in in thePETSc ordering.\<br>
<span class="m_-4930927281515451023HOEnZb"><font color="#888888"><br>
<br>
Barry<br>
</font></span>
<div class="m_-4930927281515451023HOEnZb">
<div class="m_-4930927281515451023h5"><br>
<br>
<br>
> On Oct 8, 2017, at 8:05 AM, zakaryah . <<a href="mailto:zakaryah@gmail.com" target="_blank">zakaryah@gmail.com</a>> wrote:<br>
><br>
> I'm more confused than ever. I don't understand the output of -snes_type test -snes_test_display.<br>
><br>
> For the user-defined state of the vector (where I'd like to test the Jacobian), the finite difference Jacobian at row 0 evaluates as:<br>
><br>
> row 0: (0, 10768.6) (1, 2715.33) (2, -1422.41) (3, -6121.71) (4, 287.797) (5, 744.695) (9, -1454.66) (10, 6.08793) (11, 148.172) (12, 13.1089) (13, -36.5783) (14, -9.99399) (27, -3423.49) (28, -2175.34) (29, 548.662) (30, 145.753) (31, 17.6603)
(32, -15.1079) (36, 76.8575) (37, 16.325) (38, 4.83918)<br>
><br>
> But the hand-coded Jacobian at row 0 evaluates as:<br>
><br>
> row 0: (0, 10768.6) (1, 2715.33) (2, -1422.41) (3, -6121.71) (4, 287.797) (5, 744.695) (33, -1454.66) (34, 6.08792) (35, 148.172) (36, 13.1089) (37, -36.5783) (38, -9.99399) (231, -3423.49) (232, -2175.34) (233, 548.662) (234, 145.753) (235,
17.6603) (236, -15.1079) (264, 76.8575) (265, 16.325) (266, 4.83917) (267, 0.) (268, 0.) (269, 0.)<br>
> and the difference between the Jacobians at row 0 evaluates as:<br>
><br>
> row 0: (0, 0.000189908) (1, 7.17315e-05) (2, 9.31778e-05) (3, 0.000514947) (4, 0.000178659) (5, 0.000178217) (9, -2.25457e-05) (10, -6.34278e-06) (11, -5.93241e-07) (12, 9.48544e-06) (13, 4.79709e-06) (14, 2.40016e-06) (27, -0.000335696) (28,
-0.000106734) (29, -0.000106653) (30, 2.73119e-06) (31, -7.93382e-07) (32, 1.24048e-07) (36, -4.0302e-06) (37, 3.67494e-06) (38, -2.70115e-06) (39, 0.) (40, 0.) (41, 0.)<br>
><br>
> The difference between the column numbering between the finite difference and the hand-coded Jacobians looks like a serious problem to me, but I'm probably missing something.<br>
><br>
> I am trying to use a 3D DMDA with 3 dof, a box stencil of width 1, and for this test problem the grid dimensions are 11x7x6. For a grid point x,y,z, and dof c, is the index calculated as c + 3*x + 3*11*y + 3*11*7*z? If so, then the column numbers of the
hand-coded Jacobian match those of the 27 point stencil I have in mind. However, I am then at a loss to explain the column numbers in the finite difference Jacobian.<br>
><br>
><br>
><br>
><br>
> On Sat, Oct 7, 2017 at 1:49 PM, zakaryah . <<a href="mailto:zakaryah@gmail.com" target="_blank">zakaryah@gmail.com</a>> wrote:<br>
> OK - I ran with -snes_monitor -snes_converged_reason -snes_linesearch_monitor -ksp_monitor -ksp_monitor_true_residual -ksp_converged_reason -pc_type svd -pc_svd_monitor -snes_type newtonls -snes_compare_explicit<br>
><br>
> and here is the full error message, output immediately after<br>
><br>
> Finite difference Jacobian<br>
> Mat Object: 24 MPI processes<br>
> type: mpiaij<br>
><br>
> [0]PETSC ERROR: --------------------- Error Message ------------------------------<wbr>------------------------------<wbr>--<br>
><br>
> [0]PETSC ERROR: Invalid argument<br>
><br>
> [0]PETSC ERROR: Matrix not generated from a DMDA<br>
><br>
> [0]PETSC ERROR: See <a href="http://www.mcs.anl.gov/petsc/documentation/faq.html" rel="noreferrer" target="_blank">
http://www.mcs.anl.gov/petsc/d<wbr>ocumentation/faq.html</a> for trouble shooting.<br>
><br>
> [0]PETSC ERROR: Petsc Release Version 3.7.6, Apr, 24, 2017<br>
><br>
> [0]PETSC ERROR: ./CalculateOpticalFlow on a arch-linux2-c-opt named node046.hpc.rockefeller.intern<wbr>al by zfrentz Sat Oct 7 13:44:44 2017<br>
><br>
> [0]PETSC ERROR: Configure options --prefix=/ru-auth/local/home/z<wbr>frentz/PETSc-3.7.6 --download-fblaslapack -with-debugging=0<br>
><br>
> [0]PETSC ERROR: #1 MatView_MPI_DA() line 551 in /rugpfs/fs0/home/zfrentz/PETSc<wbr>/build/petsc-3.7.6/src/dm/<wbr>impls/da/fdda.c<br>
><br>
> [0]PETSC ERROR: #2 MatView() line 901 in /rugpfs/fs0/home/zfrentz/PETSc<wbr>/build/petsc-3.7.6/src/mat/<wbr>interface/matrix.c<br>
><br>
> [0]PETSC ERROR: #3 SNESComputeJacobian() line 2371 in /rugpfs/fs0/home/zfrentz/PETSc<wbr>/build/petsc-3.7.6/src/snes/<wbr>interface/snes.c<br>
><br>
> [0]PETSC ERROR: #4 SNESSolve_NEWTONLS() line 228 in /rugpfs/fs0/home/zfrentz/PETSc<wbr>/build/petsc-3.7.6/src/snes/<wbr>impls/ls/ls.c<br>
><br>
> [0]PETSC ERROR: #5 SNESSolve() line 4005 in /rugpfs/fs0/home/zfrentz/PETSc<wbr>/build/petsc-3.7.6/src/snes/<wbr>interface/snes.c<br>
><br>
> [0]PETSC ERROR: #6 solveWarp3D() line 659 in /ru-auth/local/home/zfrentz/Co<wbr>de/OpticalFlow/working/October<wbr>6_2017/mshs.c<br>
><br>
><br>
> On Tue, Oct 3, 2017 at 5:37 PM, Jed Brown <<a href="mailto:jed@jedbrown.org" target="_blank">jed@jedbrown.org</a>> wrote:<br>
> Always always always send the whole error message.<br>
><br>
> "zakaryah ." <<a href="mailto:zakaryah@gmail.com" target="_blank">zakaryah@gmail.com</a>> writes:<br>
><br>
> > I tried -snes_compare_explicit, and got the following error:<br>
> ><br>
> > [0]PETSC ERROR: Invalid argument<br>
> ><br>
> > [0]PETSC ERROR: Matrix not generated from a DMDA<br>
> ><br>
> > What am I doing wrong?<br>
> ><br>
> > On Tue, Oct 3, 2017 at 10:08 AM, Jed Brown <<a href="mailto:jed@jedbrown.org" target="_blank">jed@jedbrown.org</a>> wrote:<br>
> ><br>
> >> Barry Smith <<a href="mailto:bsmith@mcs.anl.gov" target="_blank">bsmith@mcs.anl.gov</a>> writes:<br>
> >><br>
> >> >> On Oct 3, 2017, at 5:54 AM, zakaryah . <<a href="mailto:zakaryah@gmail.com" target="_blank">zakaryah@gmail.com</a>> wrote:<br>
> >> >><br>
> >> >> I'm still working on this. I've made some progress, and it looks like<br>
> >> the issue is with the KSP, at least for now. The Jacobian may be<br>
> >> ill-conditioned. Is it possible to use -snes_test_display during an<br>
> >> intermediate step of the analysis? I would like to inspect the Jacobian<br>
> >> after several solves have already completed,<br>
> >> ><br>
> >> > No, our currently code for testing Jacobians is poor quality and<br>
> >> poorly organized. Needs a major refactoring to do things properly. Sorry<br>
> >><br>
> >> You can use -snes_compare_explicit or -snes_compare_coloring to output<br>
> >> differences on each Newton step.<br>
> >><br>
><br>
><br>
<br>
</div>
</div>
</blockquote>
</div>
<br>
</div>
</div>
</blockquote>
</div></div></div>
<br>
</div>
</div>
</blockquote></div><br></div>