Get rid of threading and updates to form API's

[kinnala]

I propose that we simplify the code and get rid of threading since it's basically never used and doesn't provide too great improvements to execution time even if someone would use it.

[kinnala]

I'll also do some experiments with different form API's discussed in #101 but won't promise anything.

[gdmcbain]

I propose that we simplify the code and get rid of threading since it's basically never used and doesn't provide too great improvements to execution time even if someone would use it.

This would be really good. I find it very difficult to debug forms because pdb can't get into the threads. I was actually thinking of creating a local branch that eschewed threading just for local debugging purposes.

[kinnala]

These latest additions allow supporting quite easily this Fenics-type notation:

from skfem import *

m = MeshTri()
e = ElementTriP1()

b = InteriorBasis(m, e)

from skfem.assembly.form.bilinear_form import BilinearForm

@BilinearForm
def mass(u, v, w, dx):
    return u * v * dx

def inner(u, v):
    return u[0] * v[0] + u[1] * v[1]

def grad(u):
    return u[1]

@BilinearForm
def dudv(u, v, w, dx):
    return inner(grad(u), grad(v)) * dx

A=asm(mass, b)
B=asm(dudv, b)

I don't know if it makes sense or not. I don't have much experience with Fenics. Do you think it would make sense from the viewpoint of getting more users try this library?

[kinnala]

What I'm proposing in particular is to add this kind of "Fenics compatibility layer" so that the user could do

@BilinearForm
def bilinf(u, v, w, dx):
    from skfem.compatibility.fenics import grad, inner
    return inner(grad(u), grad(v)) * dx

[gdmcbain]

Very interesting.

More generally I like the idea of having multiple constructors for an internal representation of forms. (This would mean defining e.g. inner above differently, not immediately carrying out the NumPy operations. That might be a lot of work though for I'm not sure what practical benefit.)

FEniCS/ufl is a substantial effort in how to write these things down. (Other notable ones being FreeFEM and nutils.expression.) I think the UFL library is pure Python. Can any of it be used? It's LGPL. I haven't looked inside it. (nutils is MIT.)

One drawback of UFL (or anything similar) over something more indicial like einsum (or nutils.expression) is that it might be simpler for simple problems but it gets complicated quickly, e.g. the issue in UFL of ‘grad(u) vs. nabla_grad(u)’ under the Navier–Stokes equation in the tutorial, i.e. which is ∇u and which is its transpose? That's always obvious in an indicial notation.

Does anything ever get done with the dx? Could it be hidden and supplied by the decorator? I suppose one would want to keep it in a FEniCS-compatibility layer.

[kinnala]

Now dx is not passed to the form because it is not really used. Currently using @BilinearForm and @LinearForm instead of @bilinear_form and @linear_form use the new set of arguments: u, v, w and v, w.

I want to still experiment with changing w to use also DiscreteField.

[kinnala]

I tried to look at UFL docs here: https://fenics.readthedocs.io/projects/ufl/en/latest/manual/internal_representation.html but as you notice they don't help too much in understanding how to use it somewhere else than in Fenics. I'll try to install it after writing this message and play around with it. I'm not going to immediately add any dependency but just look into if and how it could be incorporated in a compatiblity layer.

@gdmcbain In DiscreteField (https://github.com/kinnala/scikit-fem/pull/317/files#diff-f84cf598dc05a4c1d028aa63625a06b9R7), any suggestions for the names of the attributes? Currently they are f for the pointwise values of the basis, df for the values of the derivative, ddf for the values of the second derivative. I was also thinking x, dx, ddx or something. I don't expect the user to touch these expect in very special circumstances. Currently (in this branch) it would be possible to define weak forms like this:

@BilinearForm
def bilinf(u, v, w):
    return u.df[0] * v.df[0] + u.df[1] + v.df[1]

or even

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

but I guess we would want to use something like

@BilinearForm
def bilinf(u, v, w):
    from skfem.helpers import d, dot
    return dot(d(u), d(v))

to abstract out the attributes and indices, so basic use cases won't need them. Note that it's also possible to use the 'classical' syntax by using @bilinear_form instead of @BilinearForm.

By design of skfem, these attributes can mean a bit different things for different types of Element classes. For example, ElementTriRT0 returns scalar dx which should be interpreted as div(phi) where phi is a vectorial Raviart-Thomas basis function. That's why names such as values, jacobian and hessian won't make sense.

[kinnala]

It seems to me that this FFC (Fenics Form Compiler) is used to turn UFL forms into C-code. I hope we don't need to do anything as excessive.

[kinnala]

Suprisingly enough, I was able to do this:

from skfem import *

m = MeshTri()
e = ElementTriP1()

b = InteriorBasis(m, e)

from ufl import *
element = FiniteElement("CG", triangle, 1)

U = TrialFunction(element)
V = TestFunction(element)

def fenics_form(form):
    
    def new_form(u, v, w):
        def eval_u(x, derivatives=None):
            if derivatives is None:
                return eval_u.c.f
            d, = derivatives
            return eval_u.c.df[d]
        eval_u.c = u
        def eval_v(x, derivatives=None):
            if derivatives is None:
                return eval_v.c.f
            d, = derivatives
            return eval_v.c.df[d]
        eval_v.c = v
        return form(U, V)(None, {U: eval_u, V: eval_v})
    
    return new_form

@BilinearForm
@fenics_form
def mass(u, v):
    return u * v

@BilinearForm
@fenics_form
def stiffness(u, v):
    return dot(grad(u), grad(v))

M=asm(mass, b)
K=asm(stiffness, b)
(M, K)

and it seems to work!

[kinnala]

Basically, when evaluating UFL forms, you just need to provide instructions how to evaluate each symbolic function in the form of a function handle. This probably incurs quite a bit of extra function evaluations but at least we have a proof-of-concept now.

[kinnala]

When porting ex04 to use the new style forms, I realized that having w as Dict[str, DiscreteField] can make the notation too tedious. In particular, I don't like writing np.abs(w['x'].f[0]) < 0.1 and such instead of np.abs(w.x[0]) < 0.1.

[kinnala]

Figured out how to bypass that using __getattr__.

[kinnala]

An overview of major internal changes:

• Should be 99% backwards compatible with the exception of Basis.interpolate return argument which is 3-tuple instead of 2-tuple.
• Introduces a new style form syntax

@BilinearForm
def bilinf(u, v, w):
    return u * v

• This won't hopefully be the final way to use the new form syntax: @bilinear_form will be deprecated only later after all the details related to the new style forms are fixed.
• skfem.models was mostly reimplemented using the new style forms.

[kinnala]

Suprisingly it seems that assembling weak Laplacian got a bit faster:

In [17]: %timeit asm(laplace, ib)                                                                        
3.75 s ± 7.95 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

In [18]: m                                                                                               
Out[18]: Tetrahedral mesh with 30481 vertices and 163840 elements.

In [11]: %timeit asm(laplace, ib)                                                                        
3.18 s ± 8.5 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

In [12]: m                                                                                               
Out[12]: Tetrahedral mesh with 30481 vertices and 163840 elements.

It's not entirely clear to me why this happened but still a positive thing I suppose.

[kinnala]

Old:

In [24]: %timeit asm(linear_elasticity(), ib)                                                            
30.2 s ± 1.5 s per loop (mean ± std. dev. of 7 runs, 1 loop each)

New:

In [17]: %timeit asm(linear_elasticity(), ib)                                                            
19.7 s ± 234 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

[gdmcbain]

Looks good! And a nifty speed-up to boot!