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 ?