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

[Patrick Sanan]

On Mon, Dec 12, 2016 at 5:58 AM, 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.
Maybe a good way to proceed is to have the "both at once computation"
explicitly available in PETSc's Newton or other solvers that "know"
that they can gain efficiency by computing both at once. That is, the
user is always required to register a callback to form the residual,
and may optionally register a function to compute (only) the Jacobian
and/or another function to compute the combined Jacobian/Residual.
Newton can look for the combined function at the times when it would
help, and if nothing was registered would fall back to the current
technique of calling separate residual and Jacobian routines. This
would preserve the benefits of the current approach, allow
optimizations of the kind mentioned, and hopefully remain
somewhat-maintainable; in particular, this doesn't require the user to
be promised anything about the order in which Jacobians and residuals
are computed.

> 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