Zero control in time-distributed control of Heat eq. with multiple subdomains

dolfin-adjoint
1

Hello everyone,

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

  1. the region, where I want to put controll (Obstacle)
  2. 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

2

This is because your mesh function is empty (as your mesh is too coarse to resolve the area you are interested in).
You can see this by doing File("mf.pvd").write(domains) and open it in paraview.
You are also redefining variables throughout your script (like u0), so you need to be careful.

1 Like
3

You are right, thank you!
But only increasing of cells number hadn’t solved the problem.
I had to rewrite the weak formulation in a simpler way. It seems that rhs() and lhs() functions can’t handle some expressions like nuinner(grad(u), grad(v))(dx(0) + dx(1)).

Could you suggest some suitable way for this problem to set state (u) and control (f) constraints?