<div dir="ltr"><div>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...</div><div><br></div>On Sun, Dec 11, 2016 at 4:02 PM Matthew Knepley <<a href="mailto:knepley@gmail.com">knepley@gmail.com</a>> wrote:<br><div class="gmail_quote"><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr" class="gmail_msg"><div class="gmail_extra gmail_msg"><div class="gmail_quote gmail_msg">I consider this bad code management more than an analytical case for the technique, but I can see the point.<br></div></div></div></blockquote><div><br></div><div>Can you expand on that?  Do you believe automatic differentiation in general to be "bad code management"?</div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr" class="gmail_msg"><div class="gmail_extra gmail_msg"><div class="gmail_quote gmail_msg"></div></div></div><div dir="ltr" class="gmail_msg"><div class="gmail_extra gmail_msg"><div class="gmail_quote gmail_msg"><div class="gmail_msg">   - Extremely complex, heavy shape function evaluation (think super high order with first and second derivatives needing to be computed)</div><div class="gmail_msg"><br class="gmail_msg"></div></div></div></div><div dir="ltr" class="gmail_msg"><div class="gmail_extra gmail_msg"><div class="gmail_quote gmail_msg"><div class="gmail_msg">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</div><div class="gmail_msg">the basis functions tabulated. I understand that for high order, you use a tensor product evaluation, but you still tabulate in 1D. What is</div><div class="gmail_msg">being recomputed here?</div></div></div></div></blockquote><div><br></div><div>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.</div><div><br></div><div>Sometimes we even have to do this step twice: once for the un-deformed mesh and once for the deformed mesh... on every element.</div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr" class="gmail_msg"><div class="gmail_extra gmail_msg"><div class="gmail_quote gmail_msg"><div class="gmail_msg"> </div><blockquote class="gmail_quote gmail_msg" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr" class="gmail_msg"><div class="gmail_msg">  - Extremely heavy material property computations</div></div></blockquote><div class="gmail_msg"><br class="gmail_msg"></div></div></div></div><div dir="ltr" class="gmail_msg"><div class="gmail_extra gmail_msg"><div class="gmail_quote gmail_msg"><div class="gmail_msg">Yes. I have to think about this one more.</div></div></div></div><div dir="ltr" class="gmail_msg"><div class="gmail_extra gmail_msg"><div class="gmail_quote gmail_msg"><div class="gmail_msg"> </div><blockquote class="gmail_quote gmail_msg" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr" class="gmail_msg"><div class="gmail_msg">  - MANY coupled variables (we've run thousands).</div></div></blockquote><div class="gmail_msg"><br class="gmail_msg"></div></div></div></div><div dir="ltr" class="gmail_msg"><div class="gmail_extra gmail_msg"><div class="gmail_quote gmail_msg"><div class="gmail_msg">Ah, so you are saying that the process of field evaluation at the quadrature points is expensive because you have so many fields.</div><div class="gmail_msg">It feels very similar to the material case, but I cannot articulate why.</div></div></div></div></blockquote><div><br></div><div>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.</div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr" class="gmail_msg"><div class="gmail_extra gmail_msg"><div class="gmail_quote gmail_msg"><div class="gmail_msg"> I guess my gut says that really expensive material properties,</div><div class="gmail_msg">much more expensive than my top level model, should be modeled by something simpler at that level. Same feeling for using</div><div class="gmail_msg">thousands of fields.</div></div></div></div></blockquote><div><br></div><div>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.</div><div><br></div><div>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.</div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr" class="gmail_msg"><div class="gmail_extra gmail_msg"><div class="gmail_quote gmail_msg"><div class="gmail_msg">  1) I compute a Jacobian with every residual. This sucks because line search and lots of other things use residuals. </div></div></div></div></blockquote><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr" class="gmail_msg"><div class="gmail_extra gmail_msg"><div class="gmail_quote gmail_msg"><div class="gmail_msg"><br class="gmail_msg"></div><div class="gmail_msg">  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</div><div class="gmail_msg">      am reusing the residual I computed to check the convergence criterion.</div><div class="gmail_msg"><br class="gmail_msg"></div><div class="gmail_msg">Can you see a nice way to express Newton for this?</div></div></div></div></blockquote><div><br></div><div>You can see my (admittedly stupidly simple) Newton code that works this way here: <a href="https://github.com/friedmud/MOOSE.jl/blob/master/src/solvers/JuliaDenseNonlinearImplicitSolver.jl#L42">https://github.com/friedmud/MOOSE.jl/blob/master/src/solvers/JuliaDenseNonlinearImplicitSolver.jl#L42</a></div><div><br></div><div>Check the assembly code here to see how both are computed simultaneously: <a href="https://github.com/friedmud/MOOSE.jl/blob/master/src/solvers/Assembly.jl#L59">https://github.com/friedmud/MOOSE.jl/blob/master/src/solvers/Assembly.jl#L59</a></div><div><br></div><div>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).</div><div><br></div><div>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.</div><div><br></div><div>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.</div><div><br></div><div>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 :-)</div><div><br></div><div>Derek</div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr" class="gmail_msg"><div class="gmail_extra gmail_msg">
</div></div></blockquote></div></div>