Functional of subfunctions projected onto a different mesh

[David]

Hi all,

I want to estimate parameters of a flow problem, given velocity observations on a coarse mesh. The flow problem needs to be solved on different, finer mesh.

So consider a Stokes problem with a mixed P2/P1 function space. After solving, I want to extract the velocity solution, project it onto the function space of the observations. The functional is the L2-norm of the difference of the projected velocity and the observations.

However, the combination of subfunctions and the projection to the second mesh does not seem to work, see the MWE below. It fails with the error:

Python traceback:
Traceback (most recent call last):
  File "/home/fenics/local/lib/python2.7/site-packages/libadjoint/libadjoint.py", line 1305, in newfunc
    func(*args, **kwargs)
  File "/home/fenics/local/lib/python2.7/site-packages/libadjoint/libadjoint.py", line 1650, in __mat_solve_callback__
    x = A.solve(Variable(var=adj_var), b)
  File "/home/fenics/local/lib/python2.7/site-packages/dolfin_adjoint/adjlinalg.py", line 427, in solve
    x = self.basic_solve(var, b)
  File "/home/fenics/local/lib/python2.7/site-packages/dolfin_adjoint/adjlinalg.py", line 364, in basic_solve
    assembled_rhs = wrap_assemble(b.data, test)
  File "/home/fenics/local/lib/python2.7/site-packages/dolfin_adjoint/adjlinalg.py", line 537, in wrap_assemble
    assert len(form.integrals()) == 0
AssertionError

For a simple non-mixed function space, omitting the subfunction, projecting to the other mesh does not seem to be a problem. Also, if for the Stokes problem, the observation function space is defined on the SAME mesh, no errors.

What can I do? Thanks!

MWE (arbitrary & nonsensical):

from dolfin import *
from dolfin_adjoint import *

mesh = RectangleMesh(Point(0, -1), Point(6, 1), 24, 8)
mesh_obs = RectangleMesh(Point(0, -1), Point(6, 1), 12, 4)

V = VectorElement('P', mesh.ufl_cell(), 2)
P = FiniteElement('P', mesh.ufl_cell(), 1)
W = FunctionSpace(mesh, V*P)

# FAILS
V_obs = VectorFunctionSpace(mesh_obs, 'P', 1)
# WORKS
# V_obs = VectorFunctionSpace(mesh, 'P', 1)

u, p = TrialFunctions(W)
v, q = TestFunctions(W)

f = Constant((1., 0.))
f = interpolate(f, W.sub(0).collapse(), name='control')

a = inner(grad(u), grad(v))*dx + p*div(v)*dx + q*div(u)*dx
L = dot(f, v)*dx

bc = DirichletBC(W.sub(0), Constant((0, 0)), 'on_boundary && (near(x[0], 0) ||'
                 ' near(x[1], -1) || near(x[1], 1))')

# create data
parameters['adjoint']['stop_annotating'] = True
w_ref = Function(W)
solve(a == L, w_ref, bc)
u_ref, p_ref = w_ref.split(True)
u_obs = interpolate(u_ref, V_obs, name='u_obs')

# anntotate
parameters['adjoint']['stop_annotating'] = False
w = Function(W, name='w')
solve(a == L, w, bc)
u, p = split(w)
u_proj = project(u, V_obs, name='u_proj')

success = replay_dolfin(tol=0.0, stop=True)
print('replay: {}'.format(success))
adj_html('tmp/forward.html', 'forward')
adj_html('tmp//adjoint.html', 'adjoint')

J = Functional(inner(u_proj - u_obs, u_proj - u_obs)*dx)
m = Control(f)
RF = ReducedFunctional(J, m)
RF.taylor_test()
# dJdm = compute_gradient(J, m, forget=False)