Stokes square cavity flow problem

dolfinx
,
1

Hello! Sorry to bother you! I am a novice. Recently, I tried to make Stokes square cavity flow, using P2-P1 yuan. I found that my speed calculation result reached the expectation, but the pressure was very bad, I looked up the relevant materials and still couldn’t solve it. Looking forward to your reply! Here’s my program.

from mpi4py import MPI
import numpy as np
from dolfinx import mesh
from dolfinx import fem
import ufl
from ufl import inner,grad,div,dx
from petsc4py import PETSc
from dolfinx.io import XDMFFile,VTKFile

import ufl
from petsc4py.PETSc import ScalarType
from basix.ufl import element, mixed_element
from dolfinx import default_real_type, fem, la, io, mesh
from dolfinx.fem import (
    Constant,
    Function,
    dirichletbc,
    extract_function_spaces,
    form,
    functionspace,
    locate_dofs_topological,
)
from dolfinx.fem.petsc import LinearProblem
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,create_unit_square
from ufl import *

msh = mesh.create_unit_square(MPI.COMM_WORLD, 40, 40, mesh.CellType.triangle)

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

# Function to mark the lid (y = 1)
def lid(x):
    return np.isclose(x[1], 1.0)

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

P2 = element("Lagrange", msh.basix_cell(), 2, shape=(msh.geometry.dim,), dtype=default_real_type)
P1 = element("Lagrange", msh.basix_cell(), 1, dtype=default_real_type)
V_el = mixed_element([P2, P1])
V, Q = functionspace(msh, P2), functionspace(msh, P1)
W = functionspace(msh, V_el)

# Create the function space
# TH = P2 * P1
# W = fem.FunctionSpace(msh, TH)
W0, _ = W.sub(0).collapse() #局部空间自由度索引   sub 表示全局自由度索引

# No slip boundary condition
noslip = fem.Function(W0)
facets = mesh.locate_entities_boundary(msh, 1, noslip_boundary)
dofs = fem.locate_dofs_topological((W.sub(0), W0), 1, facets) #返回自由度索引
bc0 = fem.dirichletbc(noslip, dofs, W.sub(0))

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

# Since for this problem the pressure is only determined up to a constant, we pin the pressure at the point (0, 0)
# zero = fem.Function(Q)
# zero.x.set(0.0)
# dofs = fem.locate_dofs_geometrical((W.sub(1), Q), lambda x: np.isclose(x.T, [0, 0, 0]).all(axis=1))
# bc2 = fem.dirichletbc(zero, dofs, W.sub(1))

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

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

# Compute solution
solver = fem.petsc.LinearProblem(
    a, L, bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu",
                              "pc_factor_mat_solver_type": "mumps"})
w_h = solver.solve()

# Split the mixed solution and collapse
u, p = w_h.split()

# Write the solution to file
with XDMFFile(MPI.COMM_WORLD, "stokesx/pressure.xdmf", "w") as pfile_xdmf:
    p.x.scatter_forward()
    pfile_xdmf.write_mesh(msh)
    pfile_xdmf.write_function(p)

with XDMFFile(MPI.COMM_WORLD, "stokesx/velocity.xdmf", "w") as ufile_xdmf:
    u.x.scatter_forward()
    P1 = element(
        "Lagrange", msh.basix_cell(), 1, shape=(msh.geometry.dim,), dtype=default_real_type
    )
    u1 = Function(functionspace(msh, P1))
    u1.interpolate(u)
    ufile_xdmf.write_mesh(msh)
    ufile_xdmf.write_function(u1)

print('Execution Complete')
2 
# Since for this problem the pressure is only determined up to a constant, we pin the pressure at the point (0, 0)

That’s a key point, you can’t just ignore it (and comment it out). If you don’t like pinning the pressure, search this forum with keyword “nullspace” for an alternative approach.

1 Like
3

Thank you very much for your reply. If I fix the pressure, the result is ok. I have found relevant documents on how to use nullspace.

I have another question about space, I would like to ask if Q in the scalar space Q = functionspace(msh, P1) is equal to W.sub (1)? In the community I found some commands (W.ub (1), Q) that map Q to W.sub (1), can you help me answer this question?

4

Yes, W.sub(1) represents the pressure subspace of the mixed function space (velocity, pressure). If you want to apply boundary conditions, the notation you used (W.sub(1), Q) is correct. If you need to get a standalone space you’ll need to use W.sub(1).collapse(), which you already used in your code.

1 Like
5

Thank you very much for your patient answer, send you the most sincere thanks!