dolfin team mailing list archive

["Garth N. Wells" <gnw20@xxxxxxxxx>]

Johan Hake wrote:
> On Monday 25 January 2010 10:01:18 Garth N. Wells wrote:
>> Johan Hake wrote:
>>> On Monday 25 January 2010 03:10:46 Garth N. Wells wrote:
>>>> Johan Hake wrote:
>>>>> 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?
>>>> Can I somehow get a 'reference' to the vector so I don't have to use
>>>> u.vector()[:] in the expressions?
>>> That is what u.vector() gives you.
>> I figured it out eventually. There's nothing quite like a '&' symbol . .
>> . . .
> 
> Good!
> 
> In Python every name in a namespace is a reference(pointer) to an instance. 
> Johan
> 
>> Garth
>>
>>> Johan
>>>
>>>> Garth
>>>>
>>>>> 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 :)
>>>>>

I have things working nicely now when working with GenericVector. I have:

  u_vec  = u.vector()
  u0_vec = u0.vector()
  v0_vec = v0.vector()
  a0_vec = a0.vector()

  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

where beta, gamma and dt are floats. I'd like now to know a little more
about what's going on behind the scenes, especially with respect to the
creation of temporaries. Do all the operations just call the relevant
functions from the GenericVector interface? Will each operation create a
temporary?

Garth


>>>>> 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