Partial derivative of a residual

dolfin-adjoint
1

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.

1 Like
2

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

3

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?

4

Consider:

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

I did this and it worked. Thanks a bunch!