FEALPy 调试

最新推荐文章于 2023-03-23 15:46:27 发布

图灵猫

最新推荐文章于 2023-03-23 15:46:27 发布

阅读量1.2k

点赞数

分类专栏： FEALPy 文章标签：偏微分方程有限元

本文链接：https://blog.csdn.net/qq_42818403/article/details/122686097

版权

FEALPy 专栏收录该内容

4 篇文章 7 订阅

订阅专栏

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

Figure_1

我的需求是根据上述 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]$ ?(此处我猜测可能是给的重心坐标)

由于未能清晰以上疑惑, 同时想进一步测试代码, 于是我根据报错信息 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() 函数在哪里, 然后也不知道改修改哪里才能更进一步了.

辛苦大哥了, 感恩的心❤️

图灵猫

关注

0
点赞
踩
4

收藏

觉得还不错? 一键收藏
打赏
2
评论
FEALPy 调试

FEALPy 调试首先放上能够正常运行的 2维 Poisson 方程的程序源代码:# 导入并创建 PDE 模型from fealpy.pde.poisson_2d import CosCosDatafrom fealpy.mesh import MeshFactory as MFfrom fealpy.functionspace import LagrangeFiniteElementSpacefrom fealpy.boundarycondition import DirichletBCfr
复制链接

扫一扫