Gradient of Drag force

dolfin-adjoint
1

Hi everyone!
I’m trying to calculate the gradient of the drag force respect the external force f, for de Navier-Stokes problem
\rho(\partial _{t}u + u\cdot \nabla u) = \nabla\cdot \sigma(p,u)+f
where \sigma(p,u)=pId + \mu\frac{1}{2}(\nabla u + \nabla u^{T})
The drag force is defined by
D=\int_{\Gamma}<\sigma(p,u)n,e_{x}>d\Gamma
Using dolfin-adjoint, I want to calculate \frac{dD}{df}.
This is my code:

from fenics import *
from mshr import *
import numpy as np
from dolfin import *
from dolfin_adjoint import *
import matplotlib.pyplot as plt
from petsc4py import PETSc
import math

#Steady Navier-Stokes

mu = 0.009 #dynamic viscosity
y_h= 11
rho = 1  # density
U0=1     #initial horizontal velocity
L=y_h*2   # height
Re=40 #Reynold's number
D=1
mu=rho*U0*D/Re
Rer=round(Re,2)

# Create mesh
channel = Rectangle(Point(-10.0, -y_h), Point(20.0, y_h))

nr=[]


cylinder = Circle(Point(0.0, 0.0), 0.5,150)
   
domain = channel - cylinder
mesh = generate_mesh(domain, 100)
cell_markers = MeshFunction("bool", mesh, mesh.topology().dim())
 
for c in cells(mesh):
    if c.midpoint().y() < 3 and c.midpoint().y() > -3 and c.midpoint().x() < 12 and c.midpoint().x() > -3:
                cell_markers[c] = True
    else:
                cell_markers[c] = False
    
for i in range (1):
    mesh = refine(mesh, cell_markers)
    nr.append(i)
    

class BoundaryCy(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and x[0]>-0.55 and x[0]<0.55 and x[1]>-0.55 and x[1]<0.55


boundary_markers = MeshFunction("size_t", mesh, mesh.topology().dim()-1)  

cy = BoundaryCy()
cy.mark(boundary_markers, 0)

dsc=Measure('ds', domain=mesh, subdomain_data=boundary_markers)   


# Define function spaces

#Product of Function Spaces
P2 = VectorElement("CG", mesh.ufl_cell(), 2)
P1 = FiniteElement("CG", mesh.ufl_cell(), 1)
TH = MixedElement([P2, P1])
W = FunctionSpace(mesh, TH)
Q = FunctionSpace(mesh, 'CG', 1)
V = VectorFunctionSpace(mesh, 'CG', 2)



# Define boundaries
inflow= 'near(x[0], -10)'
outflow= 'near(x[0], 20)'
walls= 'near(x[1], '+str(-y_h)+') || near(x[1], '+str(y_h)+' ) '
cylinder= 'on_boundary && x[0]>-0.55 && x[0]<0.55 && x[1]>-0.55 && x[1]<0.55'

# Define boundary conditions
inflow_profile = ('1', '0')

bcu_inflowp = DirichletBC(W.sub(0), Expression(inflow_profile, degree=2), inflow)
bcu_wallsp = DirichletBC(W.sub(0), Constant((0, 0)), walls)
bcu_cylinderp = DirichletBC(W.sub(0), Constant((0, 0)), cylinder)
bcp_outflowp = DirichletBC(W.sub(1), Constant(0), outflow)
bcs = [bcu_inflowp, bcu_wallsp, bcu_cylinderp, bcp_outflowp]

w = Function(W)
u,p = split(w) 
v, q=TestFunctions(W)
f = Function(V)
ix = Constant((1,0))
class bcylinder():
    def inside(self, x, on_boundary):
        return on_boundary and x[0]>-0.55 and x[0]<0.55 and x[1]>-0.55 and x[1]<0.55


bcy=bcylinder()      

def epsilon(u):
    return sym(nabla_grad(u))

def sigma(u, p):
    return 2*mu*epsilon(u) - p*Identity(len(u))

n = FacetNormal(mesh)

a= rho*dot(dot(u, nabla_grad(u)),v)*dx - dot(f,v)*dx+ inner(sigma(u,p),epsilon(v))*dx+ dot(p*n, v)*ds - dot(mu*nabla_grad(u)*n, v)*ds + q*div(u)*dx 


solve(a == 0, w, bcs)

J=assemble(inner(dot(sigma(u,p),n),ix)*dsc)
dJdf=compute_gradient(J,Control(f))
plot(dJdf)

The solution with this code is not the real solution.
Anyone has some idea?
Thank you so much!