Dolfin-adjoint: Compute adjoint of form, form is not bilinear

I’ve posted a topic but get no reply, so I want to re-post it again …
I’m come across the error “Compute adjoint of form, form is not bilinear” using dolfin-adjoint.
When I run without optimization, solve(F==0) shows no error.
When optimization or sensitivity analysis codes added, the error comes.
Here’s the mwe (the code is simplified so control actually does not influence object function J). This is a 2D scattering problem and ABC is absorbing boundary condition. I removed obstacle to create for simple mesh, and the result will be a plane wave.

import mpi4py.MPI
from fenics import *
from fenics_adjoint import *
import numpy as np
import moola
# pip install git+https://github.com/funsim/moola.git@master

# k = 2*np.pi/0.5 # 2pi/lambda
k0 = 2*np.pi/0.5
E0 = 1.  # incident wave magnitude

mesh = UnitSquareMesh(20,20)

class ABC(SubDomain):
    def inside(self,x,on_boundary):
        return on_boundary

sub_domains = MeshFunction("size_t",mesh,mesh.topology().dim()-1)
sub_domains.set_all(2)
abc=ABC()
abc.mark(sub_domains,1)
print('1')
ele_type = FiniteElement('Lagrange',triangle,2)
V2 = FunctionSpace(mesh, ele_type*ele_type)
V = FunctionSpace(mesh, ele_type)
uh = Function(V2, name='State')
(u_r, u_i) = split(uh)
(v_r, v_i) = TrialFunction(V2)
ds = Measure('ds', domain=mesh, subdomain_data=sub_domains)
u0 = Constant(0.0)

k = Expression('k0', degree=2, name='Control', k0=k0)
k = interpolate(k, V)
inc_r = Expression('E0*cos(k0*x[0])', degree=2, k0=k0, E0=E0)
inc_r = interpolate(inc_r, V)
inc_i = Expression('E0*sin(k0*x[0])', degree=2, k0=k0, E0=E0)
inc_i = interpolate(inc_i, V)

n = FacetNormal(mesh)
F = (inner(grad(u_r),grad(v_r))-inner(grad(u_i),grad(v_i)) - (k**2*u_r*v_r-k**2*u_i*v_i))*dx + \
    (-k0*u_r*v_i-k0*u_i*v_r + k0*inc_r*v_i+k0*inc_i*v_r - v_r*inner(n,grad(inc_r))+v_i*inner(n,grad(inc_i)))*ds(1)

solve(F == 0, uh, solver_parameters={"newton_solver": {"absolute_tolerance": 1.0e-7,
                                                              "maximum_iterations": 20}})
print('2')

ur, ui = uh.split(True)
File("./results/uh_r.pvd").write(ur)

x = SpatialCoordinate(mesh)
J = assemble(inner(u_r-inc_r, u_r-inc_r)*dx)
control = Control(k)
dJdc = compute_gradient(J,control, options={'riesz_representation':'L2'})import mpi4py.MPI
from fenics import *
from fenics_adjoint import *
import numpy as np
import moola
# pip install git+https://github.com/funsim/moola.git@master

# k = 2*np.pi/0.5 # 2pi/lambda
k0 = 2*np.pi/0.5
E0 = 1.  # incident wave magnitude

mesh = UnitSquareMesh(20,20)

class ABC(SubDomain):
    def inside(self,x,on_boundary):
        return on_boundary

sub_domains = MeshFunction("size_t",mesh,mesh.topology().dim()-1)
sub_domains.set_all(2)
abc=ABC()
abc.mark(sub_domains,1)
print('1')
ele_type = FiniteElement('Lagrange',triangle,2)
V2 = FunctionSpace(mesh, ele_type*ele_type)
V = FunctionSpace(mesh, ele_type)
uh = Function(V2, name='State')
(u_r, u_i) = split(uh)
(v_r, v_i) = TrialFunction(V2)
ds = Measure('ds', domain=mesh, subdomain_data=sub_domains)
u0 = Constant(0.0)

k = Expression('k0', degree=2, name='Control', k0=k0)
k = interpolate(k, V)
inc_r = Expression('E0*cos(k0*x[0])', degree=2, k0=k0, E0=E0)
inc_r = interpolate(inc_r, V)
inc_i = Expression('E0*sin(k0*x[0])', degree=2, k0=k0, E0=E0)
inc_i = interpolate(inc_i, V)

n = FacetNormal(mesh)
F = (inner(grad(u_r),grad(v_r))-inner(grad(u_i),grad(v_i)) - (k**2*u_r*v_r-k**2*u_i*v_i))*dx + \
    (-k0*u_r*v_i-k0*u_i*v_r + k0*inc_r*v_i+k0*inc_i*v_r - v_r*inner(n,grad(inc_r))+v_i*inner(n,grad(inc_i)))*ds(1)

solve(F == 0, uh, solver_parameters={"newton_solver": {"absolute_tolerance": 1.0e-7,
                                                              "maximum_iterations": 20}})
print('2')

ur, ui = uh.split(True)
File("./results/uh_r.pvd").write(ur)

x = SpatialCoordinate(mesh)
J = assemble(inner(u_r-inc_r, u_r-inc_r)*dx)
control = Control(k)
dJdc = compute_gradient(J,control, options={'riesz_representation':'L2'})

and full error info

Traceback (most recent call last):
  File "main_adjoint_mwe.py", line 59, in <module>
    dJdc = compute_gradient(J,control, options={'riesz_representation':'L2'})
  File "/usr/local/lib/python3.6/dist-packages/pyadjoint/drivers.py", line 29, in compute_gradient
    tape.evaluate_adj(markings=True)
  File "/usr/local/lib/python3.6/dist-packages/pyadjoint/tape.py", line 140, in evaluate_adj
    self._blocks[i].evaluate_adj(markings=markings)
  File "/usr/local/lib/python3.6/dist-packages/pyadjoint/tape.py", line 46, in wrapper
    return function(*args, **kwargs)
  File "/usr/local/lib/python3.6/dist-packages/pyadjoint/block.py", line 127, in evaluate_adj
    prepared = self.prepare_evaluate_adj(inputs, adj_inputs, relevant_dependencies)
  File "/usr/local/lib/python3.6/dist-packages/dolfin_adjoint_common/blocks/solving.py", line 142, in prepare_evaluate_adj
    dFdu_form = self.backend.adjoint(dFdu)
  File "/usr/local/lib/python3.6/dist-packages/dolfin/fem/formmanipulations.py", line 44, in adjoint
    raise RuntimeError("Compute adjoint of form, form is not bilinear")
RuntimeError: Compute adjoint of form, form is not bilinear

I cannot find similiar errors on the Internet. I guess this error comes from my V2 FunctionSpace ?

Problem solve! V.collapse() shall be added for subspaces.