I am trying to compute different norm types (L1, L2, Linf) and have been using the stokes flow demo as practice (link). I am using the bottom method of solving the problem (Non-Blocked Direct Solver) and when I go to compute the infinite norm of the velocity and pressure, I get this error
Invoked with types: dolfinx.cpp.fem.Function_float64, str
To compute the infinite norm I am using this line of code
norm_inf_u = la.norm(u, 'linf')
When looking at the DOLPHIN documentation they give an example with this same line (Dolphin norm documentation)
I have also found the DOLFINx documentation (DOLPHINx Norm documentation) and when I tried the line of code,
norm_inf_u = la.norm(u.x, type=2)
I get the following error
Invoked with types: dolfinx.cpp.la.Vector_float64, int
How do I compute different norms of the velocity and pressure field?
I am using FEniCSx version 0.8, installed on linux using anaconda.
Here is my code, I have tried to remove as much as possible from it
from mpi4py import MPI
from petsc4py import PETSc
import numpy as np
import ufl
from basix.ufl import element, mixed_element
from dolfinx import fem, la
from dolfinx.fem import (
Function,
dirichletbc,
form,
functionspace,
locate_dofs_topological,
)
from dolfinx.fem.petsc import assemble_matrix_block, assemble_vector_block
from dolfinx.mesh import CellType, create_rectangle, locate_entities_boundary
from ufl import div, dx, grad, inner
# Create mesh
msh = create_rectangle(
MPI.COMM_WORLD, [np.array([0, 0]), np.array([1, 1])], [128, 128], CellType.triangle
)
# Function to mark x = 0, x = 1 and y = 0
def noslip_boundary(x):
return 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])))
P2 = element("Lagrange", msh.basix_cell(), 2, shape=(msh.geometry.dim,))
P1 = element("Lagrange", msh.basix_cell(), 1)
V, Q = functionspace(msh, P2), functionspace(msh, P1)
def mixed_direct():
# Create the Taylot-Hood function space
TH = mixed_element([P2, P1])
W = functionspace(msh, TH)
# No slip boundary condition
W0 = W.sub(0)
Q, _ = W0.collapse()
noslip = Function(Q)
facets = locate_entities_boundary(msh, 1, noslip_boundary)
dofs = locate_dofs_topological((W0, Q), 1, facets)
bc0 = dirichletbc(noslip, dofs, W0)
# Driving velocity condition u = (1, 0) on top boundary (y = 1)
lid_velocity = Function(Q)
lid_velocity.interpolate(lid_velocity_expression)
facets = locate_entities_boundary(msh, 1, lid)
dofs = locate_dofs_topological((W0, Q), 1, facets)
bc1 = dirichletbc(lid_velocity, dofs, W0)
# Collect Dirichlet boundary conditions
bcs = [bc0, bc1]
# Define variational problem
(u, p) = ufl.TrialFunctions(W)
(v, q) = ufl.TestFunctions(W)
f = Function(Q)
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(msh.comm)
ksp.setOperators(A)
ksp.setType("preonly")
# Configure MUMPS to handle pressure nullspace
pc = ksp.getPC()
pc.setType("lu")
# Compute the solution
U = Function(W)
try:
ksp.solve(b, U.x.petsc_vec)
except PETSc.Error as e:
if e.ierr == 92:
print("The required PETSc solver/preconditioner is not available. Exiting.")
print(e)
exit(0)
else:
raise e
# Split the mixed solution and collapse
u, p = U.sub(0).collapse(), U.sub(1).collapse()
# Compute norms
norm_u, norm_p = la.norm(u.x), la.norm(p.x)
norm_inf_u = la.norm(u.x, type=2)
if MPI.COMM_WORLD.rank == 0:
print(f"(D) Norm of velocity coefficient vector (monolithic, direct): {norm_u}")
print(f"(D) Norm of pressure coefficient vector (monolithic, direct): {norm_p}")
return norm_u, norm_u, u, p, U
# Solve using a non-blocked matrix and an LU solver
norm_u_3, norm_p_3, u, p, U = mixed_direct()