← Back to team overview

dolfin team mailing list archive

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

 

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


Follow ups

References