Compute gradient error during adjoint assembly

[dthillaithevan]

Hi, When using compute_gradient to determine the gradient of a functional w.r.t. a control variable in a periodic linear elasticity problem that uses a mixed functionspace, I am getting a IndexError: tuple index out of range error. I've included a MWE below which replicates this issue.

**Edit I dont to have permission to assign so mentioning instead, @dham

import ufl
from firedrake import *
from firedrake_adjoint import *# Define stress/strain functions
def eps(v):
    return sym(grad(v))
def sigma(v, Eps, lmbda, mu):
    return lmbda*tr(eps(v) + Eps)*Identity(2) + 2*mu*(eps(v)+Eps)# Create periodic mesh & mixed functionspace
n = 100
mesh = PeriodicSquareMesh(n,n,1, direction='both', quadrilateral = True)
dx = Measure('dx', domain = mesh)
V = VectorFunctionSpace(mesh, 'CG', 1)
R = FunctionSpace(mesh, 'R', 0)
W = V*R*R
elem = ufl.VectorElement("DQ", mesh.ufl_cell(), 1, variant="equispaced")
DG = FunctionSpace(mesh, elem)
U = FunctionSpace(mesh, 'CG', 1)# Define & assign material distribution (circular inclusion)
E_material = 1
E_void = 1e-06
nu = 0.3x, y = SpatialCoordinate(mesh)
radius = 0.1
tol = 1e-5
material_distribution = Function(U)
material_distribution.interpolate(conditional((x-0.5)**2 + (y-0.5)**2 - radius**2 <= tol, E_material, E_void))
# Lame's parameters
lmbda_function = (material_distribution)*nu/ ((1+nu)*(1-2*nu))
mu_function = (material_distribution) / (2*(1+nu))# Create test & trial functions
dv, dlamb_0, dlamb_1 = TrialFunctions(W)
v_, lamb_0, lamb_1 = TestFunctions(W)
lamb_ = as_vector([lamb_0, lamb_1])
dlamb = as_vector([dlamb_0, dlamb_1])w = Function(W)# Assign macro strain
strain = [-0.01,0,0]
exx, eyy, exy = strain
E_xx = Constant(exx)
E_xy = Constant(exy)
E_yx = E_xy
E_yy = Constant(eyy) 
Eps = as_matrix([[E_xx, E_xy],
            [E_yx, E_yy]])# Setup periodic problem (find microscopic fluctuation) & solve
F = inner(sigma(dv, Eps, lmbda_function, mu_function), eps(v_))*dx()
a, L = lhs(F), rhs(F)
a += dot(lamb_,dv)*dx + dot(dlamb,v_)*dx
solve(a == L, w, [])# Extract solution
v, lamb0, lamb1 = split(w)# Functional for testing purposes
ww = assemble(inner(sigma(v, Eps, lmbda_function, mu_function),grad(v))*dx)# Compute influence of exx term on ww
# dW = compute_gradient(ww, Control(Eps[0,0])) # throws error
# dW = compute_gradient(ww, Control(E_xx)) # throws error

Thanks

[stephankramer]

Just had a short look at this. It appears to be the gradient of the final assembly with regards to the w that comes out of the preceding solve, that throws a wobbly.

I think the MWE can be reduced to the following:

from firedrake import *
from firedrake_adjoint import *

n = 100
mesh = PeriodicSquareMesh(n,n,1, direction='both', quadrilateral = True)
dx = Measure('dx', domain = mesh)
V = VectorFunctionSpace(mesh, 'CG', 1)
R = FunctionSpace(mesh, 'R', 0)
# R = FunctionSpace(mesh, 'CG', 1)  - works if you use a "normal" function space
W = V*R*R
w = Function(W)
v, lamb0, lamb1 = split(w)
ww = assemble(inner(grad(v), grad(v))*dx)
dW = compute_gradient(ww, Control(w))

Which works correctly if you replace R by a full function space (same for the original MWE). So I think this is just one of those cases where the R isn't fully supported.

As a workaround, I think you could probably reformulate your problem solving F==0 without the Lagrange multipliers, solving it as a singular system by providing the appropriate nullspace vectors, and then applying the appropriate constraint afterwards through projection...

[dthillaithevan]

Thanks for this @stephankramer! I have tried using nullspaces in the past but have run into convergence issues (for non-linear problems) which is why I switched to a Lagrange multiplier approach instead, which has proven to be must more stable.

Are there plans to add fixes to support R for adjoint problems/fix this issue in the near future?

[dthillaithevan]

Hi @dham, just wanted to check if there are any updates on this?

[dham]

OK here's a significantly more minimal case. This extracts the failing operation from the adjoint, so we're now down to a failing assemble with no adjoint aspect.

from firedrake import *

n = 4
mesh = UnitSquareMesh(n,n)
V = VectorFunctionSpace(mesh, 'CG', 1)
R = FunctionSpace(mesh, 'R', 0)
# R = FunctionSpace(mesh, 'CG', 1)  - works if you use a "normal" function space
W = V*R
w = Function(W)
v, lamb0 = split(w)
assemble(derivative(inner(grad(v), grad(v))*dx, w))

[dham]

OK, I think this is fixed in #2884. Try the fix_empty_real_assembly branch.

[ryan-david-murphy]

Thank you @dham, really appreciate this.

[dthillaithevan]

This is great, thanks @dham!

[LeeLizuoLiu]

OK here's a significantly more minimal case. This extracts the failing operation from the adjoint, so we're now down to a failing assemble with no adjoint aspect.

from firedrake import *

n = 4
mesh = UnitSquareMesh(n,n)
V = VectorFunctionSpace(mesh, 'CG', 1)
R = FunctionSpace(mesh, 'R', 0)
# R = FunctionSpace(mesh, 'CG', 1)  - works if you use a "normal" function space
W = V*R
w = Function(W)
v, lamb0 = split(w)
assemble(derivative(inner(grad(v), grad(v))*dx, w))

I may want to know, whether there is a way to compute the norm of the assembled derivative? I tried norm(assemble(derivative(inner(grad(v), grad(v))*dx, w)), it doesn't work.

[connorjward]

OK here's a significantly more minimal case. This extracts the failing operation from the adjoint, so we're now down to a failing assemble with no adjoint aspect.

from firedrake import *

n = 4
mesh = UnitSquareMesh(n,n)
V = VectorFunctionSpace(mesh, 'CG', 1)
R = FunctionSpace(mesh, 'R', 0)
# R = FunctionSpace(mesh, 'CG', 1)  - works if you use a "normal" function space
W = V*R
w = Function(W)
v, lamb0 = split(w)
assemble(derivative(inner(grad(v), grad(v))*dx, w))

I may want to know, whether there is a way to compute the norm of the assembled derivative? I tried norm(assemble(derivative(inner(grad(v), grad(v))*dx, w)), it doesn't work.

Could you please open a new discussion for this? This one has already been closed.