I am trying to solve the adjoint problem for poisson equation using 2 source points (i.e Dirac). I have already done this using an Expression for the source term (f), which can be found here Poisson adjoint.
But the difference is that now I need to go over an apply() function, for a b parameter that is the assemble of the RHS and the sources points (f1, f2).
This part works fine :
from mshr import *
from dolfin import *
from dolfin_adjoint import *
import matplotlib.pyplot as plt
import numpy as np
resolution = 15
domain_geometry = Circle(dolfin.Point(0.5, 0.5), 0.5+0.01)
mesh = generate_mesh(domain_geometry, resolution)
V = FunctionSpace(mesh, "Lagrange", 1)
def boundary(x, on_boundary):
return on_boundary
# Define boundary condition
u0 = Constant(0.0)
bc = DirichletBC(V, u0, boundary)
# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f1 = PointSource(V, Point(0.2, 0.2), 1)
f2 = PointSource(V, Point(0.8, 0.8), 1)
g = Constant(0.0)
a = inner(grad(u), grad(v))*dx
L = Constant(0.0)*v*dx
A = assemble(a)
bc.apply(A)
b = assemble(L)
f1.apply(b)
f2.apply(b)
bc.apply(b)
# Compute solution
u = Function(V)
solve(A, u.vector(), b)
plot(u)
plt.colorbar(plot(u))
But the error shows up when I start the adjoint steps :
ud = Function(V)
j = inner ( u - ud , u - ud ) * dx
J = assemble(j)
bc = Control(b)
dJdf = compute_gradient (J , bc)
# solution of the adjoint problem
tape = get_working_tape()
solve_block = tape.get_blocks()[1]
adj_u = solve_block.adj_sol.vector()[:]
adj_u_func = Function(V)
adj_u_func.assign(u)
adj_u_func.vector()[:]= adj_u
plot(adj_u_func)
plt.colorbar(plot(adj_u_func))
And more precisely, at this line : bc = Control(b)
With error saying : Vector object has no attribute ‘block_variable’.
If you look into my previous example, i.e : Poisson adjoint, you can see that my source term is the control of my gradient computation. In that example, the source term was an Expression, but here I am trying to use the same thinking with 2 source terms that are PointSource elements.
Is it possible to use those source points as control ? I went directly to the b parameter because it includes both of them, but you can tell me if it’s wrong and how can I consider taking directly the Pointsources (f1 & f2).
The point-source class is not supported by dolfin-adjoint.
I would probably create smoothened point-sources, using a mollified dirac delta function (Dirac delta function - Wikipedia)
D = Constant(1) # Control parameter for dirac delta
x = ufl.SpatialCoordinate(mesh)
center = Constant((c_x, c_y))
eps = Constant(1e-6)
point_source = D * (eps/pi)/(dot(x-center, x-center) + eps**2)
Please provide a minimal reproducible example.
As the code has changed Since the original post, a new post with the current code and error message is required
No, it shouldn’t. The control variable is the magnitude of the source function, not the whole function.
If you want to vary the width of the dirac delta, add epsilon to the list of control parameters.
The adjoint solution will not have a point source in it, as the point source is not a function of u.
Think about it this way, you are adjointing the system Au=b(point-sources).
The adjoint system is
A lmbda = -dJ/du,
thus no sources.
No, it shouldn’t. The control variable is the magnitude of the source function, not the whole function.
If you want to vary the width of the dirac delta, add epsilon to the list of control parameters.
I don’t understand why you are talking about Control when the f_list is just a list of my sources (f1 & f2) that I am using before the adjoint problem. There is no reason that we can’t solve u with f_list = [f1, f2], which is what I have done by showing you the u plot.
Now when we get to the adjoint problem, we have an u solution and we try to find the adjoint of this solution. Forget this line fc = [Control(f) for f in f_list], we use what you have suggested as fc = Control(D). By doing this, we get the 2nd plot I’ve showed you.
Think about it this way, you are adjointing the system Au=b(point-sources).
The adjoint system is
A lmbda = -dJ/du,
thus no sources.
Yes the point source is not a function of u, but dJ/du depends on the derivative of our solution u showed in my 1st plot, so the lmbda term (which I guess is the adj_sol we retrive) will vary according to our solution u.
If I am not clear, I just want to put 2 dirac in the medea, and reconstruct them. If you have any idea on how I can achieve that, please don’t hesitate.
Here is a minimal working example with the correct point sources (also applied to the variational form) and a manual derivation of the adjoint equation as verification
from dolfin import *
from dolfin_adjoint import *
mesh = UnitSquareMesh(100, 100)
V = FunctionSpace(mesh, "Lagrange", 1)
def boundary(x, on_boundary):
return on_boundary
# Define boundary condition
u0 = Constant(0.0)
bc = DirichletBC(V, u0, boundary)
# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
D = Constant(1) # Control parameter for dirac delta
x = SpatialCoordinate(mesh)
center1 = Constant((0.1, 0.1))
center2 = Constant((0.1, 0.2))
eps = Constant(1e-3)
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))*dx
L = sum(inner(f, v)*dx for f in f_list)
# Compute solution
u = Function(V)
solve(a == L, u, bc)
ud = Function(V)
j = inner(u - ud, u - 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]
adj_u = solve_block.adj_sol.vector()[:]
adj_u_func = Function(V)
adj_u_func.assign(u)
adj_u_func.vector()[:] = adj_u
adj_u_func.rename("adj_dolfin_adjoint", "adj_dolfin_adjoint")
# Compute adjoint solution
bc_ad = bc
L_adj = derivative(j, u)
adj = Function(V)
adj.rename("adjoint", "adjoint")
solve(a == L_adj, adj, [bc_ad])
with XDMFFile("adjoint_solution.xdmf") as file:
file.write_checkpoint(u, "u", 0.0, append=False)
file.write_checkpoint(adj, "dolfin-adjoint", 0.0, append=True)
file.write_checkpoint(adj_u_func, "manual", 0.0, append=True)
I don’t understand what you are trying to do any longer? I’ve now shown you have to setup and solve the Poisson problem with the magnitude of the point source as a control variable.
Did you look at the code I sent you, and the outputs written to XDMF? You can clearly see the dependency of the point source in the adjoint solution.
The latter approach you are trying to do seems like a boundary element approach, which would require dense matrices (all dofs have effect on eachother). For this you could maybe use BEMPP (https://bempp.com/)
Ok i figured it out. But when you change the point sources positions, you don’t really get the adjoint for each point source, right ? you get the 2nd plot I showed you where it looks like a single ellipsoid between the two point sources.
