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.