Navier-stokes using Matrix-vector multiplication: Oasis

Linear Algebra
1

I am trying to optimize a “fluid-structure interaction” solver which makes use of the fluid-momentum and pressure poisons equation. These equations are similar to navier stokes - fenics except that it is a non-dimensional solver, has slightly different discretization and has a few additional source terms. A working code of just the flow solver (excluding other terms) is provided below.

from dolfin import *
from mshr import *
import time

T = 3              # Total time
dt = 0.01          # Time step 

# Non-dimensional numbers
Re = 75    

mesh1 = RectangleMesh(Point(0.0, 0.0), import time

T = 3              # Total time
dt = 0.01          # Time step 

# Non-dimensional numbers
Re = 75    

mesh1 = RectangleMesh(Point(0.0, 0.0), Point(2.0, 1.0), 80, 40)

pv1 = File("mesh1.pvd")
pv1 << mesh1

# MPI-rank
mpi_comm = mpi_comm_world() 
MPI.barrier(mpi_comm)
my_rank = MPI.rank(mpi_comm)

# Define function spaces
V  = VectorFunctionSpace(mesh1, 'P', 2)		        # Fluid velocity on mesh1  
Q  = FunctionSpace(mesh1, 'P', 1)		            # Fluid pressure and temperature on mesh1

# Define boundaries for combined domain
bottom = 'near(x[1], 0)'
top    = 'near(x[1], 1)'
left   = 'near(x[0], 0)'
right  = 'near(x[0], 2)'

# Define boundary conditions for combined domain
bcu_left = DirichletBC(V, Constant((1, 0)), left)
bcp_right = DirichletBC(Q, Constant(0), right)
bcu_bottom = DirichletBC(V, Constant((0, 0)), bottom)
bcu_top = DirichletBC(V, Constant((0, 0)), top)
bcu = [bcu_bottom, bcu_top, bcu_left]
bcp = [bcp_right]

# Define trial and test functions
u1  = TrialFunction(V)            # Trial function for fluid velocity on mesh1
v   = TestFunction(V)             # Test function for fluid veloctiy on mesh1
p   = TrialFunction(Q)            # Trial function for fluid pressure on mesh1
q   = TestFunction(Q)             # Test function for fluid pressure on mesh1

# Define functions
u_n       = Function(V)
us_       = Function(V)
p_n       = Function(Q)
p_        = Function(Q)

# Defining integration domains
dx1 = dolfin.dx(mesh1)
ds1 = dolfin.ds(mesh1)     

# Define expressions
U    = 0.5*(u_n + u1)
n    = FacetNormal(mesh1)
tsp  = dt
dt   = Constant(dt)
Re   = Constant(Re)

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

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

# Define variational problem for step 1
F1 = dot((u1 - u_n) / dt, v)*dx1 + dot(dot(u_n, nabla_grad(u_n)), v)*dx1 + inner(sigma(U, p_n), epsilon(v))*dx1 + dot(p_n*n, v)*ds1 - dot((1/Re)*nabla_grad(U)*n, v)*ds1 
a1 = lhs(F1)
L1 = rhs(F1)

# Define variational problem for step 2
a2 = dot(nabla_grad(p), nabla_grad(q))*dx1
L2 = dot(nabla_grad(p_n), nabla_grad(q))*dx1 - (1/dt)*div(us_)*q*dx1

# Define variational problem for step 3
a3 = dot(u1, v)*dx1
L3 = dot(us_, v)*dx1 - dt*dot(nabla_grad(p_ - p_n), v)*dx1

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

# Time
t = 0
c = 0

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

file_u  = File("u.pvd")
file_p  = File("p.pvd")

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

start_time = time.time()

while t < T:

    # Update current time
    t += tsp
    
    c += 1

    if my_rank == 0:
    	print "Step 1: Tentative velocity step"
    b1 = assemble(L1)
    [bc.apply(b1) for bc in bcu]
    solve(A1, us_.vector(), b1, 'bicgstab', 'hypre_amg')

    if my_rank == 0:
    	print "Step 2: Pressure correction step"
    b2 = assemble(L2)
    [bc.apply(b2) for bc in bcp]
    solve(A2, p_.vector(), b2, 'bicgstab', 'hypre_amg')

    if my_rank == 0:
    	print "Step 3: Velocity correction step"
    b3 = assemble(L3)
    solve(A3, us_.vector(), b3, 'bicgstab', 'sor')
    
    # Print output files
    if c > 10:

	    c = 0

	    file_u  << us_
	    file_p  << p_
	    
    # Update progress bar
    if my_rank == 0:
    	progress.update(t / T)

    # Update previous solution
    u_n.assign(us_)
    p_n.assign(p_)

end_time = time.time()
print(" total time", end_time - start_time)Point(2.0, 1.0), 80, 40)

pv1 = File("mesh1.pvd")
pv1 << mesh1

# MPI-rank
mpi_comm = mpi_comm_world() 
MPI.barrier(mpi_comm)
my_rank = MPI.rank(mpi_comm)

# Define function spaces
V  = VectorFunctionSpace(mesh1, 'P', 2)		        # Fluid velocity on mesh1  
Q  = FunctionSpace(mesh1, 'P', 1)		            # Fluid pressure and temperature on mesh1

# Define boundaries for combined domain
bottom = 'near(x[1], 0)'
top    = 'near(x[1], 1)'
left   = 'near(x[0], 0)'
right  = 'near(x[0], 2)'

# Define boundary conditions for combined domain
bcu_left = DirichletBC(V, Constant((1, 0)), left)
bcp_right = DirichletBC(Q, Constant(0), right)
bcu_bottom = DirichletBC(V, Constant((0, 0)), bottom)
bcu_top = DirichletBC(V, Constant((0, 0)), top)
bcu = [bcu_bottom, bcu_top, bcu_left]
bcp = [bcp_right]

# Define trial and test functions
u1  = TrialFunction(V)            # Trial function for fluid velocity on mesh1
v   = TestFunction(V)             # Test function for fluid veloctiy on mesh1
p   = TrialFunction(Q)            # Trial function for fluid pressure on mesh1
q   = TestFunction(Q)             # Test function for fluid pressure on mesh1

# Define functions
u_n       = Function(V)
us_       = Function(V)
p_n       = Function(Q)
p_        = Function(Q)

# Defining integration domains
dx1 = dolfin.dx(mesh1)
ds1 = dolfin.ds(mesh1)     

# Define expressions
U    = 0.5*(u_n + u1)
n    = FacetNormal(mesh1)
tsp  = dt
dt   = Constant(dt)
Re   = Constant(Re)

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

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

# Define variational problem for step 1
F1 = dot((u1 - u_n) / dt, v)*dx1 + dot(dot(u_n, nabla_grad(u_n)), v)*dx1 + inner(sigma(U, p_n), epsilon(v))*dx1 + dot(p_n*n, v)*ds1 - dot((1/Re)*nabla_grad(U)*n, v)*ds1 
a1 = lhs(F1)
L1 = rhs(F1)

# Define variational problem for step 2
a2 = dot(nabla_grad(p), nabla_grad(q))*dx1
L2 = dot(nabla_grad(p_n), nabla_grad(q))*dx1 - (1/dt)*div(us_)*q*dx1

# Define variational problem for step 3
a3 = dot(u1, v)*dx1
L3 = dot(us_, v)*dx1 - dt*dot(nabla_grad(p_ - p_n), v)*dx1

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

# Time
t = 0
c = 0

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

file_u  = File("u.pvd")
file_p  = File("p.pvd")

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

start_time = time.time()

while t < T:

    # Update current time
    t += tsp
    
    c += 1

    if my_rank == 0:
    	print "Step 1: Tentative velocity step"
    b1 = assemble(L1)
    [bc.apply(b1) for bc in bcu]
    solve(A1, us_.vector(), b1, 'bicgstab', 'hypre_amg')

    if my_rank == 0:
    	print "Step 2: Pressure correction step"
    b2 = assemble(L2)
    [bc.apply(b2) for bc in bcp]
    solve(A2, p_.vector(), b2, 'bicgstab', 'hypre_amg')

    if my_rank == 0:
    	print "Step 3: Velocity correction step"
    b3 = assemble(L3)
    solve(A3, us_.vector(), b3, 'bicgstab', 'sor')
    
    # Print output files
    if c > 10:

	    c = 0

	    file_u  << us_
	    file_p  << p_
	    
    # Update progress bar
    if my_rank == 0:
    	progress.update(t / T)

    # Update previous solution
    u_n.assign(us_)
    p_n.assign(p_)

end_time = time.time()
print(" total time", end_time - start_time)

My aim is to solve complex 3D FSI problems with Temperature transport. Also Taylor hood elements are a must. To improve overall computational speed, I decompose the terms arising in navier stokes into separate pre-assembled matrices (again here the concept is exactly similar to Oasis-fenics). “Note that I cannot incorporate Oasis directly as mine is an FSI solver with delta interpolation functions and source terms”. The theory for decomposing my equations into preassembled matrices is provided below in writing.

IMG_4711
IMG_47112381×4032 2.76 MB
IMG_4712
IMG_47122821×3618 2.88 MB

The corresponding working code script is provided here.

from dolfin import *
from mshr import *
import time

T = 3              # Total time
dt = 0.01          # Time step 

# Non-dimensional numbers
Re = 75    

mesh1 = RectangleMesh(Point(0.0, 0.0), Point(2.0, 1.0), 80, 40)

pv1 = File("mesh1.pvd")
pv1 << mesh1

# MPI-rank
mpi_comm = mpi_comm_world() 
MPI.barrier(mpi_comm)
my_rank = MPI.rank(mpi_comm)

# Define function spaces
V  = VectorFunctionSpace(mesh1, 'P', 2)		        # Fluid velocity on mesh1  
Q  = FunctionSpace(mesh1, 'P', 1)		            # Fluid pressure and temperature on mesh1

# Define boundaries for combined domain
bottom = 'near(x[1], 0)'
top    = 'near(x[1], 1)'
left   = 'near(x[0], 0)'
right  = 'near(x[0], 2)'

# Define boundary conditions for combined domain
bcu_left = DirichletBC(V, Constant((1, 0)), left)
bcp_right = DirichletBC(Q, Constant(0), right)
bcu_bottom = DirichletBC(V, Constant((0, 0)), bottom)
bcu_top = DirichletBC(V, Constant((0, 0)), top)
bcu = [bcu_bottom, bcu_top, bcu_left]
bcp = [bcp_right]

# Define trial and test functions
u1  = TrialFunction(V)            # Trial function for fluid velocity on mesh1
v   = TestFunction(V)             # Test function for fluid veloctiy on mesh1
p   = TrialFunction(Q)            # Trial function for fluid pressure on mesh1
q   = TestFunction(Q)             # Test function for fluid pressure on mesh1

# Define functions
u_n       = Function(V)
us_       = Function(V)
p_n       = Function(Q)
p_        = Function(Q)

# Defining integration domains
dx1 = dolfin.dx(mesh1)
ds1 = dolfin.ds(mesh1)     

# Define expressions
U    = 0.5*(u_n + u1)
n    = FacetNormal(mesh1)
tsp  = dt
dt   = Constant(dt)
Re   = Constant(Re)

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

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

# ............ Define optimized variational problem for step 1 .............#

Mij = assemble(dot(u1/dt, v)*dx1)											# Mass matrix 
Cij = assemble(dot(dot(u_n, nabla_grad(u1)), v)*dx1)						# Convective matrix
Kij = assemble(inner((2/Re)*epsilon(0.5*u1), epsilon(v))*dx1)				# Viscous matrix
Sij = assemble(inner(p*Identity(len(u1)), epsilon(v))*dx1)					# Pressure matrix

# ............ Surface integral matrices .............
Zij = assemble(dot(p*n, v)*ds1)
Oij = assemble(dot((1/Re)*nabla_grad(0.5*u1)*n, v)*ds1)

# ............ Matrix algebra ........................
X1 = Cij.copy()
X1 *= -1.0; X1 += Mij - Kij + Oij
Y1 = Sij.copy()
Y1 -= Zij

def optimization_lhs(Mij):      
    
    A1 = Mij.copy()
    A1 += Kij - Oij

    return A1

def optimization_rhs(X1, Y1):

    b1 = X1*u_n.vector()   
    b1 += Y1*p_n.vector()
    
    return b1

# ............ Define optimized variational problem for step 1 .............#

# Define variational problem for step 2
a2 = dot(nabla_grad(p), nabla_grad(q))*dx1
L2 = dot(nabla_grad(p_n), nabla_grad(q))*dx1 - (1/dt)*div(us_)*q*dx1

# Define variational problem for step 3
a3 = dot(u1, v)*dx1
L3 = dot(us_, v)*dx1 - dt*dot(nabla_grad(p_ - p_n), v)*dx1

# Assemble matrices
A1 = optimization_lhs(Mij)
A2 = assemble(a2)
A3 = assemble(a3)

# Time
t = 0
c = 0

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

file_u  = File("u.pvd")
file_p  = File("p.pvd")

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

start_time = time.time()

while t < T:

    # Update current time
    t += tsp
    
    c += 1

    if my_rank == 0:
    	print "Step 1: Tentative velocity step"
    b1 = optimization_rhs(X1, Y1)
    [bc.apply(b1) for bc in bcu]
    solve(A1, us_.vector(), b1, 'bicgstab', 'hypre_amg')

    if my_rank == 0:
    	print "Step 2: Pressure correction step"
    b2 = assemble(L2)
    [bc.apply(b2) for bc in bcp]
    solve(A2, p_.vector(), b2, 'bicgstab', 'hypre_amg')

    if my_rank == 0:
    	print "Step 3: Velocity correction step"
    b3 = assemble(L3)
    solve(A3, us_.vector(), b3, 'bicgstab', 'sor')
    
    # Print output files
    if c > 10:

	    c = 0

	    file_u  << us_
	    file_p  << p_
	    
    # Update progress bar
    if my_rank == 0:
    	progress.update(t / T)

    # Update previous solution
    u_n.assign(us_)
    p_n.assign(p_)

end_time = time.time()
print(" total time", end_time - start_time)

The only difference in both codes can be seen in the variational form of Step 1. The first (naive) and second (optimized) code consistently take approx 38 sec and 35 sec to run. That explains the small speedup for this problem. However the results are different. Given below, the results show that code1 produces the correct result with deeper contours in the channel however code2 produces shallower contours in the channel (behaves like it has a low Reynolds number), which is not correct.

code 1 : code1 code2 :code2

I have tried a lot in code2 to figure out the issue, but have failed so far. What I concluded is that the problem lies in the viscous matrix - Kij. Any help provided to resolve this issue would be very helpful? Thanks in advance.

2 Likes
2

Some pointers:

  • I always start solving a manufactured problem, to ensure that my scheme is implemented correctly. I would suggest considering the Taylor-Green vortex problem.
  • You should minimize your example. Check if the difference is visible after the first tentative velocity step, and reduce your code example to only contain this.
3

I will switch to the Taylor green problem and will surely reduce the problem to one time step. Thanks a lot dokken.