Plot compute_gradient

dolfin-adjoint
1

Hi everyone!

I’m trying to plot the gradient of a functional of kinetic energy of a Navier-Stokes’s solution but I get some errors.
My code is

from fenics import *
from mshr import *
from dolfin import *
from dolfin_adjoint import *


T = 90
num_steps = int(T)*25
dtn = T / num_steps

y_h= 11
rho = 1  # density
U0=1     #initial horizontal velocity
Re=70 #Reynold's number
D=1
mu=rho*U0*D/Re

channel = Rectangle(Point(-10.0, -y_h), Point(20.0, y_h))


cylinder = Circle(Point(0.0, 0.1), 0.5,35)
   
domain = channel - cylinder
mesh = generate_mesh(domain, 35)

bcu_inflow = DirichletBC(V, Expression(inflow_profile, degree=2), inflow)
bcu_walls = DirichletBC(V, Constant((0, 0)), walls)
bcu_cylinder = DirichletBC(V, Constant((0, 0)), cylinder)
bcp_outflow = DirichletBC(Q, Constant(0), outflow)
bcu = [bcu_inflow, bcu_walls, bcu_cylinder]
bcp = [bcp_outflow]

u= TrialFunction(V)
v= TestFunction(V)
#v, q= TestFunctions(W)
p=TrialFunction(Q)
q= TestFunction(Q)
u_n = Function(V)
u_ = Function(V)
p_n = Function(Q)
p_ = Function(Q)
# Define expressions used in variational forms
U = 0.5*(u_n + u)
n = FacetNormal(mesh)
f = Constant((0, 0))
k = Constant(dtn)
mu = Constant(mu)

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

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

F1 = rho*dot((u - u_n) / k, v)*dx \
    + rho*dot(dot(u_n, nabla_grad(u_n)), v)*dx \
    + inner(sigma(U, p_n), epsilon(v))*dx \
    + dot(p_n*n, v)*ds - dot(mu*nabla_grad(U)*n, v)*ds \
    - dot(f, v)*dx
a1 = lhs(F1)
L1 = rhs(F1)
# Define variational problem for step 2
a2 = dot(nabla_grad(p), nabla_grad(q))*dx
L2 = dot(nabla_grad(p_n), nabla_grad(q))*dx - (1/k)*div(u_)*q*dx
# Define variational problem for step 3
a3 = dot(u, v)*dx
L3 = dot(u_, v)*dx - k*dot(nabla_grad(p_ - p_n), v)*dx
# Assemble matrices
A1 = assemble(a1)
A2 = assemble(a2)
A3 = assemble(a3)
# Apply boundary conditions to matrices
[bc.apply(A1) for bc in bcu]
[bc.apply(A2) for bc in bcp]

for n_steps in range(500):
    t += dtn
 # Step 1: Tentative velocity step
    b1 = assemble(L1)
    [bc.apply(b1) for bc in bcu]
    solve(A1, u_.vector(), b1, 'bicgstab', 'hypre_amg')
    # Step 2: Pressure correction step
    b2 = assemble(L2)
    [bc.apply(b2) for bc in bcp]
    solve(A2, p_.vector(), b2, 'bicgstab', 'hypre_amg')
    # Step 3: Velocity correction step
    b3 = assemble(L3)
    solve(A3, u_.vector(), b3, 'cg', 'sor')
    # Update previous solution
    u_n.assign(u_)
    p_n.assign(p_)


J = assemble(inner(u_, u_)*dx)
dJdnu = compute_gradient(J,Control(mu), options={"riesz_representation": "L2"})
gra = File('gradient.pvd')
gra << dJdnu

The error is

RuntimeError: Unable to cast Python instance to C++ type (compile in debug mode for details)

If I try plot(dJdnu), I’ll get

RuntimeError: Expecting a mesh as keyword argument

I tried this without “L2 representation” with same result.

Can anyone suggest any ideas?
Thank you!

2

I think you need to import mshr to use generate_mesh.

from fenics import *
from mshr import *
3

I import mshr, but I forgot put in this code

4

Your code is not complete, as you have not defined V or Q, inflow_profile etc. Please make sure your code is executable.

5

One of the things you need to do to make your code work is:

  1. Wrap the generate_mesh(domain,35) as Mesh(generate_mesh(domain, 35)).

However, the key issue is that you are trying to export a constant value (mu is a constant) to pvd. You can simply print it as print(dJdnu.values()).

For further questions regarding this topic, please provide a minimal code that is actually executable after copy-pasting

1 Like
6

Thank you so much Dokken! This works for me. I have another question about dolfin-adjoint, but I think that is better if I post another topic.