Adjoint solution using only boundary values of the forward Poisson solution

dolfin-adjoint
1

Hello,

I am trying to solve the adjoint problem using the code bellow, with only boundary values instead of the entire u solution values.

from mshr import *
from dolfin import *
from dolfin_adjoint import *
import matplotlib.pyplot as plt
import numpy as np

resolution = 185
domain_geometry = Circle(dolfin.Point(0.5, 0.5), 0.5+0.01)
mesh = Mesh(generate_mesh(domain_geometry, resolution))

# Build function space with Lagrange multiplier
P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
R = FiniteElement("Real", mesh.ufl_cell(), 0)
W = FunctionSpace(mesh, P1 * R)

# Define variational problem
(u, c) = TrialFunction(W)
(v, d) = TestFunctions(W)

D = Constant(1) # Control parameter for dirac delta
x = SpatialCoordinate(mesh)
coords_mesh = V.tabulate_dof_coordinates()
xx = [i[0] for i in coords_mesh]
yy = [i[1] for i in coords_mesh]
center1 = Constant((min(xx, key=lambda i:abs(i-0.2)), min(yy, key=lambda i:abs(i-0.2))))
center2 = Constant((min(xx, key=lambda i:abs(i-0.8)), min(yy, key=lambda i:abs(i-0.8))))
eps = Constant(1e-6) 
f1 = (D * (eps/pi)/(dot(x-center1, x-center1)**2 + eps**2))
f2 = (D * (eps/pi)/(dot(x-center2, x-center2)**2 + eps**2))
f_list = [f1, f2]
g = Constant(0.0)
a = (inner(grad(u), grad(v)) + c*v + u*d)*dx
L = sum(inner(f, v)*dx for f in f_list)

# Compute solution
w = Function(W)
solve(a == L, w)
(u, c) = w.split()

Get u values on boundary :

import numpy as np
vertex_function = MeshFunction("size_t", mesh, 0)

class BoundaryMarker(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary

BoundaryMarker().mark(vertex_function, 1)
boundary_vertices = np.asarray(vertex_function.where_equal(1))


v_to_d = np.array(W.dofmap().dofs())
dofs = v_to_d[boundary_vertices]

x = W.tabulate_dof_coordinates()
print('Boundary x y coords, with u values on each coord')
for dof in dofs:
    print(x[dof], u.vector()[dof])

Adjoint solution :

ud = Function(W)

j = inner ( u.vector()[dof] - ud , u.vector()[dof] - ud ) * dx
J = assemble(j)

fc = [Control(D)]

Jhat = ReducedFunctional(J, fc)
dJdm = Jhat.derivative()

# solution of the adjoint problem
tape = get_working_tape()
solve_block = tape.get_blocks()[1]
breakpoint()
adj_u = solve_block.adj_sol.vector()[:]

adj_u_func = Function(V)
adj_u_func.assign(u)
adj_u_func.vector()[:]= adj_u

On these lines :

j = inner ( u.vector()[dof] - ud , u.vector()[dof] - ud ) * dx
J = assemble(j)

Error is : “Can’t add expressions with different shapes”

Thank you for your help.

2

Just use the ds measure, i.e. j = inner(u -ud, u -ud) * ds, the L^2(\partial\Omega) norm squared

1 Like
4

Without a minimal reproducible example I cannot give you any further guidance.

6

I cannot reproduce your problem as mshr is long deprecated and I don’t have it on my system. Please check your output in Paraview with:

with XDMFFile(MPI.COMM_WORLD, "adjoint.xdmf") as file:
    file.write(adj_w_func)

first

8

Again, as you insist on using mshr, a deprecated meshing software, I cannot run your code and thus not reproduce the error. I would for instance suggest that you use the UnitDiscMesh for starters, which is part of DOLFIN: Bitbucket

9

@dokken I adapted this code in order to have a UnitSquareMesh. Can you run it like this ?

Thank you.

10

You are accessing the wrong block. You should access the solve block.
You can inspect the blocks with

print(tag.get_blocks())

This works for me:

solve_block = tape.get_blocks()[0]
adj_w = solve_block.adj_sol.vector()[:]

adj_w_func = Function(W)
adj_w_func.assign(w)
adj_w_func.vector()[:] = adj_w
plot(adj_w_func)
plt.colorbar(plot(adj_w_func))
1 Like
11

Ah yes the problem was indeed the block choice !

Thank you very much.