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

Jed Brown jed at jedbrown.org
Mon Oct 23 22:09:17 CDT 2017 

There is DMCoarsenHookAdd() and
DMGetNamedGlobalVector()/DMGetNamedLocalVector() that provides a
composable way to store and transfer auxilliary variables between grids.

"zakaryah ." <zakaryah at gmail.com> writes:

> Thanks Barry.  For my problem, restriction seems more natural.  I still
> don't understand how to actually introduce the field.  As I understand it,
> the multi-grid procedures (coarsening, interpolating, etc.) are performed
> on the state variables, which for my problem can be naturally represented
> by a DMDA.  I imagine for the external fields, I create a second DMDA, but
> I'm not sure if I somehow couple it to the state variable DMDA or just pass
> it through the user defined context or what.
>
> I have another question about the DMComposite.  How do I calculate the
> value of the function for the redundant DM, as it depends on the state
> variables at all grid locations?  The example ex21 is easy, because the
> function for the redundant variable depends only on a Lagrange multiplier
> which is known to belong to the first processor.  I am hoping that I can do
> something like increment to the function value for the redundant DM?  Then
> everything is safe because it happens between VecGetArray and
> VecRestoreArray?  Or, is the thread safety due to the
> DMCompositeScatter/DMCompositeGather calls?  In that case, am I forced to
> use ADD_VALUES?
>
> My last question, for now.  Matt - you said it would be tricky to
> preallocate the composite matrix, and that I should turn off allocation and
> just stick in what I need.  Does this mean I should call
> DMSetMatrixStructureOnly with PETSC_TRUE?  In that case, can I assume that
> the initial assembly of the matrix will be slow due to allocation on the
> fly, but since the matrix always has the same non-zero structure, this is
> not an issue for the many repeated uses of the matrix that will occur?
>
> Thanks so much for all the help.
>
> On Mon, Oct 23, 2017 at 3:37 PM, Barry Smith <bsmith at mcs.anl.gov> wrote:
>
>>
>> > On Oct 23, 2017, at 2:13 PM, zakaryah . <zakaryah at gmail.com> wrote:
>> >
>> > Thanks Matt and Jed, the suggestion of DMComposite is perfect for what I
>> want to do.  Thanks Lukasz, as your implementation is also illuminating.
>> >
>> > I'm working on the code now.  I've realized that this is a good
>> opportunity to set up the code so that it will work properly with
>> multigrid.  My problem has a dependence on some external fields.  In other
>> words, there is constant data at each point in the DMDA.  I don't know how
>> to implement those so that they will be scaled properly as the grid is
>> coarsened/refined.  Any hints?
>>
>>     Depending on the meaning of the fields to coarsen one can use either
>> injection (just grab values on the coarse grid points) or restriction
>> (normally this is the transpose of linear interpolation and sort of an
>> averages values near the coarse grid points. To refine one generally uses
>> linear interpolation. The DMDA can provide these operations (it also
>> provides them to the geometric multigrid solver). DMCreateInjection(),
>> DMCreateInterpolation(), DMCreateRestriction().
>>
>>    Barry
>>
>> >
>> > Thanks again - it's amazing to me how thorough the PETSc methods are,
>> and the ease with which the user can access so many powerful methods, while
>> your support is so knowledgeable and responsive.
>> >
>> > On Sun, Oct 22, 2017 at 11:41 AM, Jed Brown <jed at jedbrown.org> wrote:
>> > Alternatively, see DMComposite and src/snes/examples/tutorials/ex22.c.
>> >
>> > Lukasz Kaczmarczyk <Lukasz.Kaczmarczyk at glasgow.ac.uk> writes:
>> >
>> > > On 22 Oct 2017, at 03:16, zakaryah . <zakaryah at gmail.com<mailto:zak
>> aryah at gmail.com>> wrote:
>> > >
>> > > OK, it turns out Lukasz was exactly correct.  With whatever method I
>> try, the solver or stepper approaches a critical point, which is associated
>> with some kind of snap-through.  I have looked into the control techniques
>> and they are pretty ingenious, and I think they should work for my problem,
>> in that I hope to continue through the critical point.  I have a technical
>> question about the implementation, though.
>> > >
>> > > Following Riks 1979 for example, the control parameter is the
>> approximate arc-length in the phase space of loading intensity and
>> displacements.  It represents one additional variable in the system, and
>> there is one additional equation in the system (in Riks, this is eq. 3.9).
>> > >
>> > > In my implementation, the displacements are implemented as a DMDA with
>> 3 dof, since I'm working in 3D.  I'm not sure about the best way to add the
>> single additional variable and equality.  The way I see it, I either give
>> up on using the DMDA, in which case I'm not sure how to efficiently
>> implement the stencil I need to calculate spatial derivatives of the
>> displacements, or I have to add a rather large number of extra variables.
>> For example, if my DMDA is WxHxD, I would have to make it (W+1)xHxD, and
>> each of the extra HxD variables will have 3 dof.  Then 3xHxD-1 variables
>> are in the nullspace (because they don't represent anything, so I would
>> have to add a bunch of zeros to the function and the Jacobian), while the
>> remaining variable is used as the control parameter.  I'm aware of other
>> methods, e.g. Crisfield 1983, but I'm interested in whether there is a
>> straightforward way to implement Riks' method in PETSc.  I'm sure I'm
>> missing something so hopefully someone can give me some hints.
>> > >
>> > > Thanks for all the help!
>> > >
>> > >
>> > > Zakaryah,
>> > >
>> > > If you like to have a peek how we doing that, you can see
>> > > http://mofem.eng.gla.ac.uk/mofem/html/_arc_length_tools_
>> 8hpp_source.html
>> > > http://mofem.eng.gla.ac.uk/mofem/html/_arc_length_tools_
>> 8cpp_source.html
>> > >
>> > > The implementation is specific features related to MoFEM
>> implementation. However, you can follow the same idea; implement shell
>> matrix, which adds column and row with controlling and controlled equation,
>> respectively, This shell matrix has to have an operator for matrix-vector
>> multiplication. Then you have to add preconditioner, which is based on Riks
>> and others. In fact you can use as well FieldSplit pre-conditioner, Riks
>> method is some variant of Schur complement.
>> > >
>> > > Such implementation allows running multi-grid preconditioner and other
>> preconditions with control equation.
>> > >
>> > > Hope that this will be helpful.
>> > >
>> > > Regards,
>> > > Lukasz
>> > >
>> > >
>> > >
>> > > On Thu, Oct 12, 2017 at 2:02 PM, zakaryah . <zakaryah at gmail.com
>> <mailto:zakaryah at gmail.com>> wrote:
>> > > 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
>> <mailto:knepley at gmail.com>> wrote:
>> > > On Wed, Oct 11, 2017 at 11:33 AM, zakaryah . <zakaryah at gmail.com
>> <mailto: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
>> <mailto:knepley at gmail.com>> wrote:
>> > > On Tue, Oct 10, 2017 at 12:08 PM, zakaryah . <zakaryah at gmail.com
>> <mailto: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
>> <mailto: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<mailto:zak
>> aryah 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
>> <mailto: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/
>> documentation/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/impls/da/fdda.c
>> > >>
>> > >> [0]PETSC ERROR: #2 MatView() line 901 in /rugpfs/fs0/home/zfrentz/
>> PETSc/build/petsc-3.7.6/src/mat/interface/matrix.c
>> > >>
>> > >> [0]PETSC ERROR: #3 SNESComputeJacobian() line 2371 in
>> /rugpfs/fs0/home/zfrentz/PETSc/build/petsc-3.7.6/src/snes/interface/snes.c
>> > >>
>> > >> [0]PETSC ERROR: #4 SNESSolve_NEWTONLS() line 228 in
>> /rugpfs/fs0/home/zfrentz/PETSc/build/petsc-3.7.6/src/snes/impls/ls/ls.c
>> > >>
>> > >> [0]PETSC ERROR: #5 SNESSolve() line 4005 in /rugpfs/fs0/home/zfrentz/
>> PETSc/build/petsc-3.7.6/src/snes/interface/snes.c
>> > >>
>> > >> [0]PETSC ERROR: #6 solveWarp3D() line 659 in
>> /ru-auth/local/home/zfrentz/Code/OpticalFlow/working/October6_2017/mshs.c
>> > >>
>> > >>
>> > >> On Tue, Oct 3, 2017 at 5:37 PM, Jed Brown <jed at jedbrown.org<mailto:
>> jed at jedbrown.org>> wrote:
>> > >> Always always always send the whole error message.
>> > >>
>> > >> "zakaryah ." <zakaryah at gmail.com<mailto: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
>> <mailto:jed at jedbrown.org>> wrote:
>> > >> >
>> > >> >> Barry Smith <bsmith at mcs.anl.gov<mailto:bsmith at mcs.anl.gov>>
>> writes:
>> > >> >>
>> > >> >> >> On Oct 3, 2017, at 5:54 AM, zakaryah . <zakaryah at gmail.com
>> <mailto: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/>
>>
>>

More information about the petsc-users mailing list