NotImplementedError when using 2D Expression as Control

dolfin-adjoint
1

Hi everyone,

I am relatively new to fenics and dolfin-adjoint and I am currently facing an error and just don’t see what I am doing wrong. :face_with_monocle: I am trying to write a simple optimization routine for a 2d advection equation. After reading the post how-to-use-expression-parameters-as-a-control I also added the derivative of that expression and the dependencies (similar to the time-dependent-wave demo just that I tried to avoid defining a class but rather use Expression() directly).

Any advice on what might be the problem would be very much appreciated :slightly_smiling_face:

from fenics import *
from fenics_adjoint import *
import matplotlib.pyplot as plt

def forward(beta, u0, u_exact, mesh, delta_t = 0.01, end_time=1, save_vtks=True):
    parameters["ghost_mode"] = "shared_facet"
    time_steps = int(end_time/delta_t)
    
    V_dg = FunctionSpace(mesh, "DG", 1)
    bcs = DirichletBC(V_dg, u_exact, "on_boundary", method='geometric')
    
    phi, v = TrialFunction(V_dg), TestFunction(V_dg)
    u_new, u_old = Function(V_dg, name='u_new'), Function(V_dg, name='u_old')
    
    f = Constant(0.0)
    n = FacetNormal(mesh)
    un = (dot(beta, n) + abs(dot(beta,n)))/2.0
    a = dot(grad(v), - beta*phi)*dx + dot(jump(v), un('+')*phi('+') - un('-')*phi('-'))*dS
    a += dot(v, u_exact*phi)*ds + Constant(1/delta_t)*phi*v*dx
    L = v*f*dx + Constant(1/delta_t)*u_old*v*dx
    
    project(u0, V_dg, function=u_old)
    t = 0
    for n in range(time_steps):
            t += delta_t
            u_exact.t = t

            solve(a==L, u_new, bcs)
            u_old.assign(u_new)
    
    return u_old


# define problem
m_size = 64
mesh = UnitSquareMesh(m_size, m_size)

beta_x = Expression("0.5", degree=1) 
beta_y = Expression("0.0", degree=1)
beta = Expression(("beta_x","beta_y"), degree=2, beta_x = beta_x, beta_y=beta_y)
beta.dependencies = [beta_x, beta_y]
beta.user_defined_derivatives = {beta_x: Expression("0",degree=1), beta_y : Expression("0",degree=1)}

u0 = Expression("sin(2*pi*x[0])", degree=1, t=0, name='u0')
u_exact = Expression("sin(2*pi*(x[0] - c*t))", degree=3, t=0, c=beta_x)

# solve foward system 
u_old = forward(beta, u0, u_exact, mesh)
plot(u_old)
plt.show()
# Gradient computation 
J = assemble(inner(u_old, u_old)*dx)
dJdbeta0 = compute_gradient(J, Control(beta_x))

However, If I run my code I get the warning

WARNING:root:Adjoint value is None, is the functional independent of the control variable?

and the error

NotImplementedError  Traceback (most recent call last)
<ipython-input-2-2da1fe06dc0d> in <module>
     56 # Gradient computation
     57 J = assemble(inner(u_old, u_old)*dx)
---> 58 dJdbeta0 = compute_gradient(J, Control(beta_x))

I use 2019.1.0 for fenics and dolfin-adjoint.

2

For this particular example, there is no reason to use Expression as the control variable. Your problem can be solved with:

m_size = 64
mesh = UnitSquareMesh(m_size, m_size)
beta_x = Constant(0.5)
beta_y = Constant(0)
beta = as_vector((beta_x, beta_y))

u0 = Expression("sin(2*pi*x[0])", degree=1, t=0, name='u0')
u_exact = Expression("sin(2*pi*(x[0] - c*t))", degree=3, t=0, c=beta_x)

# solve foward system 
u_old = forward(beta, u0, u_exact, mesh)
plot(u_old)
plt.show()
# Gradient computation 
J = assemble(inner(u_old, u_old)*dx)
dJdbeta0 = compute_gradient(J, Control(beta_x))

I would be really careful with wrapping C++ expressions inside other C++ expressions.

2 Likes
3

Hi Dokken,

thanks for your fast reply. I forgot to mention that I actually want beta_x and beta_y to be non-constant in future. I only set them constant to simplify the example for the beginning. Do you have any idea on how to get around this problem in that case?

4

Hi again,

I would have a second related question. Considering only the simple example where as dokken correctly pointed out :+1: that its enough to use constants beta_x = Constant(0.5) beta_y = Constant(0) beta = as_vector((beta_x, beta_y))

Does anyone perhaps know why the taylor_test

Jhat = ReducedFunctional(J, Control(beta_x))
h = Constant(0.1)
conv_rate = taylor_test(ReducedFunctional(J, Control(beta_x)), beta_x, h)

results in the error

RuntimeError: Solve variational problem. Illegal keyword argument ‘function’ value ‘u_old’.

5

Finally found what was causing the RuntimeError in taylor_test: I created the function u_old but initialize it with u0 after defining the variational problem a, L. If I move the project(u0, V_dg, function=u_old) before the variational definition, then I do not get any errors.