FFT와 Python을 사용하여 3D 라플라 시안을 효율적으로 계산합니다

StackOverflow https://stackoverflow.com//questions/22044126

  •  21-12-2019
  •  | 
  •  

문제

PDE (Schrödinger 방정식)를 해결하기 위해 Laplace 연산자를 3 차원으로 계산해야합니다. 내 현재의 솔루션은 (가장 멀리 필요한 코드의 일부)입니다. 

for n in range(Ntstep): # loop         

 for i in range(self.Nixyz[0]): # internal levels of wavefunction
    wf.psi[i,:,:,:]=self.expu * wf.psi[i,:,:,:] # potential      

    if n < Ntstep - 1: # compute laplacian in 3d
        wf.psi[i,:,:,:]=\
            sf.ifft(self.expkx*sf.fft(wf.psi[i,:,:,:],
                axis=0,**fft_args),axis=0,**fft_args)
        wf.psi[i,:,:,:]=\
            sf.ifft(self.expky*sf.fft(wf.psi[i,:,:,:],
                axis=1,**fft_args),axis=1,**fft_args)
        wf.psi[i,:,:,:]=\
            sf.ifft(self.expkz*sf.fft(wf.psi[i,:,:,:],
                axis=2,**fft_args),axis=2,**fft_args)
.

더 나은 성능을 얻으려면 다음과 같이 시도 / 수행 / 생각했습니다.

  • 3D FFT를 직접하지 마십시오. 라플라 시안은 분리 가능하므로 3 개의 1D FFT로 분할 될 수 있으며, n^3에서 3n로 복잡성을 줄여야합니다. (위의 코드에서 수행)

  • 나는 MKL에 대해 NUMPY와 scipy를 편집하고, 특히 멀티 스레드 계산을 가능하게하기를 희망하며, 특히 멀티 스레드 계산을 가능하게하기를 희망합니다. 일부 작업의 경우 다중 스레드가 사용됩니다 (매트릭스 벡터 곱셈)이지만 numpy.fft 또는 scipy.fftpack도 여러 코어를 사용하지 않습니다.

  • LibFFTW und pyfftw를 컴파일하고 NP / SP에 대한 드롭 인 교체로 사용했습니다. 나는 인텔 코어 I7-3770K, 즉 4 개의 코어와 8 개의 스레드가 있습니다. FFTW가있는 2 개 또는 4 개의 스레드를 사용할 때 NP / SP와 비교하여 성능이 두 배나 올랐습니다. 하나의 나사 또는 4 개 이상의 스레드가 느려지 며

그래서 내 주요 질문은 기본적으로 다음과 같습니다 :

  1. 는 코어 / 스레드 수를 갖는 성능 저울이있는 방식으로 FFT (W)가 평행하게 될 수 있습니까? 그렇다면 고려해야 할 것은 무엇입니까? 현재 2 ~ 4 개의 스레드가 저에게 달콤한 자리 인 것 같습니다. CPU에서 8 개의 스레드를 사용할 수 있지만, 더 느리지 만 더 느리게합니다.

  2. Python 코드를 병렬로 교체하려고합니까? 예를 들어, 세 가지 1D FFT를 세 개의 다른 코어에 넣으십시오. 물론 다른 스레드에서 동일한 변수에서 읽지 않고 쓸 수 없도록해야하므로 위의 코드에서 추가 "temp"변수가 필요합니다. :

    • 나사 1 : tempa= fft (psi ..., 축= 0)
    • 나사산 2 : tempb= fft (psi ..., 축= 1)
    • 스레드 3 : tempc= fft (psi ..., 축= 1)
    • 최종 단계 : PSI= 템퍼 + TEMPB + TEMPC

  3. axis=0의 FFT는 다른 축에 대해 두 번 (!)됩니다. 이 차이를 없애고 모든 FFT를 똑같이 빨리 만들 수 있습니까?

  4. (new) 은 FFT가 완전히 최선의 선택을 해결할 수 있거나 적어도 성능면에서 항상 사용자 로리에 의한 유한 차이점 접근법이 항상 더 나아지고 있습니까?

    < / li>

    Laplacian을 효율적으로 계산하는 것은 확실하게 연구 된 주제이므로 일부 링크 또는 논문 등의 힌트조차도 도움이 될 수 있습니다.

도움이 되었습니까?

해결책

나는 정말로 대답이 없지만, 나의 FFT 라플라 시안은 당신보다 더 단순 해 보인다 : 

def laplacian3d(field, KX, KY, KZ):
    return ifft(-KX**2*fft(field, axis = 0), axis = 0) + 
        ifft(-KY**2*fft(field, axis = 1), axis = 1) + 
        ifft(-KZ**2*fft(field, axis = 2), axis = 2)
.

KX, KY 및 KZ는 다음에서 만든 3D 어레이입니다. KX, KY, KZ = meshgrid(kx, ky, kz, indexing='ij') 및 feild는 3D 실제 공간 필드 (웨이브 함수) 및 kx = 2*pi*fftfreq(len(x), (x[1]-x[0])) (균등하게 간격 된 위치를 포함하는 실제 공간 1D 어레이가있는 X가있는 X가 포함됨)

실제로, 나는 cython에서 시행 된 유한 차이를 발견했다. 

cimport numpy as np
cimport cython
import numpy as np

#3D laplacian of a complex function
@cython.boundscheck(False) # turn of bounds-checking for entire function
def laplacianFD3dcomplex(np.ndarray[double complex, ndim=3] f, double complex dx, double complex dy, double complex dz):
    cdef unsigned int i, j, k, ni, nj, nk
    cdef double complex ifactor, jfactor, kfactor, ijkfactor
    ni = f.shape[0]
    nj = f.shape[1]
    nk = f.shape[2]
    cdef np.ndarray[double complex, ndim=3] lapf = np.zeros((ni,nj,nk)) +0.0J

    ifactor = 1/dx**2
    jfactor = 1/dy**2
    kfactor = 1/dz**2
    ijkfactor = 2.0*(ifactor + jfactor + kfactor)

    for i in xrange(1,ni-1):
        for j in xrange(1, nj-1):
            for k in xrange(1, nk-1):
                lapf[i, j, k] = (f[i, j, k-1] + f[i, j, k+1])*kfactor + (f[i, j-1, k] + f[i, j+1, k])*jfactor + (f[i-1, j, k] + f[i+1, j, k])*ifactor - f[i,j,k]*ijkfactor
    return lapf

#3D laplacian of a real function
@cython.boundscheck(False) # turn of bounds-checking for entire function
def laplacianFD3dreal(np.ndarray[double, ndim=3] f, double dx, double dy, double dz):
    cdef unsigned int i, j, k, ni, nj, nk
    cdef double ifactor, jfactor, kfactor, ijkfactor
    ni = f.shape[0]
    nj = f.shape[1]
    nk = f.shape[2]
    cdef np.ndarray[double, ndim=3] lapf = np.zeros((ni,nj,nk))

    ifactor = 1/dx**2
    jfactor = 1/dy**2
    kfactor = 1/dz**2
    ijkfactor = 2.0*(ifactor + jfactor + kfactor)

    for i in xrange(1,ni-1):
        for j in xrange(1, nj-1):
            for k in xrange(1, nk-1):
                lapf[i, j, k] = (f[i, j, k-1] + f[i, j, k+1])*kfactor + (f[i, j-1, k] + f[i, j+1, k])*jfactor + (f[i-1, j, k] + f[i+1, j, k])*ifactor - f[i,j,k]*ijkfactor
    return lapf
.

위의 코드는 "cython_finite_diff.pyx"라는 파일에 복사되고 다음과 같이 setup.py 파일을 사용하여 컴파일 할 수 있습니다. 

#To build the cython code in the .pyx file, type in the terminal:
#"python setup.py build_ext --inplace"
from distutils.core import setup
from distutils.extension import Extension
from Cython.Build import cythonize
import numpy

extensions = [
    Extension("cython_finite_diff", ["cython_finite_diff.pyx"],
                include_dirs = [numpy.get_include()]),
]

setup(
    name = "my_cython_fd",
    ext_modules = cythonize(extensions, annotate=True),
)
.

포맷에 대해 죄송합니다. 스택 오버플로에 게시하는 멍청이입니다.또한 유한 차이 라플라피아 인은 반사 경계 조건을 제공합니다.반대쪽 경계에있는 첫 번째 행을 포함하도록 루프를 설정하여 루프를 설정하여 주기적으로 만듭니다.

라이센스 : CC-BY-SA ~와 함께 속성
제휴하지 않습니다 StackOverflow
scroll top