dolfin team mailing list archive

[Johan Hake <johan.hake@xxxxxxxxx>]

On Sunday 24 January 2010 16:39:14 Garth N. Wells wrote:
> Johan Hake wrote:
> > On Sunday 24 January 2010 16:12:37 Garth N. Wells wrote:
> >> Johan Hake wrote:
> >>> On Sunday 24 January 2010 00:03:41 Garth N. Wells wrote:
> >>>> Johan Hake wrote:
> >>>>> On Saturday 23 January 2010 14:55:08 Garth N. Wells wrote:
> >>>>>> Johan Hake wrote:
> >>>>>>> On Saturday 23 January 2010 08:42:14 Garth N. Wells wrote:
> >>>>>>>> Is it correct that behind the scenes that
> >>>>>>>>
> >>>>>>>>   U0 = Function(V)
> >>>>>>>>   U  = Function(V)
> >>>>>>>>   U0.vector()[:] = U.vector()[:]
> >>>>>>>>
> >>>>>>>> involves a GenericVector::get(..) call and a
> >>>>>>>> GenericVector::set(..) call? If so, it isn't ideal since it
> >>>>>>>> introduces unnecessary new/delete operations and unnecessary
> >>>>>>>> copying of data.
> >>>>>>>
> >>>>>>> None of GenericVector::get(..) or GenericVector::set(..) are
> >>>>>>> invoked, see __getslice__ and __setslice__ in la_post.i.
> >>>>>>>
> >>>>>>>    U0.vector()[:]
> >>>>>>>
> >>>>>>> involves
> >>>>>>>
> >>>>>>>   GenericVector::operator =(..)
> >>>>>>>
> >>>>>>> and
> >>>>>>>
> >>>>>>>   U.vector()[:]
> >>>>>>>
> >>>>>>> involves
> >>>>>>>
> >>>>>>>   GenericVector::copy()
> >>>>>>>
> >>>>>>> However the latter is unnecessary as you instead can do:
> >>>>>>>
> >>>>>>>   U0.vector()[:] = U.vector()
> >>>>>>>
> >>>>>>> invoking the assignment operator of U0's vector with U's vector.
> >>>>>>
> >>>>>> What happens if I do
> >>>>>>
> >>>>>>   x = U.vector()[:]
> >>>>>
> >>>>> It just triggers the copy method of GenericVector, which is the same
> >>>>> behavior as for other itterable Python types.
> >>>>>
> >>>>>> ? Is x a numpy array?
> >>>>>
> >>>>> No you need to call array() to accomplish that.
> >>>>
> >>>> OK. What I'm trying to do is
> >>>>
> >>>>     # Get vectors
> >>>>     u_vec  = u.vector()[:]
> >>>>     u0_vec = u0.vector()[:]
> >>>>     v0_vec = v0.vector()[:]
> >>>>     a0_vec = a0.vector()[:]
> >>>
> >>> You should not need to make a copy of the vectors here.
> >>
> >> How can I avoid it?
> >
> > a_vec and v_vec are new vectors. None of the four vectors below get
> > modified by the a_vec and v_vec expressions so no need of copying, and
> > the v0 and a0 assignment should work with GenericVectors too.
> 
> Do you mean that just
> 
>   a_vec = 1.0/(2.0*beta)*((u - u0 - v0*dt)/(0.5*dt*dt) \
>            - (1.0-2.0*beta)*a0 )
> 
> where u and u0 are GenericVectors should work?

Have you tried? 

As long as the rest (besides a0, which I assume also is a GenericVector) are 
scalars everything should just work. The Python LA interface (at least for 
GenericVector) should work more or less as the NumPy interface which I think 
is nice :)

We cannot take 1./v, where v is a GenericVector. 
 
> > Do you get any error messages?
> 
> No errors. What I have now seems to work fine.

Ok, and that is because what you do obviously works for NumPy arrays.

Johan

> Garth
> 
> >>     # Update acceleration and velocity
> >>     a_vec = (1.0/(2.0*beta))*( (u_vec - u0_vec - v0_vec*dt)/(0.5*dt*dt)\
> >>           - (1.0-2.0*beta)*a0_vec )
> >>     v_vec = dt*((1.0-gamma)*a0_vec + gamma*a_vec) + v0_vec
> >>
> >>
> >>
> >> Garth
> >>
> >>>>     a_vec = 1.0/(2.0*beta)*((u_vec-u0_vec - v0_vec*dt)/(0.5*dt*dt) \
> >>>>          - (1.0-2.0*beta)*a0_vec)
> >>>>     v_vec = dt*((1.0-gamma)*a0_vec + gamma*a_vec) + v0_vec
> >>>
> >>> The expression cause a lot of intermediate vectors, but should work, if
> >>> beta, gamma and dt are all scalars.
> >>>
> >>>>     # Update (U0 <- U0)
> >>>>     v0.vector()[:] = v_vec
> >>>>     a0.vector()[:] = a_vec
> >>>
> >>> This should work.
> >>>
> >>> Johan
> >>>
> >>>> I guess I need to call array() to make this work?
> >>>>
> >>>> Garth
> >>>>
> >>>>> Johan
> >>>>>
> >>>>>> Garth
> >>>>>>
> >>>>>>> Johan
>