using Einstein sum to replace 'for loop'
This commit is contained in:
Chuanlai Liu
parent
2d1421c09a
commit
d35c9dd431
|
|
@ -22,39 +22,26 @@ def curlFFT(geomdim,field):
|
|||
curl_fourier = np.empty(field_fourier.shape,'c16')
|
||||
|
||||
# differentiation in Fourier space
|
||||
k_s = np.zeros([3],'i')
|
||||
TWOPIIMG = 2.0j*math.pi
|
||||
for i in range(grid[2]):
|
||||
k_s[0] = i
|
||||
if grid[2]%2 == 0 and i == grid[2]//2: k_s[0] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011)
|
||||
elif i > grid[2]//2: k_s[0] -= grid[2]
|
||||
|
||||
for j in range(grid[1]):
|
||||
k_s[1] = j
|
||||
if grid[1]%2 == 0 and j == grid[1]//2: k_s[1] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011)
|
||||
elif j > grid[1]//2: k_s[1] -= grid[1]
|
||||
|
||||
for k in range(grid[0]//2+1):
|
||||
k_s[2] = k
|
||||
if grid[0]%2 == 0 and k == grid[0]//2: k_s[2] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011)
|
||||
|
||||
xi = (k_s/geomdim)[2::-1].astype('c16') # reversing the field input order
|
||||
|
||||
if dataType == 'tensor':
|
||||
for l in range(3):
|
||||
curl_fourier[i,j,k,0,l] = ( field_fourier[i,j,k,l,2]*xi[1]\
|
||||
-field_fourier[i,j,k,l,1]*xi[2]) *TWOPIIMG
|
||||
curl_fourier[i,j,k,1,l] = (-field_fourier[i,j,k,l,2]*xi[0]\
|
||||
+field_fourier[i,j,k,l,0]*xi[2]) *TWOPIIMG
|
||||
curl_fourier[i,j,k,2,l] = ( field_fourier[i,j,k,l,1]*xi[0]\
|
||||
-field_fourier[i,j,k,l,0]*xi[1]) *TWOPIIMG
|
||||
elif dataType == 'vector':
|
||||
curl_fourier[i,j,k,0] = ( field_fourier[i,j,k,2]*xi[1]\
|
||||
-field_fourier[i,j,k,1]*xi[2]) *TWOPIIMG
|
||||
curl_fourier[i,j,k,1] = (-field_fourier[i,j,k,2]*xi[0]\
|
||||
+field_fourier[i,j,k,0]*xi[2]) *TWOPIIMG
|
||||
curl_fourier[i,j,k,2] = ( field_fourier[i,j,k,1]*xi[0]\
|
||||
-field_fourier[i,j,k,0]*xi[1]) *TWOPIIMG
|
||||
k_sk = np.where(np.arange(grid[2])>grid[2]//2,np.arange(grid[2])-grid[2],np.arange(grid[2]))/geomdim[0]
|
||||
if grid[2]%2 == 0: k_sk[grid[2]//2] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011)
|
||||
|
||||
k_sj = np.where(np.arange(grid[1])>grid[1]//2,np.arange(grid[1])-grid[1],np.arange(grid[1]))/geomdim[1]
|
||||
if grid[1]%2 == 0: k_sj[grid[1]//2] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011)
|
||||
|
||||
k_si = np.arange(grid[0]//2+1)/geomdim[2]
|
||||
|
||||
kk, kj, ki = np.meshgrid(k_sk,k_sj,k_si,indexing = 'ij')
|
||||
k_s = np.concatenate((ki[:,:,:,None],kj[:,:,:,None],kk[:,:,:,None]),axis = 3).astype('c16')
|
||||
|
||||
e = np.zeros((3, 3, 3))
|
||||
e[0, 1, 2] = e[1, 2, 0] = e[2, 0, 1] = 1.0 # Levi-Civita symbols
|
||||
e[0, 2, 1] = e[2, 1, 0] = e[1, 0, 2] = -1.0
|
||||
|
||||
if dataType == 'tensor':
|
||||
curl_fourier = np.einsum('slm,ijkl,ijknm->ijksn',e,k_s,field_fourier)*TWOPIIMG
|
||||
elif dataType == 'vector':
|
||||
curl_fourier = np.einsum('slm,ijkl,ijkm->ijks',e,k_s,field_fourier)*TWOPIIMG
|
||||
|
||||
return np.fft.irfftn(curl_fourier,axes=(0,1,2),s=shapeFFT).reshape([N,n])
|
||||
|
||||
|
|
|
|||
|
|
@ -19,29 +19,22 @@ def divFFT(geomdim,field):
|
|||
div_fourier = np.empty(field_fourier.shape[0:len(np.shape(field))-1],'c16') # size depents on whether tensor or vector
|
||||
|
||||
# differentiation in Fourier space
|
||||
k_s = np.zeros([3],'i')
|
||||
TWOPIIMG = 2.0j*math.pi
|
||||
for i in range(grid[2]):
|
||||
k_s[0] = i
|
||||
if grid[2]%2 == 0 and i == grid[2]//2: k_s[0] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011)
|
||||
elif i > grid[2]//2: k_s[0] -= grid[2]
|
||||
|
||||
for j in range(grid[1]):
|
||||
k_s[1] = j
|
||||
if grid[1]%2 == 0 and j == grid[1]//2: k_s[1] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011)
|
||||
elif j > grid[1]//2: k_s[1] -= grid[1]
|
||||
|
||||
for k in range(grid[0]//2+1):
|
||||
k_s[2] = k
|
||||
if grid[0]%2 == 0 and k == grid[0]//2: k_s[2] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011)
|
||||
|
||||
xi = (k_s/geomdim)[2::-1].astype('c16') # reversing the field input order
|
||||
if n == 9: # tensor, 3x3 -> 3
|
||||
for l in range(3):
|
||||
div_fourier[i,j,k,l] = sum(field_fourier[i,j,k,l,0:3]*xi) *TWOPIIMG
|
||||
elif n == 3: # vector, 3 -> 1
|
||||
div_fourier[i,j,k] = sum(field_fourier[i,j,k,0:3]*xi) *TWOPIIMG
|
||||
k_sk = np.where(np.arange(grid[2])>grid[2]//2,np.arange(grid[2])-grid[2],np.arange(grid[2]))/geomdim[0]
|
||||
if grid[2]%2 == 0: k_sk[grid[2]//2] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011)
|
||||
|
||||
k_sj = np.where(np.arange(grid[1])>grid[1]//2,np.arange(grid[1])-grid[1],np.arange(grid[1]))/geomdim[1]
|
||||
if grid[1]%2 == 0: k_sj[grid[1]//2] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011)
|
||||
|
||||
k_si = np.arange(grid[0]//2+1)/geomdim[2]
|
||||
|
||||
kk, kj, ki = np.meshgrid(k_sk,k_sj,k_si,indexing = 'ij')
|
||||
k_s = np.concatenate((ki[:,:,:,None],kj[:,:,:,None],kk[:,:,:,None]),axis = 3).astype('c16')
|
||||
if n == 9: # tensor, 3x3 -> 3
|
||||
div_fourier = np.einsum('ijklm,ijkm->ijkl',field_fourier,k_s)*TWOPIIMG
|
||||
elif n == 3: # vector, 3 -> 1
|
||||
div_fourier = np.einsum('ijkl,ijkl->ijk',field_fourier,k_s)*TWOPIIMG
|
||||
|
||||
return np.fft.irfftn(div_fourier,axes=(0,1,2),s=shapeFFT).reshape([N,n/3])
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -22,29 +22,21 @@ def gradFFT(geomdim,field):
|
|||
grad_fourier = np.empty(field_fourier.shape+(3,),'c16')
|
||||
|
||||
# differentiation in Fourier space
|
||||
k_s = np.zeros([3],'i')
|
||||
TWOPIIMG = 2.0j*math.pi
|
||||
for i in range(grid[2]):
|
||||
k_s[0] = i
|
||||
if grid[2]%2 == 0 and i == grid[2]//2: k_s[0] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011)
|
||||
elif i > grid[2]//2: k_s[0] -= grid[2]
|
||||
|
||||
for j in range(grid[1]):
|
||||
k_s[1] = j
|
||||
if grid[1]%2 == 0 and j == grid[1]//2: k_s[1] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011)
|
||||
elif j > grid[1]//2: k_s[1] -= grid[1]
|
||||
|
||||
for k in range(grid[0]//2+1):
|
||||
k_s[2] = k
|
||||
if grid[0]%2 == 0 and k == grid[0]//2: k_s[2] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011)
|
||||
|
||||
xi = (k_s/geomdim)[2::-1].astype('c16') # reversing the field order
|
||||
|
||||
grad_fourier[i,j,k,0,:] = field_fourier[i,j,k,0]*xi *TWOPIIMG # vector field from scalar data
|
||||
|
||||
if dataType == 'vector':
|
||||
grad_fourier[i,j,k,1,:] = field_fourier[i,j,k,1]*xi *TWOPIIMG # tensor field from vector data
|
||||
grad_fourier[i,j,k,2,:] = field_fourier[i,j,k,2]*xi *TWOPIIMG
|
||||
k_sk = np.where(np.arange(grid[2])>grid[2]//2,np.arange(grid[2])-grid[2],np.arange(grid[2]))/geomdim[0]
|
||||
if grid[2]%2 == 0: k_sk[grid[2]//2] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011)
|
||||
|
||||
k_sj = np.where(np.arange(grid[1])>grid[1]//2,np.arange(grid[1])-grid[1],np.arange(grid[1]))/geomdim[1]
|
||||
if grid[1]%2 == 0: k_sj[grid[1]//2] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011)
|
||||
|
||||
k_si = np.arange(grid[0]//2+1)/geomdim[2]
|
||||
|
||||
kk, kj, ki = np.meshgrid(k_sk,k_sj,k_si,indexing = 'ij')
|
||||
k_s = np.concatenate((ki[:,:,:,None],kj[:,:,:,None],kk[:,:,:,None]),axis = 3).astype('c16')
|
||||
if n == 3: # vector, 3 -> 3x3
|
||||
grad_fourier = np.einsum('ijkl,ijkm->ijklm',field_fourier,k_s)*TWOPIIMG
|
||||
elif n == 1: # scalar, 1 -> 3
|
||||
grad_fourier = np.einsum('ijkl,ijkl->ijkl',field_fourier,k_s)*TWOPIIMG
|
||||
|
||||
return np.fft.irfftn(grad_fourier,axes=(0,1,2),s=shapeFFT).reshape([N,3*n])
|
||||
|
||||
|
|
|
|||
Loading…
Reference in New Issue