2

我正在尝试编写一段代码，以数字方式计算矢量场的卷曲到具有周期性边界条件的二阶。但是，我制作的算法非常慢，我想知道是否有人知道任何替代算法。

提供更具体的上下文：我使用 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]))

只是为了重新迭代，我正在寻找一种算法，该算法将比我拥有的周期性边界条件更快地计算向量场数组的卷曲到二阶。也许没有什么可以做到这一点，但我只是想在我继续花时间运行这个算法之前检查一下。感谢。大家提前！

4

1 回答 1

5

可能有更好的工具，但这里有一个微不足道的 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 毫秒。

于 2019-08-05T18:25:53.573 回答