Navier Stokes and Convection Diffusion

variational formulation
1

Dear all,

I am trying to solve a convection diffusion equation in the velocity field that is computed with incompressible Navier Stokes.
For this, I tried to use the FEniCS Navier Stokes cylinder tutorial and I added a convection diffusion step to it.
However, upon execution of the code I get the following error:

TypeError                                 Traceback (most recent call last)
<ipython-input-47-7430534ce1df> in <module>()
    159     b4 = assemble(L4)
    160     #[bc.apply(b4) for bc in bcs]
--> 161     solve(A4,s_,b4)
    162 
    163     # Plot solution

1 frames
/usr/local/lib/python3.7/dist-packages/dolfin/la/solver.py in solve(A, x, b, method, preconditioner)
     70     """
     71 
---> 72     return cpp.la.solve(A, x, b, method, preconditioner)

TypeError: solve(): incompatible function arguments. The following argument types are supported:
    1. (A: dolfin.cpp.la.GenericLinearOperator, x: dolfin.cpp.la.GenericVector, b: dolfin.cpp.la.GenericVector, method: str = 'lu', preconditioner: str = 'none') -> int

Invoked with: <dolfin.cpp.la.Matrix object at 0x7f193eb62e90>, Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 18166), FiniteElement('Lagrange', triangle, 1)), 18206), <dolfin.cpp.la.Vector object at 0x7f1969e922f0>, 'default', 'default

The code looks as follows:

from __future__ import print_function
from fenics import *
from mshr import *
import numpy as np
import matplotlib.pyplot as plt

T = 0.1            # final time
num_steps = 100   # number of time steps
dt = T / num_steps # time step size
mu = 0.001         # dynamic viscosity
rho = 1            # density

# Create mesh
channel = Rectangle(Point(0, 0), Point(2.2, 0.41))
cylinder = Circle(Point(0.2, 0.2), 0.05)
domain = channel - cylinder
mesh = generate_mesh(domain, 64)

# Define function spaces
V = VectorFunctionSpace(mesh, 'P', 2)
Q = FunctionSpace(mesh, 'P', 1)
S = FunctionSpace(mesh,'P',1)

# Define boundaries
inflow   = 'near(x[0], 0)'
outflow  = 'near(x[0], 2.2)'
walls    = 'near(x[1], 0) || near(x[1], 0.41)'
cylinder = 'on_boundary && x[0]>0.1 && x[0]<0.3 && x[1]>0.1 && x[1]<0.3'

# Define inflow profile
inflow_profile = ('4.0*1.5*x[1]*(0.41 - x[1]) / pow(0.41, 2)', '0')

# 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, Constant((0, 0)), cylinder)
bcp_outflow = DirichletBC(Q, Constant(0), outflow)
bcs_inflow = DirichletBC(S, Constant(1),inflow)
bcu = [bcu_inflow, bcu_walls, bcu_cylinder]
bcp = [bcp_outflow]
bcs = [bcs_inflow]

# Define trial and test functions
u = TrialFunction(V)
v = TestFunction(V)
p = TrialFunction(Q)
q = TestFunction(Q)
s = TrialFunction(S)
z = TestFunction(S)

# Define functions for solutions at previous and current time steps
u_n = Function(V)
u_  = Function(V)
p_n = Function(Q)
p_  = Function(Q)
s_n = Function(S)
s_ = Function(S)

# Define expressions used in variational forms
U  = 0.5*(u_n + u)
n  = FacetNormal(mesh)
f  = Constant((0, 0))
fs  = Constant(0)
k  = Constant(dt)
mu = Constant(mu)
rho = Constant(rho)

# 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

# Define variational problem for step 4
F4 = rho*((s - s_n) / k) *z*dx \
   + rho*dot(u_n, grad(s))*z*dx \
   - fs*z*dx
a4 = lhs(F4)
L4 = rhs(F4)

# Assemble matrices
A1 = assemble(a1)
A2 = assemble(a2)
A3 = assemble(a3)
A4 = assemble(a4)

# Apply boundary conditions to matrices
[bc.apply(A1) for bc in bcu]
[bc.apply(A2) for bc in bcp]
[bc.apply(A4) for bc in bcs]

# Create XDMF files for visualization output

xdmffile_u = XDMFFile('navier_stokes_cylinder/velocity.xdmf')
xdmffile_p = XDMFFile('navier_stokes_cylinder/pressure.xdmf')
# Create time series (for use in reaction_system.py)
timeseries_u = TimeSeries('navier_stokes_cylinder/velocity_series')
timeseries_p = TimeSeries('navier_stokes_cylinder/pressure_series')

# Save mesh to file (for use in reaction_system.py)
File('navier_stokes_cylinder/cylinder.xml.gz') << mesh

# Create progress bar
progress = Progress('Time-stepping')
set_log_level(LogLevel.PROGRESS)

# Time-stepping
t = 0

for n in range(num_steps):
    # Update current time
    t += dt

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

    # Step 4: Scalar solve step
    b4 = assemble(L4)
    [bc.apply(b4) for bc in bcs]
    solve(A4,s_,b4)

    # Save solution to file (XDMF/HDF5)
    xdmffile_u.write(u_, t)
    xdmffile_p.write(p_, t)

    # Save nodal values to file
    timeseries_u.store(u_.vector(), t)
    timeseries_p.store(p_.vector(), t)

    # Update previous solution
    u_n.assign(u_)
    p_n.assign(p_)
    s_n.assign(s_)

    # Update progress bar
    progress=(t / T)
    print('u max:', np.array(u_.vector().max()))

What am I doing wrong here? For some reason I am not using the correct arguments for the solve call. I am hoping you could provide me with some directions.

With kind regards,

2

It looks like you just forgot to pass the vector associated with s_ instead of s_ itself:

    #solve(A4,s_,b4)
    solve(A4,s_.vector(),b4)
1 Like
3

Thank you so much, after this change the code runs without errors.
However, I do not get the results I expect.

Now I get the following result at t=0.:
image

At t=3.:
image

At the inlet boundary, the Dirichlet condition of “Constant(1)” seems to be correctly applied as can be seen in the figures, but for some reason the scalar at the inlet is not convected across the domain. I probably made some other kind of error, which I am again struggling to find.

4

Could it be that my variational formulation is not correct?

5

The main problem here is that the time stepping loop is set up assuming explicit integration of advection, where the LHS form a4 is a mass matrix that only needs to be assembled once, but the formulation of F4 includes advection implicitly, such that a4 would need to be reassembled each step with the updated velocity. You can modify F4 to include advection explicitly like

F4 = rho*((s - s_n) / k) *z*dx \
   + rho*dot(u_n, grad(s_n))*z*dx \
   - fs*z*dx

or update the solver step like

    # Step 4: Scalar solve step
    A4 = assemble(a4)
    b4 = assemble(L4)
    [bc.apply(b4) for bc in bcs]
    [bc.apply(A4) for bc in bcs]
    solve(A4,s_.vector(),b4)

to keep the original implicit definition of F4. In both cases, you can see that s_ now advects into the domain over time.

However, you will also observe that the computed s_ is pretty noisy, which is typical of Galerkin solutions to advection-dominated problems. In general, you want to also incorporate some form of stabilization, e.g., SUPG, possibly augmented with a nonlinear shock-capturing viscosity.

3 Likes
6

Thank you very much again, for the explanation and the suggestion for further improvements.

Now, I would like to apply the solver to a different mesh.
To construct this mesh, I used the built-in mshr.

The resulting mesh looks as follows, supposed to represent a 2D stenosis:
image
image

However, when using this mesh and BCs, I am getting convergence issues. (at approximately t=0.05)
I tried coursening and refining the mesh, I tried decreasing the timestep and I tried varying the viscosity (and a combination of these variations), but to no improvement.

Do you have any suggestions for improving convergence for this case?

Currently, the full code looks like this:

from __future__ import print_function
from fenics import *
from mshr import *
import numpy as np
import matplotlib.pyplot as plt
import math
T = 1.            # final time
num_steps = 1000   # number of time steps
dt = T / num_steps # time step size
mu = 1e-3         # dynamic viscosity
rho = 1.           # density


# Create mesh
channel = Rectangle(Point(-0.01715, -0.001375), Point(0.01715, 0.001375))
cylinder1 = Ellipse(Point(0, 0.000975), 0.0017,0.0004)
cylinder2 = Ellipse(Point(0, -0.000975), 0.0017,0.0004)
rect1 = Rectangle(Point(-0.0034, 0.001375), Point(0.0034, 0.000975))
rect2 = Rectangle(Point(-0.0034, -0.001375), Point(0.0034, -0.000975))
cylinder3 = Ellipse(Point(-0.0034, 0.000975), 0.0017,0.0004)
cylinder4 = Ellipse(Point(-0.0034, -0.000975), 0.0017,0.0004)
cylinder5 = Ellipse(Point(0.0034, 0.000975), 0.0017,0.0004)
cylinder6 = Ellipse(Point(0.0034, -0.000975), 0.0017,0.0004)
domain = channel - cylinder1-cylinder2-rect1-rect2+cylinder3+cylinder4+cylinder5+cylinder6
mesh = generate_mesh(domain, 40)
plot(mesh)
plt.show()

# Define function spaces
V = VectorFunctionSpace(mesh, 'P', 2)
Q = FunctionSpace(mesh, 'P', 1)
S = FunctionSpace(mesh,'P',2)

# Define boundaries
inflow   = 'near(x[0], -0.01715)'
outflow  = 'near(x[0], 0.01715)'
walls    = 'near(x[1], -0.001375) ||' \
 + 'near(x[1], 0.001375) ||' \
 + 'on_boundary && x[0]>-0.0034 && x[0]<0.0034 && x[1]>-0.001375 && x[1]<-0.0005 ||' \
 + 'on_boundary && x[0]>-0.0034 && x[0]<0.0034 && x[1]>0.0005 && x[1]<0.001375'


plug_velocity=0.2e-6/(math.pi*0.001375**2)
# Define inflow profile
inflow_profile = ('2*plug_velocity*((0.001375-x[1])*(0.001375+x[1])) / pow(0.001375, 2)', '0')
vely=Expression(inflow_profile, plug_velocity=plug_velocity, degree=2)

# Define boundary conditions
bcu_inflow = DirichletBC(V, vely, inflow)
bcu_walls = DirichletBC(V, Constant((0, 0)), walls)
#bcu_cylinder = DirichletBC(V, Constant((0, 0)), cylinder)
bcp_outflow = DirichletBC(Q, Constant(0), outflow)
bcs_inflow = DirichletBC(S, Constant(0.),inflow)

bcu = [bcu_inflow, bcu_walls]
bcp = [bcp_outflow]
bcs = [bcs_inflow]

# Define trial and test functions
u = TrialFunction(V)
v = TestFunction(V)
p = TrialFunction(Q)
q = TestFunction(Q)
s = TrialFunction(S)
z = TestFunction(S)

# Define functions for solutions at previous and current time steps
u_n = Function(V)
u_  = Function(V)
p_n = Function(Q)
p_  = Function(Q)
s_n = Function(S)
s_ = Function(S)

# Define expressions used in variational forms
U  = 0.5*(u_n + u)
n  = FacetNormal(mesh)
f  = Constant((0, 0))
fs  = Constant(0.)
k  = Constant(dt)
mu = Constant(mu)
rho = Constant(rho)

# 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

# Define variational problem for step 4

F4 = rho*((s - s_n) / k)*z*dx \
   + rho*dot(u_n, grad(s))*z*dx \
   - fs*z*dx
a4 = lhs(F4)
L4 = rhs(F4)

# 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 XDMF files for visualization output
xdmffile_u = XDMFFile('navier_stokes_cylinder/velocity.xdmf')
xdmffile_p = XDMFFile('navier_stokes_cylinder/pressure.xdmf')

# Create time series (for use in reaction_system.py)
timeseries_u = TimeSeries('navier_stokes_cylinder/velocity_series')
timeseries_p = TimeSeries('navier_stokes_cylinder/pressure_series')

# Save mesh to file (for use in reaction_system.py)
File('navier_stokes_cylinder/cylinder.xml.gz') << mesh

# Create progress bar
progress = Progress('Time-stepping')
set_log_level(LogLevel.PROGRESS)

# Time-stepping
t = 0
for n in range(num_steps):

    # Update current time
    t += dt
    if n%10==0:
        print(t)
    flowsec=0.4e-6
    baseflow = 0.6e-6
    vely.plug_velocity = (baseflow+flowsec*math.sin(4*math.pi*t-0.5*math.pi))/(math.pi*0.001375**2)

    # 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')
	
    # Step 4: Scalar solve step
    A4 = assemble(a4)
    [bc.apply(A4) for bc in bcs]
    b4 = assemble(L4)
    [bc.apply(b4) for bc in bcs]
    solve(A4,s_.vector(),b4)
    if n%100==0:

    # Plot solution
        plot(u_, title='Velocity')
        plt.show()
        plot(p_, title='Pressure')
        plt.show()
        #plot(s_, title='Scalar')
        c = plot(s_, mode='color')
        plt.colorbar(c)
        plt.show()

    # Save solution to file (XDMF/HDF5)
    xdmffile_u.write(u_, t)
    xdmffile_p.write(p_, t)

    # Save nodal values to file
    timeseries_u.store(u_.vector(), t)
    timeseries_p.store(p_.vector(), t)

    # Update previous solution
    u_n.assign(u_)
    p_n.assign(p_)
    s_n.assign(s_)
    
    # Update progress bar
    progress=(t / T)
    #print('u max:', np.array(u_.vector().max()))
7

Actually, with a higher mesh density and a finer timestep I got a working code.
Thank you so much for your help.

The reason I was undertaking this exercise, is to obtain a MWE for my previous post (Strain history in Newtonian flow). My aim is to implement the frobenius norm of the strain as a source term in the convection diffusion equation.

In the previous post I referred to the Oasis solver without a MWE, so it would be difficult for anyone to answer my questions or to figure out the problem myself.
I will give this implementation another attempt and maybe resubmit my question regarding “Strain history in convection diffusion equation” if I do not succeed.