Minimeze function in Navier-Sokes's problem

dolfin-adjoint
1

Hi everyone!
I’m trying to use minimeze function of dolfin adjoint to solve this problem
\min_{g} J=\int_{\Omega}u(T)²d\Omega, where T is the flast step and g is a Dirichlet boundary condition, and u is the solution of Navier-Stokes equation.
I idid this example first Dirichlet BC control of the Stokes equations and does it. But when I try in this problem doesn’t work.
My code is (I use the Chorin’s method)

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

T = 90
num_steps = int(T)*25
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 = generate_mesh(domain, 30)

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

# 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]

for n_steps in range(600):

    # 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')
    # Plot solution
    plot(u_, title='Velocity')
    # Update previous solution
    u_n.assign(u_)
    p_n.assign(p_)


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

When I try to do this, I get “Solving linear variational problem” infinite times for hours.
Any suggestion?
Thank you!

2

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})
3

Thank you so much @dokken. This code works for me. Now I’m triying to minimize the gradient because I want it to be as uniform as possible.

J=assemble(inner(grad(u),grad(u))*dx)

I think it is a better problem and its works.
I’m doing the next problem too

J=assemble(inner(dot(sigma(u,p),n),ix)*ds)
g_opt=minimize(J,Control(g))

Minimize the the drag force around the cylinder. But I obtain the next error.

*** Error:   Unable to solve nonlinear system with NewtonSolver.
*** Reason:  Newton solver did not converge because maximum number of iterations reached.

Do you know what is the problem?
Thank you!

4

You Need to create a minimal code example that reproduces the error message for me to be of any help

5 
from fenics import *
from mshr import *
import numpy as np
from dolfin import *
from dolfin_adjoint import *
import matplotlib.pyplot as plt

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

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

cylinder = Circle(Point(0.0, 0.0), 0.5,100)
   
domain = channel - cylinder
mesh = generate_mesh(domain, 60)

#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)

g = Function(V)


# 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), g, 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 = Constant((0, 0))
ix = Constant((1,0))

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 + inner(sigma(u,p),epsilon(v))*dx+ dot(p*n, v)*ds - dot(mu*nabla_grad(u)*n, v)*ds + q*div(u)*dx 
L= 0

solve(a == L, w, bcs)

J=assemble(inner(dot(sigma(u,p),n),ix)*ds)
Jhat = ReducedFunctional(J,Control(g))
g_opt = minimize(Jhat)

*** Error:   Unable to solve nonlinear system with NewtonSolver.
*** Reason:  Newton solver did not converge because maximum number of iterations reached.