← Back to team overview

dolfin team mailing list archive

Re: Initialization of discrete functions (related to bug 24)

 



Anders Logg wrote:
On Tue, Apr 22, 2008 at 09:21:48PM +0200, Dag Lindbo wrote:
Anders Logg wrote:
On Mon, Apr 21, 2008 at 09:35:30PM +0200, Dag Lindbo wrote:
Why do I care? I'm looking at the current state of the LinearPDE class.
Here, the solution vector is a member function, and the solution
Function is initialized with a reference to this vector. I.e. the
LinearPDE (not the DiscreteFunction) owns the solution vector. This is
of course a problem if the LinearPDE goes out of scope before the
solution function (bug 24). It is also a bit counter-intuitive.

If it was possible to do an initialization like above, then the solve()
method would simply do
u.init(mesh, form, 1);
GenericVector& x = u.vector();
<solve>
There is a member local_vector in DiscreteFunction that could be used
for this. If this is nonzero, then the DiscreteFunction owns its data.

See if you can figure out a good way for the DiscreteFunction to know
that it should take responsibility for the vector here.

Right... I don't see how to do this without breaking the encapsulation of
both Function and DiscreteFunction (making LinearPDE friend in both).
Bundle attached.

Garth, does this look OK?

In essence, the member x is removed from the LinearPDE class. in
LinearPDE::solve I do

  (...)
  Vector b;
  Vector* x = new Vector();
  (... call solver etc)

  u.init(mesh, *x, a, 1);
  DiscreteFunction& uu = dynamic_cast<DiscreteFunction&>(*u.f);
  uu.local_vector = x;

/Dag
Looks like a good temporary solution to me.

I expect when we're done and happy with the linear algebra classes
(which should be soon), we will have a similar party with the function
classes... :-)

The fact that the pointer to x goes out of scope and trusts the DiscreteFunction to clean up the vector makes me a bit queasy (this can probably break in the future, causing a big leak). You should probably verify independently that nothing broke because of this revised LinearPDE (but for me it has been working nicely today).

On a more general note, I'm pretty happy with the state of the LA now that op[] is back and the down_cast<Foo> mechanism has been explained.

As a benchmark I ran a Inc. Navier-Stokes solver (based on the old module, which uses a lot of 'manual' LA) and got a modest slowdown:
(tip yesterday)    53.071u 6.708s
(0.7.2)            50.803u 1.700s
It may also be worth mentioning that I got the same numerical values (for some residual norm) even after hundreds of time steps.

/Dag

ok. It might be possible to find some optimizations when things have
stabilized.


Try accessing the underlying data type directly and operate on it. It will avoid the overhead of virtual function call and more inlining is possible.

Garth


Follow ups

References