AXPY for Dolfin functions?

-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1

On 10/14/2012 10:30 PM, Nico Schlömer wrote:
> New question #211211 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/211211
>
> Given two Dolfin Functions out of the same space, let's say
> FunctionSpace(mesh, 'Lagrange', 1), is it possible to scale and add
> those functions?
>
> w = alpha * u + beta * v

You can accomplish that by using the axpy member function of the
underlying coefficient vectors, for instance

_u.vector()->axpy(1-_omega,*_u0.vector());

where, of course, _u and _u0 are two function objects belonging to the
same function space.

HTH,

- --
Andre

>
> --Nico
>
-----BEGIN PGP SIGNATURE-----
Version: GnuPG v1.4.12 (GNU/Linux)
Comment: Using GnuPG with Mozilla - http://www.enigmail.net/

iQEcBAEBAgAGBQJQeyzIAAoJEA79ggnbq9dmxNYH/jqlH2q1x6LxadwYWBcmmSRw
Zg6UAfvHPv/gKwiT2x8wYkW9IazzR1KvBQ6ojgDlCQ5tb8Mv7PvJWRPeRq0nMz/N
4aHX/4IpXQopCObNp46o3m+Edog2Gsp5Ddpomj+SPeEn6YlCbFB2UipzOSjTbT52
JHcRQHjWraz+hd2a6QV7OgDHzyOYat8SK5dHCTujtJnOM/NQ4L0HRmqq/RxmSr47
uiREQclYNlChgi+fu7LAzI60ewiDahHszx+X/QRbGmop+Va74Z0JGZ7P/2l2HrAs
MMTQt+luJ0bfoIc6gLzRYFzmgGsBXu5kjz6IeYYGe8oo8syWnfnfAOTGGgbnFAk=
=SD1u
-----END PGP SIGNATURE-----

It seems that, as usual, solve() is your friend. An expression

solve(u*v*dx == (sin(t)*Tr*v + cos(t)*Ti*v)*dx, A)

will perform

sin(t)*Tr + cos(t)*Ti

and return the result in A.

Whoops, cross-posting here.

@Andre:
Is AXPY save? I'm don't know how the coefficients of, for example, higher-order elements are stored so I don't really know if adding the underlying vectors would actually add the functions that belong to them.

On Sun, Oct 14, 2012 at 11:31 PM, Nico Schlömer
<email address hidden> wrote:
> Question #211211 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/211211
>
> Nico Schlömer posted a new comment:
> Whoops, cross-posting here.
>
> @Andre:
> Is AXPY save? I'm don't know how the coefficients of, for example, higher-order elements are stored so I don't really know if adding the underlying vectors would actually add the functions that belong to them.

I assumed alpha and beta to be constants in your code snippet. If you
want them to be non-constant functions, then vector.axpy() cannot
be used. You could write your own small vector routine where you
interpolate your scaling functions into the the *same* space first and
then do the scaling and multiplication "pointwise". But your approach
is more elegant :)

Best
Andre

>
> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.

Oh, what I called sin(t) and cos(t) are actually constants in space, sorry
for the confusion.
I would guess that adding the underlying vectors serves you well with
lagrange-1 elements, probably of all orders. I perform this stunt for
nedelec elements, and there I'm not so sure anymore already.
On Oct 14, 2012 11:55 PM, "Andre Massing" <
<email address hidden>> wrote:

> Your question #211211 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/211211
>
> Andre Massing posted a new comment:
> On Sun, Oct 14, 2012 at 11:31 PM, Nico Schlömer
> <email address hidden> wrote:
> > Question #211211 on DOLFIN changed:
> > https://answers.launchpad.net/dolfin/+question/211211
> >
> > Nico Schlömer posted a new comment:
> > Whoops, cross-posting here.
> >
> > @Andre:
> > Is AXPY save? I'm don't know how the coefficients of, for example,
> higher-order elements are stored so I don't really know if adding the
> underlying vectors would actually add the functions that belong to them.
>
> I assumed alpha and beta to be constants in your code snippet. If you
> want them to be non-constant functions, then vector.axpy() cannot
> be used. You could write your own small vector routine where you
> interpolate your scaling functions into the the *same* space first and
> then do the scaling and multiplication "pointwise". But your approach
> is more elegant :)
>
> Best
> Andre
>
> >
> > --
> > You received this question notification because you are a member of
> > DOLFIN Team, which is an answer contact for DOLFIN.
>
> --
> You received this question notification because you asked the question.
>

If the two functions to be scaled and added are in the same function space, then scaling and adding their vectors as Andre originally suggested is safe and a good approach. This should be true for all the families of element spaces supported.

--
Marie

On 15. okt. 2012, at 00:11, Nico Schlömer <email address hidden> wrote:

> Question #211211 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/211211
>
> Nico Schlömer posted a new comment:
> Oh, what I called sin(t) and cos(t) are actually constants in space, sorry
> for the confusion.
> I would guess that adding the underlying vectors serves you well with
> lagrange-1 elements, probably of all orders. I perform this stunt for
> nedelec elements, and there I'm not so sure anymore already.
> On Oct 14, 2012 11:55 PM, "Andre Massing" <
> <email address hidden>> wrote:
>
>> Your question #211211 on DOLFIN changed:
>> https://answers.launchpad.net/dolfin/+question/211211
>>
>> Andre Massing posted a new comment:
>> On Sun, Oct 14, 2012 at 11:31 PM, Nico Schlömer
>> <email address hidden> wrote:
>>> Question #211211 on DOLFIN changed:
>>> https://answers.launchpad.net/dolfin/+question/211211
>>>
>>> Nico Schlömer posted a new comment:
>>> Whoops, cross-posting here.
>>>
>>> @Andre:
>>> Is AXPY save? I'm don't know how the coefficients of, for example,
>> higher-order elements are stored so I don't really know if adding the
>> underlying vectors would actually add the functions that belong to them.
>>
>> I assumed alpha and beta to be constants in your code snippet. If you
>> want them to be non-constant functions, then vector.axpy() cannot
>> be used. You could write your own small vector routine where you
>> interpolate your scaling functions into the the *same* space first and
>> then do the scaling and multiplication "pointwise". But your approach
>> is more elegant :)
>>
>> Best
>> Andre
>>
>>>
>>> --
>>> You received this question notification because you are a member of
>>> DOLFIN Team, which is an answer contact for DOLFIN.
>>
>> --
>> You received this question notification because you asked the question.
>>
>
> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.

Okay, sweet. Thanks for letting me know!

Gah. I'm getting bitten by the complex values against.

W = V*V
T = Function(W)
Tr, Ti = T.split()
Tr.vector()

yields

*** Error: Unable to access vector of degrees of freedom.
*** Reason: Cannot access a non-const vector from a subfunction.
*** Where: This error was encountered inside Function.cpp.
*** Process: 0

I guess I'm going to have to stick to solve(). :/

Alternatively, I could copy over (assign()) the parts:

        tmpr.assign(Tr)
        vr = tmpr.vector()
        vr *= np.sin(t)
        tmpi.assign(Ti)
        vi = tmpi.vector()
        vi *= np.cos(t)
        vr += vi

From my dolfin-adjoint perspective:

It would be really nice to have a better syntax for doing this, over accessing .vectors() or calling solve():

* Calling solve is inefficient, because it goes through the whole machinery of assembling and solving matrices, when none of it is necessary.

* Modifying the vector obscures the mathematical intent of the statement. If there was a method (I'm not proposing this interface, just trying to make things concrete) like

  w.linear_combination(zip([alpha, beta], [u, v]))

then dolfin-adjoint could annotate that as one equation, and adjoin it. But if it's written as a sequence of operations like

  w.vector().zero()
  w.vector().axpy(alpha, u.vector())
  w.vector().axpy(beta, v.vector())

then from my perspective it's much harder to see what's going on and annotate it. The mathematical intent is distributed over more lines of code. It would be possible, but it's much less clean, and so I'm reluctant to do it.

As an anecdote, the quasi-geostrophic model does a lot of such linear combinations, and because of the above when we want to use the adjoint we have to call solve(). This causes a significant inefficiency.

Would it be possible to add an expressive method to the Function class for assigning linear combinations of other functions?

On 10/15/2012 01:35 AM, Nico Schlömer wrote:
> Question #211211 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/211211
>
> Nico Schlömer posted a new comment:
> Gah. I'm getting bitten by the complex values against.
>
> W = V*V
> T = Function(W)
> Tr, Ti = T.split()
> Tr.vector()
>
> yields
>
> *** Error: Unable to access vector of degrees of freedom.
> *** Reason: Cannot access a non-const vector from a subfunction.
> *** Where: This error was encountered inside Function.cpp.
> *** Process: 0
>
> I guess I'm going to have to stick to solve(). :/
>

Nico, did you try

   T.split(deepcopy=True)

?

cf for instance https://answers.launchpad.net/dolfin/+question/199077 or
http://fenicsproject.org/documentation/dolfin/1.0.0/python/programmers-reference/functions/function/Function.html#dolfin.functions.function.Function.split

--
Marie

@Marie
deepcopy works of course. Whether I copy the values when doing split() or whether I call tmpr.assign(Tr) later on doesn't make a difference. The linear combination only happens optionally (it's for plotting), so I'd prefer not to deepcopy immediately.

@Patrick
I principally agree with you. Following FEniCS philosophy of having a clean interface to common FEM operations would certainly favor a modification here.
Something as general as adding functions like
   F2 = alpha * F + beta * G
or
  F2 *= 4.0
  F2 += G
would probably not work since F and G might be of different element families. Throwing an exception if the families are different would work.
Such functions would also have to handle the case where you add two functions of hte same element family with different order. -- Not sure if there already is functionality for interpreting a Lagrange-1 function as a Lagrange-2.

On 10/15/2012 01:01 PM, Nico Schlömer wrote:
> Question #211211 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/211211
>
> Nico Schlömer posted a new comment:
> @Marie
> deepcopy works of course. Whether I copy the values when doing split() or whether I call tmpr.assign(Tr) later on doesn't make a difference. The linear combination only happens optionally (it's for plotting), so I'd prefer not to deepcopy immediately.
>
> @Patrick
> I principally agree with you. Following FEniCS philosophy of having a clean interface to common FEM operations would certainly favor a modification here.
> Something as general as adding functions like
> F2 = alpha * F + beta * G
> or
> F2 *= 4.0
> F2 += G
> would probably not work since F and G might be of different element families. Throwing an exception if the families are different would work.
> Such functions would also have to handle the case where you add two functions of hte same element family with different order. -- Not sure if there already is functionality for interpreting a Lagrange-1 function as a Lagrange-2.

The problem is not different element families, but rather mixing the
abstractions. This is particular eminent in the Python interface of
DOLFIN where UFL and DOLFIN classes are mixed. F2 would in the above
example be a UFL expression, as a (Generic)Function in DOLFIN does not
have any linear algebra operators. I guess one therefore need to add
such functionality more explicitly like Patrick suggested.

Johan