Screenshot instructions:
Windows
Mac
Red Hat Linux
Ubuntu
Click URL instructions:
Rightclick on ad, choose "Copy Link", then paste here →
(This may not be possible with some types of ads)
From: Roy Stogner <roystgnr@ic...>  20090310 16:51:28

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. My thoughts so far: We ought to be able to put System::adjoint_solve() into the library. 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. 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. 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? 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?  Roy 
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 
From: Roy Stogner <roystgnr@ic...>  20090310 17:53:12

On Tue, 10 Mar 2009, Derek Gaston wrote: > On Mar 10, 2009, at 10:51 AM, Roy Stogner wrote: > >> 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. Great. 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 inplace 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 nontransient 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 pseudotransient problems like FINS we'll need to repeat matrix assembly anyways to get rid of the mass matrix terms. >> 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. 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 singlestep methods. For multistep methods things get trickier, but that's at least a good start. > Currently my idea is to build a finitedifference 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 mesoscale solves for material > properties.... at every quadrature point! Cool stuff. > 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 > ask. 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 adjointofthediscretization 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 convectiondiffusion should be sufficient to start.  Roy 
From: Derek Gaston <friedmud@gm...>  20090310 18:18:25

On Mar 10, 2009, at 11:52 AM, Roy Stogner wrote: > 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 inplace 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 nontransient 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 pseudotransient problems like FINS > we'll need to repeat matrix assembly anyways to get rid of the mass > matrix terms. Doing a full transpose is what I would suggest. It will be less error prone and easier to debug. We can figure out later if we can take some shortcuts. In Epetra I believe you use an Epetra_RowMatrixTransposer to create a transpose. See here: http://trilinos.sandia.gov/packages/docs/dev/packages/epetra/doc/html/classEpetra__RowMatrixTransposer.html#z431_0 At any rate... just getting it going with Petsc will suffice for now. As long as we're doing a full transpose I'm confident we can make it happen with Trilinos. > 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. I'm in the same boat myself... having grown a similar but completely incompatible nonlinear solver framework in house.... > Not ridiculous  that's exactly how I'd do it. But I was hoping > there'd be a better way. Not that I know of... and I've spent some time talking to Estep about this. If anyone would know... he would. > "You're never going to be able to" would probably be going too far, > but "never going to want to" is probably right. Hmmm.... I don't know how to take the derivative of a solve with respect to a variable... are you saying you actually think that's possible? Now that I say that... I guess it's similar to getting sensitivities from an adjoint solve! Interesting... > Sounds good. I'll commemorate the event by adding a stub > adjoint_solve() to SVN tonight. ;) Sweet... that way I can claim it's "In Development"! > 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. I also have Michael Pernice ( https://inlportal.inl.gov/portal/server.pt?open=512&objID=372&PageID=3035&cached=true&mode=2 ) ... who has done a lot of this kind of work in the past. He has funding to work with me to get this working in our code. Actually... I believe he's being funded full time... and _my_ funding is to work with _him_! Supposedly he's already down the path of building finite differenced Jacobians in our framework. I need to ping him to see exactly where he is. My first advice to him was to just do everything manually... build an implicit system fill it up with the transpose of a finite differenced derivative.... set the right hand side and solve to see if he can get the analytical solution for a simple problem. Like I say... I don't know where he's at on that. > Figuring out the linear algebra interface is job #1. So far I like you're idea for the steady state adjoint. Let's work on that. > > 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 adjointofthediscretization 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 > convectiondiffusion should be sufficient to start. Yes... some problems where we can actually analytically write out the adjoint would help a lot... that way we have exact solutions we can compare to. I'll get Mike Pernice to generate a few of these. This is sounding pretty promising! I'm already smelling the papers that will be written! Derek 
From: Roy Stogner <roystgnr@ic...>  20090310 18:47:35

On Tue, 10 Mar 2009, Derek Gaston wrote: > Doing a full transpose is what I would suggest. It will be less error prone > and easier to debug. We can figure out later if we can take some shortcuts. > > In Epetra I believe you use an Epetra_RowMatrixTransposer to create a > transpose. See here: > http://trilinos.sandia.gov/packages/docs/dev/packages/epetra/doc/html/classEpetra__RowMatrixTransposer.html#z431_0 > > At any rate... just getting it going with Petsc will suffice for now. As > long as we're doing a full transpose I'm confident we can make it happen with > Trilinos. Sounds good  full transpose it is. > I'm in the same boat myself... having grown a similar but completely > incompatible nonlinear solver framework in house.... Yeah. How'd that happen again? ;) You could have gotten some of your features for free (including those finite differenced Jacobians) just by building on top of FEMSystem, and the new features you've added seem to be things I'm going to want to shoehorn into FEMSystem myself eventually anyways. >> "You're never going to be able to" would probably be going too far, >> but "never going to want to" is probably right. > > Hmmm.... I don't know how to take the derivative of a solve with respect to a > variable... are you saying you actually think that's possible? Now that I > say that... I guess it's similar to getting sensitivities from an adjoint > solve! Interesting... Exactly. In fact, if you're doing the multiscale solves in libMesh then with any luck we might be able to get your derivatives for you "for free". If you're doing molecular dynamics or some such it's probably harder but certainly not impossible; ask Paul. > This is sounding pretty promising! I'm already smelling the papers that will > be written! Absolutely!  Roy 
From: Derek Gaston <friedmud@gm...>  20090310 21:11:23

On Mar 10, 2009, at 12:47 PM, Roy Stogner wrote: >> I'm in the same boat myself... having grown a similar but >> completely incompatible nonlinear solver framework in house.... > > Yeah. How'd that happen again? ;) You could have gotten some of > your features for free (including those finite differenced Jacobians) > just by building on top of FEMSystem, and the new features you've > added seem to be things I'm going to want to shoehorn into FEMSystem > myself eventually anyways. Well... there are a couple of reasons. The main one is that ours is engineered for matrix free solves. We can swap Residual objects in and out just by editing our input file. Also... our BC's are all specified completely separately from internal residuals. This allows us to swap BC's all the time (again from input file). We get a lot of reuse this way... for instance our DerichletBC object is used all the time by multiple codes and our Diffsion residual object is inherited from to form quite a few of our more complex physics residuals. Our Jacobian information is geared towards generating (/approximating) one block of a Jacobian for physics based preconditioning. We almost never have true jacobians. Finally.... we can couple together any arbitrary set of variables to residuals just by editing the input file. This allows us to turn off contributions from coupling or even run reduced physics simulations without having to do any coding. All in all... it just doesn't quite fit with what you have going on. Your system is geared more towards solving a _specific_ problem. Ours is geared towards solving families of multiphysics problems. Nothing wrong with either approach... but I can't see doing what we do with FEMSystem at the current time. Trust me... I wouldn't mind some convergence here... I just haven't yet seen how it would work. > Exactly. In fact, if you're doing the multiscale solves in libMesh > then with any luck we might be able to get your derivatives for you > "for free". If you're doing molecular dynamics or some such it's > probably harder but certainly not impossible; ask Paul. Unfortunately... no... the solves aren't happening in libMesh. We're actually calling out to a Fortran90 code that does phasefield simulations... I don't even want to think about doing an adjoint for that system (or calculating derivatives). Derek 
From: Roy Stogner <roystgnr@ic...>  20090310 21:26:29

On Tue, 10 Mar 2009, Derek Gaston wrote: > On Mar 10, 2009, at 12:47 PM, Roy Stogner wrote: > > Well... there are a couple of reasons. The main one is that ours is > engineered for matrix free solves. Can't be done with FEMSystem without changing the code  but it wouldn't be a huge change. > We can swap Residual objects in and out just by editing our input > file. That could have been done on top of FEMSystem. I did something very similar on top of deal.II a long time ago. > Also... our BC's are all specified completely separately from > internal residuals. That could have been done on top of FEMSystem... but this is actually one of the things on my long term list of "stuff that needs to be added *to* FEMSystem". It just seems to me like you saw my FEMSystem examples and thought "this gets subclassed to solve one specific system at a time" when you should have thought "this can get subclassed to do whatever we want it to, including assembling residual subobjects based on an input config file". You'd still have had to write much of the same code, but you would have gotten time integration for free, you would have gotten nonlinear solvers for free, and you would have gotten finite differenced Jacobians for free. That would have been something, at least.  Roy 
From: Derek Gaston <friedmud@gm...>  20090310 21:22:02

On Mar 10, 2009, at 12:47 PM, Roy Stogner wrote: > Yeah. How'd that happen again? ;) You could have gotten some of > your features for free (including those finite differenced Jacobians) > just by building on top of FEMSystem, and the new features you've > added seem to be things I'm going to want to shoehorn into FEMSystem > myself eventually anyways. Also... forgot to mention brevity. For instance... the entire residual function for the viscous flux term in the energy equation for navier stokes is only 3 lines of code. Most of our objects average around 36 lines of code... This allows nonC++ people (and dare I say it... even ENGINEERS!) to solve fairly complex PDE's.... Oh  also wanted to mention.... about finite differenced jacobians. We're actually planning on adding some structures to libMesh that will allow you to generate finite differenced jacobians. The idea is to use some of the advanced methods from Trilinos / Petsc. In particular there are some nice coloring algorithms in NOX that will give you HUGE speed boosts. We'll have to talk about what these should things look like. For now I'm thinking that they take a matrix and a pointer to a residual function (the same one that can be attached to the NonlinearImplicitSystem) and it'll fill up that matrix on command (passing in a solution vector and whatever else is necessary). Derek 
From: Roy Stogner <roystgnr@ic...>  20090310 21:33:52

On Tue, 10 Mar 2009, Derek Gaston wrote: > In particular there are some nice coloring algorithms in NOX that > will give you HUGE speed boosts. Compared to finite differencing the whole residual in one direction at a time? Sure. It's exactly what I did on top of deal.II. Compared to finite differencing each element residual? No, it's not going to be noticeably faster... in fact until we properly optimize our FE reinits it'll probably be noticeably slower. But don't let me be a wet blanket  I am definitely *not* trying to dissuade you from adding this capability! Being able to turn on finite differenced Jacobians in any NonlinearImplicitSystem would be fantastic! > We'll have to talk about what these should things look like. For > now I'm thinking that they take a matrix and a pointer to a residual > function (the same one that can be attached to the > NonlinearImplicitSystem) and it'll fill up that matrix on command > (passing in a solution vector and whatever else is necessary). That sounds like the way to go.  Roy 
From: Roy Stogner <roystgnr@ic...>  20090312 16:02:07

On Tue, 10 Mar 2009, Derek Gaston wrote: >> Sounds good. I'll commemorate the event by adding a stub >> adjoint_solve() to SVN tonight. ;) You know, for loose definitions of "tonight". > Sweet... that way I can claim it's "In Development"! adjoint_solve() stubs (and QOI function attachment, etc) are in now. I'm thinking the QOI function for nonFEMSystem stuff will just be expected to directly zero out and sum up a new rhs? At the time adjoint_solve gets called, the old rhs ought to be near zero for nonlinear problems and meaningless for linear ones anyway.  Roy 
From: Roy Stogner <roystgnr@ic...>  20090312 22:24:34

On Thu, 12 Mar 2009, Roy Stogner wrote: > adjoint_solve() stubs (and QOI function attachment, etc) are in now. And now adjoint_solve() implementations for nontransient problems are in. I may or may not have time to futz with this again next week, so Derek, if you or someone working with you is willing to be the guinea pig who determines whether or not any of this actually works, now's the time to start. There's definitely some room for debate on the APIs, both older and new, too. The NonlinearSolver, in particular, has jacobian, matvec, etc. as function pointers which probably ought to be in NonlinearImplicitSystem. And I'm not sure whether the DiffSolver ought to have separate solve() and adjoint_solve() functions; it makes some of the code simpler, and could allow for some optimizations later, but feels redundant.  Roy 
From: Kirk, Benjamin (JSCEG311) <benjamin.kirk1@na...>  20090313 12:50:55

> The NonlinearSolver, in particular, has jacobian, > matvec, etc. as function pointers... i used function pointers as the method for interfacing with user code before I understood functors. Is it time to refactor the interfaces to use function objects? I think we could provide a default functor that used a userprovided function pointer to maintain API compatibility, while also providing a base class users can inherit from to present a more C++ish interface? Ben 
From: Roy Stogner <roystgnr@ic...>  20090313 15:16:02

On Fri, 13 Mar 2009, Kirk, Benjamin (JSCEG311) wrote: >> The NonlinearSolver, in particular, has jacobian, >> matvec, etc. as function pointers... > > i used function pointers as the method for interfacing with user code before > I understood functors. Is it time to refactor the interfaces to use > function objects? I think so, but I'm still undecided as to the exact details. Actually, though, in this case it wasn't the function pointers that bothered me as much as their location. It's a little nonintuitive for Solvers to be tied to Systems (either directly as in my DiffSolver stuff or indirectly through residual/jacobian function pointers), when the System could just be an argument to Solver::solve(). On the other hand, we'd generally want to construct one Solver per System anyway, to make it easy for users to set custom persystem settings, so maybe the point is moot. > I think we could provide a default functor that used a userprovided > function pointer to maintain API compatibility, Well, for API compatibility in this case it would be sufficient to test NonlinearSolver::jacobian, etc. for NULL before using the functor; no default needed. > while also providing a base class users can inherit from to present > a more C++ish interface? Something like an "AlgebraicSystem", which passes rhs/matrix or residual/jacobian values to a solver as requested? Maybe.  Roy 
Sign up for the SourceForge newsletter:
No, thanks