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

[Derek Gaston]

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

>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20161212/9f121906/attachment.html>