External operators

[nbouziani]

This PR follows the discussion initiated here: https://github.com/firedrakeproject/ufl/issues/15, cf. the following PR https://github.com/firedrakeproject/ufl/pull/16

ExternalOperator is a subclass of Coefficient which carries operands and a derivative multi-index. Its purpose is to enable the incorporation in a form of operators that are not directly expressible in UFL. A number of important classes of function fall into this category:

• Implicit constitutive relations. I.e. nonlinear solids or fluids where the stress is computed by solving a point wise nonlinear equation every time the operator is evaluated.

• Explicit functions not easily written in UFL such as neural nets which might be used in learned parametrisations or learned regularisers for inverse problems.

• Parametrisations which are inconvenient to write in UFL because they, for example, contain a lot of switching and interpolation from lookup tables.

ufl.ExternalOperator handles the symbolic side (e.g. differentiation), a practical implementation is needed to state how does one correlate the operands to evaluate the ExternalOperator.

This practical implementation is already implemented in Firedrake: see firedrakeproject/firedrake#1674

The implementation of the adjoint capability is also already done, enabling to tackle more generally PDE-constrained optimization problems.

Examples of applications cover:

• Implicit constitutive relations (e.g. Glen's flow law in glaciology)
• Regularizer for inverse problems (e.g. seismic inversion)
• Subgrid parametrization

Mathematically, an ExternalOperator N can be defined as follows:

V_1 x ... x V_k  --> X
u1, ..., uk    |--> N(u1, ..., uk)

where u1, ..., uk are Coefficients belonging respectively to the function spaces V_1, ..., V_k, and where X is the function space associated to the external operator N.

In practice, when assembling a form containing N, we:

• Interpolate the operands (u1, ..., uk) into X
• Evaluate N(u1^{X}, ..., uk^{X})

Obviously, we need to explicitly define how to evaluate N, this is done at the Firedrake stage where we have an AbstractPointwiseOperator class that subclasses ExternalOperator and that provides a method to evaluate the ExternalOperator.

For instance, for the implicit constitutive law example the evaluate method will solve the nonlinear constitutive law pointwisely.

For the explicit function case, here is toy example illustrating how looks the pseudo code to solve the poisson's equation with a rhs f where inner(u-f, v) is replaced by inner(p2, v):

...
V = FunctionSpace(...)
...
u = Function(V)
...
p = point_expr(lambda x, y: x - y, function_space=V)
p2 = p(u, f)
...
F = inner(grad(p2), grad(v) * dx + inner(p2, v) * dx
solve(F == 0, u, ...)
...

[meg-simula]

I like this concept! Do not have time to look at the code at the moment, what do you think @michalhabera @wence- ?

[dham]

We've had an offline review from @wence- and we think that, while the interface is right, we are wrong to think of this as a terminal. In fact it should be an operator. We think that making this change will make this PR a lot less intrusive in the UFL code by removing a whole load of special casing. @nbouziani will push updates soon.

[michalhabera]

Thanks @dham , I like the concept as well as @meg-simula . I am just missing a pure UFL example code, @nbouziani provided point_expr which is a black box wrt. UFL. I will have a look at the code after announced updates.

[nbouziani]

Thanks @dham , I like the concept as well as @meg-simula . I am just missing a pure UFL example code, @nbouziani provided point_expr which is a black box wrt. UFL. I will have a look at the code after announced updates.

The updates have been pushed

[jhale]

First of all, I am generally positive about this PR, I see the need for this functionality (a generic mapping/operator/function from Expressions to Coefficients, and support for AD).

One thing I'd question is the naming ExternalOperator, particularly the External prefix. Could you explain the reasoning behind External a bit further? I understand that concretely in Firedrake you eventually want to generate a call to External code (outside UFL). But we have use cases where we might want to ufl.substitute certain terms for manually derived UFL Expressions. Possible alternatives for naming might be: AbstractOperator, NotionalOperator, Map. It would be good to keep things in UFL as generic and 'mathematical' as possible.

[dham]

Hi Jack,

It's "external" in the sense that it's external to UFL. It's a foreign function interface. I think that that's a better description than the other suggestions you make, which I think mean different things. I think I don't understand your use case - if you want to express your operator in UFL then you can already do that, right?

[jedbrown]

Is it intended that this capability could be used to implement an exact or Roe Riemann solver for DG? (These involve a Newton solve for a "star" state and an eigendecomposition, respectively.) What if the equation of state is implicit f(pressure; density, energy) = 0?

[jedbrown]

Thanks, could you point me at the Firedrake example you mentioned?

[mscroggs]

Can we close this now that #181 is open?

[nbouziani]

Can we close this now that #181 is open?

Yes but for a different reason. This PR is based on the pre-dual implementation and is now obsolete, this PR should have been closed. As for #181, it is wrong for the same reason. I also don't understand why #181 was introduced.

[jhale]

@nbouziani So where is the post dual PR?