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

Question

我有一个实现二维拉普拉斯算子的代码 finite differences integration method for partial differential equations, using the roll method of Numpy :

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 方法？

Answer 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

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

Implementing 2D Laplacian in Cython for periodic boundary counditions

python

numpy

cython

python-2.7