Stokes problem with nullspace to ensure unicity

dolfinx
2

Are you sure you have applied boundary conditions on all facets?
I will illustrate this on a massively simplified problem:


import numpy as np
from matplotlib import pyplot as plt
import ufl
import dolfinx.fem.petsc

from dolfinx import cpp as _cpp
from dolfinx import fem, plot
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,
                          locate_entities_boundary)
from ufl import as_tensor, div, dx, grad, inner, ds, sym, dot

from mpi4py import MPI
from petsc4py import PETSc


mu = 1         # dynamic viscosity
nx = 16
ny = 4*nx
# Create mesh
msh = create_rectangle(MPI.COMM_WORLD,
                       [np.array([0, 0]), np.array([0.5, 2])],
                       [nx, ny],
                       CellType.triangle)

x = ufl.SpatialCoordinate(msh)


# Function spaces
P2 = ufl.VectorElement("Lagrange", msh.ufl_cell(), 2)
P1 = ufl.FiniteElement("Lagrange", msh.ufl_cell(), 1)
TH = P2 * P1
W = FunctionSpace(msh, TH)
V, V_to_W = W.sub(0).collapse()
Q, Q_to_W = W.sub(1).collapse()


def div_cyl(v):
    return (1/x[0])*(x[0]*v[0]).dx(0) + v[1].dx(1)


def epsilon(v):
    return sym(grad(v))


# Define variational problem
(u, p) = ufl.TrialFunctions(W)
(v, q) = ufl.TestFunctions(W)

uD = dolfinx.fem.Function(V)
tdim = msh.topology.dim
fdim = tdim - 1
msh.topology.create_connectivity(fdim, tdim)
boundary_facets = dolfinx.mesh.exterior_facet_indices(msh.topology)
boundary_dofs = dolfinx.fem.locate_dofs_topological(
    (W.sub(0), V), fdim, boundary_facets)


bcs = [dolfinx.fem.dirichletbc(uD, boundary_dofs, W.sub(0))]  # works
bcs = []  # Nullspace test fails

# Normal stokes             
a = fem.form(ufl.inner(ufl.grad(u), ufl.grad(v)) * dx +
             ufl.div(v)*p*dx + q*ufl.div(u)*dx)
# Cylindrical stokes
a = form(mu * inner(epsilon(u), grad(v)) *
         x[0] * dx - inner(div_cyl(u), q) * x[0] * dx - inner(p, div_cyl(v)) * x[0] * dx)


# Assemble LHS matrix and RHS vector
A = fem.petsc.assemble_matrix(a, bcs=bcs)
A.assemble()

# assemble A and apply BCs
ns_vec = dolfinx.fem.Function(W)
ns_vec.x.array[Q_to_W] = 1
dolfinx.la.orthonormalize([ns_vec.vector])
nullspace = PETSc.NullSpace().create(vectors=[ns_vec.vector])
assert nullspace.test(A)

with empty bcs either of the nullspace tests, fails.
However, if you add in bcs, either formulation has the given nullspace.

For further posts, please think carefully about how you can remove code from your problem to illustrate the problem.

I guess your issue is that you actually set values on the pressure, so you do not have a null space issue?

2 Likes