From: Derek Gaston <friedmud@gm...>  20090310 17:07:55

On Mar 10, 2009, at 10:51 AM, Roy Stogner wrote: > > Derek, Dr. Carey tells me you told him that INL's planning to do some > adjointbased parameter sensitivity work with libMesh up there? We > ought to coordinate; I'm getting ready to try out some hacks to do the > exact same thing with libMesh/FINS. Yes... we are... and yes we should. > For nontransient nonJacobianfree systems, this would just do an > assembly, take the matrix transpose, assemble a new rhs with whatever > quantity of interest functional the user specifies (via attachment or > subclass, same as with assembly), and solve. adjoint_solve() will do > the right thing if you've got an exact Jacobian and a formulation for > which the adjoint of the discretization equals the discretization of > the adjoint; if you've only got an inexact jacobian and/or a > transposeinconsistent formulation then adjoint_solve() will probably > still be the best we can do automatically. This sounds perfect as a first step. > For transient FEMSystems, it would be easy to take the negative of the > time derivative and integrate backwards in time for one timestep. It > would be harder but doable to add automatic everytimestep restarts > and make full time integration of the adjoint problem nearly > automatic. This is fine with me. I'll comment on how I'm going to deal with this further down. > For transient nonFEMSystems, automatic adjoint solves are probably > harder: users' assembly routines have to mix mass matrix and Jacobian > terms freely. Setting a negative EquationSystems::parameters > "dt" might work for some codes? Don't try to do anything for them. Make them make the appropriate switches... just provide the infrastructure. I suppose that if you wanted to default to negating dt it wouldn't hurt. > For Jacobianfree systems, Dr. Ghattas seemed to think you were out of > luck  you can solve the forward problem by approximating a Jacobian > multiply by a finite difference, but there's no such approximation for > a Jacobiantranspose multiply, is there? There isn't... unfortunately. BUT... there is a path forward. Currently my idea is to build a finitedifference jacobian to transpose for the adjoint. This might sound ridiculous... but it does get us started. All of our codes are going to be Jacobian free... there is no other way (we're doing stuff like mesoscale solves for material properties.... at every quadrature point! You're never going to get an analytical jacobian for that!). So for now I'm ok with incurring the penalty of building a finitedifference jacobian to get some adjoint work started. For me (just as in FEMSystem) I can set a negative dt and compute a finite differenced jacobian to transpose. Right now, I'm not so worried about transient. Just being able to do some adjoint based error indication / adaptivity would be a huge gain for us. I have 4 months of funding for working on this that I need to burn before October of next year... so I'm going to be ramping up on this myself. I have quite a bit of learning to do... I understand the basic concepts but not a lot of the math surrounding adjoints... so there's definitely going to be a bit of a ramp up period for me. If there is anything you need from me in the short term don't hesitate to ask. Derek 