Hi, Everyone!
Currently, I am decoupling and solving a 3D model:
(\partial B)/(\partial t)+\curl E=0 in \Omega B\cdot n=0 on \partial \Omega
When I take the divergence of this equation, I get：
\div B=0 in \Omega
Discretize this equation and obtain the numerical scheme as follows：
（B_h^n-B_h^{n-1},C)+(\curl E_h^n,C)=0
It can be deduced from the numerical scheme：
\div B_h^n= \div B_h^{n-1}=...=\div B_h^0
When the initial condition \div B_h^0=0 is given, it can be obtained that \div B_h^n=0.
from petsc4py import PETSc
from mpi4py import MPI
import numpy as np
import ufl
from dolfinx import mesh, fem, default_real_type
from dolfinx.fem.petsc import assemble_matrix, assemble_vector, apply_lifting, create_vector, set_bc
import time
from ufl import dx, grad, dot, lhs, rhs,inner,cross,curl,div,Constant,as_vector
from basix.ufl import mixed_element,element
from numpy import sin, cos, pi, exp
from ufl import sin as ufl_sin, cos as ufl_cos, pi as ufl_pi
# 精确解
Rm=1
left =0.0
right = 1.0
bottom = 0.0
top = 1.0
near = 0.0
far = 1.0
def ExactU(x,t):
factor = exp(-t) + 1.0
temp1=factor * sin(pi * x[0]) *sin(pi * x[1]) *sin(pi * x[2])
temp2=np.zeros_like(x[0])
temp3=factor * (x[2]**2 - x[2]) *sin(pi * x[0]) * sin(pi * x[1])
return [temp1,temp2,temp3]
def ExactE(x,t):
sint = sin(pi*t)+1.5
temp1=sint*(x[1]**2-x[1])*sin(pi*x[2])
temp2=sint*(x[0]**2-x[0])*sin(pi*x[2])
temp3=sint*(x[0]**2-x[0])*sin(pi*x[1])
return [temp1,temp2,temp3]
def ExactB1(x,t):
cost=(1/ufl_pi)*cos(ufl_pi*t)-1.5*t
temp1=cost*ufl_pi*(x[0]**2-x[0])*(ufl_cos(ufl_pi*x[1])-ufl_cos(ufl_pi*x[2]))
temp2=cost*(ufl_pi*(x[1]**2-x[1])*ufl_cos(ufl_pi*x[2])-(2*x[0]-1)*ufl_sin(ufl_pi*x[1]))
temp3=cost*(2*x[0]-2*x[1])*ufl_sin(ufl_pi*x[2])
return as_vector([temp1, temp2, temp3])
def ExactB(x,t):
cost=(1/pi)*cos(pi*t)-1.5*t
temp1=cost*pi*(x[0]**2-x[0])*(cos(pi*x[1])-cos(pi*x[2]))
temp2=cost*(pi*(x[1]**2-x[1])*cos(pi*x[2])-(2*x[0]-1)*sin(pi*x[1]))
temp3=cost*(2*x[0]-2*x[1])*sin(pi*x[2])
return [temp1, temp2, temp3]
# 固定终止时间
T = 1 / 2
# 结果记录
hs = []
errors_L2 = []
errors_max = []
div_Bh=[]
start_time = time.time()
for i in range(3):
n = 2 ** (i + 2) # 空间划分数 n = 2^{i+2}
h = 1.0 / n
num_steps = 2 ** (2 * i + 3) # 时间划分数 tau = 2^{2i+4}
dt = T / num_steps
t = 0.0
# 创建 3D 网格（单位立方体）
domain = mesh.create_unit_cube(MPI.COMM_WORLD, n, n, n, mesh.CellType.tetrahedron)
x = ufl.SpatialCoordinate(domain)
# 磁场空间与电场空间的定义
VB = fem.functionspace(domain, ("RT", 1)) # degree=1, 表示 RT0元
VE = fem.functionspace(domain, ("N1curl", 1)) # degree=1, 表示 ND0元（最低阶棱单元）
E_n=fem.Function(VE)
B_n = fem.Function(VB)
B_n.interpolate(lambda x: ExactB(x, t))
B_D = fem.Function(VB)
B_D.interpolate(lambda x: ExactB(x, t))
# 施加 Dirichlet 边界条件（3D 边界）
tdim = domain.topology.dim
fdim = tdim - 1
domain.topology.create_connectivity(fdim, tdim)
boundary_facets = mesh.exterior_facet_indices(domain.topology)
bc = fem.dirichletbc(B_D, fem.locate_dofs_topological(VB, fdim, boundary_facets))
# 变分形式（3D 形式相同，但 dx 现在是 3D 积分）
B, C = ufl.TrialFunction(VB), ufl.TestFunction(VB)
G = inner(B,C)*dx - inner(B_n,C)*dx + dt*inner(curl(E_n),C)*dx
a = fem.form(lhs(G))
L = fem.form(rhs(G))
# 组装矩阵和向量
A = assemble_matrix(a, bcs=[bc])
A.assemble()
b = create_vector(L)
B_h =fem.Function(VB)
# 定义求解器（直接法：LU 分解）
solver = PETSc.KSP().create(domain.comm)
solver.setOperators(A)
# solver.setType(PETSc.KSP.Type.PREONLY)
# solver.getPC().setType(PETSc.PC.Type.LU)
solver.setType(PETSc.KSP.Type.BCGS)
pc1 = solver.getPC()
pc1.setType(PETSc.PC.Type.HYPRE)
pc1.setHYPREType("boomeramg")
# 时间步进
for n_step in range(num_steps):
t += dt
E_n.interpolate(lambda x: ExactE(x,t))
B_D.interpolate(lambda x: ExactB(x, t))
bc = fem.dirichletbc(B_D, fem.locate_dofs_topological(VB, fdim, boundary_facets))
# 更新矩阵
A.zeroEntries()
assemble_matrix(A, a, bcs=[bc])
A.assemble()
# 重新组装右端向量
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, B_h.x.petsc_vec)
B_h.x.scatter_forward()
B_n.x.array[:] = B_h.x.array
# 计算误差（在更高阶空间）
# V_dg = fem.functionspace(domain, ("Discontinuous Lagrange", 2, (3,)))
V_dg = fem.functionspace(domain, ("Lagrange", 1, (3,)))
B_ex = fem.Function(V_dg)
B_ex.interpolate(lambda x: ExactB(x,t))
B_ex1 = fem.Function(VB)
B_ex1.interpolate(lambda x:ExactB(x,t))
# L2 误差
error_L2 = np.sqrt(domain.comm.allreduce(fem.assemble_scalar(fem.form((B_h - B_ex)**2 * ufl.dx)), op=MPI.SUM))
div_Bh = np.sqrt(domain.comm.allreduce(fem.assemble_scalar(fem.form(div(B_h)**2 * dx)), op=MPI.SUM))
error_max= np.max(np.abs(B_ex1.x.array - B_h.x.array))
if MPI.COMM_WORLD.rank == 0:
hs.append(h)
errors_L2.append(error_L2)
errors_max.append(error_max)
print(f"[i={i}] n={n}, num_steps={num_steps}, dt={dt:.3e}, h={h:.3e}")
print(f" → L2 error = {error_L2:.3e}, max error={error_max:.3e}, divBh={div_Bh:.3e} ")
end_time = time.time() # 计时结束
elapsed = end_time - start_time
# 输出收敛阶
if MPI.COMM_WORLD.rank == 0:
print("\nEstimated convergence rates (L2):")
for j in range(1, len(hs)):
r1 = np.log(errors_L2[j-1]/errors_L2[j]) / np.log(hs[j-1]/hs[j])
r2 = np.log(errors_max[j-1]/errors_max[j]) / np.log(hs[j-1]/hs[j])
print(f"h={hs[j-1]:.3e} → h={hs[j]:.3e}, rate ≈ {r1:.2f},rate ≈ {r2:.2f} ")
print(f"\nTotal runtime: {elapsed:.2f} seconds")
The code corresponding to B_h^0 is：
VB = fem.functionspace(domain, ("RT", 1)) # degree=1, 表示 RT0元
VE = fem.functionspace(domain, ("N1curl", 1)) # degree=1, 表示 ND0元（最低阶棱单元）
E_n=fem.Function(VE)
B_n = fem.Function(VB)
B_n.interpolate(lambda x: ExactB(x, t)
But I calculated that \div B_h^0 is not close to 0，I hope someone can help me solve this problem. How can this initial condition be guaranteed?
Looking forward to your reply！