Lisandro, A followup to our previous discussion. It sounds to me like you are "actually" solving an n+1 unknown nonlinear problem where the "special" unknown is kept secret from SNES and managed somehow by the application code?
You can guess how I feel about this :-). PETSc/SNES is suppose to be good enough to allow you to feed it the ENTIRE nonlinear problem in a way that would allow as efficient solution as if you "handled the special unknown" specially. In particular for this problem the intended solution is to use the PETSc DMComposite object via DMCompositeCreate() and define your vector as two parts; a vector part and a special "array" part. The DMCompositeAddArray() allows you to provide the "special part". Unfortunately I only have this coded when the "vector part" comes from a DA. There is a place holder for a general vector that uses a VecScatter to communicate ghost points but that is sadly not implemented. You may want to look at this construct to see if it is suitable for your friends needs. And what we need to add (note that though DMComposite can be used with DMMG it does not have to be, it can be used directly with a SNES also. Barry On Mon, 29 Oct 2007, Barry Smith wrote: > > So SNES is actually only "seeing" a system of size n? And you are > NOT using Newton's method on that nonlinear system of size n? Instead > you are solving some other linear system (defined by your solution > process) and pretending it is the Newton solution? > > Barry > > > On Mon, 29 Oct 2007, Lisandro Dalcin wrote: > > > On 10/29/07, Barry Smith <bsmith at mcs.anl.gov> wrote: > > > This may be an error in our design, we may need a shell KSP. It would > > > be > > > nice not to need it, but ... > > > You seem to be saying you want something much like a KSP shell? > > > > Indeed. Perhaps we need one in the end, but I would prefer to avoid that. > > > > > What are the operators for this KSP (the A and B?) do you set them in > > > the > > > SNESSetJacobian()? Is b1 the same size (n+1) as b? Is the Eisenstat walker > > > test used on BOTH KSP::solve(b1,x1) and KSP::solve(b2,x2)? > > > > The operators A and B are the ones computed with the function set in > > SNESSetJacobian(). So both linear solves should be done with the > > matrices, same preconditioner, same tolerances (perhaps determined > > with EW method), etc, but different rhs's. From those two solves, I > > next determine the actual solution update for SNES, and an additional > > update for an 'internal' scalar dof managed inside my aplication code. > > > > Imagine a hypothetical formulation of a nonlinear Stokes-like problem, > > were you have exactly one 'pressure' dof in all your mesh, but this > > dof is coupled with all velocity variables (thus generating a dense > > row and column in your Jacobian matrix). Then you can eliminate this > > scalar degree of freedom, and solve your linear problem by using the > > Schur complement. Then you need two linear solves for your 'momentum' > > equations, using the same 'momentum' block of your jacobian matrix. > > > > > > >
