Dear all,
I have been stuck with an optimal control problem for a few days. I want to find the sensitivity of a cost function with respect to a forcing term of a PDE. I have solved the problem writing my own Matlab code and the Taylor test passes. However, when I implement it in Dolfinx, the Taylor test does not pass. I present here a simplified version of the problem with the code.
Consider the Poisson equation over a squared domain \Omega with homogeneous Dirichlet boundaries \Gamma:
-\Delta u = f in \Omega
u = 0 at \Gamma
I define the following cost function:
J(u) = \frac{1}{2}\int_{\Omega}||u-u_{obs}||^{2}d\Omega
where u_{obs} is the desired profile of u.
u_{obs}(x,y) = x(1-x)y(1-y)
This problem can be found in https://fenicsproject.org/pub/course/lectures/2015-08-logg-beijing/lecture_12_sensitivities.pdf solved with dolfin-adjoint.
Using the adjoint method, I compute the gradient of the cost function J with respect to the design parameter \boldsymbol{m} (i.e., f in this case) as follows:
\frac{dJ}{d\textbf{m}} = \frac{\partial{J}}{\partial\textbf{m}} + \left(\frac{\partial{\textbf{R}}}{\partial\textbf{m}}\right)^{T} \boldsymbol{\lambda}
where \boldsymbol{\lambda} is the adjoint variable, solution of the adjoint problem:
\left(\frac{\partial{\textbf{R}}}{\partial\textbf{u}}\right)^{T}\boldsymbol{\lambda} = -\frac{\partial{J}}{\partial\textbf{u}}
My code for solving this problem is below:
import dolfinx
from ufl import Measure, inner, grad, div, TestFunction, TrialFunction, adjoint, derivative
from dolfinx.mesh import CellType, create_rectangle
from dolfinx import fem, nls, io
from mpi4py import MPI
from petsc4py import PETSc
import numpy as np
import time
import dolfinx.nls.petsc
import dolfinx.fem.petsc
# mesh and function space
mesh = create_rectangle(MPI.COMM_WORLD, [[0, 0], [1, 1]], [64, 64], CellType.triangle)
V = fem.FunctionSpace(mesh, ("CG", 1))
u = fem.Function(V)
v = TestFunction(V)
# sensitivity will be calculated wrt this
f = fem.Function(V) # The forcing f is the design variable
f.x.array[:] = 0.0
dx = Measure("dx")
# Residual of the variational form of Poisson problem
R = inner(grad(u), grad(v))*dx - inner(f,v)*dx
# boundary conditions
def left(x):
return np.isclose(x[0], 0)
def right(x):
return np.isclose(x[0], 1)
def top(x):
return np.isclose(x[1], 1)
def bottom(x):
return np.isclose(x[1], 0)
fdim = mesh.topology.dim - 1
uD = fem.Function(V)
left_facets = dolfinx.mesh.locate_entities_boundary(mesh, fdim, left)
right_facets = dolfinx.mesh.locate_entities_boundary(mesh, fdim, right)
top_facets = dolfinx.mesh.locate_entities_boundary(mesh, fdim, top)
bottom_facets = dolfinx.mesh.locate_entities_boundary(mesh, fdim, bottom)
with uD.vector.localForm() as loc:
loc.set(0.0)
bc1 = fem.dirichletbc(uD, dolfinx.fem.locate_dofs_topological(V, fdim, left_facets))
bc2 = fem.dirichletbc(uD, dolfinx.fem.locate_dofs_topological(V, fdim, right_facets))
bc3 = fem.dirichletbc(uD, dolfinx.fem.locate_dofs_topological(V, fdim, top_facets))
bc4 = fem.dirichletbc(uD, dolfinx.fem.locate_dofs_topological(V, fdim, bottom_facets))
bcs = [bc1, bc2, bc3, bc4]
# Solver for the direct problem
problem_dir = fem.petsc.NonlinearProblem(R, u, bcs)
solver = nls.petsc.NewtonSolver(MPI.COMM_WORLD, problem_dir)
solver.rtol = 1e-16
# Compute solution
solver.solve(u)
# Define desired u profile
def profile(x):
return x[0]*(1-x[0])*x[1]*(1-x[1])
u_obs = fem.Function(V)
u_obs.interpolate(profile)
# Define cost function
J = (1/2)*inner(u-u_obs, u-u_obs)*dx
J_form = fem.form(J)
J_val = fem.assemble_scalar(J_form)
# Bilinear and linear form of the adjoint problem
lhs = adjoint(derivative(R, u))
rhs = -derivative(J, u)
# Create adjoint problem solver
problem = dolfinx.fem.petsc.LinearProblem(lhs, rhs, bcs=bcs,petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
# solve adjoint vector
lmbda = problem.solve()
# Compute sensitivity: dJ/dm
# Partial derivative of J w.r.t. m
dJdm_form = fem.form(derivative(J, f))
dJdm = fem.petsc.assemble_vector(dJdm_form)
dJdm.assemble()
# partial derivative of R w.r.t. m
dRdm_form = fem.form(adjoint(derivative(R, f)))
dRdm = fem.petsc.assemble_matrix(dRdm_form)
dRdm.assemble()
# Calculate the gradient of J with respect to m (the forcing f)
dJdm += dRdm*lmbda.vector
dJdm.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
Now I define the following functions to perform the Taylor test:
# Function to compute the cost
def J_cost(m_n):
f.x.array[:] = m_n
solver.solve(u)
J_val = fem.assemble_scalar(J_form)
return J_val
# Function to compute the gradient of the cost w.r.t. f
def grad_J(m_n):
f.x.array[:] = m_n
# Direct solution
solver.solve(u)
# Adjoint solution
lmbda = problem.solve()
# Reassemble dJdm
dJdm = fem.petsc.assemble_vector(dJdm_form)
dJdm.assemble()
# Reassemble dRdm
dRdm = fem.petsc.assemble_matrix(dRdm_form)
dRdm.assemble()
# Compute gradient
dJdm += dRdm*lmbda.vector
dJdm.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.REVERSE)
return dJdm
# Taylor test
def taylor_test(cost, grad, m_0, p=1e-6, n=5):
g0 = grad(m_0)
l0 = cost(m_0)
reminder = []
perturbance = []
for i in range(0, n):
l1 = cost(m_0+p)
reminder.append(l1 - l0 - sum(np.transpose(g0)*p))
perturbance.append(p)
p /= 2
conv_rate = convergence_rates(reminder, perturbance)
return reminder, perturbance, conv_rate
# Compute convergence rate (they should be 2.0 for small p)
def convergence_rates(r, p):
cr = [] # convergence rates
for i in range(1, len(p)):
cr.append(np.log(r[i] / r[i - 1])
/ np.log(p[i] / p[i - 1]))
return cr
m_0 = np.zeros(len(u.x.array))
r = taylor_test(J_cost, grad_J, m_0)
print("Convergence rate:", r[2])
The output should be 2.0, but gives
Convergence rate: [1.8214896414517914, 1.458183253415033, 0.6623141576735944, 0.33714267999803627]
Does anybody know what I am doing wrong?
I done a similar problem in this post Calculate adjoint derivatives w.r.t to multiple Constants where the control is a constant, but when the control is a function, I can’t get the expected solution in dolfinx.
If something is not clear, please let me know.
Thank you all in advance.