Impose periodic BC in pressure driven stokes flow

[Exon]

Hi Community,

I am solving the 2D Stokes flow similar to what was done in this post. However, instead of having periodic BCs on all four boundaries, I am prescribing pressures on the inlet and outlet (x=-12 and x=24 in my case) and having periodic BCs on y=0 and y=12. Here’s my MWE modified from demo_stokes.py and several dolfinx_mpc demos:

from mpi4py import MPI
from petsc4py import PETSc

import numpy as np

import ufl
from basix.ufl import element, mixed_element
from dolfinx import fem, default_scalar_type
from dolfinx.fem import (Constant, Function, dirichletbc,form, functionspace,
                        locate_dofs_topological)
import dolfinx.fem.petsc
from dolfinx.mesh import CellType, create_rectangle, locate_entities, meshtags, locate_entities_boundary
from ufl import div, dx, ds, grad, inner, Measure, FacetNormal

# Define params
x_range = [-12,24]
y_range = [0,12]
p_in = 800
p_out = 0
mu = 3.6

# Create mesh
msh = create_rectangle(MPI.COMM_WORLD, [np.array([x_range[0], y_range[0]]), np.array([x_range[1], y_range[1]])],
                    [300, 100], CellType.triangle)


# Function to mark periodic boundaries
def periodic_boundary(x):
    return np.isclose(x[1], y_range[1], atol=1e-8)

def periodic_relation(x):
    out_x = np.zeros_like(x)
    out_x[0] = x[0]
    out_x[1] = y_range[1] - x[1] 
    #out_x[2] = x[2]
    return out_x

# Function to mark the inflow
def inflow(x):
    return np.isclose(x[0], x_range[0])

# Function to mark the outflow
def outflow(x):
    return np.isclose(x[0], x_range[1])


boundaries = [(1, periodic_boundary),
            (2, inflow),
            (3, outflow)] # (boundary_marker, locator)


# We now loop through all the boundary conditions and create `MeshTags` identifying the facets for each boundary condition.
facet_indices, facet_markers = [], []
gdim = msh.topology.dim
fdim = gdim - 1
for (marker, locator) in boundaries:
    facets = locate_entities(msh, fdim, locator)
    facet_indices.append(facets)
    facet_markers.append(np.full_like(facets, marker))
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = meshtags(msh, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])
msh.topology.create_connectivity(msh.topology.dim-1, msh.topology.dim)

"""
# Verify Boundaries
from dolfinx.io import XDMFFile
with XDMFFile(msh.comm, "facet_tags_verify.xdmf", "w") as xdmf:
    xdmf.write_mesh(msh)
    xdmf.write_meshtags(facet_tag, msh.geometry)
"""

# Define Surface Integral for Nuemann Conditions
ds = Measure("ds", domain=msh, subdomain_data=facet_tag)

# `P2` corresponds to a continuous piecewise quadratic basis (vector) and 
# `P1` to a continuous piecewise linear basis (scalar).


P2 = element("Lagrange", msh.basix_cell(), 2, shape=(msh.geometry.dim,))
P1 = element("Lagrange", msh.basix_cell(), 1)


# Create the Taylor-Hood function space
TH = mixed_element([P2, P1])
W = functionspace(msh, TH)
W0, _ = W.sub(0).collapse()
# Pressures
pressure_in = Constant(msh,PETSc.ScalarType(p_in))
pressure_out = Constant(msh,PETSc.ScalarType(p_out))

# Collect Dirichlet boundary conditions
bcs = []
# Suface Normal
n = FacetNormal(msh)
# Define variational problem
(u, p) = ufl.TrialFunctions(W)
(v, q) = ufl.TestFunctions(W)
f = Function(W0)
a = form((mu*inner(grad(u), grad(v)) - inner(p, div(v)) + inner(div(u), q)) * dx)
L = form(inner(f, v) * dx - pressure_in * inner(n,v) * ds(2) - pressure_out * inner(n,v) * ds(3))

# Set BC with mpc
from dolfinx_mpc import LinearProblem, MultiPointConstraint
mpc = MultiPointConstraint(W)
mpc.create_periodic_constraint_topological(W.sub(0), facet_tag, 1, periodic_relation, bcs, default_scalar_type(1))
mpc.create_periodic_constraint_topological(W.sub(1), facet_tag, 1, periodic_relation, bcs, default_scalar_type(1))
mpc.finalize()

# Solve
problem = LinearProblem(a, L, mpc, bcs=bcs, petsc_options={
    "ksp_type": "preonly", "pc_type": "lu", "pc_factor_mat_solver_type": "mumps"})  # define linear problem
U = problem.solve()
u, p = U.sub(0).collapse(), U.sub(1).collapse()

Here are the meshtags I saved to verify boundary:

And here are the results:

that do not really make sense to me. Am I doing something wrong? Thanks in advance!

[tuderic]

First, I am wondering why you do not impose Dirichlet BCs for the given pressures at the inlet and outlet, but use them in the variational formulation.

Second, you would need to constrain each component of the velocity with periodic relations, i.e.

mpc.create_periodic_constraint_topological(W.sub(0), facet_tag, 1, periodic_relation, bcs, default_scalar_type(1))

needs to be

mpc.create_periodic_constraint_topological(W.sub(0).sub(0), facet_tag, 1, periodic_relation, bcs, default_scalar_type(1))
mpc.create_periodic_constraint_topological(W.sub(0).sub(1), facet_tag, 1, periodic_relation, bcs, default_scalar_type(1))

[Exon]

Hi @tuderic, thanks a lot for your reply! In fact your post inspired me a lot. Thanks again.

@tuderic

First, I am wondering why you do not impose Dirichlet BCs for the given pressures at the inlet and outlet, but use them in the variational formulation.

cuz I believe they are Neumann conditions which should be weakly imposed? I think

dofs_in = locate_dofs_topological(W.sub(1), 1, facet_tag.find(2))
bc_in = dirichletbc(PETSc.ScalarType(p_in), dofs_in, W.sub(1))
dofs_out = locate_dofs_topological(W.sub(1), 1, facet_tag.find(3))
bc_out = dirichletbc(PETSc.ScalarType(p_out), dofs_in, W.sub(1))

should also work if I change the variational form to

a = form((mu*inner(grad(u), grad(v)) + inner(p, div(v)) + inner(div(u), q)) * dx)
L = form(inner(f, v) * dx)

just like you did in that post but I couldnt get a mearningful result. Here’s the modified MWE:

from mpi4py import MPI
from petsc4py import PETSc

import numpy as np

import ufl
from basix.ufl import element, mixed_element
from dolfinx import fem, default_scalar_type
from dolfinx.fem import (Constant, Function, dirichletbc,form, functionspace,
                        locate_dofs_topological)
import dolfinx.fem.petsc
from dolfinx.mesh import CellType, create_rectangle, locate_entities, meshtags, locate_entities_boundary
from ufl import div, dx, ds, grad, inner, Measure, FacetNormal

x_range = [-12,24]
y_range = [0,12]
p_in = 800
p_out = 0
mu = 3.6

# Create mesh
msh = create_rectangle(MPI.COMM_WORLD, [np.array([x_range[0], y_range[0]]), np.array([x_range[1], y_range[1]])],
                    [300, 100], CellType.triangle)


# Function to mark no slip boundaries
def periodic_boundary(x):
    return np.isclose(x[1], y_range[1], atol=1e-8)

def periodic_relation(x):
    out_x = np.zeros_like(x)
    out_x[0] = x[0]
    out_x[1] = y_range[1] - x[1] 
    #out_x[2] = x[2]
    return out_x

# Function to mark the inflow
def inflow(x):
    return np.isclose(x[0], x_range[0])

# Function to mark the outflow
def outflow(x):
    return np.isclose(x[0], x_range[1])


boundaries = [(1, periodic_boundary),
            (2, inflow),
            (3, outflow)] # (boundary_marker, locator)


# We now loop through all the boundary conditions and create `MeshTags` identifying the facets for each boundary condition.
facet_indices, facet_markers = [], []
gdim = msh.topology.dim
fdim = gdim - 1
for (marker, locator) in boundaries:
    facets = locate_entities(msh, fdim, locator)
    facet_indices.append(facets)
    facet_markers.append(np.full_like(facets, marker))
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = meshtags(msh, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])
msh.topology.create_connectivity(fdim, gdim)

P2 = element("Lagrange", msh.basix_cell(), 2, shape=(msh.geometry.dim,))
P1 = element("Lagrange", msh.basix_cell(), 1)
#V, Q = functionspace(msh, P2), functionspace(msh, P1)

# Create the Taylor-Hood function space
TH = mixed_element([P2, P1])
W = functionspace(msh, TH)

W0, _ = W.sub(0).collapse()
# pressure inlet outlet conditions
dofs_in = locate_dofs_topological(W.sub(1), 1, facet_tag.find(2))
bc_in = dirichletbc(PETSc.ScalarType(p_in), dofs_in, W.sub(1))
dofs_out = locate_dofs_topological(W.sub(1), 1, facet_tag.find(3))
bc_out = dirichletbc(PETSc.ScalarType(p_out), dofs_in, W.sub(1))

# Collect Dirichlet boundary conditions
bcs = [bc_in, bc_out]

# Define variational problem
(u, p) = ufl.TrialFunctions(W)
(v, q) = ufl.TestFunctions(W)
f = Function(W0)
a = form((mu*inner(grad(u), grad(v)) + inner(p, div(v)) + inner(div(u), q)) * dx)
L = form(inner(f, v) * dx)

from dolfinx_mpc import LinearProblem, MultiPointConstraint
mpc = MultiPointConstraint(W)
mpc.create_periodic_constraint_topological(W.sub(0).sub(0), facet_tag, 1, periodic_relation, bcs, default_scalar_type(1))
mpc.create_periodic_constraint_topological(W.sub(0).sub(1), facet_tag, 1, periodic_relation, bcs, default_scalar_type(1))
mpc.create_periodic_constraint_topological(W.sub(1), facet_tag, 1, periodic_relation, bcs, default_scalar_type(1))
mpc.finalize()

problem = LinearProblem(a, L, mpc, bcs=bcs, petsc_options={
    "ksp_type": "preonly", "pc_type": "lu", "pc_factor_mat_solver_type": "mumps"})  # define linear problem
U = problem.solve()
u, p = U.sub(0).collapse(), U.sub(1).collapse()

from dolfinx.io import XDMFFile
with XDMFFile(MPI.COMM_WORLD, "velocity_periodic_test1.xdmf", "w") as ufile_xdmf:
    u.x.scatter_forward()
    P1 = element("Lagrange", msh.basix_cell(), 1, shape=(msh.geometry.dim,))
    u1 = Function(functionspace(msh, P1))
    u1.interpolate(u)
    u1.name = "Velocity"
    ufile_xdmf.write_mesh(msh)
    ufile_xdmf.write_function(u1)

with XDMFFile(MPI.COMM_WORLD, "pressure_periodic_test1.xdmf", "w") as pfile_xdmf:
    p.x.scatter_forward()
    p.name = "Pressure"
    pfile_xdmf.write_mesh(msh)
    pfile_xdmf.write_function(p)

that gives me:

@tuderic

Second, you would need to constrain each component of the velocity with periodic relations, i.e.

With my weakly imposed pressure bcs, constraining each component of the velocity did produce a more resonable pressure field, but velocities with extremely larger magnitudes (it sort of make sense cuz there were nothing to stop the flows and the fact that they are uniform?)

[dokken]

@Exon

mpc.create_periodic_constraint_topological(W.sub(1), facet_tag, 1, periodic_relation, bcs, default_scalar_type(1))

Why are you imposing periodicity in pressure?

What do you want the boundary conditions at the side walls to be?

[Exon]

Hi @dokken , thank you for your reply!

@dokken

Why are you imposing periodicity in pressure?

I did it because I want to simulate flow in a small subdomain of a larger domain that is much wider in y-directions (consider a infinitely long (in y direction) pipe), and I am gonna put some symmetric obstructions in it. I thought velocity and pressure are coupled so I might need to do it for mixed elements. Maybe the periodicity in pressure is unnecessary?
Update: Getting rid of the pressure periodicity did give the same result by a visual inspection

@dokken

What do you want the boundary conditions at the side walls to be?

If you are refering to the top and bottom boundary. I just want them to be periodic to simulate a a infinitely long (in y direction) pipe. For left and right boundary I still want them to be inlet and outlet with prescribed pressure to drive the flow.

[dokken]

By only prescribing pressure on the side walls, you are enforcing \nabla u \cdot n=0 on the left and right boundary.

@Exon

I did it because I want to simulate flow in a small subdomain of a larger domain that is much wider in y-directions (consider a infinitely long (in y direction) pipe), and I am gonna put some symmetric obstructions in it. I thought velocity and pressure are coupled so I might need to do it for mixed elements. Maybe the periodicity in pressure is unnecessary?

If pressure_inlet (left) != pressure_outlet (right), then having periodic conditions in the top/bottom direction does not make sense.

In general, having strong conditions on pressure doesn’t make sense, see for instance: Pressure boundary conditions for blood flows | Chinese Annals of Mathematics, Series B

Here is a heavily adapted example to illustrate how to solve this problem:

from mpi4py import MPI
from petsc4py import PETSc

import numpy as np

import ufl
from basix.ufl import element, mixed_element
from dolfinx import fem, default_scalar_type
from dolfinx.fem import (Constant, Function, dirichletbc,form, functionspace,
                        locate_dofs_topological)
import dolfinx.fem.petsc
from dolfinx.mesh import CellType, create_rectangle, locate_entities, meshtags, locate_entities_boundary
from ufl import div, dx, ds, grad, inner, Measure, FacetNormal
from pathlib import Path
outdir = Path("test")

x_range = [-12,24]
y_range = [0,12]
p_in = 10
p_out = 2
mu = 3.6

# Create mesh
msh = create_rectangle(MPI.COMM_WORLD, [np.array([x_range[0], y_range[0]]), np.array([x_range[1], y_range[1]])],
                    [300, 150], CellType.triangle)


# Function to mark no slip boundaries
def periodic_boundary(x):
    return np.isclose(x[1], y_range[1], atol=1e-6)

def periodic_relation(x):
    out_x = np.zeros_like(x)
    out_x[0] = x[0]
    out_x[1] = y_range[1] - x[1] 
    #out_x[2] = x[2]
    return out_x

# Function to mark the inflow
def inflow(x):
    return np.isclose(x[0], x_range[0])

# Function to mark the outflow
def outflow(x):
    return np.isclose(x[0], x_range[1])


boundaries = [(1, periodic_boundary),
            (2, inflow),
            (3, outflow)] # (boundary_marker, locator)


# We now loop through all the boundary conditions and create `MeshTags` identifying the facets for each boundary condition.
facet_indices, facet_markers = [], []
gdim = msh.topology.dim
fdim = gdim - 1
for (marker, locator) in boundaries:
    facets = locate_entities(msh, fdim, locator)
    facet_indices.append(facets)
    facet_markers.append(np.full_like(facets, marker))
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = meshtags(msh, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])
msh.topology.create_connectivity(fdim, gdim)

P2 = element("Lagrange", msh.basix_cell(), 2, shape=(msh.geometry.dim,))
P1 = element("Lagrange", msh.basix_cell(), 1)

# Create the Taylor-Hood function space
TH = mixed_element([P2, P1])


W = functionspace(msh, TH)

from dolfinx_mpc import LinearProblem, MultiPointConstraint

# pressure inlet outlet conditions
bcs = []


mpc = MultiPointConstraint(W)
mpc.create_periodic_constraint_topological(W.sub(0).sub(0), facet_tag, 1, periodic_relation, bcs, default_scalar_type(1))
mpc.create_periodic_constraint_topological(W.sub(0).sub(1), facet_tag, 1, periodic_relation, bcs, default_scalar_type(1))
mpc.finalize()


# Define variational problem
(u, p) = ufl.TrialFunctions(mpc.function_space)
(v, q) = ufl.TestFunctions(mpc.function_space)
W0, _ = mpc.function_space.sub(0).collapse()

ds = ufl.Measure("ds", domain=msh, subdomain_data=facet_tag)
x = ufl.SpatialCoordinate(msh)
f = dolfinx.fem.Constant(msh, dolfinx.default_scalar_type((0.,0.)))
a = form((mu*inner(grad(u), grad(v)) - inner(p, div(v)) - inner(q, div(u))) * dx)
L = inner(f, v) * dx
n = ufl.FacetNormal(msh)
L -= p_in*ufl.dot(v, n)*ds(2)
L -= p_out*ufl.dot(v, n)*ds(3)
L = form(L)


problem = LinearProblem(a, L, mpc, bcs=bcs, petsc_options={
    "ksp_type": "preonly", "pc_type": "lu", "pc_factor_mat_solver_type": "mumps"})  # define linear problem
U = problem.solve()
u, p = U.sub(0).collapse(), U.sub(1).collapse()

print(problem.solver.getConvergedReason())


with dolfinx.io.VTXWriter(msh.comm, outdir / "demo_stokes_u.bp", u, engine="BP4") as vtx:
    vtx.write(0.0)
with dolfinx.io.VTXWriter(msh.comm, outdir / "demo_stokes_p.bp", p, engine="BP4") as vtx:
    vtx.write(0.0)

By changing the values in the last argument of:

mpc = MultiPointConstraint(W)
mpc.create_periodic_constraint_topological(W.sub(0).sub(0), facet_tag, 1, periodic_relation, bcs, default_scalar_type(1))
mpc.create_periodic_constraint_topological(W.sub(0).sub(1), facet_tag, 1, periodic_relation, bcs, default_scalar_type(1))
mpc.finalize()

to something else than 1, you can see that the periodicity in flow is properly enforced.

[Exon]

@dokken

By changing the values in the last argument of:

mpc = MultiPointConstraint(W)
mpc.create_periodic_constraint_topological(W.sub(0).sub(0), facet_tag, 1, periodic_relation, bcs, default_scalar_type(1))
mpc.create_periodic_constraint_topological(W.sub(0).sub(1), facet_tag, 1, periodic_relation, bcs, default_scalar_type(1))
mpc.finalize()

to something else than 1, you can see that the periodicity in flow is properly enforced.

Thank you for the explanation! This does work but I’m just wondering why scale = 1 won’t work because that’s sort of exactly what I want (to simulate a uniform part of a laminar flow in a larger domain, say in a large pipe but far from any no slip walls, and I am going to put in some periodic geometries in). scale = 1 just gives me unphysical results. I know it might be more of a fluid mechanics problem but would appreciate any insight that you provides.

[dokken]

@Exon

t (to simulate a uniform part of a laminar flow in a larger domain, say in a large pipe but far from any no slip walls, and I am going to put in some periodic geometries in)

If it is laminar (Stokes flow) far away from any boundaries, and the pressure drop is known, can’t one analytically derive the fluid velocity?
As far as I can tell, these are the equations you want to solve
Find (u, p) in \Omega=[-12,24]\times[0,12] such that

\begin{align} -\mu \Delta u + \nabla p &= 0 &\qquad \text{in } \Omega\\ \nabla \cdot u &= 0 &\qquad \text{in } \partial \Omega\\ u(x, 0) &= u(x, 12)\\ -\mu\frac{\partial u}{\partial n} + p n &= p_{in}&\qquad \text{on } \partial \Omega_{in} \\ -\mu\frac{\partial u}{\partial n} + p n &= p_{out}&\qquad \text{on } \partial \Omega_{out} \end{align}

if there is no divergence of u, we can use the divergence theorem

\int_\Omega \nabla \cdot u~\mathrm{d}x=\int_{\partial \Omega} u\cdot n~\mathrm{d}s = 0

as you have periodic conditions on top and bottom meaning that the flux from left has to be the flux from right (when integrated).
As there is no internal forcing, and the pressure is uniform on the inlet and outlet, I see no way of the velocity changing direction (i.e. flow up or down).

These are just my thoughts regarding your problem, and not an analytical statement from a text-book, so take it with a pinch of salt.

[Exon]

Thanks again for your patience and detailed explaination! Sorry for the late reply. Since this is more of a fluid mechanics problem I went to one of my lecturers who suggested that periodic probably does not make sense for my specific problem. I have switched to a perfect slip BC and had some problem, but I think its better for me to create a separate post.