计算3-D立方体的最快算法

Question

我想写一段代码，用周期性边界条件，计算向量场的旋度，数值到二阶。然而，我做的算法非常慢，我想知道是否有人知道任何替代算法。你知道吗

给出更具体的上下文：我使用3xAxBxC numpy数组作为向量场，其中第一个轴表示笛卡尔方向（x，y，z），a，B，C表示笛卡尔方向上的箱子数（即分辨率）。举个例子，向量场F=np.零（（3,64,64,64）），其中Fx=F[0]本身就是一个64x64x64笛卡尔晶格。到目前为止，我的解决方案是使用以3点为中心的差分模板来计算导数，并使用嵌套循环来迭代所有不同的维度，使用模运算来强制执行周期性边界条件（参见下面的示例）。然而，随着分辨率的增加（A，B，C的大小），这开始需要很长的时间（最多2分钟，如果我在模拟中这样做几百次的话，加起来——这只是一个较大算法的一小部分）。我想知道是否有人知道做这件事的另一种方法？你知道吗

import numpy as np

F =np.array([np.ones([128,128,128]),2*np.ones([128,128,128]),
               3*np.ones([128,128,128])])


VxF =np.array([np.zeros([128,128,128]),np.zeros([128,128,128]),
               np.zeros([128,128,128])])

for i in range(0,128):
     for j in range(0,128):
          for k in range(0,128):
           VxF[0][i,j,k] = 0.5*((F[2][i,(j+1)%128,k]-
                    F[2][i,j-1,k])-(F[1][i,j,(k+1)%128]-F[1][i,j,k-1]))
           VxF[1][i,j,k] = 0.5*((F[0][i,j,(k+1)%128]-
                    F[0][i,j,k-1])-(F[2][(i+1)%128,j,k]-F[2][i-1,j,k]))
           VxF[2][i,j,k] = 0.5*((F[1][(i+1)%128,j,k]-
                    F[1][i-1,j,k])-(F[0][i,(j+1)%128,k]-F[0][i,j-1,k]))

为了再次迭代，我正在寻找一种算法，它能在给定周期边界条件下，以比我现有的算法更快的速度计算向量场数组到二阶的旋度。也许没有什么可以做到这一点，但我只想在我继续运行这个算法之前检查一下。谢谢。提前通知大家！你知道吗

Answer 1

可能有更好的工具来实现这一点，但是这里有一个numba的普通200x加速：

import numpy as np
from numba import jit

def pure_python():
    F =np.array([np.ones([128,128,128]),2*np.ones([128,128,128]),
                3*np.ones([128,128,128])])


    VxF =np.array([np.zeros([128,128,128]),np.zeros([128,128,128]),
                np.zeros([128,128,128])])

    for i in range(0,128):
        for j in range(0,128):
            for k in range(0,128):
                VxF[0][i,j,k] = 0.5*((F[2][i,(j+1)%128,k]-
                            F[2][i,j-1,k])-(F[1][i,j,(k+1)%128]-F[1][i,j,k-1]))
                VxF[1][i,j,k] = 0.5*((F[0][i,j,(k+1)%128]-
                            F[0][i,j,k-1])-(F[2][(i+1)%128,j,k]-F[2][i-1,j,k]))
                VxF[2][i,j,k] = 0.5*((F[1][(i+1)%128,j,k]-
                            F[1][i-1,j,k])-(F[0][i,(j+1)%128,k]-F[0][i,j-1,k]))

    return VxF

@jit(fastmath=True)
def with_numba():
    F =np.array([np.ones([128,128,128]),2*np.ones([128,128,128]),
                3*np.ones([128,128,128])])


    VxF =np.array([np.zeros([128,128,128]),np.zeros([128,128,128]),
                np.zeros([128,128,128])])

    for i in range(0,128):
        for j in range(0,128):
            for k in range(0,128):
                VxF[0][i,j,k] = 0.5*((F[2][i,(j+1)%128,k]-
                            F[2][i,j-1,k])-(F[1][i,j,(k+1)%128]-F[1][i,j,k-1]))
                VxF[1][i,j,k] = 0.5*((F[0][i,j,(k+1)%128]-
                            F[0][i,j,k-1])-(F[2][(i+1)%128,j,k]-F[2][i-1,j,k]))
                VxF[2][i,j,k] = 0.5*((F[1][(i+1)%128,j,k]-
                            F[1][i-1,j,k])-(F[0][i,(j+1)%128,k]-F[0][i,j-1,k]))

    return VxF

纯Python版本在我的机器上需要13秒，而numba版本需要65毫秒