Alternatively, see DMComposite and src/snes/examples/tutorials/ex22.c. Lukasz Kaczmarczyk <[email protected]> writes:
> On 22 Oct 2017, at 03:16, zakaryah . > <[email protected]<mailto:[email protected]>> 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 . > <[email protected]<mailto:[email protected]>> 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 > <[email protected]<mailto:[email protected]>> wrote: > On Wed, Oct 11, 2017 at 11:33 AM, zakaryah . > <[email protected]<mailto:[email protected]>> 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 > <[email protected]<mailto:[email protected]>> wrote: > On Tue, Oct 10, 2017 at 12:08 PM, zakaryah . > <[email protected]<mailto:[email protected]>> 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 > <[email protected]<mailto:[email protected]>> 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 . >> <[email protected]<mailto:[email protected]>> 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 . >> <[email protected]<mailto:[email protected]>> 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 >> <[email protected]<mailto:[email protected]>> wrote: >> Always always always send the whole error message. >> >> "zakaryah ." <[email protected]<mailto:[email protected]>> 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 >> > <[email protected]<mailto:[email protected]>> wrote: >> > >> >> Barry Smith <[email protected]<mailto:[email protected]>> writes: >> >> >> >> >> On Oct 3, 2017, at 5:54 AM, zakaryah . >> >> >> <[email protected]<mailto:[email protected]>> 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/>
