Renamed: Better form arguments and docstrings.

[kinnala]

When defining bilinear forms you always write

@bilinear_form
def bilin(u, du, v, dv, w):
    pass

namely, the function arguments are in the same order always.

It would be convenient to have a snippet like this and other examples available in the docstring(s) of bilinear_form (and linear_form) so it can be quickly copy-pasted from REPL (help(bilinear_form)) without remembering the order of arguments and what w contains.

[gdmcbain]

the function arguments are in the same order always

Nearly, but only in problems of order two; consider docs/examples/creeping.py

https://github.com/kinnala/scikit-fem/blob/35ec4ca7aa891216ad8b948bcbde77183c9dd7cf/docs/examples/creeping.py#L20-L21

This variability necessitates the recourse to inspect.signature in the definition of bilinear_form

https://github.com/kinnala/scikit-fem/blob/35ec4ca7aa891216ad8b948bcbde77183c9dd7cf/skfem/assembly/asm.py#L169

Feeling that this was a bit eldritch, I thought about how one might craft an order-independent signature. What about grouping the derivatives of the trial and test functions?

def bilin([u, ...]: List[ndarray], [v, ...]: List[ndarray], w: Any):  # pseudocode

Note that the length of a List isn't part of its type. Then inside the body one would access the ordinates with u[0] and for k = 1, 2, … the k-th derivatives with u[k].

Against this suggestion: it's backwards-incompatible and it increases the proliferation of indices.

[gdmcbain]

Also on docstrings and bilinear and linear forms, it's unfortunate that at present there's no point giving a particular bilinear form a docstring since it gets lost in the decorating. I think in particular of those in skfem.models.

It might be possible to get around this using functools.wraps.

[kinnala]

I was today playing around with the extended keyword argument/dict unpacking syntax explained here:

https://docs.python.org/3/whatsnew/3.5.html#pep-448-additional-unpacking-generalizations

def bilin(du, v, **w):
    return du[0]*v

then this could be called like

bilin(**{'u' : u, 'du' : du, 'v' : v, 'dv' : dv, ...})

so the unused parameters would always go to w.

What do you think?

Pros: Quicker to define as you don't need write everything, and can probably be implemented nicely as everything is in a dict defined in GlobalBasis object.

Cons: Can feel a bit like magic (like "what options do I have here?").

[kinnala]

I've been also wondering whether one should be able to call bilinear_form decorator like

basis = InteriorBasis(m, e)

@bilinear_form(basis)
def bilinf(du, dv, **w):
    return sum(du, dv)

so that it is more evident which basis it uses.

Then you could assemble simply like:

A = asm(bilinf)

or even

A = bilinf.asm()

and if parameter fields are needed,

A = bilinf.asm(w)

[gdmcbain]

Interesting. Some thoughts:

If the decorators for forms took the basis as an argument, how would one write skfem.models?

What about making asm a method of the basis?

class GlobalBasis:
    …
    def asm(self, 
                  kernel, 
                  vbasis: Optional[GlobalBasis] = None,
                  …)

then

A = basis.asm(bilinf)

It is of course possible to maintain both interfaces (basis.asm(bilinf, …) and asm(basis, bilinf, …)), as is done for many array methods in NumPy.

For the forms themselves, I think conceptually bilinear forms take two arguments (the unknown and the trial function) plus maybe some other stuff (w). And linear forms one (the trial function), again maybe plus some other stuff. Should this be reflected in the interface?

Another thought is that conceptually bilinear and linear forms are not intrinsically defined in terms of bases. If one begins with a partial differential equation and proceeds to the variational form by taking the inner product with a trial function, it's not immediately necessary to specify what the inner product space is, still let what basis is chosen. Only later, before assembly, does the basis have to be specified.

[kinnala]

I agree with your ideas on changing asm.

Feeling that this was a bit eldritch, I thought about how one might craft an order-independent signature. What about grouping the derivatives of the trial and test functions?

def bilin([u, ...]: List[ndarray], [v, ...]: List[ndarray], w: Any):  # pseudocode

Note that the length of a List isn't part of its type. Then inside the body one would access the ordinates with u[0] and for k = 1, 2, … the k-th derivatives with u[k].

Against this suggestion: it's backwards-incompatible and it increases the proliferation of indices.

I tried this out but used tuples instead of lists because it was more natural with respect to the current implementation of Element.gbasis.

It did simplify things quite a lot. I'm not however confident about suggesting people to write stuff like

@bilinear_form
def laplace(u, v, w):
    return u[1][0] * v[1][0] + u[1][1] * v[1][1]

It can look quite cumbersome when the expressions are long.

I'm now thinking about implementing some helpers to make these look nicer. Any ideas are welcome.

[gdmcbain]

Ha indeed, u[1][0] * v[1][0] + u[1][1] * v[1][1] is a proliferation of indices. It is logical and expressive but not clear. Does sum(u[1] * v[1]) work here? But yes in general some thought is required here; this is one of the core parts of the interface.

This is why FEniCS invented UFL and nutils invented their expressions. The latter found that without these, they had the same problem.

Though powerful, the resulting code is often lengthy, littered with colons and brackets, and hard to read.

An earlier (pre-Python) and very successful system is FreeFem++. There one starts with something like dx(u)*dx(v) + dy(u)*dy(v) but often lightens it by writing 'macros' (which would just be functions in Python), ending up with, e.g. scalar(grad(u), grad(v)). I guess already in Python one, if one wanted to get to the first of those ine could define little helpers.

def dx(u): u[1][0]
def dy(u): u[1][1]

I assume there are other existing systems around that would be worth studying too (sfepy, Feel++, escript-finley, Elmer, &c., &c.). Yeah, careful thought here is worthwhile, I think.

[gdmcbain]

def grad (u): return u[1]

is easy but I don't see anything equivalent for u[0] yet.

I expect that this kind of thing gets really tested in elasticity where second and higher order tensors can't be avoided.

[kinnala]

Random idea: we could also directly modify AST in the decorators to transform u into u[0] using NodeTransformer (and pretty much arbitrary substitutions).

[gdmcbain]

Oh! That's really interesting! I've never looked into the AST module but have heard of it yes. And it is in the Python Standard Library after all! Very interesting!

The nutils expression DSL is quite nice but I don't like so much how it relies on strings. I wonder whether AST would enable similar transformations without strings.

Another technique from the standard library that, used judiciously, might be handy here is operator overloading. Python 3's @-operator is a fantastic addition to the language. Whenever I've had that thought before though I started to worry that I didn't have the prerequisite knowledge of exterior algebra, geometric algebra, &c. I've long suspected that while vector notation, dyadics, tensor-indicial notation, &c., have been hugely beneficial in writing mathematics for continuum mechanics, they're unlikely to have been the last word.

[kinnala]

This was done a while ago, new forms will be released in 1.0.0.