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: http://www.dolfin-adjoint.org/en/latest/documentation/time-distributed-control/time-distributed-control.html
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?
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?
from fenics import *
from fenics_adjoint 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")
F = ( (u - u_0)/dt*v + nu*inner(grad(u), grad(v)) - f*v)*dx
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 = [0]
f = Function(V)
u = TrialFunction(V)
v = TestFunction(V)
u_0 = project(Constant(0), V)
F = ( (u - u_0)/dt*v + nu*inner(grad(u), grad(v)) - f*v)*dx
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.
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)[0]
Could you please answer couple of my questions to the following:
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").
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.
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 ?
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.
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.
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")
F = ( (u - u_0)/dt*v + nu*inner(grad(u), grad(v)) - f*v)*dx
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.
But from tadays commentary:
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.
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).
Your functional is: 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[0]*(x[0]-1)*x[1]*(x[1]-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:
