Hello, I want to use (P_1+RT0, P0) to solve the Stokes problem, which is
And what I’m wondering is how do I assemble
P_1+RT0.I’m using basix.ufl.enriched_element did not succeed.
Traceback (most recent call last):
File "/mnt/d/femcode/fem/Stokes_code/stokeseich.py", line 217, in <module>
enriched_element = basix.ufl.enriched_element(
[
...<2 lines>...
]
)
File "/home/zqaixj1234/miniconda3/envs/fenicsx-env/lib/python3.13/site-packages/basix/ufl.py", line 2091, in enriched_element
raise ValueError("Enriched elements on different map types not supported.")
ValueError: Enriched elements on different map types not supported.
ERROR conda.cli.main_run:execute(124): `conda run python /mnt/d/femcode/fem/Stokes_code/stokeseich.py` failed. (See above for error)
See the chart below for details
from mpi4py import MPI
from petsc4py import PETSc
import matplotlib.pylab as plt
import numpy as np
import basix.ufl
import dolfinx.fem.petsc
from dolfinx import fem, mesh
from ufl import (
SpatialCoordinate,
TestFunction,
TrialFunction,
as_vector,
cos,
div,
dx,
grad,
inner,
pi,
sin,
)
def u_ex(x):
sinx = sin(pi * x[0])
siny = sin(pi * x[1])
cosx = cos(pi * x[0])
cosy = cos(pi * x[1])
c_factor = 2 * pi * sinx * siny
return c_factor * as_vector((cosy * sinx, -cosx * siny))
def p_ex(x):
return sin(2 * pi * x[0]) * sin(2 * pi * x[1])
def source(x):
u, p = u_ex(x), p_ex(x)
return -div(grad(u)) + grad(p)
def create_bilinear_form(V, Q):
u, p = TrialFunction(V), TrialFunction(Q)
v, q = TestFunction(V), TestFunction(Q)
a_uu = inner(grad(u), grad(v)) * dx
a_up = inner(p, div(v)) * dx
a_pu = inner(div(u), q) * dx
return fem.form([[a_uu, a_up], [a_pu, None]])
def create_linear_form(V, Q):
v, q = TestFunction(V), TestFunction(Q)
domain = V.mesh
x = SpatialCoordinate(domain)
f = source(x)
return fem.form([inner(f, v) * dx, inner(fem.Constant(domain, 0.0), q) * dx])
def create_velocity_bc(V):
domain = V.mesh
g = fem.Constant(domain, [0.0, 0.0])
tdim = domain.topology.dim
domain.topology.create_connectivity(tdim - 1, tdim)
bdry_facets = mesh.exterior_facet_indices(domain.topology)
dofs = fem.locate_dofs_topological(V, tdim - 1, bdry_facets)
return [fem.dirichletbc(g, dofs, V)]
def create_nullspace(rhs_form):
null_vec = fem.petsc.create_vector_nest(rhs_form)
null_vecs = null_vec.getNestSubVecs()
null_vecs[0].set(0.0)
null_vecs[1].set(1.0)
null_vec.normalize()
nsp = PETSc.NullSpace().create(vectors=[null_vec])
return nsp
def create_preconditioner(Q, a, bcs):
p, q = TrialFunction(Q), TestFunction(Q)
a_p11 = fem.form(inner(p, q) * dx)
a_p = fem.form([[a[0][0], None], [None, a_p11]])
P = dolfinx.fem.petsc.assemble_matrix_nest(a_p, bcs)
P.assemble()
return P
def assemble_system(lhs_form, rhs_form, bcs):
A = fem.petsc.assemble_matrix_nest(lhs_form, bcs=bcs)
A.assemble()
b = dolfinx.fem.petsc.assemble_vector_nest(rhs_form)
dolfinx.fem.petsc.apply_lifting_nest(b, lhs_form, bcs=bcs)
for b_sub in b.getNestSubVecs():
b_sub.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
spaces = fem.extract_function_spaces(rhs_form)
bcs0 = fem.bcs_by_block(spaces, bcs)
dolfinx.fem.petsc.set_bc_nest(b, bcs0)
return A, b
def create_block_solver(A, b, P, comm):
ksp = PETSc.KSP().create(comm)
ksp.setOperators(A, P)
ksp.setType("minres")
ksp.setTolerances(rtol=1e-9)
ksp.getPC().setType("fieldsplit")
ksp.getPC().setFieldSplitType(PETSc.PC.CompositeType.ADDITIVE)
nested_IS = P.getNestISs()
ksp.getPC().setFieldSplitIS(("u", nested_IS[0][0]), ("p", nested_IS[0][1]))
# Set the preconditioners for each block
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")
# Monitor the convergence of the KSP
ksp.setFromOptions()
return ksp
def assemble_scalar(J, comm: MPI.Comm):
scalar_form = fem.form(J)
local_J = fem.assemble_scalar(scalar_form)
return comm.allreduce(local_J, op=MPI.SUM)
def compute_errors(u, p):
domain = u.function_space.mesh
x = SpatialCoordinate(domain)
error_u = u - u_ex(x)
H1_u = inner(error_u, error_u) * dx
H1_u += inner(grad(error_u), grad(error_u)) * dx
velocity_error = np.sqrt(assemble_scalar(H1_u, domain.comm))
error_p = -p - p_ex(x)
L2_p = fem.form(error_p * error_p * dx)
pressure_error = np.sqrt(assemble_scalar(L2_p, domain.comm))
return velocity_error, pressure_error
def solve_stokes(u_element, p_element, domain):
V = fem.functionspace(domain, u_element)
Q = fem.functionspace(domain, p_element)
lhs_form = create_bilinear_form(V, Q)
rhs_form = create_linear_form(V, Q)
bcs = create_velocity_bc(V)
nsp = create_nullspace(rhs_form)
A, b = assemble_system(lhs_form, rhs_form, bcs)
assert nsp.test(A)
A.setNullSpace(nsp)
P = create_preconditioner(Q, lhs_form, bcs)
ksp = create_block_solver(A, b, P, domain.comm)
u, p = fem.Function(V), fem.Function(Q)
w = PETSc.Vec().createNest([u.x.petsc_vec, p.x.petsc_vec])
ksp.solve(b, w)
assert ksp.getConvergedReason() > 0
u.x.scatter_forward()
p.x.scatter_forward()
return compute_errors(u, p)
def error_plot(element_u, element_p, convergence_u=None, convergence_p=None, refinements=5, N0=7):
hs = np.zeros(refinements)
u_errors = np.zeros(refinements)
p_errors = np.zeros(refinements)
comm = MPI.COMM_WORLD
for i in range(refinements):
N = N0 * 2**i
domain = mesh.create_unit_square(comm, N, N, mesh.CellType.triangle)
u_errors[i], p_errors[i] = solve_stokes(element_u, element_p, domain)
hs[i] = 1.0 / N
print(f"h: {hs[i]:.4e} Error: {u_errors[i]:.4e}")
legend = []
if convergence_u is not None:
y_value = u_errors[-1] * 1.4
plt.plot(
[hs[0], hs[-1]],
[y_value * (hs[0] / hs[-1]) ** convergence_u, y_value],
"k--",
)
legend.append(f"order {convergence_u}")
if convergence_p is not None:
y_value = p_errors[-1] * 1.4
plt.plot(
[hs[0], hs[-1]],
[y_value * (hs[0] / hs[-1]) ** convergence_p, y_value],
"k--",
)
legend.append(f"order {convergence_p}")
plt.plot(hs, u_errors, "bo-")
plt.plot(hs, p_errors, "ro-")
legend += [r"$H^1(\mathbf{u_h}-\mathbf{u}_{ex})$", r"$L^2(p_h-p_ex)$"]
plt.legend(legend)
plt.xscale("log")
plt.yscale("log")
plt.axis("equal")
plt.ylabel("Error in energy norm")
plt.xlabel("$h$")
plt.xlim(plt.xlim()[::-1])
plt.grid(True)
#分段恒压空间
# element_u = basix.ufl.element("Lagrange", "triangle", 1, shape=(2,))
# element_p = basix.ufl.element("DG", "triangle", 0)
# error_plot(element_u, element_p)
#P2-DG0
# element_u = basix.ufl.element("Lagrange", "triangle", 2, shape=(2,))
# element_p = basix.ufl.element("DG", "triangle", 0)
# error_plot(element_u, element_p, convergence_p=1)
#Crouzeix-Raviart
# element_u = basix.ufl.element("CR", "triangle", 1, shape=(2,))
# element_p = basix.ufl.element("DG", "triangle", 0)
# error_plot(element_u, element_p, 1)
#P1+RT0
enriched_element = basix.ufl.enriched_element(
[
basix.ufl.element("Lagrange", "triangle", 1),
basix.ufl.element("Raviart-Thomas", "triangle", 1),
]
)
element_u = basix.ufl.blocked_element(enriched_element, shape=(2,))
element_p = basix.ufl.element("DG", "triangle", 1)
error_plot(element_u, element_p, 2)
