Mixed-dimensional branch - Solving optimal control problem

1

Dear all,

I am trying to solve the following optimal control problem using Lagrange function

\min\limits_{m} \frac{1}{2}{ \parallel u-u_{d} \parallel}^{2}_{L^{2}(\Omega)} + \frac{\alpha}{2}{ \parallel m\parallel}^{2}_{L^{2}(\Gamma)}

subject to

- \Delta u + u^{3} - u = 0 in \Omega
\frac{ \partial u}{\partial \nu} = m on \Gamma

here the control variable m is defined only on the boundary, and therefore I am using a mixed dimensional branch.

However, when running it via terminal I encounter the following error

[problem.py] size F = 3
[problem] create list of residual forms OK
[problem.py] size J = 9
[problem] create list of jacobian forms OK, J_list size = 9
[MixedNonlinearVariationalProblem]::check_forms - To be implemented
Build mapping between 0 and 7
Segmentation fault (core dumped)

It seems to me that there is something wrong with the part where I take derivatives.

I tried to use extract_blocks as suggested by @cdaversin here. However, I think, in my case, I don’t have enough Arguments before differentiation (in the weak formulation all functions are basically trial functions and test functions appear only later as a direction in which I take derivatives)

I would appreciate any suggestions on how to resolve the issue.

Here is a MWE

from dolfin import *
mesh = RectangleMesh(Point(0.0, 0.0), Point(1.0, 2.0), 10, 10)

class Boundary(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary
    
mf = MeshFunction("size_t", mesh, mesh.geometry().dim()-1, 0)
Boundary().mark(mf,1)
boundary_mesh = MeshView.create(mf, 1)

U = FunctionSpace(mesh, "CG", 1) # forward solution
V = FunctionSpace(mesh, "CG", 1) # adjoint solution
M = FunctionSpace(boundary_mesh, "CG", 1) # control
Z = MixedFunctionSpace(U, V, M)

dx = Measure("dx", domain=Z.sub_space(0).mesh())
dBoundary = Measure("dx", domain=Z.sub_space(2).mesh())

z = Function(Z)
u, uT, m = z.sub(0),z.sub(1),z.sub(2)

F = inner(grad(u), grad(uT))*dx - m*uT*dBoundary + inner(u**3, uT)*dx - inner(u, uT)*dx
ud = Constant(3)
J = 0.5*inner(u-ud, u-ud)*dx + 0.5*Constant(1e-7)*inner(m, m)*dBoundary
bcs = []

# Derive optimality conditions
L = J + F

kkt = []
for zi in z.split():
        d = derivative(L, zi)
        kkt.append(d)
        
Hes = [] 
for Js_i in kkt:
    for zi in z.split():
        d = derivative(Js_i, zi)
        Hes.append(d)

problem = MixedNonlinearVariationalProblem(kkt, z._functions, bcs, Hes)
solver = MixedNonlinearVariationalSolver(problem)
solver.solve()
2

Hello,

From the error you mention it seems to me that MixedNonlinearVariationalProblem is trying to build a mapping between the meshes of index 0 (the rectangle) and 7 (the boundary). The cells of the boundary mesh are facets of the rectangle (the mapping already exists), but the cells of the rectangle are not entities of the boundary mesh (and that’s why you get the segfault).

It might be that you have a form involving an argument defined only on the boundary that you are trying to integrate over the rectangle domain. So I would double check the variational formulation given to MixedNonlinearVariationalProblem.

Especially, you say that test functions appear when you compute the derivatives, but I cannot see any TestFunctions(Z) ?

3

Thank you for your reply.

I have double-checked the weak formulation but still cannot find any issue there. As you can see only m is defined over the boundary and those terms are integrated on Z.sub_space(2) which is the boundary_mesh.

Were you able to reproduce the error with the MWE I provided? Maybe I am missing something obvious in the weak formulation there.

In the meantime, I have pulled the update from your container and now the error message is more explicit: RuntimeError: Cannot find common parent mesh

I meant that they implicitly appear in derivative(form, u, du). (As far as I know, one does not have to necessarily specify them as third parameter du and then they’re taken from the same space as u )