Stokes problem with nullspace to ensure unicity

dolfinx
1

I’m trying to solve a Stokes problem using the axisymmetrical hypothesis.
Following the solution provided in Stokes problem with only Dirichlet boundary condition - #2 by dokken I tried to apply the same nullspace to ensure unicity of the solution (since I impose only Dirichlet BC), however I have an error that the nullspace is not correctly imposed.

Here the example code with a testcase already in place:

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


def error_L2(uh, u_ex, degree_raise=3):
    # Create higher order function space
    ufl_el = uh.function_space.ufl_element()
    new_el = ufl_el.reconstruct(degree=degree_raise+ufl_el.degree())
    mesh = uh.function_space.mesh
    W = fem.FunctionSpace(mesh, new_el)
    # Interpolate approximate solution
    u_W = fem.Function(W)
    u_W.interpolate(uh)

    # Interpolate exact solution, special handling if exact solution
    # is a ufl expression or a python lambda function
    u_ex_W = fem.Function(W)
    if isinstance(u_ex, ufl.core.expr.Expr):
        u_expr = fem.Expression(u_ex, W.element.interpolation_points())
        u_ex_W.interpolate(u_expr)
    else:
        u_ex_W.interpolate(u_ex)

    # Compute the error in the higher order function space
    e_W = fem.Function(W)
    e_W.x.array[:] = u_W.x.array - u_ex_W.x.array

    # Integrate the error
    error = fem.form(ufl.inner(e_W, e_W) * ufl.dx)
    error_local = fem.assemble_scalar(error)
    error_global = mesh.comm.allreduce(error_local, op=MPI.SUM)
    return np.sqrt(error_global)


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)

# Cylindrical coordinates: x[0] --> r, x[1] --> z

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

# Functions to geometrically locate subsets
def wall1_boundary(x):
    return np.isclose(x[0], 0.5)

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

def inlet_boundary(x):
    return np.isclose(x[1], 0.0)

def outlet_boundary(x):
    return np.isclose(x[1], 2.0)

# Inflow velocity
def in_velocity_expression(x):
    return np.vstack((np.zeros(x.shape[1]), -2*(np.ones(x.shape[1])-x[0]**2)))


# Walls velocity
def wall1_velocity_expression(x):
    return np.vstack((np.zeros(x.shape[1]), -(3/2)*np.ones(x[1].shape[0])))


def wall2_velocity_expression(x):
    return np.vstack((np.zeros(x.shape[1]),-2*np.ones(x[1].shape[0])))


# Pressure out
def pressure_out_expression(x):
    return 14*np.ones(x[1].shape[0])



# Boundary conditions
# Wall1 velocity 
wall1_velocity = Function(V)
wall1_velocity.interpolate(wall1_velocity_expression)
facets = locate_entities_boundary(msh, 1, wall1_boundary)
dofs = locate_dofs_topological((W.sub(0), V), 1, facets)
bc0 = dirichletbc(wall1_velocity, dofs, W.sub(0))

# Out pressure 
out_pressure = Function(Q)
out_pressure.interpolate(pressure_out_expression)
facets = locate_entities_boundary(msh, 1, outlet_boundary)
dofs = locate_dofs_topological((W.sub(1), Q), 1, facets)
bc1 = dirichletbc(out_pressure, dofs, W.sub(1))

# Inflow velocity
in_velocity = Function(V)
in_velocity.interpolate(in_velocity_expression)
facets = locate_entities_boundary(msh, 1, inlet_boundary)
dofs = locate_dofs_topological((W.sub(0), V), 1, facets)
bc2 = dirichletbc(in_velocity, dofs, W.sub(0))

# Wall2 velocity
wall2_velocity = Function(V)
wall2_velocity.interpolate(wall2_velocity_expression)
facets = locate_entities_boundary(msh, 1, wall2_boundary)
dofs = locate_dofs_topological((W.sub(0), V), 1, facets)
bc3 = dirichletbc(wall2_velocity, dofs, W.sub(0))

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

# Cylindrical coordinates
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)

f = Function(V)
f = ufl.as_vector((0, -1))

n  = ufl.FacetNormal(msh)
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)
L = form(inner(f, v) * x[0] * 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)

# 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)
A.setNearNullSpace(nullspace)

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

# Compute the solution
U = Function(W)
ksp.solve(b, U.vector)

# Split the mixed solution and collapse
uh = U.sub(0).collapse()
ph = U.sub(1).collapse()


# Compute error L2
# Velocity
def u_e(x):
    return np.vstack((np.zeros(x.shape[1]), -2*(np.ones(x[1].shape[0])-x[0]**2)))

print('error_L2_velocity  =', error_L2(uh, u_e))

# Pressure
def p_e(x):
    return  7*x[1]

print('error_L2_pressure  =', error_L2(ph, p_e))
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?

1 Like