Hello,

I am trying to get the partial derivative of the residual of my linear elasticity PDE with respect to a scalar function (`density`

) and with respect to the solution (`u_sol`

).

I am able to get the derivative of a different function with respect to density, and I tried to follow that same format for getting the partials of the residual, but I am not sure that I’m taking the right approach since I want partial derivatives and not the total derivative. This is what I have

```
from fenics import *
from fenics_adjoint import *
import numpy as np
from ufl import nabla_div
# Define parameters
nelx = 30
nely = 15
nelz = 15
lx = float(nelx)
ly = float(nely)
lz = float(nelz)
E = 1.04
v = 0.3
mu = Constant(E / (2*(1+v)))
lmbda = Constant(E*v / ((1+v)*(1-2*v)))
# Create a mesh
mesh = BoxMesh.create(
[Point(0.0, 0.0, 0.0), Point(lx, ly, lz)], # define opposing corners
[nelx, nely, nelz], # number of elements in each direction
CellType.Type.hexahedron)
mesh = Mesh(mesh)
# FE function space and functions
V = VectorFunctionSpace(mesh, "CG", 1)
u = TrialFunction(V)
v = TestFunction(V)
u_sol = Function(V) # function for storing the solution
# Density field
D = FunctionSpace(mesh, "DG", 0)
density = Function(D)
density.vector()[:] = np.ones(nelx*nely*nelz)
# Bounday condition
class Left(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 0.)
left = Left()
bc = DirichletBC(V, Constant((0., 0., 0.)), left)
# Stress and strain functions
def sigma(u):
n = len(u) # size of u
return lmbda*nabla_div(u)*Identity(n) + 2*mu*epsilon(u)
def epsilon(u):
return 0.5*(nabla_grad(u) + nabla_grad(u).T)
# Define forms
f = Constant((0., 0., -1.))
L = dot(f, v) * dx
a = inner(density*sigma(u), nabla_grad(v)) * dx
# Solve
A, b = assemble_system(a, L, bc)
solver = KrylovSolver("cg", "ilu")
solver.solve(A, u_sol.vector(), b)
# Define control variables
control_d = Control(density)
control_u = Control(u_sol)
# Calculate gradient of compliance; this works
J = assemble(action(L, u_sol))
dJ_dd = compute_gradient(J, control_d)
# Calculate the partials of the residual; this does not work
R = assemble(action(a, u_sol) - L)
pR_pd = compute_gradient(R, control_d) # <---- error occurs here
pR_pu = compute_gradient(R, control_u)
```

And this is the error I’m getting

```
WARNING:root:Adjoint value is None, is the functional independent of the control variable?
Traceback (most recent call last):
File "mwe_partials.py", line 70, in <module>
pR_pd = compute_gradient(R, control_d) # <---- error occurs here
File "/home/carolina/anaconda3/envs/mdaoproject/lib/python3.8/site-packages/pyadjoint/drivers.py", line 31, in compute_gradient
grads = [i.get_derivative(options=options) for i in m]
File "/home/carolina/anaconda3/envs/mdaoproject/lib/python3.8/site-packages/pyadjoint/drivers.py", line 31, in <listcomp>
grads = [i.get_derivative(options=options) for i in m]
File "/home/carolina/anaconda3/envs/mdaoproject/lib/python3.8/site-packages/pyadjoint/control.py", line 51, in get_derivative
return self.control._ad_convert_type(0., options=options)
File "/home/carolina/anaconda3/envs/mdaoproject/lib/python3.8/site-packages/pyadjoint/tape.py", line 46, in wrapper
return function(*args, **kwargs)
File "/home/carolina/anaconda3/envs/mdaoproject/lib/python3.8/site-packages/fenics_adjoint/types/function.py", line 125, in _ad_convert_type
return compat.function_from_vector(self.function_space(), value, cls=Function)
File "/home/carolina/anaconda3/envs/mdaoproject/lib/python3.8/site-packages/fenics_adjoint/types/compat.py", line 206, in function_from_vector
vector = vector._cpp_object
AttributeError: 'float' object has no attribute '_cpp_object'
```

Any advice would be appreciated. Thanks in advance.