Solving for shear stress and using mixed function spaces

1

Hi all, this is a two-fold question. I am trying to solve for the wall shear stress in a basic channel flow. To do so, I believe that I need to use a MixedFunctionSpace as a vector function space for the velocity and scalar function space for the pressure. However, I get the following error when trying to solve for the velocity (line 165):

“Preconditioner number of local columns 4600 does not equal output vector size 4066”

Link for gmsh msh

I believe this issue has something to do with my mixed function space Y as a past script which didn’t used a mixed function space did not have this issue. That said, I believe I need the mixed function space as stated above. My code is below:

# Import stuff
from dolfin import *
import numpy as np
import meshio

set_log_level(1)

# Import mesh from gmsh using meshio
msh = meshio.read("../../domain/rectangle.msh")

# Extract cells and physical data
def create_mesh(mesh, cell_type, prune_z=True):
    cells = mesh.get_cells_type(cell_type)
    cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
    out_mesh = meshio.Mesh(points=mesh.points, cells={cell_type: cells}, cell_data={"name_to_read":[cell_data]})
    if prune_z:
        out_mesh.prune_z_0()
    return out_mesh

# Triangle mesh
triangle_mesh = create_mesh(msh, "triangle", prune_z=True)
meshio.write("mesh.xdmf", triangle_mesh)

mesh = Mesh()
with XDMFFile("mesh.xdmf") as infile:
    infile.read(mesh)
mvc = MeshValueCollection("size_t", mesh, 1) 

# Line mesh
line_mesh = create_mesh(msh, "line", prune_z=True)
meshio.write("mf.xdmf", line_mesh)

with XDMFFile("mf.xdmf") as infile:
    infile.read(mvc, "name_to_read")
mf = MeshFunction("size_t", mesh, mvc)

domains = MeshFunction("size_t", mesh, mesh.topology().dim())

# # # # # # # # # # # # # # # # # # # # # # # # # # # # # #  # # # # 
# # -------------------- Begin Problem ------------------------- # # 
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # #  # # # # 

# Define knowns
T         = 3.0           # time duration
num_steps = 1000          # number of time steps
dt        = T / num_steps # time step size
mu        = 0.01          # dynamic viscosity
rho       = 1             # density

# Generate Finite Element Space
vp = VectorElement('P', mesh.ufl_cell(), 2)
qp = FiniteElement('P', mesh.ufl_cell(), 1)

# Mixed elements
mex = MixedElement([vp, qp])

# Generate Function Spaces
V = FunctionSpace(mesh, vp)
Q = FunctionSpace(mesh, qp)

## Mixed Function Space
Y = FunctionSpace(mesh, mex)

# Trial Functions
(u, p) = TrialFunctions(Y)
#y   = TrialFunction(Y)
#u, p = split(y)

# Test Functions
(v, q) = TestFunctions(Y)

# Define Functions
u_n = Function(V)
u_  = Function(V)
p_n = Function(Q)
p_  = Function(Q)

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

# Define strain rate tensor/symmetric gradient
def epsilon(u):
    return 0.5*(nabla_grad(u) + nabla_grad(u).T)

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

# Surface Traction
Trac = sigma(u, p)*n

# Normal and tangential components
Tn = inner(Trac, n)
Tt = Trac - Tn*n

# Assemble functions for stress
shear_stress = Function(V)
Lt = (1 / FacetArea(mesh)) * inner(v, Tt) * ds

# Progess bar
progress = Progress('Time-stepping', num_steps)

# Inflow profile
inflow_profile = Expression(('4.0*1.5*x[1]*(0.25 - x[1]) / pow(0.25, 2)', '0'), degree=2)

# Assign Boundary Conditions using markers from gmsh
bcu_inflow  = DirichletBC(Y.sub(0), inflow_profile, mf, 5)
bcp_outflow = DirichletBC(Y.sub(1), Constant(0), mf, 6)
bcu_walls   = DirichletBC(Y.sub(0), Constant((0,0)), mf, 7)

# Collect Boundary Conditions
bcu = [bcu_inflow, bcu_walls]
bcp = [bcp_outflow]

# # # # # # # # # # # # # #
#  Variational Problems   #
# # # # # # # # # # # # # #

# Step 1: Tentative Velocity
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)

# Step 2: New Pressure, P_n
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

# Step 3: Corrected Velocity, u^{n+1}
a3 = dot(u, v)*dx
L3 = dot(u_, v)*dx - k*dot(nabla_grad(p_ - p_n), v)*dx

# Time-stepping
t = 0

for j in range(num_steps):

    # Update time
    inflow_profile.t = t

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

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

    # Update current time
    t += dt

    # Step 1: Tentative Velocity
    b1 = assemble(L1)
    [bc.apply(b1) for bc in bcu]
    solve(A1, u_.vector(), b1, "mumps", "default")

    # Step 2: New Pressure, P_n
    b2 = assemble(L2)
    [bc.apply(b2) for bc in bcp]
    solve(A2, p_.vector(), b2, "mumps", "default")

    # Step 3: Corrected Velocity
    b3 = assemble(L3)
    solve(A3, u_.vector(), b3, "mumps", "default")

    # Assemble shear stress
    assemble(Lt, tensor=shear_stress.vector())

    # Save Solutions
    if (j % 10 == 00):

      # Save solution as vtk
      ufile = File('vtk/channel_velo/velo_%04d.pvd' % j)
      ufile << (u_, t)

      # Save solution as vtk
      shearfile = File('vtk/channel_shear/shear_%04d.pvd' % j)
      shearfile << (shear_stress, t)

    # Update Previous Solution
    u_n.assign(u_)
    p_n.assign(p_)

    # Update progress bar
    set_log_level(LogLevel.PROGRESS)
    progress += 1
2

If you are using a splitting method and solving separate systems for the velocity and pressure, you should define the variational forms for those systems using TrialFunctions and TestFunctions from the separate spaces V and Q. You only need to use mixed spaces if you are setting up a monolithic system of equations to solve for the pressure and velocity simultaneously.

3

Hi kamensky,

Thank you for the response. If I do not used a mixed space, and instead define my functions spaces as follows:

# Define Function Space
V = VectorFunctionSpace(mesh, 'P', 2)
Q = FunctionSpace(mesh, 'P', 1)

# Define Trial Functions
u = TrialFunction(V)
p = TrialFunction(Q)

# Define Test Functions
v = TestFunction(V)
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)

As well as update my BC formulation as:

# Assign Boundary Conditions using markers from gmsh
bcu_inflow  = DirichletBC(V, inflow_profile, mf, 5)
bcp_outflow = DirichletBC(Q, Constant(0), mf, 6)
bcu_walls   = DirichletBC(V, Constant((0,0)), mf, 7)

Then I get an error of

Did you combine test or trial functions from different spaces?
The Arguments found are:
  v_1
  v_0
  v_1

Because the stress tensor is a function of both velocity and pressure:

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

How would I go about addressing this?

4

I didn’t try downloading the mesh and running the code, but I’m guessing that this error is coming from the assembly of Lt. Based on my interpretation of the code, I believe you want to use the actual solution Functions u_ and p_ to define the traction Trac, not the TrialFunctions u and p. Keep in mind that TrialFunctions (and TestFunctions) are just symbolic arguments to variational forms, and are not associated with vectors of degrees of freedom.

5

Sweet mother of god, I think that was it! Thank you.

Looking at the solutions in Paraview, they seem reasonable, but for the case of this simple channel flow problem, the maximum simulated shear stress is 0.16 (and stays here under steady state) while the analytic solution says that the max shear stress should be 0.24. Any ideas? I think the mesh should be fine enough. Are there other mistakes I have overlooked for my shear stress formulation?

Edit: I think it’s because the code is computing the average stress on each cell in the mesh

6

To obtain the facet-averaged traction from assembling Lt, I believe you’d want to have the test function v in the space of cell-wise constants, i.e., VectorFunctionSpace(mesh,"DG",0), instead of the space of continuous piecewise linear functions, V.