[petsc-users] A number of questions about DMDA with SNES and Quasi-Newton methods
zakaryah .
zakaryah at gmail.com
Thu Oct 12 13:02:59 CDT 2017
Thanks for the response, Matt - these are excellent questions.
On theoretical grounds, I am certain that the solution to the continuous
PDE exists. Without any serious treatment, I think this means the
discretized system should have a solution up to discretization error, but
perhaps this is indeed a bad approach.
I am not sure whether the equations are "really hard to solve". At each
point, the equations are third order polynomials of the state variable at
that point and at nearby points (i.e. in the stencil). One possible
complication is that the external forces which are applied to the interior
of the material can be fairly complex - they are smooth, but they can have
many inflection points.
I don't have a great test case for which I know a good solution. To my
thinking, there is no way that time-stepping the parabolic version of the
same PDE can fail to yield a solution at infinite time. So, I'm going to
try starting there. Converting the problem to a minimization is a bit
trickier, because the discretization has to be performed one step earlier
in the calculation, and therefore the gradient and Hessian would need to be
recalculated.
Even if there are some problems with time-stepping (speed of convergence?),
maybe I can use the solutions as better test cases for the elliptic PDE
solved via SNES.
Can you give me any additional lingo or references for the fracture problem?
Thanks, Zak
On Wed, Oct 11, 2017 at 8:53 PM, Matthew Knepley <knepley at gmail.com> wrote:
> On Wed, Oct 11, 2017 at 11:33 AM, zakaryah . <zakaryah at gmail.com> wrote:
>
>> Many thanks for the suggestions, Matt.
>>
>> I tried putting the solvers in a loop, like this:
>>
>> do {
>> NewtonLS
>> check convergence
>> if (converged) break
>> NRichardson or NGMRES
>> } while (!converged)
>>
>> The results were interesting, to me at least. With NRichardson, there
>> was indeed improvement in the residual norm, followed by improvement with
>> NewtonLS, and so on for a few iterations of this loop. In each case, after
>> a few iterations the NewtonLS appeared to be stuck in the same way as after
>> the first iteration. Eventually neither method was able to reduce the
>> residual norm, which was still significant, so this was not a total
>> success. With NGMRES, the initial behavior was similar, but eventually the
>> NGMRES progress became erratic. The minimal residual norm was a bit better
>> using NGMRES than NRichardson, but neither combination of methods fully
>> converged. For both NRichardson and NGMRES, I simply used the defaults, as
>> I have no knowledge of how to tune the options for my problem.
>>
>
> Are you certain that the equations have a solution? I become a little
> concerned when richardson stops converging. Its
> still possible you have really hard to solve equations, it just becomes
> less likely. And even if they truly are hard to solve,
> then there should be physical reasons for this. For example, it could be
> that discretizing the minimizing PDE is just the
> wrong thing to do. I believe this is the case in fracture, where you
> attack the minimization problem directly.
>
> Matt
>
>
>> On Tue, Oct 10, 2017 at 4:08 PM, Matthew Knepley <knepley at gmail.com>
>> wrote:
>>
>>> On Tue, Oct 10, 2017 at 12:08 PM, zakaryah . <zakaryah at gmail.com> wrote:
>>>
>>>> Thanks for clearing that up.
>>>>
>>>> I'd appreciate any further help. Here's a summary:
>>>>
>>>> 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.
>>>>
>>>> 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.
>>>>
>>>> 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.
>>>>
>>>
>>> The normal thing to do here (I think) is to engage solvers which do not
>>> depend on that particular point. So using
>>> NRichardson, or maybe NGMRES, to get past that. I would be interested to
>>> see if this is successful.
>>>
>>> Matt
>>>
>>>
>>>> I have checked the hand-coded Jacobian very carefully and I am
>>>> confident that it is correct.
>>>>
>>>> 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?
>>>>
>>>> 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?
>>>>
>>>> Thanks a lot in advance for any help.
>>>>
>>>> On Sun, Oct 8, 2017 at 5:57 AM, Barry Smith <bsmith at mcs.anl.gov> wrote:
>>>>
>>>>>
>>>>> 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.\
>>>>>
>>>>>
>>>>> Barry
>>>>>
>>>>>
>>>>>
>>>>> > On Oct 8, 2017, at 8:05 AM, zakaryah . <zakaryah at gmail.com> wrote:
>>>>> >
>>>>> > I'm more confused than ever. I don't understand the output of
>>>>> -snes_type test -snes_test_display.
>>>>> >
>>>>> > 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:
>>>>> >
>>>>> > 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)
>>>>> >
>>>>> > But the hand-coded Jacobian at row 0 evaluates as:
>>>>> >
>>>>> > 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.)
>>>>> > and the difference between the Jacobians at row 0 evaluates as:
>>>>> >
>>>>> > 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.)
>>>>> >
>>>>> > 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.
>>>>> >
>>>>> > 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.
>>>>> >
>>>>> >
>>>>> >
>>>>> >
>>>>> > On Sat, Oct 7, 2017 at 1:49 PM, zakaryah . <zakaryah at gmail.com>
>>>>> wrote:
>>>>> > 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
>>>>> >
>>>>> > and here is the full error message, output immediately after
>>>>> >
>>>>> > Finite difference Jacobian
>>>>> > Mat Object: 24 MPI processes
>>>>> > type: mpiaij
>>>>> >
>>>>> > [0]PETSC ERROR: --------------------- Error Message
>>>>> --------------------------------------------------------------
>>>>> >
>>>>> > [0]PETSC ERROR: Invalid argument
>>>>> >
>>>>> > [0]PETSC ERROR: Matrix not generated from a DMDA
>>>>> >
>>>>> > [0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/d
>>>>> ocumentation/faq.html for trouble shooting.
>>>>> >
>>>>> > [0]PETSC ERROR: Petsc Release Version 3.7.6, Apr, 24, 2017
>>>>> >
>>>>> > [0]PETSC ERROR: ./CalculateOpticalFlow on a arch-linux2-c-opt named
>>>>> node046.hpc.rockefeller.internal by zfrentz Sat Oct 7 13:44:44 2017
>>>>> >
>>>>> > [0]PETSC ERROR: Configure options --prefix=/ru-auth/local/home/zfrentz/PETSc-3.7.6
>>>>> --download-fblaslapack -with-debugging=0
>>>>> >
>>>>> > [0]PETSC ERROR: #1 MatView_MPI_DA() line 551 in
>>>>> /rugpfs/fs0/home/zfrentz/PETSc/build/petsc-3.7.6/src/dm/impl
>>>>> s/da/fdda.c
>>>>> >
>>>>> > [0]PETSC ERROR: #2 MatView() line 901 in
>>>>> /rugpfs/fs0/home/zfrentz/PETSc/build/petsc-3.7.6/src/mat/int
>>>>> erface/matrix.c
>>>>> >
>>>>> > [0]PETSC ERROR: #3 SNESComputeJacobian() line 2371 in
>>>>> /rugpfs/fs0/home/zfrentz/PETSc/build/petsc-3.7.6/src/snes/in
>>>>> terface/snes.c
>>>>> >
>>>>> > [0]PETSC ERROR: #4 SNESSolve_NEWTONLS() line 228 in
>>>>> /rugpfs/fs0/home/zfrentz/PETSc/build/petsc-3.7.6/src/snes/im
>>>>> pls/ls/ls.c
>>>>> >
>>>>> > [0]PETSC ERROR: #5 SNESSolve() line 4005 in
>>>>> /rugpfs/fs0/home/zfrentz/PETSc/build/petsc-3.7.6/src/snes/in
>>>>> terface/snes.c
>>>>> >
>>>>> > [0]PETSC ERROR: #6 solveWarp3D() line 659 in
>>>>> /ru-auth/local/home/zfrentz/Code/OpticalFlow/working/October
>>>>> 6_2017/mshs.c
>>>>> >
>>>>> >
>>>>> > On Tue, Oct 3, 2017 at 5:37 PM, Jed Brown <jed at jedbrown.org> wrote:
>>>>> > Always always always send the whole error message.
>>>>> >
>>>>> > "zakaryah ." <zakaryah at gmail.com> writes:
>>>>> >
>>>>> > > I tried -snes_compare_explicit, and got the following error:
>>>>> > >
>>>>> > > [0]PETSC ERROR: Invalid argument
>>>>> > >
>>>>> > > [0]PETSC ERROR: Matrix not generated from a DMDA
>>>>> > >
>>>>> > > What am I doing wrong?
>>>>> > >
>>>>> > > On Tue, Oct 3, 2017 at 10:08 AM, Jed Brown <jed at jedbrown.org>
>>>>> wrote:
>>>>> > >
>>>>> > >> Barry Smith <bsmith at mcs.anl.gov> writes:
>>>>> > >>
>>>>> > >> >> On Oct 3, 2017, at 5:54 AM, zakaryah . <zakaryah at gmail.com>
>>>>> wrote:
>>>>> > >> >>
>>>>> > >> >> 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,
>>>>> > >> >
>>>>> > >> > No, our currently code for testing Jacobians is poor quality
>>>>> and
>>>>> > >> poorly organized. Needs a major refactoring to do things
>>>>> properly. Sorry
>>>>> > >>
>>>>> > >> You can use -snes_compare_explicit or -snes_compare_coloring to
>>>>> output
>>>>> > >> differences on each Newton step.
>>>>> > >>
>>>>> >
>>>>> >
>>>>>
>>>>>
>>>>
>>>
>>>
>>> --
>>> What most experimenters take for granted before they begin their
>>> experiments is infinitely more interesting than any results to which their
>>> experiments lead.
>>> -- Norbert Wiener
>>>
>>> https://www.cse.buffalo.edu/~knepley/ <http://www.caam.rice.edu/~mk51/>
>>>
>>
>>
>
>
> --
> What most experimenters take for granted before they begin their
> experiments is infinitely more interesting than any results to which their
> experiments lead.
> -- Norbert Wiener
>
> https://www.cse.buffalo.edu/~knepley/ <http://www.caam.rice.edu/~mk51/>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20171012/66e92843/attachment-0001.html>
More information about the petsc-users
mailing list