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

[Barry Smith]

   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
>