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