I will try it and report back. Derek
On Mon, Dec 12, 2016 at 1:15 PM Barry Smith <[email protected]> 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 <[email protected]> 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 <[email protected]> wrote: > > > > > On Dec 11, 2016, at 10:58 PM, Derek Gaston <[email protected]> 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 <[email protected]> > 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 > > > >
