Navier-Stokes bc optimal conrtol

dolfin-adjoint
1

I’m tying to perform optimal control on the Dirichlet BC g. I want to minimize the recirculation area, the functional here implemented is not performing this task, it is only a tryout. The problem is that after the forward problem is solved and the minimization is called nothing appens and the program is then killed after some time. Could you suggest where the problem is?
Thank you.

from fenics import *
from mshr import *
from dolfin import *
import numpy as np
from dolfin_adjoint import *
from matplotlib import pyplot, rc

set_log_level(LogLevel.ERROR)


T = 2
dtn = 0.0005
num_steps = int(T/dtn)




mu = 1E-3 #dynamic viscosity
rho = 1  # density
Re=70 #Reynold's number


mesh = Mesh("our_mesh.xml")

cd = MeshFunction("size_t",mesh,"our_mesh_physical_region.xml")
fd = MeshFunction("size_t",mesh,"our_mesh_facet_region.xml")

         



Q = FunctionSpace(mesh, 'CG', 1)
V = VectorFunctionSpace(mesh, 'CG', 2)



g = Function(V, name="Control")

bot = mesh.coordinates().min(axis=0)[1]+0.1
top = mesh.coordinates().max(axis=0)[1]

# width of inlet channel
H = top - bot

# # 
U_in = 1.5

inflow_profile = Expression(('-4*U_in*(x[1]-bot)*(x[1]-top)/H/H' , '0'), bot=bot, top=top, H=H, U_in=U_in, degree=2)


# # Define boundary conditions
# # Inlet velocity
bcu_inlet = DirichletBC(V, inflow_profile,fd,1) 
# # No slip
bcu_wall1 = DirichletBC(V, Constant((0., 0.)),fd, 3)
bcu_wall2 = DirichletBC(V, Constant((0., 0.)),fd, 4)
bcu_wall3 = DirichletBC(V, Constant((0., 0.)),fd, 5)
    
# # Fixing outflow pressure
bcp_outflow = DirichletBC(Q, Constant(0.),fd, 2)  

bcu_jet = DirichletBC(V, g, fd ,  6)

bcu = [bcu_inlet,bcu_wall1,bcu_wall2,bcu_wall3,bcu_jet]

bcp = [bcp_outflow]

# #Define trial and test functions
u,v = TrialFunction(V), TestFunction(V)
p,q =TrialFunction(Q), TestFunction(Q)
# # Define functions for solutions at previous and current time steps
u_n = Function(V)
u_ = Function(V)
p_n = Function(Q)
p_ = Function(Q)
u_aux = Function(V)
# Define expressions used in variational forms
U = 0.5*(u_n + u)
n = FacetNormal(mesh)
f = Constant((0, 0))
ix= Constant((1,0))
k = Constant(dtn)
mu = Constant(mu)

# Define symmetric gradient
def epsilon(u):
    return sym(nabla_grad(u))

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


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]
[bc.apply(A3) for bc in bcu]

As = [A1,A2,A3]

# Create solvers
solvers = list(map(lambda x: LUSolver(), range(3)))


# Set matrices
for s, A in zip(solvers, As):
    s.set_operator(A)


import tqdm
timer = tqdm.tqdm("Solving forward problem", total=num_steps)
outfile = XDMFFile("u.xdmf")
outfile.write_checkpoint(u_, "u", 0, append=False)
outfile.write_checkpoint(p_, "p", 0, append=True)

Jtemp = assemble(inner(u_,u_)*dx)
Jlist = [Jtemp]

for n_steps in range(num_steps):
    timer.update(1)
    # Step 1: Tentative velocity step
    b1 = assemble(L1)
    [bc.apply(b1) for bc in bcu]
    solvers[0].solve(u_.vector(), b1)
    # Step 2: Pressure correction step
    b2 = assemble(L2)
    [bc.apply(b2) for bc in bcp]
    solvers[1].solve(p_.vector(), b2)
    # Step 3: Velocity correction step
    b3 = assemble(L3)
    solvers[2].solve(u_.vector(), b3)
    # Plot solution
    # Update previous solution
    u_n.assign(u_)
    p_n.assign(p_)
    outfile.write_checkpoint(u_, "u", n_steps+1, append=True)
    outfile.write_checkpoint(p_, "p", n_steps+1, append=True)

    Jtemp = assemble(inner(u_,u_)*dx)
    Jlist.append(Jtemp)


outfile.close()


print(u_.ufl_shape)


from ufl import conditional


J = 0 
for i in range(1, len(Jlist)):
    J += 0.5*(Jlist[i-1] + Jlist[i])*float(dtn)

# J = assemble(inner(conditional(u_[0] <  Constant((0.0)), u_[0], Constant((0.0))), conditional(u_[0] < Constant((0.0)), u_[0], Constant((0.0))))*dx)

# J = assemble(0.5*inner(curl(u_),curl(u_))*dx)
m = Control(g)
Jhat = ReducedFunctional(J,m)
g_opt = minimize(Jhat, options={"disp": True})
2

Without the mesh and associated markers, it is not possible to reproduce the error, and it is unlikely that anyone can give you any guidance.