Hello everyone,

I’m trying to find time-distributed control for Heat eq. with multiple subdomains.

The subdomains represent:

- the region, where I want to put controll (Obstacle)
- boundaries with NeumannBC (Top, Bottom, Left and Right)

I used this tutorial:

https://bitbucket.org/dolfin-adjoint/pyadjoint/src/ma…

But got no control. (it’s zero)

Here is my code:

`import numpy as np from fenics import * from dolfin import * from fenics_adjoint import * from collections import OrderedDict 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) class Bottom(SubDomain): def inside(self, x, on_boundary): return near(x[1], 0.0) class Top(SubDomain): def inside(self, x, on_boundary): return near(x[1], 1.0) class Obstacle(SubDomain): def inside(self, x, on_boundary): return (between(x[1], (0.4, 0.6)) and between(x[0], (0.4, 0.6))) left = Left() top = Top() right = Right() bottom = Bottom() obstacle = Obstacle() nx = 8 ny = 8 mesh = UnitSquareMesh(nx, ny) t = 0.0 Time = T = 1 num_steps = 20 dt = Time / num_steps domains = MeshFunction("size_t", mesh, mesh.topology().dim()) domains.set_all(0) obstacle.mark(domains, 1) boundaries = MeshFunction("size_t", mesh, mesh.topology().dim() - 1) boundaries.set_all(0) left.mark(boundaries, 1) top.mark(boundaries, 2) right.mark(boundaries, 3) bottom.mark(boundaries, 4) #for simplicity all constant are equal a0 = Constant(1.0) a1 = Constant(1.0) g1 = Constant("0.01") g2 = Constant("0.01") g3 = Constant("0.01") g4 = Constant("0.01") u_d = Expression('sin(pi*x[1])*exp( -pi*pi*t )*pow(sin(pi*x[0]) , 3)', degree=2, t=0) u_0 = interpolate(u_d, V) ds = ds(subdomain_data = boundaries) dx = dx(subdomain_data = domains) ctrls = OrderedDict() t = float(dt) while t <= Time: ctrls[t] = Function(V) t += float(dt) def solve_heat(ctrls): u = TrialFunction(V) v = TestFunction(V) f = Function(V, name="source") u_0 = Function(V, name="solution") d = Function(V, name="data") F = ( (u - u_0)/dt*v*(dx(0) + dx(1)) + inner(a0*grad(u), grad(v))*dx(0) + inner(a1*grad(u), grad(v))*dx(1) - f*v*dx(1) - g1*v*ds(1) - g2*v*ds(2) - g3*v*ds(3) - g4*v*ds(4) ) a, L = lhs(F), rhs(F) t = float(dt) j = 0.5*float(dt)*assemble((u_0 - d)**2*(dx(0) + dx(1))) while t <= T: # Update source term from control array f.assign(ctrls[t]) # Update data function u_d.t = t d.assign(interpolate(u_d, V)) # Solve PDE solve(a == L, u_0) # Implement a trapezoidal rule if t > T - float(dt): weight = 0.5 else: weight = 1 j += weight*float(dt)*assemble((u_0 - d)**2*(dx(0) + dx(1))) # Update time t += float(dt) return u_0, d, j u, d, j = solve_heat(ctrls) alpha = Constant(1e-2) regularisation = alpha/2*sum([1/dt*(fb-fa)**2*(dx(0) + dx(1)) for fb, fa in zip(list(ctrls.values())[1:], list(ctrls.values())[:-1])]) J = j + assemble(regularisation) m = [Control(c) for c in ctrls.values()] rf = ReducedFunctional(J, m) opt_ctrls = minimize(rf, options = {"maxiter": 50}) x = [c((0.5, 0.5)) for c in opt_ctrls] plt.plot(x) plt.show()`

Hope someone can help.

Best regards

Ruslan