Figure 2 PDE constrained framework paper

Asked by John

Hi guys, I am trying to reproduce the simple heat equation inversion example mentioned in your framework paper. Unfortunately, I can not reproduce Fig 2c from the paper. I was simply using the following:

V = FunctionSpace(mesh, "CG", 1)
u = Function(V, name="State")
m = Function(V, name="Control")
v = TestFunction(V)

F = (inner(grad(u), grad(v)) - m*v)*dx
bc = DirichletBC(V, 0.0, "on_boundary")
solve(F == 0, u, bc)

u_desired = Expression("exp(-1/(1-x[0]*x[0])-1/(1-x[1]*x[1]))")

J = Functional((0.5*inner(u-u_desired, u-u_desired))*dx)

proj_udes = project(u_desired,V)

reduced_functional = ReducedFunctional(J, SteadyParameter(m))
m_opt = minimize(reduced_functional, method = "L-BFGS-B", bounds =(0.0,0.5),
                options = {"gtol": 1e-12,"ftol": 1e-12,"disp": True})

proj_error = project(u-u_desired,V)

save_file = File("optim_heat.pvd")
# desire solution
save_file << proj_udes
# L1 error
save_file << proj_error
# optimal forcing
save_file << m_opt

Am I doing something wrong?

Question information

Language:
English Edit question
Status:
Solved
For:
dolfin-adjoint Edit question
Assignee:
No assignee Edit question
Solved by:
John
Solved:
Last query:
Last reply:
Revision history for this message
Simon Funke (simon-funke) said :
#1

Hi John,

the problem is that the state variable (in your case "u") is not updated when calling "minimize".

You can fix that by solving the equations again after the optimisation finishes and before calculating the output files.

For example:

[...]
 # desired solution
proj_udes = project(u_desired,V)
plot(u_desired, mesh, interactive=True, title="Desired")

# L1 error
solve(F == 0, u, bc) # Need to solve the problem again to get the state variable for the optimal control (the control, m, is updated automatically)
proj_error = project(u-u_desired,V)
plot(proj_error, interactive=True, title="L1-Error")
[...]

A full working example can be found here:
https://gist.github.com/funsim/5769594

Hope this helps,

Simon

Revision history for this message
John (kaleyduck-at) said :
#2

Dear Simon,

Thanks for the solution. I have been trying to convert the stationary problem (above) in to a transient problem using the timestepping library as detailed in http://dolfin-adjoint.org/documentation/timestepping.html. Unfortunately I don't seem to find that in my installation. I have also downloaded the snapshot version and I don't seem to find it there either.

What would your advice be if I the real problem that I seek to optimise using L-BFGS problems is a time-dependent coupled PDE (pre-conditioned using ILU) i.e., what would be the best way to go about them using dolfin-adjoint?

A related question is can I just plugin ODEs as well as PDEs using the system.add_solve() method?

Best regards,
John

Revision history for this message
John (kaleyduck-at) said :
#3

Dear Simon,

Thanks for the solution. I have been trying to convert the stationary problem (above) in to a transient problem using the timestepping library as detailed in http://dolfin-adjoint.org/documentation/timestepping.html. Unfortunately I don't seem to find that in my installation. I have also downloaded the snapshot version and I don't seem to find it there either.

What would your advice be if I the real problem that I seek to optimise using L-BFGS method is a time-dependent coupled PDE (pre-conditioned using ILU) i.e., what would be the best way to go about them using dolfin-adjoint?

A related question is can I just plugin ODEs as well as PDEs using the system.add_solve() method?

Best regards,
John

Revision history for this message
Simon Funke (simon-funke) said :
#4

Hi John,

have a look at the application repository of dolfin-adjoint: https://code.launchpad.net/~libadjoint/dolfin-adjoint/applications
It contains a directory "optimal_control/bump_function_transient", which solves the transient version.

Patrick, can you comment on the ODE related question, please?

Simon

Revision history for this message
John (kaleyduck-at) said :
#5

Hi Simon,

Thanks for your quick reply. Does this mean that the time stepping module is not active anymore?

Best wishes,
John

Revision history for this message
Simon Funke (simon-funke) said :
#6

The time stepping module is still actively developed, but I have not used it yet. James is the main developer of it and can hopefully help you with your question.

Revision history for this message
Patrick Farrell (pefarrell) said :
#7

Dear John,

The timestepping module is separate functionality that has only just been added to dolfin-adjoint (but not yet included in a released version). It's a different way of writing your forward model to dolfin that captures more structure. If you have already implemented your forward model in dolfin and are happy with it, you may not need it.

As Simon said, there are several examples of optimisation of transient PDEs in the tests/ directory of dolfin-adjoint (these involve writing models in the usual dolfin style).

I didn't see the ODE question; we'll discuss that separately in the other thread.

Cheerio

Patrick

Revision history for this message
James Maddison (jamesmadd) said :
#8

John,

The timestepping code is an experimental module that can often lead to performance improvements, particularly if the time discretisation is explicit. It's primarily intended for cases with complex time discretisations, or in order to maximise efficiency. It is much less thoroughly tested than the FEniCS system or the main dolfin-adjoint code, and is currently available only in source form.

However this module certainly isn't a requirement for transient modelling with FEniCS and dolfin-adjoint -- this can be achieved via approaches as described in the FEniCS documentation, and there are many examples of the standard approach both within DOLFIN demos and dolfin-adjoint tests. If you're looking at reproducing or extending previous results then I'd recommend that this would be the best way to start.

Revision history for this message
John (kaleyduck-at) said :
#9

Dear Simon, Patrick and James,

Thanks for explaining me the difference between the timestepping module and the usual manual Rothe type discretisation that is generally done for transient problems. From a more object oriented point of view, it would be really useful for the users to have the time-stepping library that has a communication interface to the adjoint code -- this becomes quite relevant for 3D problems where people (like me) try to couple reaction-diffusion type equations to fluid-flow equations and optimise them.

Perhaps, it might be useful for me to wait a little bit longer until the unit-tests on the timestepping code is mature enough and then port everything to using that instead of the manual discretisation.

Thanks a bunch to three of you for sparing time to answer my questions.

Best wishes,
John