FEniCS / ffcx

Next generation FEniCS Form Compiler for finite element forms
https://fenicsproject.org
Other
153 stars 41 forks source link

Factorization failure for certain conditionals #246

Open mrambausek opened 4 years ago

mrambausek commented 4 years ago

While probably quite a corner case, I managed to stumble into that...

import ufl
from ffcx.codegeneration import jit

elmt = ufl.FiniteElement("CG", ufl.triangle, 1)

u1 = ufl.Coefficient(elmt)
u2 = ufl.Coefficient(elmt)
ctol = ufl.Constant(ufl.triangle)

def cond(tol):
    return ufl.conditional(u2 < tol, u1, u2) * ufl.dx

# works
jit.compile_forms([ufl.derivative(cond(0), (u1, u2))])

# works
jit.compile_forms([ufl.derivative(cond(ctol), u1)])

# works
jit.compile_forms([ufl.derivative(cond(ctol), u2)])

# fails
jit.compile_forms([ufl.derivative(cond(ctol), (u1, u2))])
michalhabera commented 4 years ago

@mrambausek Can you include also an error message? Thanks!

mrambausek commented 4 years ago
~/opt/fenics/2019.1.0-fedora31-workstation/fenics-py/lib64/python3.7/site-packages/ffcx/ir/uflacs/analysis/factorization.py in compute_argument_factorization(S, rank)
    340             for o in expr.ufl_operands:
--> 341                 F.add_edge(i, F.e2i[o])
    342 

KeyError: Zero((), (), ())
mrambausek commented 4 years ago

this fails as well:

jit.compile_forms([ufl.derivative(cond(1e-6), (u1, u2))]) 

which is quite surprising since the only change is from a float "0" to a float "1e-6"

mrambausek commented 4 years ago
~/opt/fenics/2019.1.0-fedora31-workstation/fenics-py/lib64/python3.7/site-packages/ffcx/ir/uflacs/analysis/factorization.py in compute_argument_factorization(S, rank)
    340             for o in expr.ufl_operands:
--> 341                 F.add_edge(i, F.e2i[o])
    342 

KeyError: Zero((), (), ())

What seems imporant, is that at least one of the derivatives wrt. (u1, u2) has to vanish. This conditional, for example, works:

def cond(tol):
    return ufl.conditional(u2 < tol, u1*u2, 2*u1*u2) * ufl.dx
mrambausek commented 4 years ago

At least for the MW(F)E, a possible workaround seems to be

v1, v2 = ufl.TestFunctions(u1.ufl_element() * u2.ufl_element())
jit.compile_forms([ufl.derivative(cond(ctol), u2, v2) + ufl.derivative(cond(ctol), u1, v1)])
michalhabera commented 4 years ago

Just note that taking a (!distributional) derivative of a conditional is in most use cases mathematical non-sense, unless you add an extra delta measure which comes from the differentiation of a discontinuous function. What UFL and FFC do is that they differentiate piecewise, which is not always correct.

But I'll have a better look at why this issue appeared. I rewrote some factorization parts a while ago when allowing tensor-valued expressions to be factorized. So the issue originated probably at that point.

mrambausek commented 4 years ago

Thanks a lot for looking into that -- if you find the time.

I would say, a valid use case for a distributional derivative of a conditional is when you have an energy density of which the derivative exists and is smooth, but has an unstable numerical evaluation (e.g., one needs to take limits). Antiderivatives of sigmoidal functions are candidates for such potentials. As a remedy, one can -- by means of a conditional -- replace the original energy density by a quadratic approximation around certain points. Now, in case of coupled problems, e.g., magneto-mechanics, problems as described here could appear in certain situations. Of course, this is still a corner case, given the special case that triggered the problem. I am happy with my current workaround, but I would highly appreciate a fix or at least a warning if it really is an unsupported operation.