[petsc-users] Simultaneously compute Residual+Jacobian in SNES

[Derek Gaston]

I will try it and report back.

Derek

On Mon, Dec 12, 2016 at 1:15 PM Barry Smith <bsmith at mcs.anl.gov> wrote:

>
>    Derek,
>
>     Still it would be interesting to see for your application if always
> computing them together ended up being faster than doing "extra" work for
> the Jacobian. This should be trivial to try.
>
>   Barry
>
> > On Dec 12, 2016, at 12:11 PM, Derek Gaston <friedmud at gmail.com> wrote:
> >
> > I agree Barry - we got a bit sidetracked here.
> >
> > We were debating the relative merits of doing that (computing a Jacobian
> with every residual) and musing a bit about the possibility of whether or
> not SNES could be enhanced to make this more efficient.
> >
> > I think I've been convinced that it's not worth pursuing the SNES
> enhancement... the upside is limited and the complexity is not worth it.
> >
> > Derek
> >
> > On Mon, Dec 12, 2016 at 12:46 PM Barry Smith <bsmith at mcs.anl.gov> wrote:
> >
> > > On Dec 11, 2016, at 10:58 PM, Derek Gaston <friedmud at gmail.com> wrote:
> > >
> > > A quick note: I'm not hugely invested in this idea... I'm just talking
> it out since I started it.  The issues might outweigh potential gains...
> > >
> > > On Sun, Dec 11, 2016 at 4:02 PM Matthew Knepley <knepley at gmail.com>
> wrote:
> > > I consider this bad code management more than an analytical case for
> the technique, but I can see the point.
> > >
> > > Can you expand on that?  Do you believe automatic differentiation in
> general to be "bad code management"?
> > >
> > >    - Extremely complex, heavy shape function evaluation (think super
> high order with first and second derivatives needing to be computed)
> > >
> > > I honestly do not understand this one. Maybe I do not understand high
> order since I never use it. If I want to compute an integral, I have
> > > the basis functions tabulated. I understand that for high order, you
> use a tensor product evaluation, but you still tabulate in 1D. What is
> > > being recomputed here?
> > >
> > > In unstructured mesh you still have to compute the reference->physical
> map for each element and map all of the gradients/second derivatives to
> physical space.  This can be quite expensive if you have a lot of shape
> functions and a lot of quadrature points.
> > >
> > > Sometimes we even have to do this step twice: once for the un-deformed
> mesh and once for the deformed mesh... on every element.
> > >
> > >
> > >   - Extremely heavy material property computations
> > >
> > > Yes. I have to think about this one more.
> > >
> > >   - MANY coupled variables (we've run thousands).
> > >
> > > Ah, so you are saying that the process of field evaluation at the
> quadrature points is expensive because you have so many fields.
> > > It feels very similar to the material case, but I cannot articulate
> why.
> > >
> > > It is similar: it's all about how much information you have to
> recompute at each quadrature point.  I was simply giving different
> scenarios for why you could end up with heavy calculations at each
> quadrature point that feed into both the Residual and Jacobian calculations.
> > >
> > > I guess my gut says that really expensive material properties,
> > > much more expensive than my top level model, should be modeled by
> something simpler at that level. Same feeling for using
> > > thousands of fields.
> > >
> > > Even if you can find something simpler it's good to be able to solve
> the expensive one to verify your simpler model.  Sometimes the
> microstructure behavior is complicated enough that it's quite difficult to
> wrap up in a simpler model or (like you said) it won't be clear if a
> simpler model is possible without doing the more expensive model first.
> > >
> > > We really do have models that require thousands (sometimes tens of
> thousands) of coupled PDEs.  Reusing the field evaluations for both the
> residual and Jacobian could be a large win.
> > >
> > >   1) I compute a Jacobian with every residual. This sucks because line
> search and lots of other things use residuals.
> > >
> > >   2) I compute a residual with every Jacobian. This sound like it
> could work because I compute both for the Newton system, but here I
> > >       am reusing the residual I computed to check the convergence
> criterion.
> > >
> > > Can you see a nice way to express Newton for this?
> > >
> > > You can see my (admittedly stupidly simple) Newton code that works
> this way here:
> https://github.com/friedmud/MOOSE.jl/blob/master/src/solvers/JuliaDenseNonlinearImplicitSolver.jl#L42
> > >
> > > Check the assembly code here to see how both are computed
> simultaneously:
> https://github.com/friedmud/MOOSE.jl/blob/master/src/solvers/Assembly.jl#L59
> > >
> > > Lack of line search makes it pretty simple.  However, even with this
> simple code I end up wasting one extra Jacobian evaluation once the
> convergence criteria has been reached.  Whether or not that is advantageous
> depends on the relative tradeoffs of reusable element computations vs
> Jacobian calculation and how many nonlinear iterations you do (if you're
> only doing one nonlinear iteration every timestep then you're wasting 50%
> of your total Jacobian calculation time).
> > >
> > > For a full featured solver you would definitely also want to have the
> ability to compute a residual, by itself, when you want... for things like
> line search.
> > >
> > > You guys have definitely thought a lot more about this than I have...
> I'm just spitballing here... but it does seem like having an optional
> interface for computing a combined residual/Jacobian could save some codes
> a significant amount of time.
> > >
> > > This isn't a strong feeling of mine though.  I think that for now the
> way I'll do it is simply to "waste" a residual calculation when I need a
> Jacobian :-)
> >
> >    Derek,
> >
> >       I don't understand this long conversation. You can compute the
> Jacobian with the Function evaluation. There is no need for a special
> interface, just compute it! Then provide a dummy FormJacobian that does
> nothing. Please let us know if this does not work. (Yes you will need to
> keep the Jacobian matrix inside the FormFunction context so you have access
> to it but that is no big deal.)
> >
> >    Barry
> >
> > >
> > > Derek
> >
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20161212/026c359d/attachment-0001.html>