Good evening,
in the Stokes problem using P1-P1 finite elements with a stabilization term (I implemented the easiest, i.e. the Brezzi-Pitkaranta stabilization), one should expect a convergence like O(h) for the H1 error in the velocity and the pressure. In my case this does not happens and following previous advice I tried first to solve this problem using a direct solver, but still the results are not right.
The minimal code is the following:
import numpy as np
import ufl
from dolfinx import cpp as _cpp
from dolfinx import fem
from dolfinx.fem import (Constant, Function, FunctionSpace, dirichletbc,
extract_function_spaces, form,
locate_dofs_geometrical, locate_dofs_topological)
from dolfinx.io import XDMFFile
from dolfinx.mesh import (CellType, GhostMode, create_rectangle, create_unit_cube,
locate_entities_boundary, exterior_facet_indices)
from ufl import div, dx, grad, inner
from mpi4py import MPI
from petsc4py import PETSc
# Create mesh
N = 32
msh = create_unit_cube(MPI.COMM_WORLD,
N, N, N,
CellType.hexahedron, GhostMode.none)
'''# 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])))'''
def u_ex(x):
return np.vstack((4*(x[1]-1)*x[1]*(x[0]-1)*(x[0]-1)*x[0]*x[0]*(x[1]-0.5), -4*(x[0]-0.5)*x[1]*x[1]*x[0]*(x[1]-1)*(x[1]-1)*(x[0]-1), np.zeros_like(x[0])))
def p_ex(x):
return x[0]+x[1]+x[2]-1.5
def f_ex(x):
return np.vstack((1-24*(-1 + x[0])*(x[0]-1)*x[0]*x[0]*(-0.5 + x[1])- 48*(0.166667 - x[0] + x[0]*x[0])*(-1 + x[1])*(-0.5 + x[1])*x[1], 1+24*(-0.5 + x[0])*(-1 + x[1])*(x[1]-1)*x[1]*x[1] + 8*(-1 + x[0])*(-0.5 + x[0])*x[0]*(1 - 6*x[1] + 6*x[1]*x[1]), np.ones_like(x[0])))
P2 = ufl.VectorElement("Lagrange", msh.ufl_cell(), 1)
P1 = ufl.FiniteElement("Lagrange", msh.ufl_cell(), 1)
# Create the function space
TH = P2 * P1
W = FunctionSpace(msh, TH)
W0, _ = W.sub(0).collapse()
W1, _ = W.sub(1).collapse()
'''# No slip boundary condition
noslip = Function(V)
facets = locate_entities_boundary(msh, 1, noslip_boundary)
dofs = locate_dofs_topological((W.sub(0), V), 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), V), 1, facets)
bc1 = dirichletbc(lid_velocity, dofs, W.sub(0))'''
uD = Function(W0)
uD.interpolate(u_ex)
uD.vector.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
# Create facet to cell connectivity required to determine boundary facets
tdim = msh.topology.dim
fdim = tdim - 1
msh.topology.create_connectivity(fdim, tdim)
boundary_facets = exterior_facet_indices(msh.topology)
boundary_dofs = locate_dofs_topological((W.sub(0),W0), fdim, boundary_facets)
bc1 = dirichletbc(uD, boundary_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)
zero.x.set(-1.5)
zero.vector.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
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 = [bc1, bc2]
# Define variational problem
(u, p) = ufl.TrialFunctions(W)
(v, q) = ufl.TestFunctions(W)
f = Function(W0)
f.interpolate(f_ex)
h = PETSc.ScalarType(ufl.sqrt(2)*1./N)
beta = PETSc.ScalarType(0.1)
a = form((inner(grad(u), grad(v)) + inner(p, div(v)) + inner(div(u), q)) * dx - h*h*beta*ufl.inner(ufl.nabla_grad(p),ufl.nabla_grad(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")
ksp.getPC().setType("lu")
ksp.getPC().setFactorSolverType("mumps")
opts = PETSc.Options()
opts["mat_mumps_icntl_14"] = 80 # Increase MUMPS working memory
opts["mat_mumps_icntl_24"] = 0 # Option to support solving a singular matrix (pressure nullspace)
opts["mat_mumps_icntl_25"] = 0 # Option to support solving a singular matrix (pressure nullspace)
opts["ksp_error_if_not_converged"] = 1
opts["ksp_monitor"] = None
ksp.setFromOptions()
# Compute the solution
U = Function(W)
ksp.solve(b, U.vector)
# Split the mixed solution and collapse
u = U.sub(0).collapse()
p = U.sub(1).collapse()
L2_error = fem.form(ufl.inner(u - uD, u - uD) * ufl.dx)
error_local = fem.assemble_scalar(L2_error)
error_L2 = np.sqrt(msh.comm.allreduce(error_local, op=MPI.SUM))
if msh.comm.rank == 0:
print(f"Error_L2 : {error_L2:.2e}")
eh = u - uD
error_H10 = fem.form(ufl.inner(ufl.grad(eh), ufl.grad(eh)) * ufl.dx)
E_H10 = np.sqrt(msh.comm.allreduce(fem.assemble_scalar(error_H10), op=MPI.SUM))
if msh.comm.rank == 0:
print(f"H01-error: {E_H10:.2e}")
pD = Function(W1)
pD.interpolate(p_ex)
L2_errorp = fem.form(ufl.inner(p - pD, p - pD) * ufl.dx)
error_localp = fem.assemble_scalar(L2_errorp)
error_L2p = np.sqrt(msh.comm.allreduce(error_localp, op=MPI.SUM))
if msh.comm.rank == 0:
print(f"Error_L2p : {error_L2p:.2e}")
# Write the solution to file
with XDMFFile(MPI.COMM_WORLD, "out_stokes/u.xdmf", "w") as ufile_xdmf:
ufile_xdmf.write_mesh(msh)
ufile_xdmf.write_function(u)
with XDMFFile(MPI.COMM_WORLD, "out_stokes/p.xdmf", "w") as pfile_xdmf:
pfile_xdmf.write_mesh(msh)
pfile_xdmf.write_function(p)
with XDMFFile(MPI.COMM_WORLD, "out_stokes/uex.xdmf", "w") as pfile_xdmf:
pfile_xdmf.write_mesh(msh)
pfile_xdmf.write_function(uD)
The error for a mesh with respectively 8, 16, 32 elements are:
Error_L2 : 5.08e-03
H01-error: 3.39e-02
Error_L2p : 2.67e+00
Error_L2 : 1.41e-03
H01-error: 1.14e-02
Error_L2p : 2.90e+00
Error_L2 : 3.68e-04
H01-error: 3.71e-03
Error_L2p : 3.03e+00
Which are note corrects (in particular the one of the pressure??). Can someone give me an hand on this problem that is driving me crazy?