Manipulating forms for adjoint calculation

1

Hi,

I am trying to compute sensitivities using the adjoint method in DOLFIN (without using dolfin-adjoint), similar to this post.

My adjoint equation is as follows:
\frac{\partial{R}}{\partial{U}}\cdot \lambda^T = \frac{\partial{J}}{\partial{U}}, where \frac{\partial{J}}{\partial{U}} = \frac{\partial{J}}{\partial{\psi}}\frac{\partial{\psi}}{\partial{U}}

However I am struggling to assemble the RHS of the adjoint equation (\frac{\partial{J}}{\partial{\psi}}\frac{\partial{\psi}}{\partial{U}}) in a consistent manner to solve to adjoint equation, resulting in a TypeError: unsupported operand type(s) for *: 'Form' and 'Form' error.

Please see adapted MWE below:

import numpy as np
from dolfin import *
from dolfin_adjoint import *

n = 10
mesh = UnitSquareMesh(n, n)


V = VectorFunctionSpace(mesh, "CG", 1)
S = FunctionSpace(mesh, "DG", 0)
u = Function(V)  # test function (displacment field)
v = TestFunction(V)  # Test function
v_ = TrialFunction(V)  # trial function

# kinematics
lambda_ = 0.5769
mu = Function(S)
mu.vector()[:] = 0.3846
P = 2.0*mu*sym(grad(u)) + lambda_ * \
    tr(sym(grad(u)))*Identity(len(u))
dx = Measure("dx", metadata={"quadrature_degree": 4})
R = inner(grad(v), P)*dx

# def left(x): return np.isclose(x[0], 0)
# def right(x): return np.isclose(x[0], 1)

def left(x, on_boundary):
    tol = 1E-14
    return on_boundary and near(x[0], 0, tol)

def right(x, on_boundary):
    tol = 1E-14
    return on_boundary and near(x[0], 1, tol)

bc_L = DirichletBC(V, Constant((0,0)), left)
bc_R = DirichletBC(V, Constant((1,1)), right)
bcs = [bc_L, bc_R]

J = derivative(R, u, v_)

# solve the problem
problem = NonlinearVariationalProblem(R, u, bcs, J=J)
solver = NonlinearVariationalSolver(problem)
solver.solve()

# set an arbitrary functinal
psi = (mu*dot(u, u))
J = (psi**2)*dx

# Derivative of J w.r.t psi
dJ_dpsi = 2*psi*dx

# partial derivative of J w.r.t. mu
dJdmu = assemble(derivative(J, mu))[:]

# partial derivative of R w.r.t. mu
dRdmu = assemble(adjoint(derivative(R, mu))).array()  # partial derivative

lmbda = Function(V)
# solve the adjoint vector
lhs = adjoint(derivative(R, u))
# RHS using chainrule
rhs = -(dJ_dpsi * derivative(psi*dx, u)) # TypeError: unsupported operand type(s) for *: 'Form' and 'Form'

problem = LinearVariationalProblem(lhs, rhs,lmbda, bcs=bcs)
solver = LinearVariationalSolver(problem)
solver.solve()
lmbda_vec = lmbda.vector()[:]
dJdmu += np.dot(dRdmu,lmbda_vec)

Is there way to perform the \frac{\partial{J}}{\partial{\psi}}\frac{\partial{\psi}}{\partial{U}} calculation such that the output (rhs) can be used in the solver?

I have tried assembling the forms and then performing the multiplication (assemble(dJ_dpsi) * assemble(derivative(psi*dx, u)), however this is not compatible with the solver inputs since the rhs is now a vector rather than form.

Thanks!