# Initial condition and desired solution in Opt control of heat eq

Hello,

I want to set initial state condition and desired solution as constants in Opt. control of heat eq in 2D.
As I can see from this example:
desired behavior defined through expression for given data d.
Is there an easier way to set Initial condition and desiered solution as constants, than to define such function in data expression, that approximates their values?
And if there is, how should it look like for demo above?

For the data you can just use the Constant wrapper. For initial conditions, you should project a constant to the appropriate function so space.

Thank you, Mr. Dokken!

I tried this before and got odd result. But now knowing, this is the right way, I found out, that the problem was in getting optial solution of PDE (with optimal control as forcing term) and not in solving optimal control problem.
So I fixed some things, and the result got better. But there is still some issue.
Speaking about time-distributed control demo, I defined 2 as desired solution and 0 as initial cond.:

data = Constant(2)

def solve_heat(ctrls):

u_0 = interpolate(Constant(0), V)

I got this solution (in the middle of domain for each time step): As you can see, at the last time step u = 2.25. So it exceeds desired 2.0. The same with other values of desired solution. data = 4: u = 4.5, data = 8: u = 9.

It’s probably still something wrong with getting optimal solution of PDE. What I do for it, is just solving forward heat eq. with opt. control as forcing term.

Is there way to get optimal solution of PDE (at each time step) directly from dolfin adjoint, by which optimal control was found?

First of all, interpolation is not annotated in dolfin-adjoint, so I would use a projection.

It doesn’t solve the issue. Everything remains the same.

You need to add a minimal working example for anyone to help you. See Read before posting: How do I get my question answered? on how to set up an MWE

Code is almost identical to the demo:

from fenics import *
from collections import OrderedDict
import matplotlib.pyplot as plt

data = Constant(2)
nu = Constant(1e-5)

mesh = UnitSquareMesh(8, 8)
V = FunctionSpace(mesh, "CG", 1)

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):
u = TrialFunction(V)
v = TestFunction(V)
f = Function(V, name="source")
u_0 = project(Constant(0), V)
d = Function(V, name="data")
a, L = lhs(F), rhs(F)
bc = DirichletBC(V, 0, "on_boundary")
t = float(dt)
j = 0.5*float(dt)*assemble((u_0 - d)**2*dx)
while t <= T:
f.assign(ctrls[t])
data.t = t
d.assign(project(data, V))
solve(a == L, u_0, bc)
if t > T - float(dt):
weight = 0.5
else:
weight = 1
j += weight*float(dt)*assemble((u_0 - d)**2*dx)
t += float(dt)
return u_0, d, j

u, d, j = solve_heat(ctrls)

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})

# Then I solve forward problem

steps = 20
nx = ny = 8

def solve_heat_forward_cut(ctrls):
t = 0.0
global solution_graf
solution_graf = 
f = Function(V)
u = TrialFunction(V)
v = TestFunction(V)
u_0 = project(Constant(0), V)
a, L = lhs(F), rhs(F)
u = Function(V)
for n in range(steps-1):
f.assign(ctrls[n])
u = Function(V)
solve(a == L, u)
solution_graf.append(u(0.5,0.5))
u_0.assign(u)
t += dt

solve_heat_forward_cut(opt_ctrls)
print("Solution at the last time step = ", solution_graf[-1])

plt.plot(solution_graf)
plt.show()

So, as far as I can understand, you are solving the following problem:
Find the source f such that u=0 at t=dt, u=2 at t>dt.
As this is not a continuous heat profile in time, you cannot expect the optimal solution.
Other things that you could check is how many iterations does the optimization algorithm use (are you hitting max number of iterations?), and what is the value of the functional.

I haven’t done it yet, because I was actually more concerned about represinting solution of PDE. It seems to me, like the solution of the optimization problem might be right one, but issue arises when finding a solution of PDE with given optimal control.
So what I have done to represent the solution of PDE (function solve_heat_forward_cut), is this the right way? Is’nt there any more suitable/more direct way, which maybe was thought out by the developers?

You can re-evaluate your functional for any set of controls by rf(ctrls)
and get the solution at the last time step with:
u.block_variable.saved_output.
However, I do not get why you dont just call solve_heat(opt_ctrls) instead of creating a new version of it.
From a developer perspective, adding overloading of output is rather memory intensive, especially for time dependent problems, and the user must manually define how they would like to post-process the optimal solution. Dolfin-adjoint only saves the latest time-step, as you in your loop keep on re-assigning data to the variational problem.

Since your output figure is a point evaluation of a function, that is not something you use in your forward problem, there is no way for dolfin-adjoint to represent that operation.

1 Like

Thank you again, Mr. Dokken.
u.block_variable.saved_output is almost what I was looking for.
I tried solve_heat(opt_ctrls) before. It gives only the last timestep, meanwhile I’d like to investigate the solution at each time step. So thats why I wrote forward solver, full version of which saves not only point evaluation but full surface at each timestep.
Also I wasn’t quite shure, whether solve_heat(opt_ctrls) gives the right solution. Now it’s clear, cause u.block_variable.saved_output and solve_heat(opt_ctrls) give identical output.
Also solve_heat(opt_ctrls) can’t be applied diractly. opt_ctrls must have float range of indices between 0 and T:

Dict_Opt_ctrls = OrderedDict()
t = float(dt)
t_float = []
while t <= T:
t_float.append(t)
t += float(dt)
for i in range(len(t_float)):
Dict_Opt_ctrls[t_float[i]] = opt_ctrls[i]
solution = solve_heat(Dict_Opt_ctrls)

Dear Mr. Dokken,

So, as far as I can understand, you are solving the following problem:
Find the source f such that u=0 at t=dt, u=2 at t>dt.

line  u_0 = project(Constant(0), V) in the MWE above sets u=0 at t=dt. I can see why, I didn’t need it and changed it to how it was in demo u_0 = Function(V, name="solution").

1. But according to the demo: u = 0 \quad \textrm{for } \Omega \times \{0\} i.e. u=0 at t=0. I don’t undertand how and where we set that in the code.

2. How can I set another initial condition u = u_initial, t=0, where u_initial is a fenics solution of the given PDE. I.e. where do I have to put u_initial ?

Best regards

In the demo, for the first time step u(0)=u0, Which is
u0=Function(V). Dolfin default every function to have a zero value. To set another initial condition, you have to change u0.

But if I need to change u_0 to set initial condition why does u_0 = project(Constant(0),V) set u=0 at t=dt, but not u=0 at t=0?

And for another initial cond. u_initial should it be just like the following ?:
u_0 = Function(V)
u_0.assign(u_initial)

I do not follow your argumentation. Please supply a minimal example explaining what you do not get.
It is almost half a year since this post was active, and this is not very fresh in my memory.

Ok, sorry.

The MWE is basically identical to the demo, except one line u_0 = project(Constant(0), V)

So then MWE:

data = Constant(2)
nu = Constant(1e-5)

mesh = UnitSquareMesh(8, 8)
V = FunctionSpace(mesh, "CG", 1)

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):
u = TrialFunction(V)
v = TestFunction(V)
f = Function(V, name="source")
u_0 = project(Constant(0), V)
d = Function(V, name="data")
a, L = lhs(F), rhs(F)
bc = DirichletBC(V, 0, "on_boundary")
t = float(dt)
j = 0.5*float(dt)*assemble((u_0 - d)**2*dx)
while t <= T:
f.assign(ctrls[t])
data.t = t
d.assign(project(data, V))
solve(a == L, u_0, bc)
if t > T - float(dt):
weight = 0.5
else:
weight = 1
j += weight*float(dt)*assemble((u_0 - d)**2*dx)
t += float(dt)
return u_0, d, j

u, d, j = solve_heat(ctrls)

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})


You gave the following commentary:

Which I understood as u_0 = project(Constant(0), V) makes u=0 at t=dt and also makes discontinuous heat profile in time.

I see that u_0 sets initial condition at t=0 (not t=dt), and in case u_0 = Function(V) the initial condition is u=0 at t=0 and it doesn’t make discontinuous heat profile in time.

So I’m a little bit confused.

Let me try to rephrase.
Your forward model, starts with u = 0 at t=0,
which can either be achieved by u_0=Function(V) or u_0=project(Constant(0), V).
J = \int_{0}^T\int_{Omega}(u-2)^2\mathrm{d}x\mathrm{d}t + \text{regularization of source}.
As the equation you are solving is continuous in time, u(0.1) wil be quite close to zero, and far away from the value 2, that you would like to obtain.
In the demo that you started with, the data one wants to match is:

data = Expression("16*x*(x-1)*x*(x-1)*sin(pi*t)", t=0, degree=4)


which is zero at t = 0, and thus you will be able to get close to the analytical solution.
To assign a non-zero value as the initial guess (for instance 1), you can use:

    u_0 = project(Constant(1), V)

1 Like

Ok, got it.

Thank you for the explanation!