在 Cython 中为周期性边界条件实现 2D 拉普拉斯算子

【问题标题】：Implementing 2D Laplacian in Cython for periodic boundary counditions在 Cython 中为周期性边界条件实现 2D 拉普拉斯算子
【发布时间】：2015-05-29 02:44:47
【问题描述】：

我有一个代码，使用roll method of Numpy 为偏微分方程的finite differences 积分方法实现二维拉普拉斯算子：

def lapOp(u):
    """
    This is the laplacian operator on 2D array
    of stencil of 4th accuracy terms
    """
    lap = ((4.0/3.0)*np.roll(u,1,axis=0) + (4.0/3.0)*np.roll(u,-1,axis=0) + (4.0/3.0)*np.roll(u,1,axis=1) + (4.0/3.0)*np.roll(u,-1,axis=1) -5.0*u)
    lap -= ((1.0/12.0)*np.roll(u,2,axis=0) + (1.0/12.0)*np.roll(u,-2,axis=0) + (1.0/12.0)*np.roll(u,2,axis=1) + (1.0/12.0)*np.roll(u,-2,axis=1))
    lap = lap / hh
    return lap

我想对我的代码进行 cythonize - roll 方法可以在我的 pyx 代码中工作，还是应该使用 C 实现 roll 方法？

【问题讨论】：

标签： python python-2.7 numpy cython

【解决方案1】：

简短的回答是：roll 可以在 Cython 中工作，但不会更快（任何？）。

如果你想要速度，你可能应该完全避免使用类似roll 的东西（它很慢，因为每次调用它都会创建一个完整的副本），而是使用索引来获取大块 numpy 数组 u 的视图。您不应该需要 Cython，并且可能不会从中受益。

下面是一个不完整的示例（足以说明要点）：

def lapOp(u):
    lap = np.empty_like(u)
    # this bit is equivalent to (4.0/3)*np.roll(u,1,axis=0)
    lap[1:,:] = (4.0/3.0)*u[:-1,:]
    lap[0,:] = (4.0/3.0)*u[-1,:]

    # add (4.0/3)*np.roll(u,-1,axis=0)
    lap[:-1,:] += (4.0/3.0)*u[1:,:]
    lap[-1,:] += (4.0/3.0)*u[0,:]

    # add (4.0/3)*np.roll(u,1,axis=1)
    lap[:,1:] += (4.0/3.0)*u[:,:-1]
    lap[:,0] += (4.0/3.0)*u[:,-1]

    # the remainder is left as a rather tedious exercise for the reader

    return lap/hh

【讨论】：