Minimeze function in Navier-Sokes's problem

One question that makes me wonder is why you have num_steps

when your loop is not related

Secondly, I would suggest you start using fewer time steps.

You are also using an unstable finite element pair for solving the fluid problem. Here V should be “CG”, 2.

Starting with a fixed inlet profile (from an initial velocity in the whole domain being zero) can also cause instabilities.

I am also uncertain about your choice of functional. You are trying to minimize the flow in the whole domain (L^2-norm of u given a fixed inlet condition), by changing the boundary condition on the object. Will this result in any reasonable physical solution.

You have also not specified at what point of the code it gets stuck (is it in forward solve, or is it during the optimization process)?

Consider the following code (with quite a few modifications for speedup) , including the tracking of the progress of the first forward solve (using tqdm, which can be installed with pip3 install tqdm --user), and writing the first forward solution to file. The code still struggles with converging to a good solution in the optimization process, which I believe is due to the choice of functional.

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

T = 1
num_steps = int(T)*50
dtn = T / num_steps

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


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


cylinder = Circle(Point(0.0, 0.1), 0.5,25)
   
domain = channel - cylinder
mesh = Mesh(generate_mesh(domain, 30))

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.45 && x[1]<0.65'

# Define boundary conditions
inflow_profile = ('0.95', '0.31')

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

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

#Define trial and test functions
u= TrialFunction(V)
v= TestFunction(V)
p=TrialFunction(Q)
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]

# Create solvers
solver1 = LUSolver(A1,'bicgstab', 'hypre_amg')
solver2 = LUSolver(A2, 'bicgstab', 'hypre_amg')
solver3 = LUSolver(A3, 'cg', 'sor')



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)
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]
    solver1.solve(u_.vector(), b1)
    # Step 2: Pressure correction step
    b2 = assemble(L2)
    [bc.apply(b2) for bc in bcp]
    solver2.solve(p_.vector(), b2)
    # Step 3: Velocity correction step
    b3 = assemble(L3)
    solver3.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)

outfile.close()
J = assemble(inner(u_, u_)*dx)
Jhat = ReducedFunctional(J, Control(g))
g_opt = minimize(Jhat, options={"disp": True})