On Tue, 10 Mar 2009, Derek Gaston wrote:
> On Mar 10, 2009, at 10:51 AM, Roy Stogner wrote:
>> For non-transient non-Jacobian-free 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
>> transpose-inconsistent formulation then adjoint_solve() will probably
>> still be the best we can do automatically.
> This sounds perfect as a first step.
The only question left there is: what do we do for the linear algebra
interface? PETSc's "KSPSolveTranspose" seems to be an internal
routine that doesn't work with their real KSPs. Their
MatCreateTranspose just creates a proxy matrix, and if all it
implements is MatMult then it may not work with their preconditioners.
It looks to me like calling MatTranspose in-place is what we want to
do... if only we can support the same behavior with our other solvers.
I don't see anything in Epetra better than SetUseTranspose, which
looks like it just sets up another proxy matrix with the same
limitations. You're the Trilinos guy; any ideas?
We could have a setting in System which tells it to take DenseMatrix
transposes during global matrix assembly, and then we wouldn't ened
*any* linear solver support... but I don't like that idea. Ideally
for non-transient problems we shouldn't need to repeat matrix assembly
at all, just flip around what we calculated for the last Newton step.
Hmmm... on the other hand, for pseudo-transient problems like FINS
we'll need to repeat matrix assembly anyways to get rid of the mass
>> For transient non-FEMSystems, 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.
Unfortunately for me, Ben wrote the core of FINS before my brilliant
FEMSystem ideas revolutionized libMesh (or before I added a big wad of
dubious code, depending on your perspective). We're just interested
in steady solutions right now, but I don't know if we'll be able to
(or if we'll want to) rewrite in FEMSystem before we become interested
in transient results.
> 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.
That should do for single-step methods. For multistep methods things
get trickier, but that's at least a good start.
> Currently my idea is to build a finite-difference jacobian to transpose for
> the adjoint. This might sound ridiculous... but it does get us started.
Not ridiculous - that's exactly how I'd do it. But I was hoping
there'd be a better way.
> All of our codes are going to be Jacobian free... there is no other
> way (we're doing stuff like meso-scale solves for material
> properties.... at every quadrature point!
> You're never going to get an analytical jacobian for that!).
"You're never going to be able to" would probably be going too far,
but "never going to want to" is probably right.
> 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.
Sounds good. I'll commemorate the event by adding a stub
adjoint_solve() to SVN tonight. ;-)
> 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.
The key here is cooperation: You have a supply of time and a demand
for adjoint math expertise - Paul Bauman has a supply of adjoint math
expertise and a demand for an excuse to get away from DPLR this
summer. We can probably work something out.
> If there is anything you need from me in the short term don't hesitate to
Figuring out the linear algebra interface is job #1.
Some help building up test cases is probably the next step. There's
a number of reasons why complex formulations like our apps might not
work with an adjoint-of-the-discretization method, so if it doesn't
work, to rule out "bad software" as one of those reasons we'll want to
set up some basic example codes too. Diffusion and
convection-diffusion should be sufficient to start.