I/O of `ufl` forms

[francesco-ballarin]

It would be great if the I/O of forms was introduced as part of newfl.

In the context of my work in reduced order modelling, in the former implementation in RBniCS I had several routines that took a ufl form, processed it in some way, and gave back another ufl form as output. The preprocessing is so complex that for problems of moderate mesh size it takes more CPU time than the assembly of the resulting vectors/matrices themselves. I have not reimplemented yet this feature in RBniCSx.

The implementation of I/O is surely somewhat complicated, in the sense that you wouldn't want to save for instance the mesh itself, but only the expressions and the integrals. This is why, for instance, simply pickling the form doesn't work. Furthermore, also simply storing the __repr__ of the form does not work either, since a string like Mesh #1 will typically appear in the __repr__ of a form , and when loading the form there is no guarantee that your actual mesh is still the first one.

[jorgensd]

It would be good to have some examples of this complex processing to understand what is taking time as well (maybe there is a way of speeding it up).

[francesco-ballarin]

There are at least three uses cases within reduced order modelling that I have in mind:

1. empirical interpolation method: introduced in this work, it is a method to approximate the assembly of a form to make it amenable to model reduction. Say that you have a form f1(mu) * v * dx + f2 * v * dx, where mu represents the input parameter, f1(mu) a parameter dependent forcing term, f2 a parameter independent forcing term, v a test function, the preprocessing needed at the ufl is roughly the following:
   1. loop over all addends, and determine which of them containing a parameter dependent term. For instance, in the example above f1(mu) * v * dx will be selected for modification since f1 depends on mu, while f2 * v * dx will be left unchanged.
   2. out of f1(mu) * v * dx, extract which part of the expression depend on mu, namely f1(mu).
   3. replace f1(mu) with a suitable approximation, obtained through an application of the proposed method. The approximation will replace f1 with a linear combination of several functions, where by several it may typically happen tens or hundreds functions. This is part of the reason why the preprocessing is expensive in this case.
   4. restrict integration to a subset of the mesh, again specified by the method itself. In particular, this requires replacing measures dx over the entire domain with a measure dx(i), where i is an appropriate marker of a subset of the cells, and regenerate the integral accordingly. This will ensure that the evaluation of the final form is compatible with the requirements of the reduced order model.
2. shape parametrization: this is a case in which the parameter mu affects the geometry of the domain itself. Since it is challenging to integrate mesh motion methods into reduced order models, the standard way of handing parameter dependent domain is to take a reference mesh, a map between the deformed (parameter dependent) domain and the reference one, and use the map to pull back the weak form of the PDE from the deformed domain to the reference one. For efficiency of the resulting reduced model, the map is typically shared among all cells in a subdomain (rather than being defined cell by cell, as in FE). Given the expression of the map in each subdomain and a form on the deformed domain, the preprocessing I used to have in RBniCS is to carry out such a pull back symbolically. Part of the reason why the implementation is currently slow is that I used sympy in some of the calculations: a complete reimplementation using only ufl would probably be a benefit the overall performance.
3. separation of variables: some model reduction methods operate using separation of variables. For instance, if the domain is a unit square and dxy is a volume measure on the square, they aim to approximate integrals f * v * dxy with an approximation given by the combination of the product of fx * vx * dx and fy * vy * dy, where f, fx, fy are functions on the square, unit interval in the x direction, unit interval in the y direction, and similarly for v, vx and vy. I don't have yet an implementation for such methods, but the goal would be to have the user type the f * v * dxy and then have ufl write out the corresponding x and y forms by implementing the appropriate separation of variables employed by the method.

ufl is the ideal setting for carrying out this preprocessing. For instance, one could argue that 2 is done (slightly differently) at the level of the form compiler for the FE assembly, so why do this at the ufl level? The rationale is that 1, 2 and 3 often come together, for instance if the problem has both shape parametrization and parameter dependent terms then it is critical that 2 (first) and 1 (afterwards) need both to be carried out.

This feature request has anyway a potentially broader audience that the specific usage 1, 2, 3 above.