I have issues with convergence and bounadry recognition

[quinso]

Hello,everyone
I’m trying to solve the stokes problem on the unite square \Omega={]0,1[}^2 (with only Dirichlet conditions) :

-\Delta u+ \nabla p=f~~~~ in ~\Omega
\nabla.u=0~~~~ in ~\Omega
u=0~~~ in \partial \Omega

such that:
f(x,y)= -\Delta u_{ex}+\nabla p_{ex}
u_{ex}=(-256y(y − 1)(2y − 1)x^2(x − 1)^2,256x(x − 1)(2x − 1)y^2(y − 1)^2)
p_{ex}= (x-0.5)(y-0.5)

The code doesn’t return any errors, but from visualizing the solutions (velocity and pressure), i noticed that :

• the pressure doesn’t seem to converge to the exact solution.
• the Dirichlet conditions are not reconized, since the velocity is not nul in the boundaries of the domain (as shown in the last figure).

Could you help me please to fixe these problems. Thank you in advance.

Here is my code :

from dolfinx import plot
from mpi4py import MPI
from petsc4py import PETSc
import numpy as np
import ufl
from basix.ufl import element, mixed_element
from dolfinx import mesh, fem, plot, io
from dolfinx.fem import (Function, dirichletbc, form, functionspace,
                         locate_dofs_topological)
from dolfinx.fem.petsc import assemble_matrix_block, assemble_vector_block
from dolfinx.mesh import CellType, create_rectangle, locate_entities_boundary
from ufl import div, dx, grad, inner
import matplotlib.pyplot as plt
import pyvista
from dolfinx.plot import vtk_mesh
import warnings
warnings.filterwarnings("ignore")

# Define boundary conditions
def noslip_boundary(x):
    return np.logical_or(np.logical_or(np.isclose(x[0], 0.0), np.isclose(x[0], 1.0)),np.logical_or(np.isclose(x[1], 0.0), np.isclose(x[1], 1.0)))
def u_expression(x):
    return np.stack((np.zeros(x.shape[1]), np.zeros(x.shape[1])))

def mixed_direct(TH,msh):

    # Define function spaces
    W = functionspace(msh, TH)
    W0, _ = W.sub(0).collapse()
    W1,_ = W.sub(1).collapse()
    u_d=Function(W0)

    # Dirichlet Conditions
    u_d.interpolate(u_expression)
    facets = locate_entities_boundary(msh, 1, noslip_boundary)
    dofs = locate_dofs_topological((W.sub(0), W0), 1, facets)
    bc0 = dirichletbc(u_d, dofs, W.sub(0))
    bcs = [bc0]

    # Define variational problem
    (u, p) = ufl.TrialFunctions(W)
    (v, q) = ufl.TestFunctions(W)
    u_exact = Function(W0)
    p_exact=Function(W1)
    p_exact.interpolate(lambda x: (x[0]-0.5)*(x[1]-0.5))
    u_exact.interpolate(lambda x: (-256*x[1]*(x[1]-1)*(2*x[1]-1)*(x[0]**2)*(x[0]-1)**2, 256*x[0]*(x[0]-1)* 
   (2*x[0]-1)*(x[1]**2)*(x[1]-1)**2))
    f=-div(grad(u_exact))+grad(p_exact)
    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
    ksp.solve(b, U.vector)
    
    # Split the mixed solution and collapse
    u, p = U.sub(0).collapse(), U.sub(1).collapse()
    u.x.scatter_forward()
    p.x.scatter_forward()
    return u,p

def compute_errors(u,p,msh,TH):
    W = functionspace(msh, TH)
    W0, _ = W.sub(0).collapse()
    W1,_ = W.sub(1).collapse()
    # Define velocity exact solution
    u_exact = Function(W0)
    u_exact.interpolate(lambda x: (-256*x[1]*(x[1]-1)*(2*x[1]-1)*(x[0]**2)*(x[0]-1)**2, 256*x[0]*(x[0]-1)* 
    (2*x[0]-1)*(x[1]**2)*(x[1]-1)**2))
    #Define pressure exact solution
    p_exact=Function(W1)
    p_exact.interpolate(lambda x: (x[0]-0.5)*(x[1]-0.5))

    # L2 error for Velocity
    L2_error = fem.form(ufl.inner(u - u_exact, u- u_exact) * ufl.dx)
    error_local = fem.assemble_scalar(L2_error)
        error_L2 = np.sqrt(msh.comm.allreduce(error_local, op=MPI.SUM))

     # H1 error for Velocity
    H1_u_1 = inner(u - u_exact, u- u_exact) * ufl.dx
    H1_u_2 = inner(grad(u - u_exact), grad(u- u_exact)) * ufl.dx
    H1_u=fem.form(H1_u_1+H1_u_2)
    H1_error=fem.assemble_scalar(H1_u)
    error_H1= np.sqrt(msh.comm.allreduce(H1_error, op=MPI.SUM))

     # L2 error for presure
    L2_error_2 = fem.form(ufl.inner(-p - p_exact, -p- p_exact) * ufl.dx)
    error_local_2 = fem.assemble_scalar(L2_error_2)
    error_L2_2 = np.sqrt(msh.comm.allreduce(error_local_2, op=MPI.SUM))
    return error_L2,error_H1,error_L2_2

#Convergence curve for velocity

errors_L2_u = [ ]
errors_L2_p=[ ]
errors_H1_u=[ ]
h=[]
num_elements_list = [4,8,16, 32,64]
for i in num_elements_list:
    mesh = create_rectangle(MPI.COMM_WORLD, ((0.0, 0.0), (1.0, 1.0)), (i,i), CellType.triangle)
    P2 = element("Lagrange", mesh.basix_cell(), 2, shape=(mesh.geometry.dim,))
    P1 = element("Lagrange", mesh.basix_cell(), 1)
    TH = mixed_element([P2, P1])
    u,p=mixed_direct(TH,mesh)
    error_u,error_H1,error_p=compute_errors(u,p,mesh,TH)
    errors_L2_u.append(error_u)
    errors_L2_p.append(error_p)
        errors_H1_u.append(error_H1)
    h.append(1./i)
    V = fem.FunctionSpace(mesh, P2)
    num_dofs_local = (V.dofmap.index_map.size_local) * V.dofmap.index_map_bs
    num_dofs_global = V.dofmap.index_map.size_global * V.dofmap.index_map_bs
fig = plt.figure(figsize=(14, 6))
plt.subplot(1, 2, 1)
legend=[]
plt.plot(h, errors_L2_u, "bo-")
plt.plot(h, errors_H1_u, "ro-")
legend += [r"$L^2(\mathbf{u_h}-\mathbf{u}_{ex})$",r"$H^1(\mathbf{u_h}-\mathbf{u}_{ex})$"]
plt.legend(legend)
plt.xlim(plt.xlim()[::-1])
plt.xlabel("h")
plt.ylabel("error")
plt.title("convergence curve for velocity")
plt.subplot(1, 2, 2)
legend=[]
plt.plot(h, errors_L2_u, "bo-")
plt.plot(h, errors_H1_u, "ro-")
legend += [r"$L^2(\mathbf{u_h}-\mathbf{u}_{ex})$",r"$H^1(\mathbf{u_h}-\mathbf{u}_{ex})$"]
plt.legend(legend)
plt.xscale("log")
plt.yscale("log")
plt.xlabel("log(h)")
plt.ylabel("log(error)")
plt.axis("equal")
plt.title("conergence curve for velocity")
plt.show()

# Convergence curve for pressure
fig = plt.figure(figsize=(14, 6))
plt.subplot(1, 2, 1)
legend=[]
plt.plot(h, errors_L2_p, "bo-")
legend += [r"$L^2(\mathbf{p_h}-\mathbf{p}_{ex})$"]
plt.legend(legend)
plt.xlim(plt.xlim()[::-1])
plt.xlabel("h")
plt.ylabel("error")
plt.title("convergence curve for pressure")
plt.subplot(1, 2, 2)
legend=[]
plt.plot(h, errors_L2_p, "ro-")
legend += [r"$L^2(\mathbf{p_h}-\mathbf{p}_{ex})$"]
plt.legend(legend)
plt.xscale("log")
plt.yscale("log")
plt.xlabel("log(h)")
plt.ylabel("log(error)")
plt.axis("equal")
plt.title("convergence curve for pressure")
plt.show()

#Visualizing the numeric solution

mesh = create_rectangle(MPI.COMM_WORLD, ((0.0, 0.0), (1.0, 1.0)), (8,8), CellType.triangle)
P2 = element("Lagrange", mesh.basix_cell(), 2, shape=(mesh.geometry.dim,))
P1 = element("Lagrange", mesh.basix_cell(), 1)
V = fem.FunctionSpace(mesh, P2)
TH = mixed_element([P2, P1])
u_h,p_h=mixed_direct(TH,mesh)
pyvista.start_xvfb()
topology, cell_types, geometry = vtk_mesh(V)
values = np.zeros((geometry.shape[0], 3), dtype=np.float64)
values[:, :len(u_h)] = u_h.x.array.real.reshape((geometry.shape[0], len(u_h)))

# Create a point cloud of glyphs
function_grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)
function_grid["u_h"] = values
glyphs = function_grid.glyph(orient="u_h", factor=0.2)

# Create a pyvista-grid for the mesh
grid = pyvista.UnstructuredGrid(*vtk_mesh(mesh, mesh.topology.dim))

# Create plotter
plotter = pyvista.Plotter()
plotter.add_mesh(grid, style="wireframe", color="k")
plotter.add_mesh(glyphs)
plotter.view_xy()
if not pyvista.OFF_SCREEN:
plotter.show()
else:
fig_as_array = plotter.screenshot("glyphs.png")

The convergence_curve for velocity was good ( i got a line, the convergence order is close to 2)
as you can see the velocity is not nul in the bounadries.

[dokken]

Please note that you need to set the pressure null space, as the pressure is only uniquely determined up to a constant when you prescribe pure Dirichlet conditions to the velocity on the whole boundary.

Here is a minimal working example (using DOLFINx v0.7.x):

from mpi4py import MPI
from petsc4py import PETSc
import numpy as np
import ufl
import dolfinx.fem.petsc
import dolfinx.io
from basix.ufl import element, mixed_element
from dolfinx import mesh, fem
from dolfinx.fem import (Function, dirichletbc, form, functionspace,
                         locate_dofs_topological)
from dolfinx.mesh import CellType, create_rectangle, locate_entities_boundary
from ufl import div, dx, grad, inner
import matplotlib.pyplot as plt
import warnings
warnings.filterwarnings("ignore")

# Define boundary conditions


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


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


def mixed_direct(TH, msh):

    # Define function spaces
    W = functionspace(msh, TH)
    W0, _ = W.sub(0).collapse()
    W1, W1_to_W = W.sub(1).collapse()
    u_d = Function(W0)

    # Dirichlet Conditions
    u_d.interpolate(u_expression)
    facets = locate_entities_boundary(msh, 1, noslip_boundary)
    dofs = locate_dofs_topological((W.sub(0), W0), 1, facets)
    bc0 = dirichletbc(u_d, dofs, W.sub(0))
    bcs = [bc0]

    # Define variational problem
    (u, p) = ufl.TrialFunctions(W)
    (v, q) = ufl.TestFunctions(W)
    u_exact = Function(W0)
    p_exact = Function(W1)
    p_exact.interpolate(lambda x: (x[0]-0.5)*(x[1]-0.5))
    u_exact.interpolate(lambda x: (-256*x[1]*(x[1]-1)*(2*x[1]-1)*(x[0]**2)*(x[0]-1)**2, 256*x[0]*(x[0]-1) *
                                   (2*x[0]-1)*(x[1]**2)*(x[1]-1)**2))
    f = -div(grad(u_exact))+grad(p_exact)
    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")

    null_vec = dolfinx.fem.Function(W)
    null_vec.x.array[W1_to_W] = 1.0
    dolfinx.la.orthonormalize([null_vec.vector])
    nsp = PETSc.NullSpace().create(vectors=[null_vec.vector])
    A.setNearNullSpace(nsp)
    A.setNullSpace(nsp)
    assert nsp.test(A)

    # 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
    ksp.solve(b, U.vector)

    # Split the mixed solution and collapse
    u, p = U.sub(0).collapse(), U.sub(1).collapse()
    u.x.scatter_forward()
    p.x.scatter_forward()
    # Normalize p_ex over unit square
    p_avg = W.mesh.comm.allreduce(
        fem.assemble_scalar(fem.form(p*dx)), op=MPI.SUM)
    p.x.array[:] -= p_avg

    return u, p


def compute_errors(u, p, msh, TH):
    W = functionspace(msh, TH)
    W0, _ = W.sub(0).collapse()
    W1, _ = W.sub(1).collapse()
    # Define velocity exact solution
    u_exact = Function(W0)
    u_exact.interpolate(lambda x: (-256*x[1]*(x[1]-1)*(2*x[1]-1)*(x[0]**2)*(x[0]-1)**2, 256*x[0]*(x[0]-1) *
                                   (2*x[0]-1)*(x[1]**2)*(x[1]-1)**2))
    # Define pressure exact solution
    p_exact = Function(W1)
    p_exact.interpolate(lambda x: (x[0]-0.5)*(x[1]-0.5))

    # L2 error for Velocity
    L2_error = fem.form(ufl.inner(u - u_exact, u - u_exact) * ufl.dx)
    error_local = fem.assemble_scalar(L2_error)
    error_L2 = np.sqrt(msh.comm.allreduce(error_local, op=MPI.SUM))

    # H1 error for Velocity
    H1_u_1 = inner(u - u_exact, u - u_exact) * ufl.dx
    H1_u_2 = inner(grad(u - u_exact), grad(u - u_exact)) * ufl.dx
    H1_u = fem.form(H1_u_1+H1_u_2)
    H1_error = fem.assemble_scalar(H1_u)
    error_H1 = np.sqrt(msh.comm.allreduce(H1_error, op=MPI.SUM))

    # L2 error for presure
    L2_error_2 = fem.form(ufl.inner(-p - p_exact, -p - p_exact) * ufl.dx)
    error_local_2 = fem.assemble_scalar(L2_error_2)
    error_L2_2 = np.sqrt(msh.comm.allreduce(error_local_2, op=MPI.SUM))
    return error_L2, error_H1, error_L2_2

# Convergence curve for velocity


errors_L2_u = []
errors_L2_p = []
errors_H1_u = []
h = []
num_elements_list = [8, 16, 32, 64, 128, 256]
for i in num_elements_list:
    mesh = create_rectangle(
        MPI.COMM_WORLD, ((0.0, 0.0), (1.0, 1.0)), (i, i), CellType.triangle)
    P2 = element("Lagrange", mesh.basix_cell(), 2, shape=(mesh.geometry.dim,))
    P1 = element("Lagrange", mesh.basix_cell(), 1)
    TH = mixed_element([P2, P1])
    u, p = mixed_direct(TH, mesh)
    with dolfinx.io.VTXWriter(mesh.comm, f"u_{i}.bp", [u]) as bp:
        bp.write(0.0)
    with dolfinx.io.VTXWriter(mesh.comm, f"p_{i}.bp", [p]) as bp:
        bp.write(0.0)

    error_u, error_H1, error_p = compute_errors(u, p, mesh, TH)
    errors_L2_u.append(error_u)
    errors_L2_p.append(error_p)
    errors_H1_u.append(error_H1)
    h.append(1./i)
    V = fem.functionspace(mesh, P2)
    num_dofs_local = (V.dofmap.index_map.size_local) * V.dofmap.index_map_bs
    num_dofs_global = V.dofmap.index_map.size_global * V.dofmap.index_map_bs
fig = plt.figure(figsize=(14, 6))
plt.subplot(1, 2, 1)
legend = []
plt.plot(h, errors_L2_u, "bo-")
plt.plot(h, errors_H1_u, "ro-")
legend += [r"$L^2(\mathbf{u_h}-\mathbf{u}_{ex})$",
           r"$H^1(\mathbf{u_h}-\mathbf{u}_{ex})$"]
plt.legend(legend)
plt.xlim(plt.xlim()[::-1])
plt.xlabel("h")
plt.ylabel("error")
plt.title("convergence curve for velocity")
plt.subplot(1, 2, 2)
legend = []
plt.plot(h, errors_L2_u, "bo-")
plt.plot(h, errors_H1_u, "ro-")
legend += [r"$L^2(\mathbf{u_h}-\mathbf{u}_{ex})$",
           r"$H^1(\mathbf{u_h}-\mathbf{u}_{ex})$"]
plt.legend(legend)
plt.xscale("log")
plt.yscale("log")
plt.xlabel("log(h)")
plt.ylabel("log(error)")
plt.axis("equal")
plt.title("conergence curve for velocity")
plt.savefig("u_rate.png")

# Convergence curve for pressure
fig = plt.figure(figsize=(14, 6))
plt.subplot(1, 2, 1)
legend = []
plt.plot(h, errors_L2_p, "bo-")
legend += [r"$L^2(\mathbf{p_h}-\mathbf{p}_{ex})$"]
plt.legend(legend)
plt.xlim(plt.xlim()[::-1])
plt.xlabel("h")
plt.ylabel("error")
plt.title("convergence curve for pressure")
plt.subplot(1, 2, 2)
legend = []
plt.plot(h, errors_L2_p, "ro-")
legend += [r"$L^2(\mathbf{p_h}-\mathbf{p}_{ex})$"]
plt.legend(legend)
plt.xscale("log")
plt.yscale("log")
plt.xlabel("log(h)")
plt.ylabel("log(error)")
plt.axis("equal")
plt.title("convergence curve for pressure")
plt.savefig("p_rate.png")

The visualization is another issue, you should use the collapsed space of u_h, not V, as they might be ordered differently.
i.e.

@quinso

topology, cell_types, geometry = vtk_mesh(V)

should be
topology, cell_types, geometry = vtk_mesh(u_h.function_space)

[quinso]

Thank you so much Mr. dokken for your help, Now every thing is working very well.
Have a good day,

[quinso]

hi, i have another question, what can i do to visualize the exact solution with the numeric solution to compare between them?.
Should i interpolate the exact solution and do the same as i did with the numeric solution?
Thank you in advance.

[dokken]

Yes, that would be appropriate

[quinso]

Hello, I’m currently attempting to address the same stokes problem mentioned above. This time, I’m taking Neumann condition (\nabla u-pI).n=0 (the variationnal forme will not change) at the outlet and Dirichlet condition u=0 elsewhere. I used the code you previously supplied. However, I’m encountering an issue. Upon visualizing the numerical solution, I’ve observed that the velocity at the outlet does not match the expected values. Could you please help me in identifying the source of this problem?
Below is my code with the modifications:

from mpi4py import MPI
from petsc4py import PETSc
import numpy as np
import ufl
import dolfinx.fem.petsc
import dolfinx.io
from basix.ufl import element, mixed_element
from dolfinx import mesh, fem
from dolfinx.fem import (Function, dirichletbc, form, functionspace,
                         locate_dofs_topological)
from dolfinx.mesh import CellType, create_rectangle, locate_entities_boundary
from ufl import div, dx, grad, inner
import matplotlib.pyplot as plt
import warnings
warnings.filterwarnings("ignore")

# Define boundary conditions


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


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


def mixed_direct(TH, msh):

    # Define function spaces
    W = functionspace(msh, TH)
    W0, _ = W.sub(0).collapse()
    W1, W1_to_W = W.sub(1).collapse()
    u_d = Function(W0)

    # Dirichlet Conditions
    u_d.interpolate(u_expression)
    facets = locate_entities_boundary(msh, 1, noslip_boundary)
    dofs = locate_dofs_topological((W.sub(0), W0), 1, facets)
    bc0 = dirichletbc(u_d, dofs, W.sub(0))
    bcs = [bc0]

    # Define variational problem
    (u, p) = ufl.TrialFunctions(W)
    (v, q) = ufl.TestFunctions(W)
    u_exact = Function(W0)
    p_exact = Function(W1)
    p_exact.interpolate(lambda x: (x[0]-0.5)*(x[1]-0.5))
    u_exact.interpolate(lambda x: (-256*x[1]*(x[1]-1)*(2*x[1]-1)*(x[0]**2)*(x[0]-1)**2, 256*x[0]*(x[0]-1) *
                                   (2*x[0]-1)*(x[1]**2)*(x[1]-1)**2))
    f = -div(grad(u_exact))+grad(p_exact)
    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")

    null_vec = dolfinx.fem.Function(W)
    null_vec.x.array[W1_to_W] = 1.0
    dolfinx.la.orthonormalize([null_vec.vector])
    nsp = PETSc.NullSpace().create(vectors=[null_vec.vector])
    A.setNearNullSpace(nsp)
    A.setNullSpace(nsp)
    '''assert nsp.test(A)'''

    # 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
    ksp.solve(b, U.vector)

    # Split the mixed solution and collapse
    u, p = U.sub(0).collapse(), U.sub(1).collapse()
    u.x.scatter_forward()
    p.x.scatter_forward()
    # Normalize p_ex over unit square
    p_avg = W.mesh.comm.allreduce(
        fem.assemble_scalar(fem.form(p*dx)), op=MPI.SUM)
    p.x.array[:] -= p_avg

    return u, p

the graph of numeric solution is :

[dokken]

I am not sure what you expect as a solution to this problem.
Having the natural boundary condition does not seem appropriate in this case.
You are generating a flow by an internal force, and it is not clear to me that \frac{\partial u}{\partial n} - pn = 0 on the outlet with the chosen u_exact, p_exact, as just by doing quick maths on this, you would get something like

grad({-256 y (y - 1) (2 y - 1) x^2 (x - 1)^2, 256 x ((x - 1) (2 x - 1) y^2 (y - 1)^2)}.{1, 0}) = (512 (1 - 2 x) (-1 + x) x (-1 + y) y (-1 + 2 y), -256 (-1 + x)^2 x^2 (1 - 6 y + 6 y^2))

for the \frac{\partial u}{\partial n}, while pn is just ((x-0.5)*(y-0.5), 0)

[quinso]

Ok, i think i understand the problem, so if i use \frac{\partial u}{\partial n}-pn=g and calulate g that satisfy this boundary condition then the problem will be solved. However, this adjustment will change the variational form as well.
think you so much for you help,

[quinso]

"Hello, I’ve addressed the calculation problem as you suggested, but I’m still having the same graph as before, showing random flow on the outlet. My questions are:

• Is it normal for the flow to appear random on the outlet?
• Is it expected to have a lack of convergence in this scenario, particularly when applying Neumann conditions on the outlet?"
• Is there alternative way to properly handle Neumann condition for this Stokes problem?

Thank you in advance,
Here is my code:

from mpi4py import MPI
from petsc4py import PETSc
import numpy as np
import ufl
import dolfinx.fem.petsc
import dolfinx.io
from basix.ufl import element, mixed_element
from dolfinx import mesh, fem
from dolfinx.fem import (Function, dirichletbc, form, functionspace,
                         locate_dofs_topological)
from dolfinx.mesh import CellType, create_rectangle, locate_entities_boundary
from ufl import div, dx, grad, inner, ds
import matplotlib.pyplot as plt
import warnings
warnings.filterwarnings("ignore")

# Define boundary conditions


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


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


def mixed_direct(TH, msh):

    # Define function spaces
    W = functionspace(msh, TH)
    W0, _ = W.sub(0).collapse()
    W1, W1_to_W = W.sub(1).collapse()
    u_d = Function(W0)

    # Dirichlet Conditions
    u_d.interpolate(u_expression)
    facets = locate_entities_boundary(msh, 1, noslip_boundary)
    dofs = locate_dofs_topological((W.sub(0), W0), 1, facets)
    bc0 = dirichletbc(u_d, dofs, W.sub(0))
    bcs = [bc0]

    # Define variational problem
    (u, p) = ufl.TrialFunctions(W)
    (v, q) = ufl.TestFunctions(W)
    u_exact = Function(W0)
    p_exact = Function(W1)
    p_exact.interpolate(lambda x: (x[0]-0.5)*(x[1]-0.5))
    u_exact.interpolate(lambda x: (-256*x[1]*(x[1]-1)*(2*x[1]-1)*(x[0]**2)*(x[0]-1)**2, 256*x[0]*(x[0]-1) *
                                   (2*x[0]-1)*(x[1]**2)*(x[1]-1)**2))
    f = -div(grad(u_exact))+grad(p_exact)
    g = Function(W0)
    g.interpolate(lambda x: (-512*(x[1]-1)*(2*x[1]-1)*x[0]*(x[0]-1)*(2*x[0]-1)*x[1]-(x[0]-0.5)*(x[1]-0.5), 
    -256*(x[0]**2)*((x[0]-1)**2)*(1-6*x[1]+6*x[1]**2)))

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

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

    null_vec = dolfinx.fem.Function(W)
    null_vec.x.array[W1_to_W] = 1.0
    dolfinx.la.orthonormalize([null_vec.vector])
    nsp = PETSc.NullSpace().create(vectors=[null_vec.vector])
    A.setNearNullSpace(nsp)
    A.setNullSpace(nsp)
    '''assert nsp.test(A)'''

    # 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
    ksp.solve(b, U.vector)

    # Split the mixed solution and collapse
    u, p = U.sub(0).collapse(), U.sub(1).collapse()
    u.x.scatter_forward()
    p.x.scatter_forward()
    # Normalize p_ex over unit square
    p_avg = W.mesh.comm.allreduce(
        fem.assemble_scalar(fem.form(p*dx)), op=MPI.SUM)
    p.x.array[:] -= p_avg

    return u, p

[dokken]

You have messed up the signs in your variational form.
Here is a working example, where I rely on UFL to derive f and g, instead of hard-coding them:

from mpi4py import MPI
from petsc4py import PETSc
import numpy as np
import ufl
import dolfinx.fem.petsc
import dolfinx.io
from basix.ufl import element, mixed_element
from dolfinx import mesh, fem
from dolfinx.fem import (Function, dirichletbc, form, functionspace,
                         locate_dofs_topological)
from dolfinx.mesh import CellType, create_rectangle, locate_entities_boundary
from ufl import div, dx, grad, inner, ds
import matplotlib.pyplot as plt
import warnings
warnings.filterwarnings("ignore")

# Define boundary conditions


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


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


def mixed_direct(TH, msh):

    # Define function spaces
    W = functionspace(msh, TH)
    W0, _ = W.sub(0).collapse()
    W1, W1_to_W = W.sub(1).collapse()
    u_d = Function(W0)

    # Dirichlet Conditions
    u_d.interpolate(u_expression)
    facets = locate_entities_boundary(msh, 1, noslip_boundary)
    dofs = locate_dofs_topological((W.sub(0), W0), 1, facets)
    bc0 = dirichletbc(u_d, dofs, W.sub(0))
    bcs = [bc0]

    def u_ex(x):
        return (-256*x[1]*(x[1]-1)*(2*x[1]-1)*(x[0]**2)*(x[0]-1)**2, 256*x[0]*(x[0]-1) *
                                   (2*x[0]-1)*(x[1]**2)*(x[1]-1)**2)
    
    def p_ex(x):
        return (x[0]-0.5)*(x[1]-0.5)
    # Define variational problem
    (u, p) = ufl.TrialFunctions(W)
    (v, q) = ufl.TestFunctions(W)
    x = ufl.SpatialCoordinate(mesh)
    n = ufl.FacetNormal(mesh)
    mu = 1
    u_exact = ufl.as_vector((u_ex(x)))
    f = -mu*div(grad(u_exact)) + grad(p_ex(x))
    g = -mu * ufl.dot(ufl.grad(u_exact), n) + p_ex(x) * n

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

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

    null_vec = dolfinx.fem.Function(W)
    null_vec.x.array[W1_to_W] = 1.0
    dolfinx.la.orthonormalize([null_vec.vector])
    nsp = PETSc.NullSpace().create(vectors=[null_vec.vector])
    A.setNearNullSpace(nsp)
    A.setNullSpace(nsp)
    '''assert nsp.test(A)'''

    # 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
    ksp.solve(b, U.vector)

    # Split the mixed solution and collapse
    u, p = U.sub(0).collapse(), U.sub(1).collapse()
    u.x.scatter_forward()
    p.x.scatter_forward()
    # Normalize p_ex over unit square
    p_avg = W.mesh.comm.allreduce(
        fem.assemble_scalar(fem.form(p*dx)), op=MPI.SUM)
    p.x.array[:] -= p_avg

    return u, p

if __name__ == "__main__":
    mesh = create_rectangle(MPI.COMM_WORLD, ((0.0, 0.0), (1.0, 1.0)), (8,8), CellType.triangle)
    P2 = element("Lagrange", mesh.basix_cell(), 2, shape=(mesh.geometry.dim,))
    P1 = element("Lagrange", mesh.basix_cell(), 1)
    V = fem.functionspace(mesh, P2)
    TH = mixed_element([P2, P1])
    u_h,p_h=mixed_direct(TH,mesh)
    with dolfinx.io.VTXWriter(mesh.comm, "u.bp", [u_h], engine="BP4") as bp:
        bp.write(0.0)

Which returns

[quinso]

Thank you so much Sir for your help and you support, every thing works perfectly now. i just have a question: why do we have to take g with negative sign. in my problem i defined g= \frac{\partial u}{\partial n}-pn , shouldn’t we define g with g=\frac{\partial u}{\partial n}+pn (with only the sign of ‘p’ reversed)?. i dont unserstand why we should also change the sign of velocity.
Have a good day.

[dokken]

In the Stokes equation, the diffusive term and pressure gradient terms have opposite signs, which is also how the natural outlet bc is defined.