Stokes equations for poiseuille flow

dolfinx
1

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.

2

Please add the code you have used to attempt to add these boundary conditions, as it increases the likelihood of anyone trying to help you.

Use 3x` encapsulation for any code, ie

```python
#add code here 

```
3

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

We create a {py:class}Mesh <dolfinx.mesh.Mesh>, define functions for

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).

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]

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

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

-

A block-diagonal preconditioner will be used with the iterative

solvers for this problem:

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

### Nested matrix solver

We assemble the bilinear form into a nested matrix A, and call the

assemble() method to communicate shared entries in parallel. Rows

and columns in A that correspond to degrees-of-freedom with

Dirichlet boundary conditions wil be zeroed by the assembler, and a

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():
    b_sub.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)

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

solution process.

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])

… [rest of your code] …

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

### Monolithic block iterative solver

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)
‘’’

4

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

5

‘’'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 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])))

-

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)

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

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

a = form([[inner(grad(u), grad(v)) * dx, inner(p, div(v)) * dx],
[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)
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])
b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)

# 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.

6

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.

7

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

a = form([[inner(grad(u), grad(v)) * dx, inner(p, div(v)) * dx],
          [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])
    b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)

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

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