# Partial derivative of a residual

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 *
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.)
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):

# 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))

# 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
``````

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
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
File "/home/carolina/anaconda3/envs/mdaoproject/lib/python3.8/site-packages/pyadjoint/tape.py", line 46, in wrapper
return function(*args, **kwargs)
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'
``````

Why don’t you simply call `derivative(a, u_sol)`?

No reason besides ignorance. So I did as you suggested and replaced the last three lines with this

``````pR_pd = derivative(a - L, density)
pR_pu = derivative(a - L, u_sol)
``````

And I no longer get an error, which is great. However, I’m not exactly sure what to do with this.

What I’m trying to get is the Jacobian of the residual with respect to `density` and `u_sol` in matrix form. In vector form, the residual `R = A*u_sol - b` should be of size n, the number of degrees of freedom and also the size of `u_sol.vector()` and `b`. The size of `density` is m, the number of elements. I’d expect the partial of `R` wrt `density` to be a matrix of size (n x m) and the partial of `R` wrt `u_sol` to be a matrix of size (n x n).

When I call `derivative(a - L, density)` this gives me a ufl form, which I tried to assemble into a matrix using `assemble(pR_pd)` but then I get this error:

``````Calling FFC just-in-time (JIT) compiler, this may take some time.
Traceback (most recent call last):
File "<stdin>", line 1, in <module>
File "/home/carolina/anaconda3/envs/mdaoproject/lib/python3.8/site-packages/fenics_adjoint/assembly.py", line 17, in assemble
output = backend.assemble(*args, **kwargs)
File "/home/carolina/anaconda3/envs/mdaoproject/lib/python3.8/site-packages/dolfin/fem/assembling.py", line 202, in assemble
tensor = _create_tensor(comm, form, dolfin_form.rank(), backend, tensor)
File "/home/carolina/anaconda3/envs/mdaoproject/lib/python3.8/site-packages/dolfin/fem/assembling.py", line 426, in _create_tensor
raise RuntimeError("Unable to create tensors of rank %d." % rank)
RuntimeError: Unable to create tensors of rank 3.
``````

So I’m not sure how to get the matrices that I want. Do you have any suggestions?

Consider:

``````R = action(a, u_sol) - L
dRdu = assemble(derivative(R, u_sol))
dRdrho = assemble(derivative(R, density))
``````
1 Like

I did this and it worked. Thanks a bunch!