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',
degree=1)
f_2 = Expression('pow(x[0]-0.1,2)+pow(x[1]-0.3,2)<0.05*0.05 ? 0.1 : 0',
degree=1)
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
u_n.assign(u_)
J = assemble(inner(u_n3, u_n3)*dx)
dJdK = compute_gradient(J, control)
print(dJdK)
output of the above code is
f_4161