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

Question

我正在尝试编写一段代码，以数值方式计算具有周期性边界条件的二阶矢量场的旋度。但是，我做的算法很慢，我想知道是否有人知道任何替代算法。

为了提供更具体的上下文：我使用 3xAxBxC numpy 数组作为我的矢量场，其中第一个轴指的是笛卡尔方向 (x,y,z)，A、B、C 指的是在那个笛卡尔方向上的垃圾箱（即分辨率）。因此，例如，我可能有一个向量场 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 毫秒。

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

Fastest algorithm for computing 3-D curl

python

algorithm

physics

numerical-computing

numerical-methods