Sensitivity analysis with boundary conditions as parameters

[dokken]

@JLorenteMacias

# 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

You could use:

V_c = FunctionSpace(mesh, "Real", 0)
u_R = Function(V_c)
# ...
dF_m = assemble(derivative(F, u_R))

or as @Calva

@Calva

grad = assemble(-action(adjoint(derivative(F, u_R)), lmbda))

suggested

lmbda = Function(V)
grad = assemble(-action(adjoint(derivative(F, u_R)), lmbda))
print(grad.get_local())

which would give you a single number