平流扩散模型二维阵列的迭代时间过长

Question

我试图用Python来模拟二维平流扩散，结果遇到了麻烦。

我使用指示符和嵌套for循环来操作包含特定粒子浓度的2d数组中的值，现在要花费很长时间来迭代。是否有更有效的迭代方法？

事先非常感谢，

# Initialize u
    u_UPW = np.zeros((Nx+1,Ny+1))

    # Initial conditions  
    u_UPW[:] = u0

    # Time loop
    for n in range(num_timesteps):

        # Impose boundary conditions
        # Dirichlet
        u_UPW[range(0, int(Ny/3)),0] = u_dir[0]
        u_UPW[range(int(Ny/3), int(Ny*2/3)+1),0] = u_dir[1]
        u_UPW[range(int(Ny*2/3)+1, Ny+1),0] = u_dir[2]
        
        #Neumann
        u_UPW[Ny-1,:] = u_UPW[Ny,:] + dy*u_neu[0]
        u_UPW[1,:] = u_UPW[0,:] - dy*u_neu[1]
        u_UPW[:,Nx-1] = u_UPW[:,Nx] - dx*u_neu[2]

        # Compute u @ n+1
        l,j = np.indices(u_UPW.shape)
        for l in range(1,Nx):
            for j in range(1,Ny):
                u_UPW[l,j] = u_UPW[l,j]
                + dt*(mu[0]*((u_UPW[l+1,j]-2*u_UPW[l,j]+u_UPW[l-1,j])/dx**2
                            +(u_UPW[l,j+1]-2*u_UPW[l,j]+u_UPW[l,j-1])/dy**2)
                        - u*(u_UPW[l,j]-u_UPW[l-1,j])/dx 
                        - v*(u_UPW[l,j]-u_UPW[l,j-1])/dy)

Answer 1

我假设dt，mu[0]，dx，dy，u，v都是标量。您可以使用相同的公式来消除内部的两个循环，根据所需的索引，u_UPW数组向上或向下移动，或左或右移动1。

例如，代替u_UPW[l,j]，只使用u_UPW[1:Nx-1,1:Ny-1]，而不是u_UPW[l+1,j]使用u_UPW[2:Nx,1:Ny-1]。表达式中的所有元素都是大小为Nx-2 x Ny-2的数组，用于每次外部循环的迭代更新u_UPW[1:Nx-1,1:Ny-1]。

顺便说一句，您的标量代码可能有一个错误，因为您正在就地更新u_UPW。因此，以后的计算将使用以上或左边单元格的新值。您可以通过将新计算的值放在不同的数组中，然后在双内循环之后分配u_UPW来解决这个问题。当你一次做数组计算时，这个问题不会发生。