Subject: [Ffc] FFC code generation for inline coefficient evaluation (i.e. coefficients defined as expressions in UFL)

From: Florian Rathgeber <[email-masked]>
Date: Fri, 23 Jul 2010 13:40:15 +0200

Hi,

Profiling DOLFIN assembly shows that for forms with many coefficients or
computationally expensive expressions assigned to coefficients, a great
deal of the whole assembly is spent in expression evaluation.

An idea I wanted to investigate was giving an expresssion in UFL and
have FFC generate corresponding code to evaluate it inline during tensor
tabulation. That should have the potential of eliminating some of the
overhead.

However, currently that is only possible within very limited bounds. FFC
does not yet support UFL conditionals, which essentially limits possible
expressions to depend on coordinates only and be the same over the whole
domain (which is not very useful for most form coefficients, e.g. source
terms).

This brings me to my questions:
1) Do you see a reasonable chance for speeding up assembly in the way
described?
2) Is my assessement correct? Or do I miss some capabilities of FFC that
support my use case?
3) Since UFL basically support arbitrary python syntax, would it be
reasonable to support translation for user-defined Expressions given in
python (either overloading eval in a class defived from Expression or
passing C++ code to the SWIG constructor of Expression)? FFC would
generate code to evaluate these Expressions inline during tensor
tabulation. That would of course require the Expression to be stateless
and have not input other than the coordinates. But at least it would
provide much more flexibility and also consistency with the DOLFIN
interface in general
(sidenote: it is possible to specify Expressions in this way in UFL, but
they are treated just like Coefficients by FFC)

Illustration example (UFL file):
  from dolfin import Expression

  code_expr = """
  class Source : public Expression
  {
  public:
    void eval(Array<double>& values, const Array<double>& x) const
    {
      values[0] = ((x[0]-0.25)*(x[0]-0.25)+(x[1]-0.5)*(x[1]-0.5) <
0.025) ? 1.0 : 0.0;
    }
  };
  """

  scalar = FiniteElement("Lagrange", triangle, 1)

  v  = TestFunction(scalar)
  u  = TrialFunction(scalar)
  u0 = Coefficient(scalar)
  f = Expressions(code_expr)

  c = 0.05
  k = 0.1
  alpha = 1.0

  a = v*u*dx + k*alpha*u*v*dx + k*c*dot(grad(v), grad(u))*dx
  L = v*u0*dx + k*v*f*dx

Florian

---

From: Kristian Ølgaard <[email-masked]>
Date: Fri, 23 Jul 2010 12:51:07 +0100

On 23 July 2010 12:40, Florian Rathgeber
<[email-masked]> wrote:
> -----BEGIN PGP SIGNED MESSAGE-----
> Hash: SHA1
>
> Hi,
>

Did you have a look at the dolfin/demo/pde/spatial-coordinates/cpp
demo? It is essentially the dolfin/demo/pde/poisson/cpp demo but it
uses source and boundary terms defined using coordinates and is not
the same over the whole domain (I don't really get your point here).
The coordinates are defined as the coordinates of the current
integration point located on the current cell which is being
integrated. Have a look at the generated code, then it might be
clearer what's going on.

Kristian

> -----BEGIN PGP SIGNATURE-----
> Version: GnuPG v2.0.12 (MingW32)
>
> iEYEARECAAYFAkxJf58ACgkQ8Z6llsctAxYHqACgtzcgezJn1VrCFp9N7UcE40kq
> 9/oAoLJPMhaosQvthrOy0dGAgOkNiYh5
> =KB3/
> -----END PGP SIGNATURE-----
>
>
> _______________________________________________
> Mailing list: https://launchpad.net/~ffc
> Post to     : [email-masked]
> Unsubscribe : https://launchpad.net/~ffc
> More help   : https://help.launchpad.net/ListHelp
>
>

---

From: Florian Rathgeber <[email-masked]>
Date: Fri, 23 Jul 2010 14:02:34 +0200

I guess I haven't made my point very clear. The demo you gave is a good
example of what's currently possible and where the limit is I think.

Say you take dolfin/demo/pde/spatial-coordinates/cpp and want to
restrict the source term f to only a circle in the center of the domain.
How would you do that? I wouldn't know another way than a conditional
and I couldn't find a way to get that working with FFC.

I was not talking about constants, but about expressions that are not
uniform over the whole domain, should have made that more clear.

Florian

---

From: Kristian Ølgaard <[email-masked]>
Date: Fri, 23 Jul 2010 13:22:34 +0100

> and I couldn't find a way to get that working with FFC.

Now it makes more sense, I think you would need a conditional indeed.

>> provide much more flexibility and also consistency with the DOLFIN

I think it sounds complicated, much better to defined the function
directly in terms of coordinates.

>> they are treated just like Coefficients by FFC)

I don't think this is the case, it is possible in the Python interface
to DOLFIN, but not in UFL.

>>      values[0] = ((x[0]-0.25)*(x[0]-0.25)+(x[1]-0.5)*(x[1]-0.5) <

In UFL you can do:

f = conditional( (x[0]-0.25)**2 + (x[1]-0.5)**2 <  0.025), 1.0, 0.0  )

although no output is generated yet because Conditionals are not
supported yet, but the code which will be generated should look an
awful lot like what you have i.e.

f = ((x[0]-0.25)*(x[0]-0.25)+(x[1]-0.5)*(x[1]-0.5) < 0.025) ? 1.0 : 0.0;

Kristian

> iEYEARECAAYFAkxJhNoACgkQ8Z6llsctAxYaOgCeNN1gpBh2h3QzTrnevHMcKcXn
> gXAAn0SaoqOou75hAayf8+Y8Y5cVZurn
> =tNNo

---

From: Florian Rathgeber <[email-masked]>
Date: Fri, 23 Jul 2010 14:30:53 +0200

>> to DOLFIN, but not in UFL.

You are right, it is not in the UFL specification. It is technically
possible with FFC though. I passed the example UFL file I gave to FFC
and it was accepted, however the generated code was the very same as if
I had just written

  f = Coefficient(scalar)

instead of

  f = Expression(code_expr)

Florian

---

From: "Garth N. Wells" <[email-masked]>
Date: Fri, 23 Jul 2010 13:36:10 +0100

> tabulation. That should have the potential of eliminating some of the

There is some scope for this, but its potential will be limited for
complicated problems where the coefficient lies is another solution
vector. An issue at the moment is that with the DOLFIN Function design,
a number of virtual function calls, and call-back function call from the
generated code, are required to just get a scalar coefficient. 

Off the top of my head, some areas that I think can be addressed are:

- Explore use of templates to limit virtual calls

- Tabulating an arbitrary number of element tensors at each
ufc::tabulate_tensor call. This might make the cost of the virtual
function calls negligible.

Garth

---

From: Kristian Ølgaard <[email-masked]>
Date: Fri, 23 Jul 2010 19:33:27 +0100

> and it was accepted, however the generated code was the very same as if

No, it is technically possible in PyDolfin. Not UFL or FFC
Anyway, I've implemented support for conditionals see if it works.
(dolfin/demo/pde/conditional)
I only tested '<' so far but the other operators should work the same way.

Kristian

> iEYEARECAAYFAkxJi3wACgkQ8Z6llsctAxZoNwCgg0mSVE6b7jupVWLKyoQoFrbZ
> e/wAoIXtwpNyC+40gF1DVG9UGmvXBnCa
> =Z4Co

---

From: Florian Rathgeber <[email-masked]>
Date: Mon, 02 Aug 2010 14:08:06 +0200

>> I only tested '<' so far but the other operators should work the same way.

Thank you for implementing that and sorry for not coming back to this
earlier. As far as I can see, it is only implemented for quadrature. Is
there a fundamental problem that would inhibit implementing it for
tensor representation as well?

Florian

---

From: "Garth N. Wells" <[email-masked]>
Date: Mon, 02 Aug 2010 14:17:09 +0100

> tensor representation as well?

Yes. The tensor contraction relies on computing everything which is not
related to cell geometry a priori, so it won't be able to handle
conditionals.

Garth 

> iEYEARECAAYFAkxWtSYACgkQ8Z6llsctAxZwpACfaAKeewvmEknig7ySe9nxGaM1
> 2ZwAn07Mkka+BMcvbuVlWMMV1iVRF0QP
> =0E1w