Hi, everyone.
I would like to solve a mixed poisson problem with homogeneous flux boundary conditions：
\sigma =-\grad u,
\div \sigma = -f. in \Omaga
The boundary condition is：
\sigma \cdot n=0 on \partial \Omega.
And the exact solution of \sigma that I set already satisfies \sigma\cdot n=0 on the boundary.
This is my code：
from mpi4py import MPI # type: ignore
from petsc4py import PETSc
import numpy as np
import ufl
from dolfinx.fem.petsc import assemble_matrix, assemble_vector, apply_lifting, create_vector, set_bc
from basix.ufl import element, mixed_element
from dolfinx import default_real_type, fem, io, mesh
from ufl import TestFunctions, TrialFunctions, div, inner,dx
def Exactu(x):
return (1/3)*x[0]**3 - (1/2)*x[0]**2 + (1/3)*x[1]**3 - (1/2)*x[1]**2
def Exact_sigma(x):
temp1= x[0]**2 - x[0]
temp2= x[1]**2 - x[1]
return [temp1,temp2]
def source_term(x):
return -2*x[0] - 2*x[1] + 2
def clamped_boundary(x):
return np.isclose(x[0], 0.0) | np.isclose(x[0], 1.0) | np.isclose(x[1], 0.0) | np.isclose(x[1], 1.0)
hs = []
errors_sigma_L2 = []
errors_u_L2= []
for i in range(4):
n = 2 ** (i + 3)
h = 1.0 / n
msh = mesh.create_unit_square(MPI.COMM_WORLD, n, n, mesh.CellType.triangle)
k = 1
#Q_el = element("BDMCF", msh.basix_cell(), k, dtype=default_real_type) mesh.CellType.triangle
Q_el = element("RT", msh.basix_cell(), k)
P_el = element("DG", msh.basix_cell(), k - 1)
V = fem.functionspace(msh, mixed_element([Q_el, P_el]))
(sigma, u) = TrialFunctions(V)
(tau, v) = TestFunctions(V)
V0 = V.sub(0) # RT元
V1 = V.sub(1) # P_k 元
Q, _ = V0.collapse()
P, _ = V1.collapse()
f = fem.Function(P)
f.interpolate(source_term)
sigma_D = fem.Function(Q)
sigma_D.interpolate(Exact_sigma)
tdim = msh.topology.dim
fdim = tdim - 1
facets = mesh.locate_entities_boundary(msh,fdim,clamped_boundary)
dofs = fem.locate_dofs_topological((V0, Q), fdim, facets)
bc = fem.dirichletbc(sigma_D, dofs, V0)
G1 = inner(sigma, tau) * dx +inner(u, div(tau)) * dx + inner(div(sigma), v) * dx
G2 = -inner(f, v) * dx
a = fem.form(G1)
L = fem.form(G2)
A = assemble_matrix(a, bcs=[bc])
A.assemble()
b = create_vector(L)
w_h = fem.Function(V)
solver = PETSc.KSP().create(msh.comm)
solver.setOperators(A)
solver.setType(PETSc.KSP.Type.PREONLY)
solver.getPC().setType("lu")
solver.getPC().setFactorSolverType("mumps")
# 设置MUMPS专门选项
opts = PETSc.Options()
opts["mat_mumps_icntl_14"] = 80 # 增加工作内存
opts["mat_mumps_icntl_24"] = 1 # 支持奇异矩阵
opts["mat_mumps_icntl_25"] = 0 # 支持奇异矩阵
opts["ksp_error_if_not_converged"] = 1 # 错误检查
solver.setFromOptions()
with b.localForm() as loc_b:
loc_b.set(0)
assemble_vector(b, L)
# 施加边界条件
apply_lifting(b, [a], [[bc]])
b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
set_bc(b, [bc])
# 求解线性系统
solver.solve(b, w_h.x.petsc_vec)
w_h.x.scatter_forward()
sigma_h = w_h.sub(0).collapse()
u_h = w_h.sub(1).collapse()
# 计算误差
V_ex = fem.functionspace(msh, ("DG", k+2))
Vv_ex = fem.functionspace(msh, ("RT", k+2))
# u 的误差
u_ex = fem.Function(V_ex)
u_ex.interpolate(Exactu)
sigma_ex = fem.Function(Vv_ex)
sigma_ex.interpolate(Exact_sigma)
# L2 误差
error_u_L2 = np.sqrt(msh.comm.allreduce(fem.assemble_scalar(fem.form(inner(u_h - u_ex, u_h - u_ex) * dx)), op=MPI.SUM))
error_sigma_L2 = np.sqrt(msh.comm.allreduce(fem.assemble_scalar(fem.form(inner(sigma_h - sigma_ex, sigma_h - sigma_ex) * dx)), op=MPI.SUM))
# 存储误差和网格尺寸
hs.append(h)
errors_u_L2.append(error_u_L2)
errors_sigma_L2.append(error_sigma_L2)
if MPI.COMM_WORLD.rank == 0:
print(f"网格尺寸: {n}x{n}, h = {h:.3f}")
print(f"原始变量u的L2误差: {error_u_L2:.2e}")
print(f"通量sigma的L2误差: {error_sigma_L2:.2e}")
print("-" * 50)
# 计算并输出收敛阶
if MPI.COMM_WORLD.rank == 0 and len(hs) > 1:
print("\n" + "="*60)
print("收敛阶分析:")
print("="*60)
print(f"{'网格尺寸':<12} {'h':<8} {'u误差':<12} {'u收敛阶':<10} {'σ误差':<12} {'σ收敛阶':<10}")
print("-" * 60)
for i in range(len(hs)):
if i == 0:
order_u = 0.0
order_sigma = 0.0
else:
order_u = np.log(errors_u_L2[i-1] / errors_u_L2[i]) / np.log(hs[i-1] / hs[i])
order_sigma = np.log(errors_sigma_L2[i-1] / errors_sigma_L2[i]) / np.log(hs[i-1] / hs[i])
n = 2 ** (i + 3)
print(f"{n}x{n:<8} {hs[i]:<8.3f} {errors_u_L2[i]:<12.2e} {order_u:<10.2f} {errors_sigma_L2[i]:<12.2e} {order_sigma:<10.2f}")
The following is my running result：
h error_u convergence of u error_sigma convergence of sigma
8x8 0.125 3.29e-01 0.00 2.95e-02 0.00
16x16 0.062 2.11e-01 0.64 1.47e-02 1.00
32x32 0.031 1.53e-01 0.46 7.37e-03 1.00
64x64 0.016 1.70e-01 -0.15 3.68e-03 1.00
I don’t know what the problem is. I’m looking forward to your reply and would be very grateful.