Subject: [UFL-dev] control statements in UFL

From: Kristian Oelgaard <[email-masked]>
Date: Wed, 20 Aug 2008 09:26:28 +0200

Hi,

Garth, Anders and I have discussed the possibility of defining control
statements in FFC/UFL. This can be useful for many things but if we
take DG advection diffusion as an example we have the following on an
exterior facet:

scalar = FiniteElement("Discontinuous Lagrange", "triangle", 1)
vector = VectorElement("Lagrange", "triangle", 2)

v = TestFunction(scalar)
u = TrialFunction(scalar)

b   = Function(vector)
n   = FacetNormal("triangle")
# outflow facet switch, computed in DOLFIN (has values 1.0 or 0.0)
of  = Function(constant)

a = dot(mult(v, n), mult(b, of*u))*ds

This works, but means that one has to compute the value of 'of' in DOLFIN.
This is currently handled by OutflowFacet in SpecialFunctions.h but to me
it seems silly to add an additional argument to a form the value of which
depends on the value of other arguments to the form. Wouldn't it be much
nicer to let FFC/UFL produce code that can decide which values to use in
tabulate_tensor()?

So the question is how we can do this. I think the following could work
but I don't know for sure and I definitely do not know if it's the best
solution:

A new FFC/UFL type with possible interface like this:

Condition( left, logical_operator, right, return_value, else (optional) )

that can work as a placeholder for an expression or alternatively just
represent a real number (a switch/bool).

The syntax in case of a real number for advection-diffusion:

condition = Condition(dot(n, b), ">", 0.0, 1.0)
The logical operator ">" should be in Python syntax, then we can always
map it in ufcformat.py (or any other format)

Then

a = dot(mult(v, n), mult(b, condition*u))*ds

and in tabulate_tensor() we would have code that looks something like:

cond0 = 0.0;
if (w[0][0]*w[1][0] + w[0][1]*w[1][1] > 0.0)
   cond0 = 1.0

A[0] = .... cond0* .....

Alternatively, condition could work like a placeholder

condition = Condition(dot(n, b), ">", 0.0, mult(b,u))

a = dot(mult(v, n), condition))*ds

The code would now be:

cond0[0] = 0.0;
cond0[1] = 0.0;

if (w[0][0]*w[1][0] + w[0][1]*w[1][1] > 0.0)
{
   cond0[0] = w[1][0]*...
   cond0[1] = w[1][1]*...
}

For a given tensor we can keep track of the number of conditions through
its own index and check if any conditions are present like we do for
coefficients, transforms, derivatives etc.

How would this fit into the UFL framework? Could it even work, or does
anyone see a cleaner way of doing it?

Kristian

---

From: "Martin Sandve Alnæs" <[email-masked]>
Date: Wed, 20 Aug 2008 10:02:37 +0200

> Condition( left, logical_operator, right, return_value, else (optional) )

I'm not a fan of the "else (optional)". I didn't understand
the example below intuitively before I read this line.

Maybe we can have named conditional operators (mirroring fortran):

  eq(left, right) # ==
  ne(left, right) # !=
  le(left, right) # <=
  lt(left, right) # <
  ge(left, right) # >=
  gt(left, right) # >

And the conditional can then be more "structured" (shorter argument
list, in a way):

  c = Conditional(a, lt(left, right), b)

mirroring "c = a if (left<right) else b" in python,
or:

  c = Conditional(lt(left, right), a, b)

mirroring "c = (left<right) ? a : b" in C++.

Which do you prefer?

> represent a real number (a switch/bool).

a and b can easily be any UFL expression, while left and right must of
course be scalar.

> condition = Condition(dot(n, b), ">", 0.0, 1.0)

I my syntax suggestions this will be:
  condition = Conditional(gt(dot(n, b), 0.0), 1.0, 0.0)
or
  condition = Conditional(1.0, gt(dot(n, b), 0.0), 0.0)

or written differently:

  cond = gt(dot(n, b), 0.0)
  condition = Conditional(cond, 1.0, 0.0)
or
  condition = Conditional(1.0, cond, 0.0)

> a = dot(mult(v, n), mult(b, condition*u))*ds

In UFL:

  a = dot(v*n, b*condition*u)*ds

> condition = Condition(dot(n, b), ">", 0.0, mult(b,u))

What's the difference?

> if (w[0][0]*w[1][0] + w[0][1]*w[1][1] > 0.0)

Ok.

> coefficients, transforms, derivatives etc.

What do you mean? Sounds like FFC details, maybe I don't need to know.

> anyone see a cleaner way of doing it?

I think this should be ok in UFL, although of course it causes
some mathematical issues w.r.t. differentiation. I guess that
often the derivative of a Conditinal can simply be another
Conditinal with the same condition and the derivatives of
the other arguments. We just have to be careful.

I'll implement this in UFL, and see if I spot any difficulties.

--
Martin

---

From: Anders Logg <[email-masked]>
Date: Wed, 20 Aug 2008 10:14:53 +0200

>   c = Conditional(lt(left, right), a, b)

I think this second option looks better.

But wouldn't it be possible to just overload <, >, = etc for UFLObject
so that one may do

  c = Conditional(left < right, a, b)

or even

  a = dot(v*n, (left < right)*u + (1 - left < right)*sin(u))*ds

We could have left < right return a Conditional object which would
return either 1 or 0 (in the generated code) depending on whether the
comparison is true or false.

-- 
Anders

---

From: "Martin Sandve Alnæs" <[email-masked]>
Date: Wed, 20 Aug 2008 10:30:09 +0200

>  c = Conditional(left < right, a, b)

The reason we can't do that is that those operators
(at least ==) are used for other things, like putting
expressions in sets and dicts which we cannot avoid.
If we allow it for <, >, <=, >=, but not ==, people
will be confused and make errors.

> comparison is true or false.

1 or 0 is a very special case, we should have the Conditional(condition, a, b).

-- 
Martin

---

From: Kristian Oelgaard <[email-masked]>
Date: Wed, 20 Aug 2008 10:48:52 +0200

> >>   c = Conditional(lt(left, right), a, b)

I think that this looks like a good solution if it's possible to implement.

Kristian

> _______________________________________________
> UFL-dev mailing list
> [email-masked]
> http://fenics.org/mailman/listinfo/ufl-dev
>

---

From: Anders Logg <[email-masked]>
Date: Wed, 20 Aug 2008 10:53:44 +0200

> will be confused and make errors.

ok.

-- 
Anders

---

From: "Martin Sandve Alnæs" <[email-masked]>
Date: Wed, 20 Aug 2008 15:56:43 +0200

This is now implemented in UFL, that is, the syntax is possible:
 a = conditional(gt(f, 0.0), 1.0, 0.0)

The left and right expressions are required to have the same shape
and free indices.

In the automatic differentiation I issue a warning if the condition
depends on what you differentiate w.r.t., but don't touch the condition
and simply differentiate the left and right expressions separately.

--
Martin