Bad convergence of stokes problem

Errors
1

Hello everyone, it’s been months since I’ve been trying to solve poiseuille Stokes problem without succeeding. The code seems to be correct; I obtain the solution and the corresponding graph. However, when attempting to study convergence, I noticed that errors are increasing instead of decreasing as the mesh is refined.

Here are the L2 errors for velocity i got:

  • 8x8 : 6.7252891620447405e-16
  • 16x16 : 1.4505225737978914e-15
  • 32x32: 3.4758634833870467e-15
  • 64x64: 9.849036242099696e-15

My Stokes problem is:
-\Delta u+\nabla p=0 in \Omega=]0,2.5[ \times ]-0.5,0.5[
\nabla.u=0, in \Omega=]0,2.5[ \times ]-0.5,0.5[
u=(0.25-y^2,0) , on x=0 or x=2.5
u=(0,0), on y=-0.5 or y=0.5

My exact solutions are: u(x,y)=(0.25-y^2,0) and p(x,y)=-2x

My questions are:

  • am i taking the wrong exact solutions for velocity and pessure?
  • For the variational form, I included the pressure term with a negative sign. Is it better to put a positive sign and adjust the pressure sign after?
  • how can i solve the bad convergence of my problem?

I really need your help. Have a good day.

here is my code:

import dolfinx
from dolfinx import plot
from dolfinx import mesh, fem
import numpy as np
import ufl
from dolfinx.fem import Function, functionspace, dirichletbc, locate_dofs_topological, form
from dolfinx.cpp.mesh import CellType
from dolfinx.fem import locate_dofs_geometrical, locate_dofs_topological
from dolfinx.io import XDMFFile
from dolfinx.mesh import locate_entities_boundary,create_rectangle
from mpi4py import MPI
import dolfinx.fem.petsc
from petsc4py import PETSc
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")

# Create mesh
mesh = create_rectangle(MPI.COMM_WORLD,
                     [np.array([0, -0.5]), np.array([2.5, 0.5])],
                     [16, 16], CellType.triangle)


def wall_boundary(x):
    return np.logical_or(np.isclose(x[1], -0.5), np.isclose(x[1],0.5 ))

def in_out_boundary(x):
    return np.logical_or(np.isclose(x[0], 0.0), np.isclose(x[0], 2.5))

# Lid velocity
def lid_velocity_expression(x):
    return np.vstack((0.25-x[1]**2,np.zeros(x.shape[1])))

P2 = ufl.VectorElement("Lagrange", mesh.ufl_cell(), 2)
P1 = ufl.FiniteElement("Lagrange", mesh.ufl_cell(), 1)
V, Q = functionspace(mesh, P2), functionspace(mesh, P1)

# Create the function space

TH = P2 * P1
W = functionspace(mesh, TH)
W0,_ = W.sub(0).collapse()
W1, W1_to_W = W.sub(1).collapse()

# No slip boundary condition
noslip = Function(W0)
facets = locate_entities_boundary(mesh, 1, wall_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(mesh, 1, in_out_boundary)
dofs = locate_dofs_topological((W.sub(0), W0), 1, facets)
bc1 = 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 = Function(W1)
with zero.vector.localForm() as zero_local:
    zero_local.set(0.0)
dofs = locate_dofs_geometrical((W.sub(1), W1),
                               lambda x: np.isclose(x.T, [0, 0, 0]).all(axis=1))
bc2 = dirichletbc(zero, dofs, W.sub(1))

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

# 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(mesh.comm)
ksp.setOperators(A)
ksp.setType("preonly")
ksp.getPC().setType("lu")
ksp.getPC().setFactorSolverType("superlu_dist")

# Compute the solution
U = Function(W)
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()

# Define velocity exact solution
u_exact = Function(W0)
u_exact.interpolate(lid_velocity_expression)
# Define pressure exact solution
p_exact = Function(W1)
p_exact.interpolate(lambda x: -2*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(mesh.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(mesh.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(mesh.comm.allreduce(error_local_2, op=MPI.SUM))

pyvista.start_xvfb()
topology, cell_types, geometry = vtk_mesh(u.function_space)
values = np.zeros((geometry.shape[0], 3), dtype=np.float64)
values[:, :len(u)] = u.x.array.real.reshape((geometry.shape[0], len(u)))

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

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

The errors are of the magnitude of Machine precision, thus you cannot expect a decrease. You have errors of order 10e-30, as you have taken the square root of the error

3

Thank you for your answer,
I understand perfectly what you mean, what recommendations would you suggest to address this issue effectively?

4

Any reason why you couldn’t continue the discussion in Convergence problem for poiseuille flow of stokes equations ?
Purposely opening duplicate posts without an appropriate justification may lead to us closing down the second one, i.e. this one.
(And, no, the fact that the domain is not a square anymore but a rectangle is not an appropriate justification, in general).

what recommendations would you suggest to address this issue effectively

Pick a u and p which are not quadratic and linear, respectively, but still solution of the Stokes problem.

5

I apologize, i believed that due to encountering another issues, it would be better to create a separate post.
The main issue I’m facing is the lack of convergence, even when attempting to change the exact solutions.

I completely understand if this post needs to be deleted.

Have a good day;

6

No worries. We’ll keep this open since we’ve already started the discussion here, but please keep this in mind in future.

You haven’t changed them in the right way. One way could be as suggested above

Pick a u and p which are not quadratic and linear, respectively, but still solution of the Stokes problem.

If you don’t have any solution in mind that satisfies that, a possible choice could be the Taylor-Green vortex solution, i.e. something like u = (\sin(2 \pi x) \cos(2 \pi y), - \cos(2 \pi x) \sin(2 \pi y)) on the unit square. Please double check that it is a divergence free field, and determine the corresponding pressure expression from the momentum equation.

1 Like
7

Thank you very much for your help. I will verify the velocity solution you suggested and calculate the corresponding pressure. I will then update you with my findings.
Your help is greatly appreciated.

8

Excuse me, sir. It seems there is an issue with the example you provided. The velocity
u is not compatible with the Poiseuille flow I am looking for. That’s why i took u=(1-y^2,0) as exact solution for velocity.

9

But since you are using a finite element space that contains the solution, then you get machine precision errors, i.e. the solution is as close to the exact solution as it can be represented on your machine.

1 Like
10

as @dokken says, and other people have already said in previous posts, you must NOT be using a Poiseuille flow if your goal is to do a convergence study. The expression I sent you contains sin and cos on purpose, since they cannot be represented exactly by the FE space.

1 Like
11

Ah ok, i inderstand now. thank so much Mr @dokken and @francesco-ballarin for your answers.

12

Poisseuille flow in a pipe or annulus is typically good for a convergence check since the analytical solutions do not belong to a polynomial space.
However, this is typically solved in a cylindrical coordinate system. See, e.g.,

and the code developed by @kamilova-maths

1 Like