Determining time dependent source function

Hi,
I am trying to determine time dependent source function and the solution. But, when I plot the solution it is irrelevant. Error also not small enough. And is there any way for ploting the source function?
Can you please check the code…
Thanks in advance

from fenics import *

from fenics_adjoint import *
from collections import OrderedDict
data = Expression("16*x[0]*(x[0]-1)*x[1]*(x[1]-1)*sin(pi*t)", t=0, degree=4)
mesh = UnitSquareMesh(8, 8)
V = FunctionSpace(mesh, "CG", 1)
nu= interpolate(Expression("1e-5",degree=2),V)
dt = Constant(0.1)
T = 2
ctrls = OrderedDict()
t = float(dt)
while t <= T:
    ctrls[t] = Function(V)
    t += float(dt)
def solve_heat(ctrls,annotate=False):
    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 + inner(nu*grad(u), grad(v)) - f*v)*dx
    a, L = lhs(F), rhs(F)
    bc = DirichletBC(V, 0, "on_boundary")
    t = float(dt)
    j = assemble((u_0 - d)**2*dx)
    while t <= T:
        # Update source term from control array
        f.assign(ctrls[t])
        # Update data function
        data.t = t
        d.assign(interpolate(data, V))
        # Solve PDE
        solve(a == L, u_0, bc)
        j=assemble((u_0 - d)**2*dx)
        # Update time
        t += float(dt)
    return u_0, d, j
u, d, j = solve_heat(ctrls,annotate=True)
alpha = Constant(1e-1)
regularisation = alpha/2*sum([1/dt*(fb-fa)**2*dx 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})
error_u= errornorm(u,d)
print ("error in state function   :", error_u)
from matplotlib import pyplot
x = [c((0.5, 0.5)) for c in opt_ctrls]
pyplot.plot(x)
pyplot.ylim([-3, 3])

Hi,

At least one source of error is that you do not sum up the functional contributions on each time step.
I guess you should do j+=assemble((u_0-d)**2*dx) if you are interested in a functional that is not only dependent on the last time step.

To plot your controls, you can either plot each source function using matplotlib:

for c in opt_ctrls:
    plt.figure()
    plot(c)
    plt.show()

Alternatively, since you have many time-steps, I would consider storing the solution to pvd or xdmf-format, visualizing in paraview:

opt_f = File("opt_f.pvd")
out_func = Function(V)
for c in opt_ctrls:
    out_func.assign(c)
    opt_f << out_func 

There are two mistakes in your error comparison:

  1. You do not compare with the optimal solution, but with the initial solution.
  2. Your alpha parameter is quite high, if you loosen up on your regularization, the error will decrease.
    I suggest the following
# Creating a variable on the tape for the last solution of the PDE
from pyadjoint.placeholder import Placeholder
u_last = Placeholder(u)

alpha = Constant(1e-2)
regularisation = alpha/2*sum([1/dt*(fb-fa)**2*dx 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})
# Evaluate functional at optimum
rf(opt_ctrls)
# Error of u at last time-step
error_u= errornorm(u_last.saved_output,d)