FEALPy 调试
首先放上能够正常运行的 2维 Poisson 方程的程序源代码:
# 导入并创建 PDE 模型
from fealpy.pde.poisson_2d import CosCosData
from fealpy.mesh import MeshFactory as MF
from fealpy.functionspace import LagrangeFiniteElementSpace
from fealpy.boundarycondition import DirichletBC
from scipy.sparse.linalg import spsolve
import matplotlib.pyplot as plt
pde = CosCosData()
# 创建网格
domain = [0, 1, 0, 1]
mesh = MF.boxmesh2d(domain, nx=15, ny=15, meshtype='tri')
# mesh.uniform_refine(n=1)
# 建立 Lagrange 有限元函数空间
space = LagrangeFiniteElementSpace(mesh, p=2)
# 创建 Dirichlet 边界条件对象
bc = DirichletBC(space, pde.dirichlet)
# 创建离散代数系统
node = mesh.entity('node')
uh = space.function()
A = space.stiff_matrix()
F = space.source_vector(pde.source)
# 处理边界条件
A, F = bc.apply(A, F, uh)
# 求解离散代数系统
uh[:] = spsolve(A, F)
print('uh[0:10]:', uh[0:10])
# 计算 L2 和 H1 误差
L2Error = space.integralalg.L2_error(pde.solution, uh)
H1Error = space.integralalg.L2_error(pde.gradient, uh.grad_value)
print('L2Error:', L2Error, '\tH1Error:', H1Error)
# 绘制数值解图像和网格图像
fig = plt.figure()
axes = fig.add_subplot(1, 2, 1, projection='3d')
uh.add_plot(axes, cmap='rainbow')
axes = fig.add_subplot(1, 2, 2)
mesh.add_plot(axes)
plt.show()
运行结果:
uh[0:10]: [ 1. 0.9781476 0.91354546 0.80901699 0.66913061 0.5
0.30901699 0.10452846 -0.10452846 -0.30901699]
L2Error: 8.352207463706391e-05 H1Error: 0.009575809898399913
我的需求是根据上述 Poisson 方程求解程序修改成一个求解抛物型方程的例子
我修改的程序如下(不能成功运行):
# 导入并创建 PDE 模型
from fealpy.pde.parabolic_model_2d import SinSinExpData # 改动
from fealpy.mesh import MeshFactory as MF
from fealpy.functionspace import LagrangeFiniteElementSpace
from fealpy.boundarycondition import DirichletBC
from scipy.sparse.linalg import spsolve
import matplotlib.pyplot as plt
pde = SinSinExpData()
# 创建网格
domain = [0, 1, 0, 1]
mesh = MF.boxmesh2d(domain, nx=15, ny=15, meshtype='tri')
timeline = pde.time_mesh(0, 1, 100) # 添加(好像是 fealpy 时间离散的函数)
# 建立 Lagrange 有限元函数空间
space = LagrangeFiniteElementSpace(mesh, p=2)
# 创建 Dirichlet 边界条件对象
bc = DirichletBC(space, pde.dirichlet)
# 创建离散代数系统
uh = space.function()
A = space.stiff_matrix()
F = space.source_vector(pde.source)
# 处理边界条件
A, F = bc.apply(A, F, uh)
# 求解离散代数系统
uh[:] = spsolve(A, F)
print('uh[0:10]:', uh[0:10])
# 计算 L2 和 H1 误差
L2Error = space.integralalg.L2_error(pde.solution, uh)
H1Error = space.integralalg.L2_error(pde.gradient, uh.grad_value)
print('L2Error:', L2Error, '\tH1Error:', H1Error)
# 绘制数值解图像和网格图像
fig = plt.figure()
axes = fig.add_subplot(1, 2, 1, projection='3d')
uh.add_plot(axes, cmap='rainbow')
axes = fig.add_subplot(1, 2, 2)
mesh.add_plot(axes)
plt.show()
报错信息:
run boxmesh2d with time: 0.000916731000000004
run serial_construct_matrix with time: 0.018699856999999986
Traceback (most recent call last):
File "/Users/dasha/Desktop/Developer/Python之旅/项目/微分方程数值解测试/Poisson方程/demo1_20.py", line 22, in <module>
F = space.source_vector(pde.source)
File "/Users/dasha/Library/Python/3.7/lib/python/site-packages/fealpy/decorator/coordinates.py", line 13, in add_attribute
return func(*args, **kwargs)
File "/Users/dasha/Library/Python/3.7/lib/python/site-packages/fealpy/functionspace/LagrangeFiniteElementSpace.py", line 975, in source_vector
fval = f(pp)
File "/Users/dasha/Library/Python/3.7/lib/python/site-packages/fealpy/decorator/coordinates.py", line 13, in add_attribute
return func(*args, **kwargs)
TypeError: source() missing 1 required positional argument: 't'
根据报错提示第 22 行 F = space.source_vector(pde.source)中source() 缺少 1 个必需的位置参数：'t'. 定位到 source 函数位置:
@cartesian
def source(self, p, t):
""" The right hand side of Possion equation
INPUT:
p: array object, N*2
"""
return (2*np.pi**2-1)*self.solution(p, t)
疑惑 1 : 发现 source() 需要 p,t 两个参数, 但报错提示信息却只提示“缺少 1 个必需的位置参数：'t'”?
于是我去查看 Poisson 方程的对应位置
sourc()源代码如下:@cartesian def source(self, p): """ The right hand side of Possion equation INPUT: p: array object, """ x = p[..., 0] y = p[..., 1] val = -2*(x**2 + y**2) return val
疑惑 2 : 为什么 Poisson 方程求解程序中没有传入参数 p 直接运行 F = space.source_vector(pde.source) 就能不报错成功运行? 我进一步 Dbug 发现 p 在 source() 中被调用时已经被赋值
Dbug 截图
p的详细信息p的部分内容
疑惑 3 : p 到底是在什么时候被赋值的, 并且按照我的理解, p 应该代表的是空间方向的离散求解区域(即此处的 x, y), 而求解区域’domain = [0, 1, 0, 1]’, 如此的话为什么 p 的两个轴的取值不是
[
0
,
1
]
[0,1]
[0,1] ?(此处我猜测可能是给的重心坐标)
由于未能清晰以上疑惑, 同时想进一步测试代码, 于是我根据报错信息 TypeError: source() missing 1 required positional argument: 't', 直接将 F = space.source_vector(pde.source) 修改成 F = space.source_vector(pde.source(t = 0.5)),再次运行抛物型方程求解程序
报错信息:
run boxmesh2d with time: 0.0008686629999999917
Traceback (most recent call last):
File "/Users/dasha/Desktop/Developer/Python之旅/项目/微分方程数值解测试/Poisson方程/demo1_20.py", line 22, in <module>
F = space.source_vector(pde.source(t = 0.5))
File "/Users/dasha/Library/Python/3.7/lib/python/site-packages/fealpy/decorator/coordinates.py", line 13, in add_attribute
return func(*args, **kwargs)
TypeError: source() missing 1 required positional argument: 'p'
run serial_construct_matrix with time: 0.018582748000000038
疑惑 4 : 为什么现在 p 又成了“必需的位置参数”?
好吧, 那我大不了再把 p 赋值给你罢了, 于是我使用获取节点坐标的代码 node = mesh.entity('node'), 并将 node 传入函数,修改细节如下:
# 导入并创建 PDE 模型
from fealpy.pde.parabolic_model_2d import SinSinExpData # 改动
from fealpy.mesh import MeshFactory as MF
from fealpy.functionspace import LagrangeFiniteElementSpace
from fealpy.boundarycondition import DirichletBC
from scipy.sparse.linalg import spsolve
import matplotlib.pyplot as plt
pde = SinSinExpData()
# 创建网格
domain = [0, 1, 0, 1]
mesh = MF.boxmesh2d(domain, nx=15, ny=15, meshtype='tri') # 添加
timeline = pde.time_mesh(0, 1, 100) # 添加(好像是 fealpy 时间离散的函数)
# 建立 Lagrange 有限元函数空间
space = LagrangeFiniteElementSpace(mesh, p=2)
# 创建 Dirichlet 边界条件对象
bc = DirichletBC(space, pde.dirichlet)
# 创建离散代数系统
uh = space.function()
A = space.stiff_matrix()
F = space.source_vector(pde.source(node, t = 0.5)) # 修改
# 处理边界条件
A, F = bc.apply(A, F, uh)
# 求解离散代数系统
uh[:] = spsolve(A, F)
print('uh[0:10]:', uh[0:10])
# 计算 L2 和 H1 误差
L2Error = space.integralalg.L2_error(pde.solution, uh)
H1Error = space.integralalg.L2_error(pde.gradient, uh.grad_value)
print('L2Error:', L2Error, '\tH1Error:', H1Error)
# 绘制数值解图像和网格图像
fig = plt.figure()
axes = fig.add_subplot(1, 2, 1, projection='3d')
uh.add_plot(axes, cmap='rainbow')
axes = fig.add_subplot(1, 2, 2)
mesh.add_plot(axes)
plt.show()
再次运行抛物型方程求解程序.
报错信息:
run boxmesh2d with time: 0.0009375699999999432
run serial_construct_matrix with time: 0.017745607999999913
Traceback (most recent call last):
File "/Users/dasha/Desktop/Developer/Python之旅/项目/微分方程数值解测试/时间分数阶Fisher方程/demo1_20.py", line 22, in <module>
F = space.source_vector(pde.source(node, t = 0.5))
File "/Users/dasha/Library/Python/3.7/lib/python/site-packages/fealpy/decorator/coordinates.py", line 13, in add_attribute
return func(*args, **kwargs)
File "/Users/dasha/Library/Python/3.7/lib/python/site-packages/fealpy/functionspace/LagrangeFiniteElementSpace.py", line 973, in source_vector
if f.coordtype == 'cartesian':
AttributeError: 'numpy.ndarray' object has no attribute 'coordtype'
定位到 source_vector() 函数:
@cartesian
def source_vector(self, f, dim=None, q=None):
p = self.p
cellmeasure = self.cellmeasure
bcs, ws = self.integrator.get_quadrature_points_and_weights()
if f.coordtype == 'cartesian':
pp = self.mesh.bc_to_point(bcs)
fval = f(pp)
elif f.coordtype == 'barycentric':
fval = f(bcs)
gdof = self.number_of_global_dofs()
shape = gdof if dim is None else (gdof, dim)
b = np.zeros(shape, dtype=self.ftype)
if p > 0:
if type(fval) in {float, int}:
if fval == 0.0:
return b
else:
phi = self.basis(bcs)
bb = np.einsum('m, mik, i->ik...',
ws, phi, self.cellmeasure)
bb *= fval
else:
phi = self.basis(bcs)
bb = np.einsum('m, mi..., mik, i->ik...',
ws, fval, phi, self.cellmeasure)
cell2dof = self.cell_to_dof() #(NC, ldof)
if dim is None:
np.add.at(b, cell2dof, bb)
else:
np.add.at(b, (cell2dof, np.s_[:]), bb)
else:
b = np.einsum('i, ik..., k->k...', ws, fval, cellmeasure)
return b
报错 AttributeError：“numpy.ndarray”对象没有属性“coordtype” 我更懵了, 这个装饰子报错我并不知道从何处下手修改.
于是考虑换一种修改方法, 根据之前的报错信息 TypeError: source() missing 1 required positional argument: 't', 我放弃将 F = space.source_vector(pde.source) 修改成 F = space.source_vector(pde.source(t = 0.5)) 这种修改方式, 改变成直接在 source() 中初始化 t =0.5,修改细节如下:
@cartesian
def source(self, p, t = 0.5): # 修改
""" The right hand side of Possion equation
INPUT:
p: array object, N*2
"""
return (2*np.pi**2-1)*self.solution(p, t)
再次运行抛物型方程求解程序.
报错信息:
run boxmesh2d with time: 0.0010344709999999813
run serial_construct_matrix with time: 0.02001516999999997
Traceback (most recent call last):
File "/Users/dasha/Desktop/Developer/Python之旅/项目/微分方程数值解测试/Poisson方程/demo1_20.py", line 24, in <module>
A, F = bc.apply(A, F, uh)
File "/Users/dasha/Library/Python/3.7/lib/python/site-packages/fealpy/boundarycondition/BoundaryCondition.py", line 27, in apply
isDDof = space.set_dirichlet_bc(uh, gD, threshold=threshold)
File "/Users/dasha/Library/Python/3.7/lib/python/site-packages/fealpy/functionspace/LagrangeFiniteElementSpace.py", line 1013, in set_dirichlet_bc
uh[isDDof] = gD(ipoints[isDDof])
File "/Users/dasha/Library/Python/3.7/lib/python/site-packages/fealpy/decorator/coordinates.py", line 13, in add_attribute
return func(*args, **kwargs)
TypeError: dirichlet() missing 1 required positional argument: 't'
根据错误信息 TypeError: dirichlet() missing 1 required positional argument: 't' 发现仍然是参数缺失, 但仔细一看发现原来的报错位置 F = space.source_vector(pde.source) 经过修改已经能够成功运行, 而报错位置变成了 A, F = bc.apply(A, F, uh), 于是定位到错误位置cartesian()
def cartesian(func):
@wraps(func)
def add_attribute(*args, **kwargs):
return func(*args, **kwargs)
add_attribute.__dict__['coordtype'] = 'cartesian'
return add_attribute
疑惑 5: 报错信息提示 dirichlet() missing 1 required positional argument: 't', 但我并找到 dirichlet() 函数在哪里, 然后也不知道改修改哪里才能更进一步了.
辛苦大哥了, 感恩的心❤️