Questions regarding output of compute_gradient()

Hi, I’m trying to compute a gradient over certain initial condition. I expected a Functional derivation which I can plot gradient over the mesh. However, the compute_gradient() only returns one single float value. Attached is a demonstration code. I found some old demo code uses Funtional() instead of assemble(), but Functional() is not defined in my fenics 2019.1.0 version.

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

T = 5.0            # final time
num_steps = 50    # number of time steps
dt = T / num_steps # time step size
eps = 0.01         # diffusion coefficient
K = 10.0           # reaction rate

# Read mesh from file
mesh = Mesh('navier_stokes_cylinder/cylinder.xml.gz')

# Define function space for velocity
W = VectorFunctionSpace(mesh, 'P', 2)
Q = FunctionSpace(mesh, 'P', 2)

# Define function space for system of concentrations
P1 = FiniteElement('P', triangle, 1)
element = MixedElement([P1, P1, P1])
V = FunctionSpace(mesh, element)

# Define test functions
v_1, v_2, v_3 = TestFunctions(V)
u_1, u_2, u_3 = TrialFunctions(V)

# Define functions for velocity and concentrations
w = Function(W)
u_n = Function(V)

# Split system functions to access components
u_n1, u_n2, u_n3 = split(u_n)

# Define source terms
f_1 = Expression('pow(x[0]-0.1,2)+pow(x[1]-0.1,2)<0.05*0.05 ? 0.1 : 0',
f_2 = Expression('pow(x[0]-0.1,2)+pow(x[1]-0.3,2)<0.05*0.05 ? 0.1 : 0',
f_3 = project(Expression('x[0]', degree=1), Q)
control = Control(f_3)

# Define expressions used in variational forms
k = Constant(dt)
K = Constant(K)
eps = Constant(eps)

# Define variational problem
F = ((u_1 - u_n1) / k)*v_1*dx + dot(w, grad(u_1))*v_1*dx \
  + eps*dot(grad(u_1), grad(v_1))*dx + K*u_1*u_n2*v_1*dx  \
  + ((u_2 - u_n2) / k)*v_2*dx + dot(w, grad(u_2))*v_2*dx \
  + eps*dot(grad(u_2), grad(v_2))*dx + K*u_n1*u_2*v_2*dx  \
  + ((u_3 - u_n3) / k)*v_3*dx + dot(w, grad(u_3))*v_3*dx \
  + eps*dot(grad(u_3), grad(v_3))*dx - K*u_n1*u_n2*v_3*dx + K*u_3*v_3*dx \
  - f_1*v_1*dx - f_2*v_2*dx - f_3*v_3*dx

a = lhs(F)
L = rhs(F)

# Create time series for reading velocity data
timeseries_w = TimeSeries('navier_stokes_cylinder/velocity_series')

# Time-stepping
t = 0
for n in range(num_steps):

    # Update current time
    t += dt

    # Read velocity from file
    timeseries_w.retrieve(w.vector(), t)

    # Solve variational problem for time step

    u_ = Function(V)
    solve(a == L, u_)

    # Update previous solution

J = assemble(inner(u_n3, u_n3)*dx)
dJdK = compute_gradient(J, control)

output of the above code is

The output is not a single float value but an internal identifier for some Function object I suppose. Since I’m missing the mesh and time-series files (and since the code seems to be really old) I can’t run your code to check it out.

What happens if you plot(dJdK)? What is the output of type(dJdK)?

Functional is deprecated, assemble is replacing that command to make dolfin-adjointmore flexible. Newer examples can be found at:

As @klunkean suggest, you should plot your function, or save it to file and visualize with paraview.

Thanks klunkean and @dokken. Sorry for those dumb questions, I’m still trying to get familiar with the software. I’m trying to solve a similar problem about get the gradient of an iterative process by mimic the above code. However, the following code did not work the same way. In between of iterations, I tried to have a relaxation factor, I wonder did that break the chain which causes problem for gradient computation.

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

mesh = UnitIntervalMesh(200)
n = FacetNormal(mesh)
def boundary(x, on_boundary):
    return on_boundary

P = FiniteElement('CG', interval, 2)
element = MixedElement([P, P, P])
T = FunctionSpace(mesh, element)

u_test, k_test, w_test = TestFunctions(T)
u, k, w = TrialFunctions(T)

T_n = Function(T)
u_n, k_n, w_n = split(T_n)
T_ = project(Expression(("0.0", "1e-7", "32000"), degree = 1), T)
u_, k_, w_ = split(T_)

gradP = Constant(10.0)
mu = Constant(1e-4)

V = FunctionSpace(mesh, 'CG', 2)
f = project(Expression("1.0", name='Control', degree=1), V)
control = Control(f)

bc_u = DirichletBC(T.sub(0), Constant(0.0), boundary)
bc_k = DirichletBC(T.sub(1), Constant(0.0), boundary)
bc_w = DirichletBC(T.sub(2), Constant(32000), boundary)
bc = [bc_u, bc_k, bc_w]

F1 = dot(grad(u), grad(u_test))*dx - gradP*u_test*dx
F2 = dot(dot(grad(u_), grad(u_)), k_test)*dx - mu*dot(grad(k), grad(k_test))*dx + mu*dot(grad(k), (k_test*n))*ds
F3 = f*w_test*dx + dot(dot(grad(u_), grad(u_)), w_test)*dx - mu*dot(grad(w), grad(w_test))*dx + mu*dot(grad(w), (w_test*n))*ds

F = F1 + F2 + F3

a = lhs(F)
L = rhs(F)

for i in range(10):

    solve(a == L, T_n, bc)
    T_.assign(project(Expression(("0.6*u_+0.4*u_n","0.6*k_+0.4*k_n","0.6*w_+0.4*w_n"), degree=1, u_n=T_n.sub(0), k_n=T_n.sub(1), w_n=T_n.sub(2), u_=T_.sub(0), k_=T_.sub(1), w_=T_.sub(2)), T))


J = assemble(inner(u_, u_)*dx)
dJdEpsilon = compute_gradient(J, control)

but I got the following error. I’m not too sure what is the difference makes this code not working but the demo works. The WARNING shows that the connection between the final output u_ and f was not identified.

WARNING:root:Adjoint value is None, is the functional independent of the control variable?
Traceback (most recent call last):
  File "/home/juntao/.vscode/extensions/ms-python.python-2019.8.30787/pythonFiles/", line 43, in <module>
  File "/home/juntao/.vscode/extensions/ms-python.python-2019.8.30787/pythonFiles/lib/python/ptvsd/", line 432, in main
  File "/home/juntao/.vscode/extensions/ms-python.python-2019.8.30787/pythonFiles/lib/python/ptvsd/", line 316, in run_file
    runpy.run_path(target, run_name='__main__')
  File "/home/juntao/anaconda3/envs/fenics/lib/python3.6/", line 263, in run_path
    pkg_name=pkg_name, script_name=fname)
  File "/home/juntao/anaconda3/envs/fenics/lib/python3.6/", line 96, in _run_module_code
    mod_name, mod_spec, pkg_name, script_name)
  File "/home/juntao/anaconda3/envs/fenics/lib/python3.6/", line 85, in _run_code
    exec(code, run_globals)
  File "/home/juntao/Desktop/FEnics_learning/", line 53, in <module>
    dJdEpsilon = compute_gradient(J, control)
  File "/home/juntao/anaconda3/envs/fenics/lib/python3.6/site-packages/pyadjoint/", line 31, in compute_gradient
    grads = [i.get_derivative(options=options) for i in m]
  File "/home/juntao/anaconda3/envs/fenics/lib/python3.6/site-packages/pyadjoint/", line 31, in <listcomp>
    grads = [i.get_derivative(options=options) for i in m]
  File "/home/juntao/anaconda3/envs/fenics/lib/python3.6/site-packages/pyadjoint/", line 51, in get_derivative
    return self.control._ad_convert_type(0., options=options)
  File "/home/juntao/anaconda3/envs/fenics/lib/python3.6/site-packages/pyadjoint/", line 46, in wrapper
    return function(*args, **kwargs)
  File "/home/juntao/anaconda3/envs/fenics/lib/python3.6/site-packages/fenics_adjoint/types/", line 125, in _ad_convert_type
    return compat.function_from_vector(self.function_space(), value, cls=Function)
  File "/home/juntao/anaconda3/envs/fenics/lib/python3.6/site-packages/fenics_adjoint/types/", line 206, in function_from_vector
    vector = vector._cpp_object
AttributeError: 'float' object has no attribute '_cpp_object'

Function.sub is not supported in the current version of dolfin-adjoint (use Function.split instead).
To differentiate through dolfin Expressions, you must provide a derivative Expression, look here for an example

You can avoid this hassle by switching to as_vector:

    tu_n, tk_n, tw_n = T_n.split()
    tu_, tk_, tw_ = T_.split()
    T_.assign(project(as_vector((0.6*tu_+0.4*tu_n,0.6*tk_+0.4*tk_n,0.6*tw_+0.4*tw_n)), T))

Note that your gradient will be 0 still because f has no effect on the functional.

Hi, Thanks for your supports. I am new to FEniCs am currently using this for a project. I am wondering can we get gradient if we use “Function.split”, I receive error message:

vector = vector._cpp_object
AttributeError: 'float' object has no attribute '_cpp_object'

Please supply a minimal example that reproduces your error message, as described in: