计算 3-D 旋度的最快算法

Question

我正在尝试编写一段代码，用于在周期性边界条件下以二阶数值计算向量场的旋度。然而，我制作的算法非常慢，我想知道是否有人知道任何替代算法。

给出更具体的上下文：我使用 3xAxBxC numpy 数组作为向量场，其中第一个轴指的是笛卡尔方向 (x,y,z)，A、B、C 指的是该笛卡尔方向上的 bin 数量（即决议）。例如，我可能有一个向量场 F = np.zeros((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

可能有更好的工具可以实现这一点，但这里有一个微不足道的 200 倍加速numba：

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 毫秒。