有限元算例-椭圆方程混合元求解#
本节将引导您使用 FEALPy 完成椭圆方程的混合有限元求解，现定义如下数学模型：
\[\begin{split}
\begin{cases} A\Delta u + c u = f & \text{in } \Omega= [0,1]^2, \\ \nabla u \cdot \boldsymbol{n} = 0 & \text{on } \partial\Omega. \end{cases}
\end{split}\]
引入通量\(\boldsymbol{p}=-A\nabla u\)，将原方程转化为一阶系统：
\[\begin{split}
\begin{cases} \boldsymbol{p} + A \nabla u = 0 & \text{in } \Omega, \\ \operatorname{div} \boldsymbol{p} + c u = f & \text{in } \Omega, \\ \boldsymbol{p} \cdot \boldsymbol{n} = 0 & \text{on } \partial\Omega. \end{cases}
\end{split}\]
其中我们定义真解为：
\[
u = \cos(2\pi x)\cos(2\pi y)
\]
扩散系数为：
\[\begin{split}
A = \begin{bmatrix} 10 & 0 \\ 0 & 10 \end{bmatrix}
\end{split}\]
反应系数为：
\[
c(x, y) = 2
\]
源项为：
\[
f(x,y)=-8\pi^2 \sin(2\pi x)\sin(2\pi y) + 2\cos(2\pi x)\cos(2\pi y) + 80\pi^2 \cos(2\pi x)\cos(2\pi y)
\]
混合变分形式为：
\[\begin{split}
\left\{
\begin{aligned}
&\int_\Omega A^{-1} \boldsymbol{p_h} \cdot \boldsymbol{v_h} \, dx - \int_\Omega u_h \, \operatorname{div} \boldsymbol{v_h} \, dx = 0, \quad \forall \boldsymbol{v_h} \in \boldsymbol{V_h}, \\
&\int_\Omega \operatorname{div} \boldsymbol{p_h} \, w_h \, dx + \int_\Omega c u_h w_h \, dx = \int_\Omega f w_h \, dx, \quad \forall w_h \in W_h.
\end{aligned}
\right.
\end{split}\]
1. 定义 PDE 模型#
from fealpy.backend import backend_manager as bm
from fealpy.decorator import cartesian
# 定义域
def domain():
return [0, 1, 0, 1]
# 扩散系数
@cartesian
def diffusion_coef(p):
val = bm.array([[10.0, 0.0], [0.0, 10.0]])
shape = p.shape[:-1] + val.shape
return bm.broadcast_to(val, shape)
# 扩散系数的逆
@cartesian
def diffusion_coef_inv(p):
val = bm.array([[0.1, 0.0], [0.0, 0.1]])
shape = p.shape[:-1] + val.shape
return bm.broadcast_to(val, shape)
# 反应系数
@cartesian
def reaction_coef(p):
val = bm.array([2.0])
shape = p.shape[:-1] + val.shape
return bm.broadcast_to(val, shape)
# 真解
@cartesian
def solution(p):
x = p[..., 0]
y = p[..., 1]
pi = bm.pi
val = bm.cos(2*pi*x)*bm.cos(2*pi*y)
return val
# 梯度
@cartesian
def gradient(p):
x = p[..., 0]
y = p[..., 1]
pi = bm.pi
val = bm.stack((
-2*pi*bm.sin(2*pi*x)*bm.cos(2*pi*y),
-2*pi*bm.cos(2*pi*x)*bm.sin(2*pi*y)), axis=-1)
return val
# 通量
@cartesian
def flux(p):
grad = gradient(p)
val = diffusion_coef(p)
val = bm.einsum('...ij, ...j->...i', val, -grad)
return val
# 源项
@cartesian
def source(p):
x = p[..., 0]
y = p[..., 1]
pi = bm.pi
val = 2*bm.cos(2*pi*x)*bm.cos(2*pi*y) + 80*pi**2*bm.cos(2*pi*x)*bm.cos(2*pi*y)
return val
2. 进行参数配置和初始化#
设置后端
from fealpy.backend import backend_manager as bm
backend = 'numpy'
device = 'cpu'
bm.set_backend(backend)
导入日志工具
from fealpy.utils import timer
from fealpy import logger
logger.setLevel('WARNING')
tmr = timer()
next(tmr)
设置初始网格和加密次数
from fealpy.mesh import TriangleMesh
mesh = TriangleMesh.from_box(domain(), nx=4, ny=4)
maxit = 4
定义误差存储矩阵
errorMatrix = bm.zeros((1, maxit), dtype=bm.float64)
3. 有限元求解#
流程包含：
(1)构建Lagrange有限元空间
(2)组装刚度矩阵 \(A\) 和载荷向量 \(F\)
(3)处理Dirichlet边界条件
(4)求解线性系统 \(A x_h = b\)
(5)计算L2误差 \(\|u - u_h\|_{L^2(\Omega)}\) \(\|\boldsymbol p - \boldsymbol p_h\|_{L^2(\Omega)}\)
(6)网格均匀加密
from fealpy.functionspace import RaviartThomasFESpace, LagrangeFESpace
from fealpy.fem import BilinearForm
from fealpy.fem import DivIntegrator, ScalarMassIntegrator
from fealpy.fem import LinearForm, ScalarSourceIntegrator
from fealpy.fem import BlockForm
from fealpy.solver import cg
for i in range(maxit):
uspace = LagrangeFESpace(mesh, p=1, ctype='D')
pspace = RaviartThomasFESpace(mesh, p=1)
tmr.send(f'第{i}次空间时间')
uh = uspace.function()
ph = pspace.function()
uGDOF = uspace.number_of_global_dofs()
pGDOF = pspace.number_of_global_dofs()
xh = bm.zeros((uGDOF+pGDOF,), dtype=bm.float64)
bform1 = BilinearForm(pspace)
bform1.add_integrator(ScalarMassIntegrator(coef=diffusion_coef_inv, q=3))
bform2 = BilinearForm((uspace,pspace))
bform2.add_integrator(DivIntegrator(coef=-1, q=3))
bform3 = BilinearForm((uspace,pspace))
bform3.add_integrator(DivIntegrator(coef=1, q=3))
bform4 = BilinearForm(uspace)
bform4.add_integrator(ScalarMassIntegrator(coef=reaction_coef, q=3))
M = BlockForm([[bform1,bform2],
[bform3.T,bform4]])
A = M.assembly()
lform = LinearForm(uspace)
lform.add_integrator(ScalarSourceIntegrator(source=source))
F = lform.assembly()
G = bm.zeros(pGDOF)
b = bm.concatenate([G,F],axis=0)
tmr.send(f'第{i}次矩阵组装时间')
G_apply = pspace.set_neumann_bc(solution)
F = bm.zeros(uGDOF, dtype=bm.float64)
b_apply = bm.concatenate([G_apply,F],axis=0)
b = b - b_apply
tmr.send(f'第{i}次边界处理时间')
from fealpy.solver import spsolve
xh[:] = spsolve(A, b, solver= 'scipy')
ph[:] = xh[:pGDOF]
uh[:] = xh[pGDOF:]
tmr.send(f'第{i}次求解器时间')
errorMatrix[0, i] = mesh.error(solution, uh)
if i < maxit-1:
mesh.uniform_refine(n=1)
tmr.send(f'第{i}次误差计算及网格加密时间')
4. 误差分析和收敛阶计算#
next(tmr)
print("最终误差",errorMatrix)
print("order : ", bm.log2(errorMatrix[0,:-1]/errorMatrix[0,1:]))
Timer received None and paused.
=================================================
ID Time Proportion(%) Label
-------------------------------------------------
1 811.782 [ms] 75.004 第0次空间时间
2 9.002 [ms] 0.832 第0次矩阵组装时间
3 0.000 [us] 0.000 第0次边界处理时间
4 1.000 [ms] 0.092 第0次求解器时间
5 1.000 [ms] 0.092 第0次误差计算及网格加密时间
6 0.000 [us] 0.000 第1次空间时间
7 10.116 [ms] 0.935 第1次矩阵组装时间
8 0.000 [us] 0.000 第1次边界处理时间
9 1.999 [ms] 0.185 第1次求解器时间
10 1.001 [ms] 0.093 第1次误差计算及网格加密时间
11 0.000 [us] 0.000 第2次空间时间
12 27.983 [ms] 2.585 第2次矩阵组装时间
13 0.000 [us] 0.000 第2次边界处理时间
14 11.003 [ms] 1.017 第2次求解器时间
15 2.002 [ms] 0.185 第2次误差计算及网格加密时间
16 1.001 [ms] 0.092 第3次空间时间
17 124.813 [ms] 11.532 第3次矩阵组装时间
18 1.964 [ms] 0.181 第3次边界处理时间
19 76.053 [ms] 7.027 第3次求解器时间
20 1.598 [ms] 0.148 第3次误差计算及网格加密时间
=================================================
最终误差 [[0.07434676 0.0195544 0.00495466 0.00124288]]
order : [1.9267765 1.98063669 1.99509364]
5. 结果可视化#
在单元重心处计算真解和数值解，并进行可视化比较
from matplotlib import pyplot as plt
bc = bm.array([[1/3, 1/3, 1/3]], dtype=bm.float64)
ps = mesh.bc_to_point(bc)
u = solution(ps)
uh = uh(bc)
fig, axes = plt.subplots(1, 2)
mesh.add_plot(axes[0], cellcolor=u, linewidths=0)
axes[0].set_title('真解', fontname='Microsoft YaHei')
mesh.add_plot(axes[1], cellcolor=uh, linewidths=0)
axes[1].set_title('数值解', fontname='Microsoft YaHei')
plt.show()
