Dear all,
I am working on adjoint sensitivity analysis.
Does anybody know any examples on how to compute sensitivities when the control parameters are Dirichlet BCs?
I can do it when the parameters appear explicitly in the equations, but I find the sensitivity analysis a bit complicated when the parameters are in the boundary conditions.
I know that fenics-adjoint deals with it, but I am looking for doing the same in dolfinx.
Thanks everyone!
As stated in
you should check out chapter 4.4 of
as it goes through what equations is required for strong enforcement of Dirichlet conditions.
I have a side question about this. I read through the above linked thesis, and, as far as I understand it implementation wise, the \partial F / \partial m matrix should be assembled separately because from a UFL perspective F does not directly depend on m.
Wouldnât it be easier, from an implementation perspective, to impose boundary conditions weakly, as in the tutorial, so that the implementation could leverage UFLâs AD for deriving the full model?
You can impose them weakly and let ufl do the job for you.
Historically people tend to like strong enforcement of boundary conditions, as it does not require a parameter to ensure coercivity.
I successfully solved an adjoint-based optimization problem using this approach (i.e., by imposing the boundary conditions weakly and then computing \frac{\partial F}{\partial m} using ufl.
The problem that I have encountered is that, because I need to define the BC as a dolfin function to use ufl.derivative, \frac{\partial F}{\partial m} is a matrix, although I will only use one of its columns (the one associated with the BC I am using as control). Therefore, for simulations in which I have many elemens (for example, 3D systems), this leads to memory problems.
Does any of you know how to solve this problem?
I have written a minimal version of the problem, with annotations, below:
# Example to get (partial_F)/(partial_m) by imposition of Dirichlet BC weakly
from dolfin import *
import numpy as np
# Problem parameters
c0 = 1 # sound speed
nu = 0.01 # kinematic viscosity
k = Constant(1e-2)
# Mesh
mesh = UnitIntervalMesh(50)
n = FacetNormal(mesh)
# Boundaries
class Left(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 0.0)
class Right(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 1.0)
# Define boundaries
boundaries = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundaries.set_all(0)
left = Left()
left.mark(boundaries, 1)
right = Right()
right.mark(boundaries, 2)
ds = Measure('ds')[boundaries]
# Function Space
V = FunctionSpace(mesh, "Lagrange", 1)
v = TestFunction(V)
# Boundary conditions
u_L = Expression("sin(omega*t)",omega=2*np.pi, t=0, degree=1)
u_D = np.loadtxt('u_R2.txt')
u_R = interpolate(Constant(float(u_D[0])), V)
# Weak BC parameters
h = 2*Circumradius(mesh)
alpha = Constant(100)
# Equation
u = Function(V)
u1 = Function(V)
u0 = Function(V)
u_dot = u -2*u1 + u0
F = (u_dot*v*dx + nu*k*dot(grad(u),grad(v))*dx + (c0*c0*k*k - nu*k)*dot(grad(u1),grad(v))*dx)
F += - inner(n, grad(v))*u*ds(1) + alpha/h*inner(u, v)*ds(1) + inner(n, grad(v))*u_L*ds(1) - alpha/h*inner(u_L, v)*ds(1)
F += - inner(n, grad(v))*u*ds(2) + alpha/h*inner(u, v)*ds(2) + inner(n, grad(v))*u_R*ds(2) - alpha/h*inner(u_R, v)*ds(2)
# Derivative of "F" w.r.t. parameters "m"
dF_m = assemble(derivative(F, u_R))
dF_m_mat = dF_m.array() # this matrix has zeros everywhere except for the first column
dF_m_vec = dF_m_mat[:,0] # this is the real (partial_F)/(partial_m) that I need
Many thanks!
You donât need to assemble \partial F / \partial m, the gradient update formula involves \partial F^* / \partial m\ \lambda.
So you first need to solve the adjoint equation, and then use the following (or something similar, I havenât checked precisely)
grad = assemble(-action(adjoint(derivative(F, u_R)), lmbda))
which is a vector and not a matrix.
Thank you for your quick response!
Firstly, let me clarify that I am solving an unsteady problem. Please consider the following discretised domain:
I define the control variable to be the time dependent Dirichlet BC at the right end. I then parametrise the control vector as the value of the BC at each time step: m = \left[u_{R_{1}}, u_{R_{2}}, \dots, u_{R_{N}}\right]^{T}, where N is the number of time steps. For simplicity, I want u_{R_{k}} to be constant across the boundary (i.e., the value of u at all the DOFs associated to \partial\Omega_{c} is the same), so that at each time step I just have one single parameter.
Therefore, as you say, I first solve the direct problem by forward integration in time, and then I solve the adjoint problem by backward integration in time. During the solution of the adjoint problem, I compute the gradient as you stated:
for k = Num\_steps : 0
\hspace{15pt} \left(\nabla J\right)_{k} = -\left(\partial F/\partial m_{k}\right)^{\ast} \cdot \lambda_{k}
end
If we focus on one single time step, I need to compute \left(\partial F/\partial m_{k}\right)^{\ast} \cdot \lambda_{k} as you suggest:
assemble(-action(adjoint(derivative(F, u_R)), lmbda))
But note that m_{k} is a scalar, so \left(\partial F/\partial m_{k}\right)^{\ast} is a vector and after all the time integration I will have a gradient vector whose size is equal to the number of time steps (which is consistent with my problem).
Consequently, my problem is that, in order to use derivative(F, u_R), I need to define u_R as a dolfin function, which is actually a vector defined at all the DOFs of the mesh. Then, derivative(F, u_R) is a matrix whose columns contain the partial derivative of F with respect to each u_R defined at every DOF of the mesh. However, I just need the column associated to the control boundary.
Is there a way of isolating the column of derivative(F, u_R) associated to the Boundary \partial\Omega_{c}?
Please, let me know if something is unclear. Also, I can supply the code if you want.
Thank you very much again!
You are stating that u_R is a constant, but you are interpolating it into a function space, allowing for variation that is non constant when differentiating. What happens if you let u_R = Constant(float(u_D[0]) and then in turn call
dF_m = assemble(derivative(F, u_R))
Thank you for your response.
That is what I tried to do in the begining, but when I do so, this is what I get:
File â/home/javierlm/Documentos/FEniCs/Optimisation_3D/Gradient.pyâ, line 104, in
dF_m = assemble(derivative(F, u_R))
File â/home/javierlm/anaconda3/envs/fenicsproject/lib/python3.9/site-packages/dolfin/fem/formmanipulations.pyâ, line 80, in derivative
raise RuntimeError(âComputing derivative of form w.r.t. â{}â. Supply Function as a Coefficientâ.format(u))
RuntimeError: Computing derivative of form w.r.t. âf_8â. Supply Function as a Coefficient
This error fixes when u_R is <class 'dolfin.function.function.Function'>.
You could use:
V_c = FunctionSpace(mesh, "Real", 0)
u_R = Function(V_c)
# ...
dF_m = assemble(derivative(F, u_R))
or as @Calva
suggested
lmbda = Function(V)
grad = assemble(-action(adjoint(derivative(F, u_R)), lmbda))
print(grad.get_local())
which would give you a single number
Thank you all for the help. @dokken approach worked for me.
Just another question, is there a similar function to FunctionSpace(mesh, "Real", 0) in dolfinx?
It is currently being worked on by @jackhale in Add Real Space by jhale ¡ Pull Request #2266 ¡ FEniCS/dolfinx ¡ GitHub
I see.
Is there an alternative in dolfinx to compute the sensitivity as in Sensitivity analysis with boundary conditions as parameters - #10 by dokken ?
Hi @dokken,
Iâve seen that the issue to create a Real function space (Add Real Space by jhale ¡ Pull Request #2266 ¡ FEniCS/dolfinx ¡ GitHub) has been closed. I wonder if there is available an alternative to get derivative(F, u_R) in dolfinx, where u_R is a Dirichlet BC. (And so that it is possible ti compute sensitivities with respect to boundary conditions as in Sensitivity analysis with boundary conditions as parameters - #7 by JLorenteMacias).