Shape Derivatives in UFL

dolfin-adjoint
1

Hi,

I’ve got a quick question regarding dolfin-adjoint, in particular the computation of shape derivatives. Assuming that I try to solve the optimization problem of the form

\min J(\Omega) = \int_\Omega (y - y_d)^2 dx

where y_d is some desired state and y(\Omega) is the solution of a PDE on \Omega.
The shape derivative for this is then given by

dJ(\Omega)[V] = \int_\Omega (y - y_d)^2 \nabla \cdot (V) - 2*(y - y_d)* \nabla y_d [V] dx

As far as I understand, if y_d is an ufl expression only depending on a SpatialCoordinate, then the shape derivative will be correct. However, if y_d is a Function or Expression, then the above code will not include the second term.

Is there any way to treat such a problem currectly when using a Function / Expression?

In the following I created a MWE highlighting the problem

from fenics import *
import numpy as np



mesh = UnitSquareMesh(25, 25)
dx = Measure('dx', mesh)
x = SpatialCoordinate(mesh)

space = FunctionSpace(mesh, 'CG', 1)
expr = Expression('sin(pi*x[0])*sin(pi*x[1])', degree=2)
y = interpolate(expr, space)

def_space = VectorFunctionSpace(mesh, 'CG', 1)
V = TestFunction(def_space)

y_d_coor = pow(x[0], 2) + pow(x[1], 2) - pow(0.5, 2)
y_d_expr = Expression('pow(x[0], 2) + pow(x[1], 2) - pow(0.5, 2)', degree=2, domain=mesh)
y_d_func = interpolate(y_d_expr, space)
y_d_list = [y_d_coor, y_d_expr, y_d_func]

J_list = [pow(y - y_d, 2)*dx for y_d in y_d_list]
dJ_ex_list = [pow(y - y_d, 2)*div(V)*dx - 2*(y - y_d)*inner(grad(y_d), V)*dx for y_d in y_d_list]
dJ_mod_list = [pow(y - y_d, 2)*div(V)*dx for y_d in y_d_list]
dJ_ad_list = []

for idx in range(len(J_list)):

    vec_ex = assemble(dJ_ex_list[idx])[:]
    vec_mod = assemble(dJ_mod_list[idx])[:]
    vec_ad = assemble(derivative(J_list[idx], x, V))[:]

    error_ex = np.max(np.abs(vec_ex - vec_ad)) / np.max(np.abs(vec_ex)) * 100
    error_mod = np.max(np.abs(vec_mod - vec_ad)) / np.max(np.abs(vec_mod)) * 100

    print('Relative error: ' + str(error_ex) + ' % (exact)')
    print('Relative error: ' + str(error_mod) + ' % (modified)')
    print('')
2 Likes
2

I suggest having a look at the following papers:

As you can see here, to get the shape derivative of a function, you need to solve the corresponding adjoint equation. This can be done by hand (as in the first paper) or by using dolfin-adjoint (as implemented in the second paper). You cannot use Expressions in either case, and you need to use SpatialCoordinate for purely geometrical dependencies

1 Like
3

I have already read the two papers, but I was wondering if there is any possibility to do it anyway, using a function. The corresponding derivatives needed can be computed for this (as in my example above).

Think of a shape identification problem, where I want to identify some kind of shape in a box, by a transmission problem. I could solve the PDE with the desired shape to get the solution, which is then used as a desired state, but there is no analytical formula for this, so that I can only use a Function, and not an UFL expression with SpatialCoordinate…

4

As I said above, to Get the shape derivative of a PDE dependent variable, you need to solve the corresponding adjoint problem. In your problem above, you have not included the solution of the PDE, and your shape derivative will be incorrect.

As I also said in the post above, dolfin-adjoint can compute these automatically for you (even for functions).

5

Could you give me any example for where this is done (like for a tracking type cost functional), where one function depends on the PDE constraint and the other does not? In the Papers you linked there is no such example… Basically I just need the code and can figure the rest out for myself.

6

I would claim you haven’t looked carefully at the dolfin-adjoint paper, as the two last examples there use volume and barycenter constraints weakly, which is not a PDE constraint.
You can even find a documented example for this at the pyadjoint examples page

7

Yes, I’ve seen that, but it does not apply, as the Constraints for Volume and Barycenter are implemented with “SpatialCoordinate”, and not using a “Function”, which is what I want to do…

8

It would be nice to see what kind of expression you are going to handle, as the problem above does not highlight that.

Below I’ve used dolfin-adjoint, via projection, to have u_d as a function.

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

mesh = UnitSquareMesh(10, 10)
V = FunctionSpace(mesh, "CG", 1)
x = SpatialCoordinate(mesh)

S = VectorFunctionSpace(mesh, "CG", 1)
s = Function(S)
ALE.move(mesh, s)
u_d = project(x[0]-x[1], V)

J = assemble(inner(u_d, u_d)*dx)

Jhat = ReducedFunctional(J, Control(s))

dJds = Jhat.derivative()

v = TestFunction(S)
dJ_exact = assemble(div(v)*inner(u_d, u_d)*dx + 2*inner(u_d, dot(grad(u_d), v))*dx)
assert(np.allclose(dJ_exact.get_local(), dJds.vector().get_local()))
9

I see, and now I get that my question was not well-formulated, and actually it is not so much related to dolfin-adjoint, but to its implementation…

What I want to do is basically to reproduce the result, i.e., the shape derivative, but actually without using dolfin-adjoint, but only the “derivative” function of UFL, as this can handle shape derivatives as well. I’ve modified your example to use the derivative function below, but here you can see that the second term of the shape derivative which comes from the pull back of u_d (that does not depend on the domain) is missing. And my question is, how to get this term automatically…

import numpy as np
from dolfin import *

mesh = UnitSquareMesh(10, 10)
dx = Measure('dx', mesh)
V = FunctionSpace(mesh, "CG", 1)
x = SpatialCoordinate(mesh)

S = VectorFunctionSpace(mesh, "CG", 1)
u_d = project(x[0]-x[1], V)

J = inner(u_d, u_d)*dx

v = TestFunction(S)

dJds = assemble(derivative(J, x, v))

dJ_exact = assemble(div(v)*inner(u_d, u_d)*dx + 2*inner(u_d, dot(grad(u_d), v))*dx)
dJ_modified = assemble(div(v)*inner(u_d, u_d)*dx)
if np.allclose(dJ_exact.get_local(), dJds.get_local()):
    print('Exact shape derivative')
elif np.allclose(dJ_modified.get_local(), dJds.get_local()):
    print('Modified shape derivative')
else:
    print('Test failed')

Again, I’m very sorry for the poor formulation of the question previously, I hope you can understand me now. And thanks a lot for your time.

10

There is an alternative package for shape derivative called FEmorph, which is the predecessor to the ufl derivative and the dolfin-adjoint implementation.

I worked on cleaning up that software a couple of years back. I do not recall if it does the pull back or not, I guess have a look at the implementation paper.

The code, as i linked to above, is not actively maintained, so it is not optimal to use.

Anyhow, concluding, i would say that what you are asking for might not be implemented. I am still curious as to what the content of your function/expression is going to be (as in many cases, you could frame it as a projection of a mixture of ufl expressions/constants)

11

Thanks a lot, I will take a look at this.
So I guess dolfin-adjoint is able to do this due to “real” AD functionality needed for it.

The problems would be, e.g., of the type discussed here: https://epubs.siam.org/doi/abs/10.1137/130930807
(Section 4), or for a particular application see https://arxiv.org/abs/1911.06819 (Section 2.3)

To generate the desired state for the tracking type cost functionals it might be necessary to numerically solve a PDE, so that the solution might not be given explicitely.

12

But as long the desired state is a solution of a PDE, you can use dolfin adjoint as above (replace the projection of x[0]*x[1] with solving any other PDE).

Is there any specific reason for not wanting to use dolfin-adjoint?

13

Okay, well I might consider this. There is no reason in particular, maybe only that I want to use my own code for the optimization to be able to better adapt it for my needs (I know that dolfin-adjoint is open source, but I guess its quicker for me to write my own code than to familiarize myself with all of the code of dolfin-adjoint)

Thanks a lot for your help and for clearing this up.

14

You do not have to use the built in optimization packages in dolfin-adjoint. As you can obtain the derivative (using Jhat.derivarive(), as I did in the example above) and you can use this in any custom optimization scheme.
This is for instance highlighted in the paper I linked to above on shape derivatives in dolfin-adjoint (Stokes example).

15

I have a follow up question in that case. Lets assume I follow along your example, but do a H^1_0 projection of the function instead (after moving the mesh, so that we have the pull back, the other case as state equation works fine). I have created this example

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

mesh = UnitSquareMesh(10, 10)
V = FunctionSpace(mesh, "CG", 1)
x = SpatialCoordinate(mesh)

S = VectorFunctionSpace(mesh, "CG", 1)
s = Function(S)

u_d_expr = pow(x[0], 2) + pow(x[1], 2) - pow(0.5, 2)

a = inner(grad(TrialFunction(V)), grad(TestFunction(V)))*dx
L = u_d_expr*TestFunction(V)*dx

def bdry(x, on_boundary):
    return on_boundary

bc = DirichletBC(V, Constant(0), bdry)

u_d = Function(V)

ALE.move(mesh, s)
solve(a==L, u_d, bc)


J = assemble(inner(u_d, u_d)*dx)

Jhat = ReducedFunctional(J, Control(s))

dJds = Jhat.derivative()

v = TestFunction(S)
dJ_exact = assemble(div(v)*inner(u_d, u_d)*dx + 2*inner(u_d, dot(grad(u_d), v))*dx)

error = np.max(np.abs(dJ_exact[:] - dJds.vector()[:])) / np.max(np.abs(dJ_exact[:])) * 100
print('Relative Error: ' + str(error) + ' %')

This does not seem to give the correct shape derivative… (as I said, if I exchange the ALE.move and solve commands and correct the formula for the exact shape derivative, it works)

16

This is due to the following:

  • If you move your mesh, then do your projection, you are projecting the spatial coordinates of the transformed mesh, and there will be no pull-back.
  • If you project the function, and then move your mesh, the nodal values (values of the function of the dof coordinate) is moved along with the mesh, and you will get the pull-back effect.

These ideas are nicely written in Berggren’s unified shape analysis paper, Example 1 and 2 (The links should take you directly to that page).

17

This is clear to me.
If I manually do an L2 projection, then everything works as expected just as when I use “project”, but running the very same code with the H^1_0 inner product does only give the correct derivative for the case that first the PDE is solved, and then the mesh is moved, but not for the other case.

18

If you move the mesh, then solve the PDE, your problem becomes more complex. Then you have to solve the corresponding adjoint problem created by the project function.
I have created a verification of this below (where Im using the H1 projection as I do not want to go into the Dirichlet BC boundary adjoint issue now).
In this code, I also verify what the shape derivatives produced by the UFL derivative is for this specific case.

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

def compute_manually():
    from ufl import replace
    mesh = UnitSquareMesh(5,5)
    V = FunctionSpace(mesh, "CG", 1)
    x = SpatialCoordinate(mesh)

    u_d_expr = x[0] + x[1]
    u,v = Function(V), TestFunction(V)
    F = inner(grad(u), grad(v))*dx + inner(u, v)*dx - inner(u_d_expr, v)*dx
    bc = DirichletBC(V, Constant(0), "on_boundary")
    solve(F==0, u)#, bc)
    j = inner(u, u)*dx
    J = assemble(j)

    # Solve adjoint
    lmbd = Function(V)
    ut = TrialFunction(V)
    a_adj = adjoint(derivative(F, u, ut))
    L_adj = -derivative(j, u, TestFunction(V))
    solve(a_adj==L_adj, lmbd)#, bc)

    # Solve boundary adjoint equation
    # FIXME: not complete, but something along these lines
    # adj_bdy = assemble(derivative(j, u, TestFunction(V))-action(a_adj,lmbd))
    # lmbd.vector()[:] += adj_bdy

    S = VectorFunctionSpace(mesh, "CG", 1)
    s = TestFunction(S)

    # Verify shape derivative
    L_tmp = action(F, lmbd)
    dFdX = derivative(L_tmp, x)
    dFex = assemble(div(s)*(inner(grad(u), grad(lmbd)) + inner(u-u_d_expr,lmbd) )*dx\
                    -inner(dot(grad(s), grad(u)), grad(lmbd)) * dx\
                    -inner(grad(u),dot(grad(s), grad(lmbd))) * dx\
                    - dot(s, grad(u_d_expr))*lmbd*dx)
    assert np.allclose(assemble(dFdX).get_local(), dFex.get_local())

    dJdX = derivative(j, x)
    dJex = assemble(div(s)*u*u*dx)
    assert np.allclose(assemble(dJdX).get_local(), dJex.get_local())
    print("Adjoint", lmbd.vector().get_local())

    return assemble(dFdX) + assemble(dJdX)


with stop_annotating():
    dJ_exact = compute_manually()


mesh = UnitSquareMesh(5, 5)
V = FunctionSpace(mesh, "CG", 1)
x = SpatialCoordinate(mesh)

S = VectorFunctionSpace(mesh, "CG", 1)
s = Function(S)

bc = DirichletBC(V, Constant(0), "on_boundary")

u_d_expr = x[0] + x[1]
u,v = Function(V), TestFunction(V)
F = inner(u,v)*dx + inner(grad(u), grad(v))*dx - inner(u_d_expr, v)*dx

ALE.move(mesh, s)
solve(F==0, u)#, bc)


J = assemble(inner(u, u)*dx)

Jhat = ReducedFunctional(J, Control(s))

dJds = Jhat.derivative()

with stop_annotating():
    dJds = Jhat.derivative()
    s1 = project(as_vector((x[0],x[1])),S)
    result = taylor_to_dict(Jhat, s, s1)
    print(result)

print(dJds.vector().get_local())
print(dJ_exact.get_local())
assert(np.allclose(dJ_exact.get_local(), dJds.vector().get_local()))
19

Here is the H_0^1 manual implementation with your source term:

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

def compute_manually():
    from ufl import replace
    mesh = UnitSquareMesh(5,5)
    V = FunctionSpace(mesh, "CG", 1)
    x = SpatialCoordinate(mesh)

    u_d_expr = pow(x[0], 2) + pow(x[1], 2) - pow(0.5, 2)
    u,v = Function(V), TestFunction(V)
    F = inner(grad(u), grad(v))*dx - inner(u_d_expr, v)*dx
    bc = DirichletBC(V, Constant(0), "on_boundary")
    solve(F==0, u, bc)
    j = inner(u, u)*dx
    J = assemble(j)

    # Solve adjoint
    lmbd = Function(V)
    ut = TrialFunction(V)
    a_adj = adjoint(derivative(F, u, ut))
    L_adj = -derivative(j, u, TestFunction(V))
    solve(a_adj==L_adj, lmbd, bc)

    S = VectorFunctionSpace(mesh, "CG", 1)
    s = TestFunction(S)

    # Verify shape derivative
    L_tmp = action(F, lmbd)
    dFdX = derivative(L_tmp, x)
    dFex = assemble(div(s)*(inner(grad(u), grad(lmbd)) + inner(-u_d_expr,lmbd) )*dx\
                    -inner(dot(grad(s), grad(u)), grad(lmbd)) * dx\
                    -inner(grad(u),dot(grad(s), grad(lmbd))) * dx\
                    - dot(s, grad(u_d_expr))*lmbd*dx)
    assert np.allclose(assemble(dFdX).get_local(), dFex.get_local())

    dJdX = derivative(j, x)
    dJex = assemble(div(s)*u*u*dx)
    assert np.allclose(assemble(dJdX).get_local(), dJex.get_local())

    return assemble(dFdX) + assemble(dJdX)


with stop_annotating():
    dJ_exact = compute_manually()


mesh = UnitSquareMesh(5, 5)
V = FunctionSpace(mesh, "CG", 1)
x = SpatialCoordinate(mesh)

S = VectorFunctionSpace(mesh, "CG", 1)
s = Function(S)

bc = DirichletBC(V, Constant(0), "on_boundary")

u_d_expr = pow(x[0], 2) + pow(x[1], 2) - pow(0.5, 2)
u,v = Function(V), TestFunction(V)
F = inner(grad(u), grad(v))*dx - inner(u_d_expr, v)*dx

ALE.move(mesh, s)
solve(F==0, u, bc)


J = assemble(inner(u, u)*dx)

Jhat = ReducedFunctional(J, Control(s))


with stop_annotating():
    dJds = Jhat.derivative()
    s1 = project(as_vector((x[0],x[1])),S)
    result = taylor_to_dict(Jhat, s, s1)
    print(result)

print(dJds.vector().get_local())
print(dJ_exact.get_local())
assert(np.allclose(dJ_exact.get_local(), dJds.vector().get_local()))
20

Dear @dokken,

thanks a lot for all this work, this already cleared up a lot of things. Now, I’ve modified your H^1_0 example so that we now actually have a tracking type example, which would be actually what I’m looking for.

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

def compute_manually():
    from ufl import replace
    mesh = UnitSquareMesh(5,5)
    V = FunctionSpace(mesh, "CG", 1)
    x = SpatialCoordinate(mesh)

    u_d_expr = pow(x[0], 2) + pow(x[1], 2) - pow(0.5, 2)
    u_d_expr = project(u_d_expr, V)
    u,v = Function(V), TestFunction(V)
    F = inner(grad(u), grad(v))*dx - Constant(1)*v*dx
    bc = DirichletBC(V, Constant(0), "on_boundary")
    solve(F==0, u, bc)
    j = inner(u - u_d_expr, u - u_d_expr)*dx
    J = assemble(j)

    # Solve adjoint
    lmbd = Function(V)
    ut = TrialFunction(V)
    a_adj = adjoint(derivative(F, u, ut))
    L_adj = -derivative(j, u, TestFunction(V))
    solve(a_adj==L_adj, lmbd, bc)

    S = VectorFunctionSpace(mesh, "CG", 1)
    s = TestFunction(S)

    # Verify shape derivative
    L_tmp = action(F, lmbd)
    dFdX = derivative(L_tmp, x)
    dFex = assemble(div(s)*(inner(grad(u), grad(lmbd)) - Constant(1)*lmbd)*dx\
                    -inner(dot(grad(s), grad(u)), grad(lmbd)) * dx\
                    -inner(grad(u),dot(grad(s), grad(lmbd))) * dx)

    assert np.allclose(assemble(dFdX).get_local(), dFex.get_local())

    dJdX = derivative(j, x)
    dJex = assemble(div(s)*inner(u - u_d_expr, u - u_d_expr)*dx - 2*(u - u_d_expr)*inner(grad(u_d_expr), s)*dx)
    assert np.allclose(assemble(dJdX).get_local(), dJex.get_local())

    return assemble(dFdX) + assemble(dJdX)



with stop_annotating():
    dJ_exact = compute_manually()


mesh = UnitSquareMesh(5, 5)
V = FunctionSpace(mesh, "CG", 1)
x = SpatialCoordinate(mesh)

S = VectorFunctionSpace(mesh, "CG", 1)
s = Function(S)

bc = DirichletBC(V, Constant(0), "on_boundary")

u_d_expr = pow(x[0], 2) + pow(x[1], 2) - pow(0.5, 2)
u_d_expr = project(u_d_expr, V)
u,v = Function(V), TestFunction(V)
F = inner(grad(u), grad(v))*dx - Constant(1)*v*dx

ALE.move(mesh, s)
solve(F==0, u, bc)


J = assemble(inner(u - u_d_expr, u - u_d_expr)*dx)

Jhat = ReducedFunctional(J, Control(s))


with stop_annotating():
    dJds = Jhat.derivative()
    s1 = project(as_vector((x[0],x[1])),S)
    result = taylor_to_dict(Jhat, s, s1)
    print(result)

print(dJds.vector().get_local())
print(dJ_exact.get_local())
assert(np.allclose(dJ_exact.get_local(), dJds.vector().get_local()))

Note, that everything works fine and as intended, if I uncomment the lines

u_d_expr = project(u_d_expr, V)

which would be the auxilliary PDE solve, and could be respresented by solving one particular equation, saving the dofs and then reinserting them to this function. Again, it might not be possible to give u_d_expr analytically. So when the line is active, the assertion fails. The critical term is just the one coming from the pull back in J

- 2*(u - u_d_expr)*inner(grad(u_d_expr), s)*dx

so when I leave this out and do the projection, everything works fine. So I guess it is just not possible to apply the pull back to arbitrary functions, but only to ufl expressions using SpatialVariable, or am I wrong?

If it helps, I can send you the examples where the corresponding parts are changed.