From: Roy Stogner <roystgnr@ic...>  20090506 16:50:25

Copying this to the devel list since it's something a few of us might want to talk about. In particular I'd be happy to hear any better ideas for automatically computing \partial R / \partial p. A finitedifferenced residual would work for FEMSystem and NonlinearImplicitSystem, but LinearImplicitSystem doesn't even let you assemble the residual without getting a big fat matrix along with it.  Roy  Forwarded message  Date: Wed, 6 May 2009 11:42:45 0500 (CDT) From: Roy Stogner <roystgnr@...> ReplyTo: Roy Stogner <roy@...> To: Andrea Hawkins <andjhawkins@...> Subject: Re: libMesh Adjoint Capabilities On Wed, 6 May 2009, Andrea Hawkins wrote: > Dr. Oden told me that in the PECOS meetings earlier this week you > had mentioned something about libMesh having some remarkable adjoint > problem capabilities. Brand new development, but starting to get there. > From what he understood, if you had the forward problem R(u,v,p) = > 0, libMesh was somehow able to formulate \partial R/ \partial u and > \partial R / \partial p to set up the adjoint problem and the > sensitivities. Is this correct? That's the whole goal, but it's only about 25% correct, so far. The good news: In the FEMSystem framework you often don't even to write all of R(u,v,p) yourself  you can just give the weighted residual parts from du/dt=F(u,v) and G(u,v)=0, and if you're using conservation variables and unstabilized Galerkin, the library can do theta method Euler or trapezoidal time integration to give you R. Nonlinear or stabilized mass matrices are a little trickier but still doable. This is what we've started testing; Vikram claims that FEMSystem::adjoint_solve() gave him a verified solution for his quantity of interest. FEMSystem will also automatically do elementwise finitedifferenced Jacobians to give you \partial R / \partial u; that's been pretty well tested for years. For any ImplicitSystem based code, libMesh will do the Jacobian transpose automatically  take a System for which you've already got a good matrix/rhs assembly routine, attach a quantityofinterest assembly routine which replaces the rhs , and call System::adjoint_solve(), and it should give you the right output. The bad news: If you want to solve a discretizationoftheadjoint instead of adjointofthediscretization problem, libMesh can't help; best to try to make those things equivalent in the first place. The adjoint solution is done on your current mesh  if you need a finer adjoint solution for goaloriented error estimation then you still have to do the refinement yourself after your solve() and before your adjoint_solve(). That may not change soon; there's no telling which users will want to do p refinement vs. h refinement, use lower tolerances on a fine mesh, etc. We don't do adjoints for transient or decoupled problems yet. I probably won't add the former this year or the latter ever. Adjoint solve semantics may change soon, and frequently  we'll probably soon pick an "official" additional vector to store adjoint solutions in instead of overwriting System::solution, for an immediate example. This stuff is bleeding edge. For instance there's separate adjoint_solve() implementations for LinearImplicitSystem, NonlinearImplicitSystem, and FEMSystem, but only the latter has been tested. The worst news: Right now users have to assemble \partial R / \partial p by hand; I haven't even yet decided what the best way is to give the library access to p. My best guess at a first cut would be a finitedifferenced residual based on automatically adjusted EquationSystems::parameters values.  Roy 