## libmesh-devel

 Re: [Libmesh-devel] libMesh Adjoint Capabilities (fwd) From: Roy Stogner - 2009-05-06 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 finite-differenced 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 Reply-To: Roy Stogner To: Andrea Hawkins 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 element-wise finite-differenced 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 quantity-of-interest 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 discretization-of-the-adjoint instead of adjoint-of-the-discretization 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 goal-oriented 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 finite-differenced residual based on automatically adjusted EquationSystems::parameters values. --- Roy 
 Re: [Libmesh-devel] libMesh Adjoint Capabilities From: Roy Stogner - 2009-05-06 20:11:52 Oh, I missed something in all that: not only do we have yet to automate \partial R / \partial p, but we leave it up to the user to evaluate the quantity of interest derivative \partial Q / \partial u. Trying to automate that is uglier because, unlike R, in general we can't necessarily express Q as a sum of Q_e on each element. That means finite differencing or automatic differencing is probably out of the question, and that means there's probably nothing we can do to make the users' task any simpler. I'd be happy to be contradicted, though. --- Roy On Wed, 6 May 2009, Roy Stogner wrote: > 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 > finite-differenced 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 > Reply-To: Roy Stogner > To: Andrea Hawkins > 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 element-wise finite-differenced > 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 quantity-of-interest > 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 discretization-of-the-adjoint instead of > adjoint-of-the-discretization 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 goal-oriented 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 > finite-differenced residual based on automatically adjusted > EquationSystems::parameters values. > --- > Roy > > ------------------------------------------------------------------------------ > The NEW KODAK i700 Series Scanners deliver under ANY circumstances! Your > production scanning environment may not be a perfect world - but thanks to > Kodak, there's a perfect scanner to get the job done! With the NEW KODAK i700 > Series Scanner you'll get full speed at 300 dpi even with all image > processing features enabled. http://p.sf.net/sfu/kodak-com > _______________________________________________ > Libmesh-devel mailing list > Libmesh-devel@... > https://lists.sourceforge.net/lists/listinfo/libmesh-devel > 
 Re: [Libmesh-devel] libMesh Adjoint Capabilities From: Derek Gaston - 2009-05-11 16:14:43 Attachments: Message as HTML On Wed, May 6, 2009 at 2:11 PM, Roy Stogner wrote: > > Oh, I missed something in all that: not only do we have yet to > automate \partial R / \partial p, but we leave it up to the user to > evaluate the quantity of interest derivative \partial Q / \partial u. > Trying to automate that is uglier because, unlike R, in general we > can't necessarily express Q as a sum of Q_e on each element. That > means finite differencing or automatic differencing is probably out of > the question, and that means there's probably nothing we can do to > make the users' task any simpler. I'd be happy to be contradicted, > though. The one thing you can do is provide some standard QuantitiyOfInterest objects that do things like integrate (or average) a variable over a subdomain or a boundary (which are pretty standard quantities of interest). Also, the "point" QoI (where you are just interested in the solution at one point) can be easily generalized. This is the direction we took in Encore at Sandia. We had a small suite of QoI objects that people could choose from... and if you needed something more specific you wrote your own QoI object.... In fact... we made things even more general than that. We had a general concept of a Postprocessor... any operation that you do to compute a quantity based on a solution. Each postprocessor object could compute it's value... but could also optionally provide quantity of interest calculations for an adjoint solve. All of this might be too specific to go into libMesh itself... just thought I would give some insight into some options. My personal feeling is that I'll engineer something that looks like this system inside our applications... but I'm not sure everyone in the world would want a system like this. Just providing a callback function for filling in the RHS of an adjoint solve is probably enough... Derek