# AssertError in solving Stokes equation

Good evening.
I am trying to solve the Stokes equation on a 2D domain that I created through gmsh. In particular, I am referring to Stokes equations using Taylor-Hood elements — DOLFINx 0.7.3 documentation as my starting point, and I am trying to use this tutorial on my geometry. I already tried to solve the Poisson equation on my geometry and everything works well, for completeness I attach my .geo file through which i generated my .msh file:

``````lc=1;
//+
SetFactory("Built-in");
//+
Point(1) = {0, 0, 0, lc};
//+
Point(2) = {22.09, 0, 0, lc};
//+
Point(3) = {13.2684, 23.2987, 0, lc};
//+
Line(1) = {2, 1};
//+
Physical Curve("inlet", 1) = {1};
//+
Point(4) = {0.0810365,2.0938, 0.,lc};
//+
Point(5) = {0.323662, 4.17507 , 0.,lc};
//+
Point(6) = {0.726423, 6.23136 , 0., lc};
//+
Point(7) = {1.28691, 8.25038, 0., lc};
//+
Point(8) = {2.00178, 10.22,0.,lc};
//+
Point(9) = {2.86673,12.1285 , 0.,lc};
//+
Point(10) = {3.87662, 13.9645 ,  0., lc};
//+
Point(11) = {5.02538, 15.7169 ,  0., lc};
//+
Point(12) = {6.30615, 17.3753 ,  0., lc};
//+
Point(13) = {7.71127, 18.9297 ,  0., lc};
//+
Point(14) = { 9.20661, 20.3483 ,  0., lc};
//+
Point(15) = {10.5994, 21.4925,0., lc};
//+
Point(16) = {11.8011, 22.3633,  0., lc};
//+
Point(17) = {12.8204, 23.0271, lc};
//+
Spline(2) = {1,4,5,6,7,8,9,10,11,12,13,14,15,16,17,3};
//+
Point(18) = {7.09, 40.  ,0.,lc};
//+
Point(19) = {14.6913,40. ,0.,lc};
//+
Point(20) = {19.8334, 26.1,0.,lc};
//+
Line(3) = {3, 18};
//+
Line(4) = {18, 19};
//+
Line(5) = {19,20};
//+
Point(21) = {20.9561, 26.3864,  0.,lc};
//+
Point(22) = {22.09,26.6246,0.,lc};
//+
Spline(6) = {20,21,22};
//+
Point(23) = {22.09, 40.,  0.,lc};
//+
Point(24) = {27.09, 40.  ,  0.,lc};
//+
Point(25) = {27.09, 27.09, 0.,lc};
//+
Line(7) = {22,23};
//+
Line(8) = {23,24};
//+
Line(9) = {24,25};
//+
Point(26) = {27.5761, 27.0856,  0.,lc};
//+
Point(27) = {28.6936, 27.0425,  0.,lc};
//+
Point(28) = {30.0321, 26.9298,  0.,lc};
//+
Point(29) = {31.4608, 26.7351,  0.,lc};
//+
Point(30) = {32.6669, 26.5097,  0.,lc};
//+
Point(31) = {33.6603, 26.2812,  0.,lc};
//+
Spline(10) = {25,26,27,28,29,30,31};
//+
Point(32) = {37.09, 40.  ,  0.,lc};
//+
Point(33) = {42.09, 40.  ,  0.,lc};
//+
Point(34) = {38.7524, 24.4511,  0.,lc};
//+
Line(11) = {31,32};
//+
Line(12) = {32,33};
//+
Line(13) = {33,34};
//+
Point(35) = {22.1893 ,  2.02761,  0.,lc};
//+
Point(36) = {22.4915 ,  4.03495,  0.,lc};
//+
Point(37) = {23.0096 , 5.99743, 0.,lc};
//+
Point(38) = {23.7679 , 7.87976,  0.,lc};
//+
Point(39) = {24.8037 ,  9.62365,  0.,lc};
//+
Point(40) = {26.1672, 11.1223,  0.,lc};
//+
Point(41) = {27.8575, 12.1619,  0.,lc};
//+
Point(42) = {29.59, 12.5 ,  0. ,lc};
//+
Point(43) = {31.3208, 12.1626,  0.,lc};
//+
Point(44) = {33.0113, 11.1236,  0.,lc};
//+
Point(45) = {34.3753 ,  9.62498,  0.,lc};
//+
Point(46) = {35.4116 ,  7.88094,  0.,lc};
//+
Point(47) = {36.17   ,  5.99837,  0.,lc};
//+
Point(48) = {36.6884 ,  4.03559,  0. ,lc};
//+
Point(49) = {36.9906 ,  2.02794,  0.,lc};
//+
Point(50) = {37.09, 0.  ,  0.,lc};
//*
Point(51) = {36.669 ,  -4.12918,  0. ,lc};
//+
Point(52) = {36.1146 ,  -6.13024,  0.,lc};
//+
Point(53) = {35.2404 ,  -8.016389999999,  0.,lc};
//+
Point(54) = {34.4245 ,  -9.9285,  0. ,lc};
//+
Point(55) = {33.6611 ,  -11.8622,  0.,lc};
//+
Point(56) = {35.2404 ,  -8.016389999999,  0.,lc};
//+
Point(57) = {32.9455,  -13.8141 ,lc};
//+
Point(58) = {32.274 ,  -15.7815,  0.,lc};
//+
Point(59) = {31.644 ,  -17.7627,  0.,lc};
//+
Point(60) = {29.496 ,  -25.7946,  0.,lc};
//+
Point(61) = {28.0614, -32.9273, 0, lc};
//+
Point(62) = {27.2791, -38.0651, 0, lc};
//+
Point(63) = {26.7022, -43.7478, 0, lc};
//+
Recursive Delete {
Point{63};
}
//+
Point(63) = {26.7022, -43.2301, 0, lc};
//+
Point(64) = {26.3437, -48.4147, 0, lc};
//+
Point(65) = {26.2375, -54.6498, 0, lc};
//+
Point(66) = {26.4723, -59.841, 0, lc};
//+
Point(67) = {27.09, -65, 0, lc};
//+
Point(68) = {40.9924, -65, 0, lc};
//+
Spline(14) = {2, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50};
//+
Spline(15) = {50, 51, 52, 53, 54, 55, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67};
//+
Line(16) = {67, 68};
//+
Point(69) = {45.2112, 20.1368, 0, lc};
//+
Point(70) = {49.0032, 15.9273, 0, lc};
//+
Point(71) = {50.6766, 13.3244, 0, lc};
//+
Point(72) = {51.8364, 11.0219, 0, lc};
//+
Point(73) = {52.6035, 9.10643, 0, lc};
//+
Point(74) = {53.4755, 6.13776, 0, lc};
//+
Point(75) = {54.0033, 3.08903, -3.9, lc};
//+
Recursive Delete {
Point{75};
}
//+
Point(75) = {54.1014, 2.06184, 0., lc};
//+
Point(76) = {54.18, 0., 0., lc};
//+
Point(77) = {54.0559, -2.58973, 0., lc};
//+
Point(78) = {54.4679, -6.17023, 0., lc};
//+
Recursive Delete {
Point{78};
}
//+
Point(78) = {53.6848, -5.15575, 0., lc};
//+
Point(79) = {52.5867, -9.1535, 0., lc};
//+
Point(80) = {51.3692, -12.0161, 0., lc};
//+
Point(81) = {50.3786, -13.8387, 0., lc};
//+
Point(82) = {48.6386, -16.4173, 0., lc};
//+
Point(83) = {46.9089, -20.181, 0., lc};
//+
Point(84) = {45.1548, -25.0617, 0., lc};
//+
Point(85) = {44.256, -28.0414, 0., lc};
//+
Point(86) = {43.7134, -30.0441, 0., lc};
//+
Point(87) = {42.9787, -33.0684, 0., lc};
//+
Point(88) = {42.1397, -37.1324, 0., lc};
//+
Point(89) = {41.456, -41.2253, 0., lc};
//+
Point(90) = {40.9274, -45.3412, 0., lc};
//+
Point(91) = {40.4947, -50.51, 0., lc};
//+
Point(92) = {40.3483, -54.657, 0., lc};
//+
Point(93) = {40.4596, -59.8421, 0., lc};
//+
Spline(17) = {68, 93, 92, 91, 90, 89, 88, 87, 86, 85, 84, 83, 82, 81, 80, 79, 78, 77, 76, 75, 74, 73, 72, 71, 70, 69, 34};

//+
Curve Loop(1) = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, -17, -16, -15, -14};
//+
Plane Surface(1) = {1};
//+
Physical Curve("wall", 18) = {2, 3, 5, 6, 7, 9, 10, 11, 13, 17, 14, 15};
//+
Physical Curve("outlet", 19) = {16};
//+
Physical Curve("C1", 20) = {4};
//+
Physical Curve("C2", 21) = {8};
//+
Physical Curve("C3", 22) = {12};
//+
Physical Surface("surface", 23) = {1};

``````

Therefore, starting from this file I generated a .msh file through the graphic interface of gmsh, so opening the .geo file and generating a 2D mesh with Frontal-Delaunay algorithm, setting the order equal to 2 and imported it as:

``````msh, cell_tags, facet_tags = gmshio.read_from_msh(
"Mesh_prova.msh", MPI.COMM_WORLD, 0, gdim=2)
facet_tags.name = "Facet markers"
``````

Just to try if everything works good, I am implementing a no-slip conditions on all the Dirichlet boundaries, and a homogeneous Neumann condition for the outlet boundary, and setting the forcing term = 0, so that I should get a null vector for my velocity and p=0 for my pressure as my solution.
After importing the mesh, I set my function spaces:

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

and my boundary condition:

``````ft = facet_tags
fdim = msh.topology.dim - 1
#Walls
noslip = np.zeros(msh.geometry.dim, dtype=PETSc.ScalarType)
bcu_walls = dirichletbc(noslip, locate_dofs_topological(V, fdim, facet_tags.find(18)), V)
# Inlet
bcu_inflow = dirichletbc(noslip, locate_dofs_topological(V, fdim, ft.find(1)),V)
#Dirichlet outlets
bcu_C1 = dirichletbc(noslip, locate_dofs_topological(V, fdim, ft.find(20)), V)
bcu_C2 = dirichletbc(noslip, locate_dofs_topological(V, fdim, ft.find(21)), V)
bcu_C3 = dirichletbc(noslip, locate_dofs_topological(V, fdim, ft.find(22)), V)
bcs = [bcu_walls,bcu_inflow,bcu_C1, bcu_C2, bcu_C3]
``````

At this point, I define the variational problem as in the tutorial:

``````(u, p) = ufl.TrialFunction(V), ufl.TrialFunction(Q)
(v, q) = ufl.TestFunction(V), ufl.TestFunction(Q)
f = Constant(msh, (PETSc.ScalarType(0), PETSc.ScalarType(0)))  # type: ignore
a = form([[inner(grad(u), grad(v)) * dx, inner(p, div(v)) * dx],
[inner(div(u), q) * dx, None]])
L = form([inner(f, v) * dx, inner(Constant(msh, PETSc.ScalarType(0)), q) * dx])  # type: ignore
# -
a_p11 = form(inner(p, q) * dx)
a_p = [[a[0][0], None],
[None, a_p11]]
``````

Then I want to solve my problem through the Nested matrix solver defined as:

``````def nested_iterative_solver():
"""Solve the Stokes problem using nest matrices and an iterative solver."""

# Assemble nested matrix operators
A = fem.petsc.assemble_matrix_nest(a, bcs=bcu_walls)
A.assemble()

# Create a nested matrix P to use as the preconditioner. The
# top-left block of P is shared with the top-left block of A. The
# bottom-right diagonal entry is assembled from the form a_p11:
P11 = fem.petsc.assemble_matrix(a_p11, [])
P = PETSc.Mat().createNest([[A.getNestSubMatrix(0, 0), None], [None, P11]])
P.assemble()
A00 = A.getNestSubMatrix(0, 0)
A00.setOption(PETSc.Mat.Option.SPD, True)
P00, P11 = P.getNestSubMatrix(0, 0), P.getNestSubMatrix(1, 1)
P00.setOption(PETSc.Mat.Option.SPD, True)
P11.setOption(PETSc.Mat.Option.SPD, True)

# Assemble right-hand side vector
b = fem.petsc.assemble_vector_nest(L)

# Modify ('lift') the RHS for Dirichlet boundary conditions
fem.petsc.apply_lifting_nest(b, a, bcs=bcs)

# Sum contributions for vector entries that are share across
# parallel processes
for b_sub in b.getNestSubVecs():
b_sub.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)

# Set Dirichlet boundary condition values in the RHS vector
bcs0 = fem.bcs_by_block(extract_function_spaces(L), bcs)
fem.petsc.set_bc_nest(b, bcs0)

# The pressure field is determined only up to a constant. We supply
# a vector that spans the nullspace to the solver, and any component
# of the solution in this direction will be eliminated during the
# solution process.
null_vec = fem.petsc.create_vector_nest(L)

# Set velocity part to zero and the pressure part to a non-zero
# constant
null_vecs = null_vec.getNestSubVecs()
null_vecs[0].set(0.0), null_vecs[1].set(1.0)

# Normalize the vector that spans the nullspace, create a nullspace
# object, and attach it to the matrix
null_vec.normalize()
nsp = PETSc.NullSpace().create(vectors=[null_vec])
assert nsp.test(A)
A.setNullSpace(nsp)

# Create a MINRES Krylov solver and a block-diagonal preconditioner
# using PETSc's additive fieldsplit preconditioner
ksp = PETSc.KSP().create(msh.comm)
ksp.setOperators(A, P)
ksp.setType("minres")
ksp.setTolerances(rtol=1e-9)
ksp.getPC().setType("fieldsplit")
ksp.getPC().setFieldSplitType(PETSc.PC.CompositeType.ADDITIVE)

# Define the matrix blocks in the preconditioner with the velocity
# and pressure matrix index sets
nested_IS = P.getNestISs()
ksp.getPC().setFieldSplitIS(("u", nested_IS[0][0]), ("p", nested_IS[0][1]))

# Set the preconditioners for each block. For the top-left
# Laplace-type operator we use algebraic multigrid. For the
# lower-right block we use a Jacobi preconditioner. By default, GAMG
# will infer the correct near-nullspace from the matrix block size.
ksp_u, ksp_p = ksp.getPC().getFieldSplitSubKSP()
ksp_u.setType("preonly")
ksp_u.getPC().setType("gamg")
ksp_p.setType("preonly")
ksp_p.getPC().setType("jacobi")

# Create finite element {py:class}`Function <dolfinx.fem.Function>`s
# for the velocity (on the space `V`) and for the pressure (on the
# space `Q`). The vectors for `u` and `p` are combined to form a
# nested vector and the system is solved.
u, p = Function(V), Function(Q)
x = PETSc.Vec().createNest([la.create_petsc_vector_wrap(u.x),
la.create_petsc_vector_wrap(p.x)])
ksp.solve(b, x)

# Save solution to file in XDMF format for visualization, e.g. with
# ParaView. Before writing to file, ghost values are updated using
# `scatter_forward`.
with XDMFFile(MPI.COMM_WORLD, "out_stokes/velocity.xdmf", "w") as ufile_xdmf:
u.x.scatter_forward()
P1 = element("Lagrange", msh.basix_cell(), 1, shape=(msh.geometry.dim,))
u1 = Function(functionspace(msh, P1))
u1.interpolate(u)
ufile_xdmf.write_mesh(msh)
ufile_xdmf.write_function(u1)
with XDMFFile(MPI.COMM_WORLD, "out_stokes/pressure.xdmf", "w") as pfile_xdmf:
p.x.scatter_forward()
pfile_xdmf.write_mesh(msh)
pfile_xdmf.write_function(p)

# Compute norms of the solution vectors
norm_u = u.x.norm()
norm_p = p.x.norm()
if MPI.COMM_WORLD.rank == 0:
print(f"(A) Norm of velocity coefficient vector (nested, iterative): {norm_u}")
print(f"(A) Norm of pressure coefficient vector (nested, iterative): {norm_p}")
return norm_u, norm_p
``````

The problem arises now, when I call the solver through:

``````norm_u_0, norm_p_0 = nested_iterative_solver()
``````

I get the following error:

``````AssertionError                            Traceback (most recent call last)
Cell In[21], line 1
----> 1 norm_u_0, norm_p_0 = nested_iterative_solver()

Cell In[20], line 50, in nested_iterative_solver()
48 null_vec.normalize()
49 nsp = PETSc.NullSpace().create(vectors=[null_vec])
---> 50 assert nsp.test(A)
51 A.setNullSpace(nsp)
53 # Create a MINRES Krylov solver and a block-diagonal preconditioner
54 # using PETSc's additive fieldsplit preconditioner

AssertionError:
``````

I tried to change boundary conditions, I also tried different solvers but I always get this error related to this null space. I tried to look on the documentation of PETSc.NullSpace but I cannot fully understand why I get this error.
Sorry for the long and maybe silly question but being a beginner with FEniCSx it is not easy for me to get to the bottom of this error.
Thank you very much for your help.

I have at tutorial_create_nullspace a tutorial about singular problems, it’s not about Stokes equations but maybe it can help you understand what is going on in the part of the code where you set up the null space.

About your case: I don’t think you should impose the null space, since you don’t have a singular problem. The problem would have been singular if Dirichlet boundary conditions were imposed on the entire boundary, but you say you have a Neumann condition at the outflow.

Finally, make sure to read https://fenicsproject.discourse.group/t/read-before-posting-how-do-i-get-my-question-answered/21: the code you post is not fully reproducible, in the sense that people cannot just copy and paste to run it.

Good evening,
thank you very much for your answer, I will look at the tutorial and see if I can find any useful information.
I understand your observation about the reproducibility of the code, but in order to pass from the .geo file to the .msh file, I used the graphic interface of gmsh. I updated my post explaining the settings I used on gmsh, but being a complex geometry the file .msh is very long to be copied and paste here, but if you have any suggestion on how to make this code reproducible, I am very happy to hear that.
I see your point about the null space, actually at first I was using Dirichlet boundary conditions everywhere and this is why I kept that part of the code, but with the actual version I do not need the null space. Indeed, by eliminating the part of the solver with the null space and by interpolating the functions correctly in order to be saved through XDMFFile, I managed to get reasonable solutions also with different boundary values.
Thank you for your precious help.