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