# Stokes equations for poiseuille flow

I want to solve Stokes’s equations by using poiseuille flow’s boundary conditions. I made many efforts but I could not succeed.
https://docs.fenicsproject.org/dolfinx/v0.7.1/python/demos/demo_stokes.html
The Above link is for Stokes equations with lid boundary conditions.
.
https://jsdokken.com/dolfinx-tutorial/chapter2/ns_code1.html#test-problem-1-channel-flow-poiseuille-flow
This link is for Navier-Stokes equations but it is for Planar Poiseuille Flow. Use these boundary conditions in Stokes equations.

Use 3x` encapsulation for any code, ie

`````````python

```
``````

import numpy as np

import ufl
from basix.ufl import element, mixed_element
from dolfinx import fem, la
from dolfinx.fem import (Constant, Function, dirichletbc,
extract_function_spaces, form, functionspace,
locate_dofs_topological, locate_dofs_geometrical)
from dolfinx.fem.petsc import assemble_matrix_block, assemble_vector_block
from dolfinx.io import XDMFFile
from dolfinx.mesh import CellType, create_rectangle, locate_entities_boundary
from ufl import div, dx, grad, inner

from mpi4py import MPI
from petsc4py import PETSc

# Create mesh

msh = create_rectangle(MPI.COMM_WORLD, [np.array([0, 0]), np.array([1, 1])],
[32, 32], CellType.triangle)

# locating geometrically subsets of the boundary, and define a function

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)

# Function to mark x = 0, x = 1 and y = 0

def noslip_boundary(x):
return np.logical_or(np.isclose(x[1], 0.0), np.isclose(x[1], 1.0))

# No-slip condition on boundaries where u = 0, y = 1, and y = 0

noslip_dofs = locate_dofs_geometrical(V, noslip_boundary)
noslip = np.zeros(msh.geometry.dim, dtype=PETSc.ScalarType) # type: ignore
bc_noslip = dirichletbc(noslip, noslip_dofs, V)

# Plannar velocity condition for outflow (x = 0)

def inflow(x):
return np.isclose(x[0], 0)

inflow_dofs = locate_dofs_geometrical(Q , inflow)
bc_inflow = dirichletbc(PETSc.ScalarType(8), inflow_dofs, Q)

# Plannar velocity condition for outflow (x = 1)

def outflow(x):
return np.isclose(x[0], 1)

outflow_dofs = locate_dofs_geometrical(Q, outflow)
bc_outflow = dirichletbc(PETSc.ScalarType(0), outflow_dofs, Q)
bcu = [bc_noslip]
bcp = [bc_inflow, bc_outflow]

# Define variational problem

(u, p) = ufl.TrialFunction(V), ufl.TrialFunction(Q)
(v, q) = ufl.TestFunction(V), ufl.TestFunction(Q)
f = Constant(msh, (PETSc.ScalarType(0),PETSc.ScalarType(0))) # type: ignore

[inner(div(u), q) * dx, None]])
L = form([inner(f, v) * dx]) # type: ignore

# solvers for this problem:

a_p11 = form(inner(p, q) * dx)
a_p = [[a[0][0], None],
[None, a_p11]]

# value of 1 will be set on the diagonal for these rows.

def nested_iterative_solver():
“”“Solve the Stokes problem using nest matrices and an iterative solver.”“”

``````# Assemble nested matrix operators
A = fem.petsc.assemble_matrix_nest(a, bcu)
A.assemble()

u_ = Function(V)

# Create a nested matrix P to use as the preconditioner. The
# top-left block of P is shared with the top-left block of A. The
# bottom-right diagonal entry is assembled from the form a_p11:
P11 = fem.petsc.assemble_matrix(a_p11, [])
P = PETSc.Mat().createNest([[A.getNestSubMatrix(0, 0), None], [None, P11]])
P.assemble()

A00 = A.getNestSubMatrix(0, 0)
A00.setOption(PETSc.Mat.Option.SPD, True)

P00, P11 = P.getNestSubMatrix(0, 0), P.getNestSubMatrix(1, 1)
P00.setOption(PETSc.Mat.Option.SPD, True)
P11.setOption(PETSc.Mat.Option.SPD, True)

# Assemble right-hand side vector
b = fem.petsc.assemble_vector_nest(L)

# Modify ('lift') the RHS for Dirichlet boundary conditions
fem.petsc.apply_lifting_nest(b, a, bcs=bcu)

# Sum contributions for vector entries that are share across
# parallel processes
for b_sub in b.getNestSubVecs():

# Set Dirichlet boundary condition values in the RHS vector
bcs0 = fem.bcs_by_block(extract_function_spaces(L), bcs=bcu)
fem.petsc.set_bc_nest(b, bcs0)
``````

# Create a null vector for the velocity and pressure space

null_vec_u = fem.petsc.create_vector(V.dofmap.index_map, V.dofmap.index_map_bs)
null_vec_p = fem.petsc.create_vector(Q.dofmap.index_map, Q.dofmap.index_map_bs)

# Set velocity part to zero and the pressure part to a non-zero constant

null_vec_u.set(0.0)
null_vec_p.set(1.0)

# Combine the null vectors into a nested vector

null_vec = PETSc.Vec().createNest([null_vec_u, null_vec_p])

# When setting up the KSP solver, add the null space

ksp.setNullSpace(PETSc.NullSpace().create(vectors=[null_vec], comm=PETSc.COMM_WORLD))

``````# nested vector and the system is solved.
``````

u, p = Function(V), Function(Q)
x = PETSc.Vec().createNest([la.create_petsc_vector_wrap(u.x),
la.create_petsc_vector_wrap(p.x)])
ksp.solve(b, x)

``````# Save solution to file in XDMF format for visualization, e.g. with
# ParaView. Before writing to file, ghost values are updated using
# scatter_forward.
``````

with XDMFFile(MPI.COMM_WORLD, “out_stokes/velocity7.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)
ufile_xdmf.write_mesh(msh)
ufile_xdmf.write_function(u1)

with XDMFFile(MPI.COMM_WORLD, “out_stokes/pressure7.xdmf”, “w”) as pfile_xdmf:
p.x.scatter_forward()
pfile_xdmf.write_mesh(msh)
pfile_xdmf.write_function(p)

``````# Compute norms of the solution vectors
``````

norm_u = u.x.norm()
norm_p = p.x.norm()
if MPI.COMM_WORLD.rank == 0:
print(f"(A) Norm of velocity coefficient vector (nested, iterative): {norm_u}“)
print(f”(A) Norm of pressure coefficient vector (nested, iterative): {norm_p}")

``````return norm_u, norm_p
``````

# Solve using PETSc MatNest

norm_u_0, norm_p_0 = nested_iterative_solver()

# Solve using PETSc block matrices and an iterative solver

norm_u_1, norm_p_1 = block_iterative_solver()
assert np.isclose(norm_u_1, norm_u_0)
assert np.isclose(norm_p_1, norm_p_0)

# Solve using PETSc block matrices and an LU solver

norm_u_2, norm_p_2 = block_direct_solver()
assert np.isclose(norm_u_2, norm_u_0)
assert np.isclose(norm_p_2, norm_p_0)

# Solve using a non-blocked matrix and an LU solver

norm_u_3, norm_p_3 = mixed_direct()
assert np.isclose(norm_u_3, norm_u_0)
‘’’

First of all, you have not used ` for the encapsulation of the code.
Secondly, please reduce it to just use one of the four methods from the demo.

Thirdly, please provide an error message or a visualization of the unexpected result

‘’'python

# The required modules are first imported:

import numpy as np

import ufl
from basix.ufl import element, mixed_element
from dolfinx import fem, la
from dolfinx.fem import (Constant, Function, dirichletbc,
extract_function_spaces, form, functionspace,
locate_dofs_topological)
from dolfinx.fem.petsc import assemble_matrix_block, assemble_vector_block
from dolfinx.io import XDMFFile
from dolfinx.mesh import CellType, create_rectangle, locate_entities_boundary
from ufl import div, dx, grad, inner

from mpi4py import MPI
from petsc4py import PETSc

# Create mesh

msh = create_rectangle(MPI.COMM_WORLD, [np.array([0, 0]), np.array([1, 1])],
[32, 32], CellType.triangle)

# Function to mark x = 0, x = 1 and y = 0

def noslip_boundary(x):
return np.logical_or(np.isclose(x[1], 0.0), np.isclose(x[1], 1))

# Function to mark the lid (y = 1)

def lid(x):
U_max = 1
y = x[1]
return 4U_maxy*(1-y)

# Lid velocity

def lid_velocity_expression(x):
return np.stack((np.ones(x.shape[1]), np.zeros(x.shape[1])))

# piecewise linear basis (scalar).

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)

# No-slip condition on boundaries where x = 0, x = 1, and y = 0

noslip = np.zeros(msh.geometry.dim, dtype=PETSc.ScalarType) # type: ignore
facets = locate_entities_boundary(msh, 1, noslip_boundary)
bc0 = dirichletbc(noslip, locate_dofs_topological(V, 1, facets), V)

# Driving (lid) velocity condition on top boundary (y = 1)

lid_velocity = Function(V)
lid_velocity.interpolate(lid_velocity_expression)
facets = locate_entities_boundary(msh, 1, lid)
bc1 = dirichletbc(lid_velocity, locate_dofs_topological(V, 1, facets))

bcs = [bc0, bc1]

# Define variational problem

(u, p) = ufl.TrialFunction(V), ufl.TrialFunction(Q)
(v, q) = ufl.TestFunction(V), ufl.TestFunction(Q)
f = Constant(msh, (PETSc.ScalarType(0), PETSc.ScalarType(0))) # type: ignore

[inner(div(u), q) * dx, None]])
L = form([inner(f, v) * dx, inner(Constant(msh, PETSc.ScalarType(0)), q) * dx]) # type: ignore

# (non-nested) matrices and a direct (LU) solver.

def mixed_direct():

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

# No slip boundary condition
W0, _ = W.sub(0).collapse()
noslip = Function(W0)
facets = locate_entities_boundary(msh, 1, noslip_boundary)
dofs = locate_dofs_topological((W.sub(0), W0), 1, facets)
bc0 = dirichletbc(noslip, dofs, W.sub(0))

# Driving velocity condition u = (1, 0) on top boundary (y = 1)
lid_velocity = Function(W0)
lid_velocity.interpolate(lid_velocity_expression)
facets = locate_entities_boundary(msh, 1, lid)
dofs = locate_dofs_topological((W.sub(0), W0), 1, facets)
bc1 = dirichletbc(lid_velocity, dofs, W.sub(0))

# Collect Dirichlet boundary conditions
bcs = [bc0, bc1]

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

# Assemble LHS matrix and RHS vector
A = fem.petsc.assemble_matrix(a, bcs=bcs)
A.assemble()
b = fem.petsc.assemble_vector(L)

fem.petsc.apply_lifting(b, [a], bcs=[bcs])

# Set Dirichlet boundary condition values in the RHS
fem.petsc.set_bc(b, bcs)

# Create and configure solver
ksp = PETSc.KSP().create(msh.comm)
ksp.setOperators(A)
ksp.setType("preonly")

# Configure MUMPS to handle pressure nullspace
pc = ksp.getPC()
pc.setType("lu")
pc.setFactorSolverType("mumps")
pc.setFactorSetUpSolverType()
pc.getFactorMatrix().setMumpsIcntl(icntl=24, ival=1)
pc.getFactorMatrix().setMumpsIcntl(icntl=25, ival=0)

# Compute the solution
U = Function(W)
try:
ksp.solve(b, U.vector)
except PETSc.Error as e:
if e.ierr == 92:
print("The required PETSc solver/preconditioner is not available. Exiting.")
print(e)
exit(0)
else:
raise e

# Split the mixed solution and collapse
u, p = U.sub(0).collapse(), U.sub(1).collapse()

# Compute norms
norm_u, norm_p = u.x.norm(), p.x.norm()
if MPI.COMM_WORLD.rank == 0:
print(f"(D) Norm of velocity coefficient vector (monolithic, direct): {norm_u}")
print(f"(D) Norm of pressure coefficient vector (monolithic, direct): {norm_p}")

return norm_u, norm_u
``````

# Solve using a non-blocked matrix and an LU solver

norm_u_3, norm_p_3 = mixed_direct()
assert np.isclose(norm_u_3, norm_u_3)

u, p = Function(V), Function(Q)
x = PETSc.Vec().createNest([la.create_petsc_vector_wrap(u.x),
la.create_petsc_vector_wrap(p.x)])

with XDMFFile(MPI.COMM_WORLD, “out_stokes/velocity81.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)
ufile_xdmf.write_mesh(msh)
ufile_xdmf.write_function(u1)

with XDMFFile(MPI.COMM_WORLD, “out_stokes/pressure81.xdmf”, “w”) as pfile_xdmf:
p.x.scatter_forward()
pfile_xdmf.write_mesh(msh)
pfile_xdmf.write_function(p)
‘’’
I want to save XDMF files for Poiseuille flow for Stokes Equation. Now the Code is running smoothly, But need suggestions to do some changes in this CODE.

Again, you are using the wrong back-ticks. Please copy those that I used above.
You also did not reduce it to a single version of the solving.

Thirdly, you need to make an effort to impose the bcs yourself.

I want to save XDMF files for Poiseuille flow for Stokes Equation. Now the Code is running smoothly, But need suggestions to make some changes in this CODE.
I made a lot of effort but I thing I am making some mistakes.

`````````python
# The required modules are first imported:

import numpy as np

import ufl
from basix.ufl import element, mixed_element
from dolfinx import fem, la
from dolfinx.fem import (Constant, Function, dirichletbc,
extract_function_spaces, form, functionspace,
locate_dofs_topological)
from dolfinx.fem.petsc import assemble_matrix_block, assemble_vector_block
from dolfinx.io import XDMFFile
from dolfinx.mesh import CellType, create_rectangle, locate_entities_boundary
from ufl import div, dx, grad, inner

from mpi4py import MPI
from petsc4py import PETSc

# We create a {py:class}Mesh <dolfinx.mesh.Mesh>, define functions for
# locating geometrically subsets of the boundary, and define a function
# for the  velocity on the lid:

# +
# Create mesh
msh = create_rectangle(MPI.COMM_WORLD, [np.array([0, 0]), np.array([1, 1])],
[32, 32], CellType.triangle)

# Function to mark x = 0, x = 1 and y = 0
def noslip_boundary(x):
return np.logical_or(np.isclose(x[1], 0.0), np.isclose(x[1], 1))

# Function to mark the inflow ()
def inflow(x):
U_max = 1
y = x[1]
return 4*U_max*y*(1-y)

# Lid velocity
def parabolic_inflow_velocity(x):
U_max = 1
H = 1
return np.stack((np.stack((4*U_max * x[1] * (H - x[1])/H**2, np.zeros(x.shape[1])))))
# -

# Two {py:class}function spaces <dolfinx.fem.FunctionSpaceBase> are
# defined using different finite elements. 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)
V, Q = functionspace(msh, P2), functionspace(msh, P1)

# Boundary conditions for the velocity field are defined:

# No-slip condition on boundaries where x = 0, x = 1, and y = 0
noslip = np.zeros(msh.geometry.dim, dtype=PETSc.ScalarType)   # type: ignore
facets = locate_entities_boundary(msh, 1, noslip_boundary)
bc0 = dirichletbc(noslip, locate_dofs_topological(V, 1, facets), V)

# Driving (lid) velocity condition on top boundary (y = 1)
inflow_velocity = Function(V)
inflow_velocity.interpolate(parabolic_inflow_velocity)
facets = locate_entities_boundary(msh, 1, inflow)
bc1 = dirichletbc(inflow_velocity, locate_dofs_topological(V, 1, facets))

# Collect Dirichlet boundary conditions
bcs = [bc0, bc1]
# -

# The bilinear and linear forms for the Stokes equations are defined
# using a a blocked structure:

# +
# Define variational problem
(u, p) = ufl.TrialFunction(V), ufl.TrialFunction(Q)
(v, q) = ufl.TestFunction(V), ufl.TestFunction(Q)
f = Constant(msh, (PETSc.ScalarType(0), PETSc.ScalarType(0)))  # type: ignore

[inner(div(u), q) * dx, None]])
L = form([inner(f, v) * dx, inner(Constant(msh, PETSc.ScalarType(0)), q) * dx])  # type: ignore
# -

# We now solve the same Stokes problem again, but using monolithic
# (non-nested) matrices and a direct (LU) solver.

def mixed_direct():

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

# No slip boundary condition
W0, _ = W.sub(0).collapse()
noslip = Function(W0)
facets = locate_entities_boundary(msh, 1, noslip_boundary)
dofs = locate_dofs_topological((W.sub(0), W0), 1, facets)
bc0 = dirichletbc(noslip, dofs, W.sub(0))

# Driving velocity condition u = (1, 0) on top boundary (y = 1)
inflow_velocity = Function(W0)
inflow_velocity.interpolate(parabolic_inflow_velocity)
facets = locate_entities_boundary(msh, 1, inflow)
dofs = locate_dofs_topological((W.sub(0), W0), 1, facets)
bc1 = dirichletbc(inflow_velocity, dofs, W.sub(0))

# Collect Dirichlet boundary conditions
bcs = [bc0, bc1]

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

# Assemble LHS matrix and RHS vector
A = fem.petsc.assemble_matrix(a, bcs=bcs)
A.assemble()
b = fem.petsc.assemble_vector(L)

fem.petsc.apply_lifting(b, [a], bcs=[bcs])

# Set Dirichlet boundary condition values in the RHS
fem.petsc.set_bc(b, bcs)

# Create and configure solver
ksp = PETSc.KSP().create(msh.comm)
ksp.setOperators(A)
ksp.setType("preonly")

# Configure MUMPS to handle pressure nullspace
pc = ksp.getPC()
pc.setType("lu")
pc.setFactorSolverType("mumps")
pc.setFactorSetUpSolverType()
pc.getFactorMatrix().setMumpsIcntl(icntl=24, ival=1)
pc.getFactorMatrix().setMumpsIcntl(icntl=25, ival=0)

# Compute the solution
U = Function(W)
try:
ksp.solve(b, U.vector)
except PETSc.Error as e:
if e.ierr == 92:
print("The required PETSc solver/preconditioner is not available. Exiting.")
print(e)
exit(0)
else:
raise e

# Split the mixed solution and collapse
u, p = U.sub(0).collapse(), U.sub(1).collapse()

# Compute norms
norm_u, norm_p = u.x.norm(), p.x.norm()
if MPI.COMM_WORLD.rank == 0:
print(f"(D) Norm of velocity coefficient vector (monolithic, direct): {norm_u}")
print(f"(D) Norm of pressure coefficient vector (monolithic, direct): {norm_p}")

return norm_u, norm_u

# Solve using a non-blocked matrix and an LU solver
norm_u_3, norm_p_3 = mixed_direct()
assert np.isclose(norm_u_3, norm_u_3)

u, p = Function(V), Function(Q)
x = PETSc.Vec().createNest([la.create_petsc_vector_wrap(u.x),
la.create_petsc_vector_wrap(p.x)])

with XDMFFile(MPI.COMM_WORLD, "out_stokes/velocity81.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)
ufile_xdmf.write_mesh(msh)
ufile_xdmf.write_function(u1)

with XDMFFile(MPI.COMM_WORLD, "out_stokes/pressure81.xdmf", "w") as pfile_xdmf:
p.x.scatter_forward()
pfile_xdmf.write_mesh(msh)
pfile_xdmf.write_function(p)

``````

@dokken I need your help. Please give me some suggestions and tell me where am I making a mistake.