dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

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.

-- 
Anders