Actions and adjoints.

[ah3557]

I'm getting an error when I try to apply actions and adjoints. Could you please help me with this?
The followiing is a minimal example and the error.

from firedrake import *
from firedrake.__future__ import interpolate
mesh = UnitSquareMesh(1, 1)
V = FunctionSpace(mesh, "CG", 1)
K = interpolate(TestFunction(V), V)
w = Function(V)
derivative(action(adjoint(K), inner(TestFunction(V), w)*dx), w)

The error:

Traceback (most recent call last):
File "/opt/firedrake/firedrake/src/ufl/ufl/domain.py", line 189, in as_domain
return extract_unique_domain(domain)
File "/opt/firedrake/firedrake/src/ufl/ufl/domain.py", line 255, in extract_unique_domain
domains = extract_domains(expr)
File "/opt/firedrake/firedrake/src/ufl/ufl/domain.py", line 248, in extract_domains
for t in traverse_unique_terminals(expr):
File "/opt/firedrake/firedrake/src/ufl/ufl/corealg/traversal.py", line 137, in traverse_unique_terminals
if op.ufl_is_terminal:
AttributeError: 'Action' object has no attribute 'ufl_is_terminal'

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
File "/opt/firedrake/codes/rfem.py", line 7, in
derivative(action(adjoint(K), inner(TestFunction(V), w)*dx), w)
File "petsc4py/PETSc/Log.pyx", line 188, in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
File "petsc4py/PETSc/Log.pyx", line 189, in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
File "/opt/firedrake/firedrake/src/firedrake/firedrake/ufl_expr.py", line 239, in derivative
mesh = as_domain(form)
File "/opt/firedrake/firedrake/src/ufl/ufl/domain.py", line 191, in as_domain
return domain.ufl_domain()
File "/opt/firedrake/firedrake/src/ufl/ufl/form.py", line 112, in ufl_domain
self._analyze_domains()
File "/opt/firedrake/firedrake/src/ufl/ufl/action.py", line 127, in _analyze_domains
self._domains = join_domains(chain.from_iterable(e.ufl_domain() for e in self.ufl_operands))
File "/opt/firedrake/firedrake/src/ufl/ufl/domain.py", line 205, in join_domains
domains = set(domains) - set((None,))
TypeError: 'MeshGeometry' object is not iterable

Thanks!

[dham]

Something is not quite right in the logic of extracting domains when applied to nested BaseForms. @nbouziani ?

[nbouziani]

I think the logic should be that we collect all the domains we find as we traverse the DAG, and then check that we only get one domain at the end of the traversal, which is what we currently do.

What I don’t understand though is why Action has this special way of chaining the domains of its operands, which assumes that domains are iterable and that is causing the above error. @pbrubeck I can see you introduced this here, am I missing something ?

@ah3557 note that the differentiation of Action with a left slot that is not a 1-form is currently not supported, so the above derivative cannot be handled at the moment. What are you trying to do ? There might be ways to avoid taking this derivative.

[ah3557]

I'm trying to solve a linear system of the form
A = action(adjoint(K), action(S, action(K, w)))

where S is a stiffness matrix, K is an interpolation matrix and w is a function.

I thought I can solve the system without assembling A.

[dham]

I think I don't understand where the solve comes into this (or the derivative, for that matter). Can you show us mathematically what you are trying to do?

[pbrubeck]

You could avoid this issue by simply providing a dummy Jacobian, so derivative is never called. But this issue is more subtle. The commit that @nbouziani refers to calls e.ufl_domains() which returns an iterable, as opposed to e.ufl_domain(), which is deprecated. But if you look at the current the FEniCS/ufl main branch, this commit seems to have been dropped last week.

[dham]

Well if you build a LVP or call solve on a linear system you also don't get a derivative, but I don't actually see where the system is here.

[pbrubeck]

Apart from that, Adjoint needs to implement ufl_domains()

[ah3557]

@dham: I'm trying to apply recovered finite element methods. Attached is the weak form I'm solving. the operator, epsilon, maps DG elements to Nedelec ones.

[dham]

OK, I think if you set up a LinearVariationalProblem then you will avoid this derivative. There is really no need for a derivative to be happening in this case. I think Firedrake will let you solve this problem either matrix-explicit or matrix-free depending on how you set the mat_type.

[pbrubeck]

There's also another bug in how we extract arguments of a complicated linear form

from firedrake import *
from firedrake.__future__ import interpolate
mesh = UnitSquareMesh(1, 1)
V = FunctionSpace(mesh, "CG", 1)
Q = FunctionSpace(mesh, "DG", 1)

# Declare a single argument for Q
q = TestFunction(Q)

K = interpolate(q, V)
S = inner(TrialFunction(V), TestFunction(V)) * dx
z = Function(Q)

# L and R are both linear forms
L = Cofunction(Q.dual())
R = action(adjoint(K), action(S, action(K, z)))
assert len(R.arguments()) == 1 # Passes
assert len(L.arguments()) == 1 # Passes

# L can be combined with a Form
F1 = L - inner(q, z)*dx
assert len(F1.arguments()) == 1 # Passes

# L cannot be combined with a complicated Action
F2 = L - R
assert len(F2.arguments()) == 1 # Fails