libMesh: A C++ Finite Element Library

On Fri, 30 Apr 2010, Derek Gaston wrote:

> 1.  I change the signature of the residual and jacobian callbacks so
> that they pass a pointer to the System currently being solved for.

We're currently passing back a pointer to the solver, ctx, and the
solver contains a pointer to the System, _system.
---
Roy

On Apr 30, 2010, at 12:09 PM, Roy Stogner wrote:

> We're currently passing back a pointer to the solver, ctx, and the
> solver contains a pointer to the System, _system.

Yes, this is what we get back from the callback from Petsc... I'm talking about the callback to user code.  As in the function pointer that you attach to nonlinear_solver->residual.  You know... the function that users write.  That one doesn't have any state information coming through to it.

Derek

On Fri, 30 Apr 2010, Derek Gaston wrote:

> On Apr 30, 2010, at 12:09 PM, Roy Stogner wrote:
>
>> We're currently passing back a pointer to the solver, ctx, and the
>> solver contains a pointer to the System, _system.
>
> Yes, this is what we get back from the callback from Petsc... I'm
> talking about the callback to user code.  As in the function pointer
> that you attach to nonlinear_solver->residual.  You know... the
> function that users write.  That one doesn't have any state
> information coming through to it.

Ah, I see.  I guess Ben just assumed that nobody would be reusing the
same residual function for multiple systems, but then
dial-an-operator type frameworks break that assumption.

Personally I'd prefer turning residual/jacobian/matvec into functor
objects over adding a System* or void* state argument.

One other possibility that Ben pointed out the other day, though: 
The DiffSolver framework originally was written to only work for
DiffSystem because it relied on being able to call assembly(bool,bool)
to request calculation of residuals, jacobians, or both at once, and
the other existing System subclasses couldn't do that.  But when all
the adjoint work required assembly(bool,bool) too, I moved the
definition up to ImplicitSystem and added implementations (sometimes
suboptimally efficient for backwards compatibility...) for
LinearImplicitSystem and NonlinearImplicitSystem.  So now we could
make DiffSolver work with any ImplicitSystem, you could use
PetscDiffSolver (or NewtonSolver) instead of PetscNonlinearSolver, and
since the assembly calls from there are class methods you'd have all
the state you need.

Thoughts?
---
Roy

On Apr 30, 2010, at 12:29 PM, Roy Stogner wrote:

> Ah, I see.  I guess Ben just assumed that nobody would be reusing the
> same residual function for multiple systems, but then
> dial-an-operator type frameworks break that assumption.

Exactly.

> Personally I'd prefer turning residual/jacobian/matvec into functor
> objects over adding a System* or void* state argument.

Hmmm.... I think I prefer the object-oriented approach over this... but I'm not married to it.

> One other possibility that Ben pointed out the other day, though: 
> The DiffSolver framework originally was written to only work for
> DiffSystem because it relied on being able to call assembly(bool,bool)
> to request calculation of residuals, jacobians, or both at once, and
> the other existing System subclasses couldn't do that.  But when all
> the adjoint work required assembly(bool,bool) too, I moved the
> definition up to ImplicitSystem and added implementations (sometimes
> suboptimally efficient for backwards compatibility...) for
> LinearImplicitSystem and NonlinearImplicitSystem.  So now we could
> make DiffSolver work with any ImplicitSystem, you could use
> PetscDiffSolver (or NewtonSolver) instead of PetscNonlinearSolver, and
> since the assembly calls from there are class methods you'd have all
> the state you need.

Well... I think I just need to move over to DiffSystem.... but that's not happening soon enough.

I think for now I'm just going to modify the signatures of the callback to provide the System being solved... and think about grander ideas later.

Objections?

Derek

On Fri, 30 Apr 2010, Derek Gaston wrote:

> On Apr 30, 2010, at 12:29 PM, Roy Stogner wrote:
>
>> Personally I'd prefer turning residual/jacobian/matvec into functor
>> objects over adding a System* or void* state argument.
>
> Hmmm.... I think I prefer the object-oriented approach over this...
> but I'm not married to it.

I went a little *too* OO with FEMSystem - in hindsight there should
have been a separate subclassable Physics object to attach which
included the residual evaluation functions, and requiring subclassing
of FEMSystem itself was more work and less flexible.

> I think for now I'm just going to modify the signatures of the
> callback to provide the System being solved... and think about
> grander ideas later.
>
> Objections?

Yes - doing so with the existing functions would break backwards
compatibility for other NonlinearSystem users; doing so by adding new
functions would make the API redundant and confusing, especially if we
had to do so again later for whatever grander idea we settle on.

If you're up to being a guinea pig on the DiffSolver(ImplicitSystem)
combination then let me know and I'll do that now; if not then I'd
suggest a private fork until we've got time to put the "grander ideas"
in place.
---
Roy

On Fri, Apr 30, 2010 at 2:14 PM, Roy Stogner <roystgnr@...> wrote:
>
> On Fri, 30 Apr 2010, Derek Gaston wrote:
>
>> Objections?
>
> Yes - doing so with the existing functions would break backwards
> compatibility for other NonlinearSystem users; doing so by adding new
> functions would make the API redundant and confusing, especially if we
> had to do so again later for whatever grander idea we settle on.

Any way to make the extra argument be a default NULL pointer or is
this something that can't really be overloaded in that way?

-- 
John

On Apr 30, 2010, at 1:14 PM, Roy Stogner wrote:

> I went a little *too* OO with FEMSystem - in hindsight there should
> have been a separate subclassable Physics object to attach which
> included the residual evaluation functions, and requiring subclassing
> of FEMSystem itself was more work and less flexible.

Yeah - but all we're talking about here is an abstract base class that consists of:

class ProblemInterface
{
public:
  virtual void residual() = 0;
  virtual void jacobian() = 0;
};

It's not going to be a huge burden...

> Yes - doing so with the existing functions would break backwards
> compatibility for other NonlinearSystem users; doing so by adding new
> functions would make the API redundant and confusing, especially if we
> had to do so again later for whatever grander idea we settle on.

I'm not certain we need a grander idea.  After thinking about it more... I really think just adding an argument might be good enough.  It really should have been there since the beginning.

About backward compatibility: come on!  This is adding _one_ argument to _two_ functions.  Most codes will be able to fix this in 10 seconds!  We can change it and announce it on libmesh-users and it will be fine.  If users are upset about that small of a change... then they should be using a libMesh release instead of svn....

> If you're up to being a guinea pig on the DiffSolver(ImplicitSystem)
> combination then let me know and I'll do that now; if not then I'd
> suggest a private fork until we've got time to put the "grander ideas"
> in place.

Private forks are the path to the darkside.  All my users pull sources direct from sourceforge... I really don't want to change that.

Derek

On Fri, 30 Apr 2010, Derek Gaston wrote:

> On Apr 30, 2010, at 1:14 PM, Roy Stogner wrote:
>
>       I went a little *too* OO with FEMSystem - in hindsight there should
>       have been a separate subclassable Physics object to attach which
>       included the residual evaluation functions, and requiring subclassing
>       of FEMSystem itself was more work and less flexible.
> 
> 
> Yeah - but all we're talking about here is an abstract base class that consists of:
> 
> class ProblemInterface
> {
> public:
>   virtual void residual() = 0;
>   virtual void jacobian() = 0;
> };
> 
> It's not going to be a huge burden...

Okay, I misunderstood what we were talking about then.  This is just
fine - it effectively *is* what I was calling a "functor" approach,
despite it not being a plain "overload operator()".  I thought you
were talking about requiring users to overload NonlinearSystem.

> About backward compatibility: come on!  This is adding _one_
> argument to _two_ functions.  Most codes will be able to fix this in
> 10 seconds!  We can change it and announce it on libmesh-users and
> it will be fine.  If users are upset about that small of a change...
> then they should be using a libMesh release instead of svn....

Fair enough.  I'll leave the decision up to Ben; he's the first person
I thought of who would be affected.
---
Roy

I thought about this... but I don't think there is a way.  I believe the signatures have to match exactly or the function pointer won't be accepted.

Could be wrong though...

Derek

On Apr 30, 2010, at 1:36 PM, John Peterson wrote:

> On Fri, Apr 30, 2010 at 2:14 PM, Roy Stogner <roystgnr@...> wrote:
>> 
>> On Fri, 30 Apr 2010, Derek Gaston wrote:
>> 
>>> Objections?
>> 
>> Yes - doing so with the existing functions would break backwards
>> compatibility for other NonlinearSystem users; doing so by adding new
>> functions would make the API redundant and confusing, especially if we
>> had to do so again later for whatever grander idea we settle on.
> 
> Any way to make the extra argument be a default NULL pointer or is
> this something that can't really be overloaded in that way?
> 
> -- 
> John
> 
> ------------------------------------------------------------------------------
> _______________________________________________
> Libmesh-devel mailing list
> Libmesh-devel@...
> https://lists.sourceforge.net/lists/listinfo/libmesh-devel

On Fri, 30 Apr 2010, Derek Gaston wrote:

> I thought about this... but I don't think there is a way.  I believe
> the signatures have to match exactly or the function pointer won't
> be accepted.

Ah, but that's what typecasts are for!

(joking!...)

> Could be wrong though...

I think the best you could do is overload the attachment methods, so
that one takes the first kind of function pointer and the other takes
the second, you store the two in different internal pointers, and then
at dispatch time you see which of the internal pointers are non-NULL.

Basically not worth it to avoid a tiny API compatibility breakage.
Especially when that breakage has the side benefit of punishing users
for not using DiffSystem.
(still joking! mostly! ...)
---
Roy

On Apr 30, 2010, at 1:44 PM, Roy Stogner wrote:

> Okay, I misunderstood what we were talking about then.  This is just
> fine - it effectively *is* what I was calling a "functor" approach,
> despite it not being a plain "overload operator()".  I thought you
> were talking about requiring users to overload NonlinearSystem.

Ah - gotcha.  Nope... not that.  But I'm not sure this new object is necessary at all.  I believe that Ben had talked about doing this before... possibly making all the callbacks look like this...

> Fair enough.  I'll leave the decision up to Ben; he's the first person
> I thought of who would be affected.

K - glad my stubbornness is overriding your good judgement!  ;-)

Derek

On Apr 30, 2010, at 1:59 PM, Roy Stogner wrote:

> Ah, but that's what typecasts are for!

Ewwww!

> (joking!...)

You'd better be!

> I think the best you could do is overload the attachment methods, so
> that one takes the first kind of function pointer and the other takes
> the second, you store the two in different internal pointers, and then
> at dispatch time you see which of the internal pointers are non-NULL.

Double ewwww!

Not to mention that we don't actually have setters for these things... we literally just do nonlinear_solver->residual = func_ptr!  So to go down this path we would do:

nonlinear_solver->new_residual_function_that_passes_more_info = func_ptr

;-)

> Basically not worth it to avoid a tiny API compatibility breakage.
> Especially when that breakage has the side benefit of punishing users
> for not using DiffSystem.
> (still joking! mostly! ...)

Heh.

Derek

On Fri, 30 Apr 2010, Derek Gaston wrote:

> On Apr 30, 2010, at 1:59 PM, Roy Stogner wrote:
>
>> Basically not worth it to avoid a tiny API compatibility breakage.

Okay, Ben says 'yes'.  Shoot us an email when you make the commit?  In
hindsight, it'll take him 10 seconds to fix his code, and it'll take
me 10 minutes to svn upgrade all our local libMesh modules to be
compatible with the fix and avoid BuildBot regression test failures,
so maybe I should have fought harder on this one...  ;-)
---
Roy

Ok - I actually got around to doing this... and just committed the changes.

Who wants to send an official looking email about the change to the user
list?

Derek

On Fri, Apr 30, 2010 at 4:44 PM, Roy Stogner <roystgnr@...:

>
>
> On Fri, 30 Apr 2010, Derek Gaston wrote:
>
>  On Apr 30, 2010, at 1:59 PM, Roy Stogner wrote:
>>
>>  Basically not worth it to avoid a tiny API compatibility breakage.
>>>
>>
> Okay, Ben says 'yes'.  Shoot us an email when you make the commit?  In
> hindsight, it'll take him 10 seconds to fix his code, and it'll take
> me 10 minutes to svn upgrade all our local libMesh modules to be
> compatible with the fix and avoid BuildBot regression test failures,
> so maybe I should have fought harder on this one...  ;-)
> ---
> Roy
>

Doh - are there not any Parameters on a System?  I was planning on stuffing
stuff into my NonlinearImplicitSystem so that I can get it back out when the
NonlinearImplicitSystem get's passed into compute_residual / jacobian.

Sigh.

For now I can key off of it's name... does someone see a better way to
recognize EXACTLY which NonlinearImplicitSystem is being solved in
compute_residual / jacobian?  I have a bunch of NonlinearImplicitSystems all
using the same compute_residual function...

Maybe there is something I'm missing.

Derek

On Tue, Jun 29, 2010 at 8:28 AM, Derek Gaston <friedmud@...> wrote:

> Ok - I actually got around to doing this... and just committed the changes.
>
> Who wants to send an official looking email about the change to the user
> list?
>
> Derek
>
>
> On Fri, Apr 30, 2010 at 4:44 PM, Roy Stogner <roystgnr@...:
>
>>
>>
>> On Fri, 30 Apr 2010, Derek Gaston wrote:
>>
>>  On Apr 30, 2010, at 1:59 PM, Roy Stogner wrote:
>>>
>>>  Basically not worth it to avoid a tiny API compatibility breakage.
>>>>
>>>
>> Okay, Ben says 'yes'.  Shoot us an email when you make the commit?  In
>> hindsight, it'll take him 10 seconds to fix his code, and it'll take
>> me 10 minutes to svn upgrade all our local libMesh modules to be
>> compatible with the fix and avoid BuildBot regression test failures,
>> so maybe I should have fought harder on this one...  ;-)
>> ---
>> Roy
>>
>
>

On Tue, 29 Jun 2010, Derek Gaston wrote:

> Doh - are there not any Parameters on a System?  I was planning on
> stuffing stuff into my NonlinearImplicitSystem so that I can get it
> back out when the NonlinearImplicitSystem get's passed into
> compute_residual / jacobian.  Sigh.

No, the Parameters object is just in the EquationSystems.

> Maybe there is something I'm missing.

*cough* FEMSystem *cough*

:-)

But seriously: it would be quite reasonable to add a Parameters object
to each System as well, if you want.  What I'd really like is a simple
API for querying parameters that would check the System first, then
check the EquationSystems second... maybe give a Parameters object an
optional Parameters* parent?
---
Roy

On Tue, Jun 29, 2010 at 11:05 AM, Derek Gaston <friedmud@...> wrote:
> Doh - are there not any Parameters on a System?  I was planning on stuffing
> stuff into my NonlinearImplicitSystem so that I can get it back out when
> theNonlinearImplicitSystem get's passed into compute_residual / jacobian.

Your Systems have an EquationSystems parent, right?  If you do have an
ES you can keep your parameters there, prefix them with "systemName_"
or something?

-- 
John

This is close to what I've done now.

If you guys think it would be useful to others I will add Parameters to System... and can even do the "parent" parameters thing Roy was talking about.

Derek

On Jun 29, 2010, at 10:14 AM, John Peterson wrote:

> On Tue, Jun 29, 2010 at 11:05 AM, Derek Gaston <friedmud@...> wrote:
>> Doh - are there not any Parameters on a System?  I was planning on stuffing
>> stuff into my NonlinearImplicitSystem so that I can get it back out when
>> theNonlinearImplicitSystem get's passed into compute_residual / jacobian.
> 
> Your Systems have an EquationSystems parent, right?  If you do have an
> ES you can keep your parameters there, prefix them with "systemName_"
> or something?
> 
> -- 
> John

On Tue, Jun 29, 2010 at 11:13 AM, Roy Stogner <roystgnr@...> wrote:
>
> But seriously: it would be quite reasonable to add a Parameters object
> to each System as well, if you want.  What I'd really like is a simple
> API for querying parameters that would check the System first, then
> check the EquationSystems second... maybe give a Parameters object an
> optional Parameters* parent?

Only thing I'd worry about is users wanting to redefine the 'built-in'
parameters libmesh uses ("linear solver tolerance", "nonlinear solver
tolerance", "linear solver maximum iterations", etc) on a per-System
basis while the internals of libmesh keep using the parent ES value
for those parameters.  Might lead to confusion/misunderstanding but
not bugs per se.

-- 
John

Actually... this would be a good reason to use Roy's suggestion of "fall-through" parameters.... it would allow you to _optionally_ override the values at the System level.  Of course that would mean changes to all of the solvers to look at the System level for parameters.

What we're really talking about here is Parameter polymorphism.  Has the idea that we store things like "linear_tolerance" as parameters outgrown itself?  If we had a linearTolerance() virtual method defined on Systems then it wouldn't be a problem... you would just ask a system for it's linearTolerance() and it would do the "right thing".

I'm not advocating this BTW... just tossing it out there for discussion.

Oh - I got my stuff working with what's current checked in.... just had to think about it a bit more.  I'm still up for changing things though.

Derek

On Jun 29, 2010, at 10:32 AM, John Peterson wrote:

> On Tue, Jun 29, 2010 at 11:13 AM, Roy Stogner <roystgnr@...> wrote:
>> 
>> But seriously: it would be quite reasonable to add a Parameters object
>> to each System as well, if you want.  What I'd really like is a simple
>> API for querying parameters that would check the System first, then
>> check the EquationSystems second... maybe give a Parameters object an
>> optional Parameters* parent?
> 
> Only thing I'd worry about is users wanting to redefine the 'built-in'
> parameters libmesh uses ("linear solver tolerance", "nonlinear solver
> tolerance", "linear solver maximum iterations", etc) on a per-System
> basis while the internals of libmesh keep using the parent ES value
> for those parameters.  Might lead to confusion/misunderstanding but
> not bugs per se.
> 
> -- 
> John
> 
> ------------------------------------------------------------------------------
> This SF.net email is sponsored by Sprint
> What will you do first with EVO, the first 4G phone?
> Visit sprint.com/first -- http://p.sf.net/sfu/sprint-com-first
> _______________________________________________
> Libmesh-devel mailing list
> Libmesh-devel@...
> https://lists.sourceforge.net/lists/listinfo/libmesh-devel

On Tue, 29 Jun 2010, Derek Gaston wrote:

> Actually... this would be a good reason to use Roy's suggestion of
> "fall-through" parameters.... it would allow you to _optionally_
> override the values at the System level.  Of course that would mean
> changes to all of the solvers to look at the System level for
> parameters.

Exactly.

> What we're really talking about here is Parameter polymorphism.

More like non-static parameters here; I can't think of any place
where we'd actually want to implement a tolerance lookup as anything
more sophisticated than "check what the user last set for our
tolerance".

> Has the idea that we store things like "linear_tolerance" as
> parameters outgrown itself?  If we had a linearTolerance() virtual
> method defined on Systems then it wouldn't be a problem... you would
> just ask a system for it's linearTolerance() and it would do the
> "right thing".

Yeah; I'd have done that long ago if it wasn't for the backwards
compatibility issue.  For flexibility there's no better solution than
Parameters, but for things we know we need we might as well avoid the
string lookup.
---
Roy

On Tue, 29 Jun 2010, Roy Stogner wrote:

> Yeah; I'd have done that long ago if it wasn't for the backwards
> compatibility issue.  For flexibility there's no better solution than
> Parameters, but for things we know we need we might as well avoid the
> string lookup.

Actually, for backwards compatibility, it would be possible to use the
same solution as in ParameterVector: add the option to add a parameter
which is actually located elsewhere in memory; then that parameter
could be accessed directly as a System member variable or indirectly
via a Parameters lookup.

Probably not worth the effort though.  A string lookup per solve is no
big deal.  If anyone is in a string lookup per quadrature point
situation, they could be caching that already without library support.
---
Roy.