Support linear combinations of fields (coefficients) in UFL

[kynan]

Currently MCFC only supports solving a linear system of equations (corresponding to a bilinear and linear form) for a coefficient and writing this coefficient back to Fluidity state. There is an implicit assumption that a solve is present in the UFL equation for the generated code to be meaningful at this point.

Other arithmetic expressions (linear combinations of fields) that can be expressed in UFL are not supported i.e.

f = state.scalar_fields["Tracer"]

# Blow up the coefficient
f2 = 2*f

state.scalar_fields["Tracer"]=f2

Relates to issue #4

[kynan]

Thinking about that again maybe it's not all that useful/necessary. DOLFIN doesn't support it either (in this way).

DOLFIN supports that through the Function class. Linear combinations of Functions are supported

[gmarkall]

I just discovered a case where you want to be able to write a Sum back to state:

c1 = 0.5
c2 = 0.5

t1 = solve(...)
t2 = solve(...)

ta = t + c1*t1 + c2*t2

as if you were implementing Runge-Kutta timestepping. I imagine as in your example, you also want to be able to return a Product into state.

[gmarkall]

Documenting a discussion with Florian:

We have establised that one way to tackle this would be to fake a coefficient when we want to return a sum or product into state.

Subclassing Coefficient and redefining __add__, __radd__, __mul__ etc might break the way it works in regular forms. We could avoid it by ensuring that the result of the linear combination is computed in time for the form kernel to use the result.

Alternatively, we could get the UFL state to fake the coefficient inside the __setitem__ method. However, this has the disadvantage that it wouldn't work if the linear combination were to be an argument of action(). The other advantages/disadvantages of this method are not clear.

[kynan]

At this point it's not clear to me how this fake coefficient representing a sum / product / linear combination could look like and how it could possibly work with UFL / fit in the UFL design.

Were you thinking about adding constructors to the Coefficient subclass taking a Sum or a Product and storing it in some attribute of the class? Or even evaluating it straight away? If so, what would evaluating mean?

I have to admit I'm rather vague still about what happens on reconstruction of a form and its components (which is afaik what is done when you call preprocess on it).

Also, this linear combination will need to be translated into a (device) kernel to actually yield a field, correct? Where is the logical place for this to happen? At the point where this resulting field is used I suppose?

So far I'm completely unclear what the more straightforward implementation would be. I'll try to come up with some pros and cons of either method:

1. Subclassing Coefficient
   • (+) transparent to use for UFL operators currently accepting a Coefficient
   • (-) need to overload __add__, __radd__, __mul__ etc. to return a Coefficient instead of a Sum, Product etc.
   • (-) for this to work probably need to add a constructor accepting a ufl.AlgebraOperator and add an attribute to store it (or its evaluation?)
   • (-) needs special case treatment when used in a form
   • (-) not clear when the evaluation (i.e. a call to the kernel evaluating the expression and returning a field) is to be done:
     • at Coefficient construction: this is likely to be horribly inefficient for nested expressions (i.e. we only want to do it at the outermost operator, but that is non-trivial to detect
     • at Coefficient assignment to state or action: then we're back to special case treatment
2. allowing ufl.AlgebraOperator to be passed back into state / action
   • (+) no special case treatment needed for UFL forms, no overloading of Coefficient
   • (+) evaluation (i.e. a call to the kernel evaluating the expression and returning a field) of the entire expression done at once on assignment
   • (-) possible duplication of evaluation if expression is assigned multiple times (could be avoided by a cache or sth. similar)

At this point I would lean towards option 2 since it seems much simpler. In fact, the evaluation could be done by a dedicated function of (functor) class which would be need to called by the entity accepting the expression (i.e. state.__setitem__ and action). This evaluation facility could also take care of caching previous evaluations and performing sanity checks i.e. making sure it's actually a linear combination. This would also allow fine-grained control over the places where we want to accept expressions rather than "plain" coefficients (although we probably want to accept expressions everywhere and would rather have the evaluation happen automatically).

[dham]

I like 2. I really don't want to subclass UFL objects unnecessarily. I also strongly suspect that multiple evaluation is a corner case which we should ignore until it becomes a big issue.

[kynan]

I think this issue gives rise to a much more general question / design decision we need to make:

• Do we want to stay really "close to the metal" (being UFL) and add as little to the tool set as possible? This will probably mean more analysis of the UFL code AST.
• Do we want to have the python interpreter do the work for us? This will probably mean adding more and more functionality / convenience wrappers to the MCFC namespace (where the UFL is executed in) with the danger of "overengineering" / re-inventing much of the DOLFIN functionality.

I looked more into the way DOLFIN deals with that and it clearly follows 2. The clever thing they do is deriving the PyDOLFIN Function class from ufl.Coefficient and cpp.Function (SWIG-wrapped). Therefore it appears to UFL as a Coefficient and to (SWIG) DOLFIN as a Function. Also the Python __repr__ is that of a ufl.Coefficient, so it is indistinguishable from the equivalent pure UFL form. I couldn't figure out yet how they deal with linear algebra operations in the end.

[kynan]

Regarding the implementation of 2.:

Whenever a linear combination (LC) of fields/coefficients (as opposed to a "plain" field/coefficient) in form of a ufl.AlgebraOperator is passed to state or solve, the following needs to happen:

1. check validity of the linear combination i.e. check it is actually a linear combination and check that element types match
   • Do they have to be the exact same type? Or are there "compatible" element types that can be summed?
   • Are we restricted to actual linear combinations? Could there ever be a use case for e.g. (point-wise) multiplication of fields?
   • Could be implemented with a tree visitor (ufl.Transfomer?)
2. generate a kernel evaluating the linear combination point-wise
   • gets all the fields as parameters (i.e. one for each coefficient in the LC)
   • constants are hard coded, since we know them at compile time
   • do we support float variables as multiplicative factors?
   • this is not a million miles from evaluating a form -> can we reuse code?
3. generate kernel call with correct parameters and pass result into state/solve
   • for the CUDA backend this will probably mean adding a temporary field

As David pointed out, we would initially ignore the problem of multiple evaluations of this kernel until we actually face it (and it affects performance).

Caveat: if you write

ta = t + c1*t1 + c2*t2
state.scalar_field["Tracer"] = ta

there will be an entry of type ufl.AlgebraOperator in the dictionary of UFL objects, whereas if you write it back to state immediately, i.e.

state.scalar_field["Tracer"] = t + c1*t1 + c2*t2

there won't. To make both of these inputs perfectly equal (as they should) the implementation should not depend on the presence of ta in the dictionary of UFL objects

[dham]

Do they have to be the exact same type? Or are there "compatible" element types that can be summed?

There are compatible types but I think it's a much less important case than the matching one (the matching one turns up in timestepping schemes all the time). Let's just support exact matches until someone complains.

Are we restricted to actual linear combinations? Could there ever be a use case for e.g. (point-wise) multiplication of fields?

Point-wise multiplication is mathematically very dodgy. Function spaces are vector spaces so the only closed operations on them which are guaranteed to be legal are linear combinations.

do we support float variables as multiplicative factors?

I think we need to. We need to support dt as a factor, for example, and dt in general changes during execution.

[kynan]

It may be more straightforward (both implementation-wise and for the user) to actually introduce an evaluation operator (like evaluate) that evaluates a linear combination of coefficients and returns a new coefficient.

• Eliminates any ambiguity (also for the user) where the expression is evaluated
• Eliminates the need for any special treatment in solve/state/action
• Helps avoiding a problem you run into for explicit time stepping with PyDOLFIN: forms containing a linear combination of coefficients from a previous time step successively "blow up" because this linear combination is never evaluated, hence making the form more and more complex with each time step.

[gmarkall]

As I understand it, isn't forcing the user to make a decision about what gets evaluated when exactly the thing we're trying to avoid?

Or is the suggestion that the user should have to give the compiler a helping hand with the analysis, by saying something like

a = evaluate(linearCombination)

so that the uflObjects dict gets forced to be filled with a new coefficient that the evaluate function returned? (if we are going to do this, I think evaluate would be a bad name - maybe something that sounds more like baking or freezing or crystallising - somehow solidifying the expression, but avoiding alluding to actually computing it at that location)

[kynan]

Yes, I think what you're describing is exactly what we want to do. The problem in PyDOLFIN is that such a liner combination is never evaluated, i.e. code like

# Time-stepping
while time < T:

    # Solve for T1-T4.
    rhs = dt*(dot(grad(q),u)*t-q*t*div(u))*dx
    t1 = VariationalProblem(m, rhs).solve()

    rhs = dt*(dot(grad(q),u)*(t+0.5*t1)-q*(t+0.5*t1)*div(u))*dx
    t2 = VariationalProblem(m, rhs).solve()

    rhs = dt*(dot(grad(q),u)*(t+0.5*t2)-q*(t+0.5*t2)*div(u))*dx
    t3 = VariationalProblem(m, rhs).solve()

    rhs = dt*(dot(grad(q),u)*(t+t3)-q*(t+t3)*div(u))*dx
    t4 = VariationalProblem(m, rhs).solve()

    t = t + 1.0/6.0*t1 * 1.0/3.0*t2 * 1.0/3.0*t3 * 1.0/6.0*t4

    time += dt

blows up t in the right-hand-side forms with every time step and JIT-compiling the forms takes longer and longer etc.

I agree that evaluate is a poor choice of name, but baking or freezing is maybe a bit too casual.

But I'm not so sure what other options we'd have than evaluating the linear operator before it gets passed to state/solve. Can you think of some?

[gmarkall]

In dolfin you would get around this by writing

t.assign(t.value() + 1.0/1.0*t1.value() + ... )

However, that particular problem shouldn't affect us anyway since we only compile the code for one timestep, as we have no iteration within the UFL.

I still feel that we should be able to get the effect of evaluate without actually having to force people to write it, but I've not had chance to think about it enough.

[kynan]

Thanks, I wasn't aware how to get this effect in PyDOLFIN.

David pointed out we might well have explicit (time stepping) loops in UFL equations for e.g. subcycling or nonlinear iterations.

I'm looking forward to any suggestions on how to avoid an explicit call to something like evaluate.

[dham]

I think it's going to be hard to avoid an explicit evaluate as it's very hard to diagnose why an expression is being put together. For example:

ta = t + c1*t1 + c2*t2

Now is ta a coefficient or an expression? If we now type:

form=ta*dx

then it clearly has to be an expression.

OTOH if ta was a Runge-Kutta sum then we get the explosion previously noted.

It's not clear to me that the system would ever be in a position to work out which of those two cases it was in.

Regards,

David

[kynan]

I don't see any option to achieve this either in python, since there's no way to overload assignment.

[gmarkall]

OK, so is the conclusion that we need to use something like evaluate?

[dham]

OK, so is the conclusion that we need to use something like evaluate?

I think so.