Convergence problem for poiseuille flow of stokes equations

[quinso]

Hello,
I’m having an issue with convergence when trying to solve the folowing stokes problem:
-\Delta u+ \nabla p=0~~ in \Omega
\nabla.u=0~~ in \Omega
u=(4y(1-y),0)~~ on \{0,1\} \times ]0,1[
u=(0,0)~~ on ]0,1[ \times\{0,1\}

I take as exact solutions:
u_{ex}=(4y(1-y),0) and p_{ex}=8(1-x)
I don’t know if i should also set the pressure in the boundary. I’m a begginer, Any help will be highly appreciated.
Thank you in advance,

Here is the code:

from dolfinx import plot

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

# Define boundary conditions

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

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


def ud_expression(x):
    return np.stack((4 * x[1] * (1.0 - x[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()

    #  Dirichlet Conditions
    u_d_1 = Function(W0)
    u_d_1.interpolate(ud_expression)
    facets = locate_entities_boundary(msh, 1, inflow_boundary)
    dofs = locate_dofs_topological((W.sub(0), W0), 1, facets)
    bc0 = dirichletbc(u_d_1, dofs, W.sub(0))
    u_d_2=Function(W0)
    facets_2 = locate_entities_boundary(msh, 1, wall_boundary)
    dofs_2 = locate_dofs_topological((W.sub(0), W0), 1, facets_2)
    bc1 = dirichletbc(u_d_2, dofs_2, W.sub(0))
    bcs = [bc0,bc1]
    
    # Define variational problem
    (u, p) = ufl.TrialFunctions(W)
    (v, q) = ufl.TestFunctions(W)
    f = dolfinx.fem.Constant(mesh, (0., 0.))
    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
    ArithmeticError
    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: (4*x[1]*(1-x[1]), np.zeros(x.shape[1])))
    # Define pressure exact solution
    p_exact = Function(W1)
    p_exact.interpolate(lambda x: 8*(1-x[0]))
    # 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")

[dokken]

You have not stated what the issue of convergence for your problem is. How does the convergence rates look like?

Note that for Poiseuille flow, it is usually common to specify the pressure difference on the in and outlet instead of prescribing velocity over the whole boundary.

[quinso]

the issue is that the error increases as as the mesh is refined. I should have the opposite.

[quinso]

Goog morning Mr. @dokken , could you please review my code?, as you can see in the graph i have attached, there is no convergence for velocity and pressure, even when i try to prescribe pressure and velocity in boundaries instead of only prescribing velocity in the whole boundary.
I feel like I should be missing something important here, but I can’t quite figure out what it is.
I appreciate your assistance, have a good day.

you find here the graph of my solution (velocity), wich is quite resonable:

[dokken]

I would suggest you have a look at the Navier-Stokes example available at: Test problem 1: Channel flow (Poiseuille flow) — FEniCSx tutorial

1. Prescribe pressure on inlet and outlet and let that determine the flow
2. Prescribe inlet, and wall conditions, but let the outlet have a natural condition.

Please also supply actual codes for the cases you have considered.
I would suggest you remove the error/convergence study for now, and consider the feasibility of a single solution, i.e. look at how the pressure and velocity field look like for a single solve at a given resolution (say 32x32 unit square).

[quinso]

OK, thank you so much for your help.

[quinso]

thank you for taking the time to answer my questions, i tried to consider pressure on the inlet and autlet boundaries as you suggested me. I’m still having a non sense solutions. comparing the exact solutions and numeric solutions i’m getting a huge difference.
here is my code:

from dolfinx import plot

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 pyvista
from dolfinx.plot import vtk_mesh
import warnings
warnings.filterwarnings("ignore")

# Define boundary conditions


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

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

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

def exact_velocity_expression(x):
    return np.vstack((4*x[1]*(1-x[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()

    #  Dirichlet Conditions

    # set p=8 at the inflow x=0
    facets2 = locate_entities_boundary(msh, 1, in_boundary)
    dofs2 = locate_dofs_topological(W1, 1, facets2)
    bc0 = dirichletbc(fem.Constant(msh, PETSc.ScalarType(8)), dofs2, W.sub(1))

    # set p=0 at the inflow x=1
    facets3 = locate_entities_boundary(msh, 1, out_boundary)
    dofs3= locate_dofs_topological(W1, 1, facets3)
    bc1= dirichletbc(fem.Constant(msh, PETSc.ScalarType(0)), dofs3, W.sub(1))
    
    # Set u=0 at the walls (y=0 and y=1)
    noslip =Function(W0)
    facets3 = locate_entities_boundary(msh, 1, wall_boundary)
    dofs3 = locate_dofs_topological((W.sub(0), W0), 1, facets3)
    bc2 = dirichletbc(noslip, dofs3, W.sub(0))

    bcs = [bc0,bc1,bc2]
    
    # Define variational problem
    (u, p) = ufl.TrialFunctions(W)
    (v, q) = ufl.TestFunctions(W)
    
    f = fem.Constant(msh, (PETSc.ScalarType(0), PETSc.ScalarType(0)))
    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()
    # Normalize p_ex over unit square
    ArithmeticError
    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(exact_velocity_expression)
    # Define pressure exact solution
    p_exact = Function(W1)
    p_exact.interpolate(lambda x: 8*(1-x[0]))

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

errors_L2_u = []
errors_L2_p = []
errors_H1_u = []
h = []
num_elements_list = [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)
    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
    print("dofs:",num_dofs_global,",","h=",1./i," ,", "Error_L2_volocity=",error_u,",","Error_L2_pressure=",error_p)

errors i got are :

graph of velocity exact_solution is:

graph of numeric solution is:

[dokken]

@quinso

a = form((inner(grad(u), grad(v)) + inner(p, div(v)) + inner(div(u), q)) * dx)

Your pressure terms have the wrong sign?
You should have both + on the viscous term and the pressure term, as you have integrated both by parts, and they start with opposite signs

[jackhale]

The pressure lives in L^2/R so you cannot set the pressure solution on the boundary. MUMPS should handle the null space with the options you are using. I don’t know what it does - you will need to check the MUMPS documentation.

On first glance I noticed two things:

1. Your pressure solution is only determined up to a constant. Perhaps you are comparing two pressure solutions in L^2 with two different additional constants?
2. Notice that your exact velocity solution is a P2 polynomial, and so it can be represented exactly in your P2 velocity space (this is why velocity error is ~1E-13). Same applies for the pressure in P1. So, you are not really testing the convergence rate with this choice of exact solution.