223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

cleaning 

trying to get rid of strange re-ordering related to ASCII table data
layout
This commit is contained in:
Martin Diehl 2019-11-28 10:09:22 +01:00
parent 4c4ccfe72e
commit 80b50f460e
1 changed files with 19 additions and 28 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

47
python/damask/grid_filters.py
View File

@ -3,21 +3,18 @@ import numpy as np
def curl(size,field): def curl(size,field):
"""Calculate curl of a vector or tensor field in Fourier space.""" """Calculate curl of a vector or tensor field in Fourier space."""
shapeFFT = np.array(np.shape(field))[0:3] grid = np.array(np.shape(field)[0:3])
grid = np.array(np.shape(field)[2::-1])
N = grid.prod() # field size
n = np.array(np.shape(field)[3:]).prod() # data size n = np.array(np.shape(field)[3:]).prod() # data size
field_fourier = np.fft.rfftn(field,axes=(0,1,2),s=shapeFFT) field_fourier = np.fft.rfftn(field,axes=(0,1,2))
curl_fourier = np.empty(field_fourier.shape,'c16')
k_sk = np.where(np.arange(grid[2])>grid[2]//2,np.arange(grid[2])-grid[2],np.arange(grid[2]))/size[0] k_sk = np.where(np.arange(grid[0])>grid[0]//2,np.arange(grid[0])-grid[0],np.arange(grid[0]))/size[0]
if grid[2]%2 == 0: k_sk[grid[2]//2] = 0 # Nyquist freq=0 for even grid (Johnson, MIT, 2011) if grid[0]%2 == 0: k_sk[grid[0]//2] = 0 # Nyquist freq=0 for even grid (Johnson, MIT, 2011)
k_sj = np.where(np.arange(grid[1])>grid[1]//2,np.arange(grid[1])-grid[1],np.arange(grid[1]))/size[1] k_sj = np.where(np.arange(grid[1])>grid[1]//2,np.arange(grid[1])-grid[1],np.arange(grid[1]))/size[1]
if grid[1]%2 == 0: k_sj[grid[1]//2] = 0 # Nyquist freq=0 for even grid (Johnson, MIT, 2011) if grid[1]%2 == 0: k_sj[grid[1]//2] = 0 # Nyquist freq=0 for even grid (Johnson, MIT, 2011)
k_si = np.arange(grid[0]//2+1)/size[2] k_si = np.arange(grid[2]//2+1)/size[2]
kk, kj, ki = np.meshgrid(k_sk,k_sj,k_si,indexing = 'ij') 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') k_s = np.concatenate((ki[:,:,:,None],kj[:,:,:,None],kk[:,:,:,None]),axis = 3).astype('c16')
@ -26,29 +23,26 @@ def curl(size,field):
e[0, 1, 2] = e[1, 2, 0] = e[2, 0, 1] = +1.0 # Levi-Civita symbol e[0, 1, 2] = e[1, 2, 0] = e[2, 0, 1] = +1.0 # Levi-Civita symbol
e[0, 2, 1] = e[2, 1, 0] = e[1, 0, 2] = -1.0 e[0, 2, 1] = e[2, 1, 0] = e[1, 0, 2] = -1.0
curl = (np.einsum('slm,ijkl,ijkm, ->ijks', e,k_s,field_fourier)*2.0j*np.pi if n == 3 else # vector, 3 -> 3 curl = (np.einsum('slm,ijkl,ijkm ->ijks', e,k_s,field_fourier)*2.0j*np.pi if n == 3 else # vector, 3 -> 3
np.einsum('slm,ijkl,ijknm,->ijksn',e,k_s,field_fourier)*2.0j*np.pi) # tensor, 3x3 -> 3x3 np.einsum('slm,ijkl,ijknm->ijksn',e,k_s,field_fourier)*2.0j*np.pi) # tensor, 3x3 -> 3x3
return np.fft.irfftn(curl,axes=(0,1,2),s=shapeFFT).reshape([N,n]) return np.fft.irfftn(curl,axes=(0,1,2))
def divergence(size,field): def divergence(size,field):
"""Calculate divergence of a vector or tensor field in Fourier space.""" """Calculate divergence of a vector or tensor field in Fourier space."""
shapeFFT = np.array(np.shape(field))[0:3] grid = np.array(np.shape(field)[0:3])
grid = np.array(np.shape(field)[2::-1])
N = grid.prod() # field size
n = np.array(np.shape(field)[3:]).prod() # data size n = np.array(np.shape(field)[3:]).prod() # data size
field_fourier = np.fft.rfftn(field,axes=(0,1,2),s=shapeFFT) field_fourier = np.fft.rfftn(field,axes=(0,1,2),s=grid)
div_fourier = np.empty(field_fourier.shape[0:len(np.shape(field))-1],'c16')
k_sk = np.where(np.arange(grid[2])>grid[2]//2,np.arange(grid[2])-grid[2],np.arange(grid[2]))/size[0] k_sk = np.where(np.arange(grid[0])>grid[0]//2,np.arange(grid[0])-grid[0],np.arange(grid[0]))/size[0]
if grid[2]%2 == 0: k_sk[grid[2]//2] = 0 # Nyquist freq=0 for even grid (Johnson, MIT, 2011) if grid[0]%2 == 0: k_sk[grid[0]//2] = 0 # Nyquist freq=0 for even grid (Johnson, MIT, 2011)
k_sj = np.where(np.arange(grid[1])>grid[1]//2,np.arange(grid[1])-grid[1],np.arange(grid[1]))/size[1] k_sj = np.where(np.arange(grid[1])>grid[1]//2,np.arange(grid[1])-grid[1],np.arange(grid[1]))/size[1]
if grid[1]%2 == 0: k_sj[grid[1]//2] = 0 # Nyquist freq=0 for even grid (Johnson, MIT, 2011) if grid[1]%2 == 0: k_sj[grid[1]//2] = 0 # Nyquist freq=0 for even grid (Johnson, MIT, 2011)
k_si = np.arange(grid[0]//2+1)/size[2] k_si = np.arange(grid[2]//2+1)/size[2]
kk, kj, ki = np.meshgrid(k_sk,k_sj,k_si,indexing = 'ij') 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') k_s = np.concatenate((ki[:,:,:,None],kj[:,:,:,None],kk[:,:,:,None]),axis = 3).astype('c16')
@ -56,26 +50,23 @@ def divergence(size,field):
divergence = (np.einsum('ijkl,ijkl ->ijk', k_s,field_fourier)*2.0j*np.pi if n == 3 else # vector, 3 -> 1 divergence = (np.einsum('ijkl,ijkl ->ijk', k_s,field_fourier)*2.0j*np.pi if n == 3 else # vector, 3 -> 1
np.einsum('ijkm,ijklm->ijkl',k_s,field_fourier)*2.0j*np.pi) # tensor, 3x3 -> 3 np.einsum('ijkm,ijklm->ijkl',k_s,field_fourier)*2.0j*np.pi) # tensor, 3x3 -> 3
return np.fft.irfftn(div_fourier,axes=(0,1,2),s=shapeFFT).reshape([N,n//3]) return np.fft.irfftn(div_fourier,axes=(0,1,2),s=grid)
def gradient(size,field): def gradient(size,field):
"""Calculate gradient of a vector or scalar field in Fourier space.""" """Calculate gradient of a vector or scalar field in Fourier space."""
shapeFFT = np.array(np.shape(field))[0:3]
grid = np.array(np.shape(field)[2::-1]) grid = np.array(np.shape(field)[2::-1])
N = grid.prod() # field size
n = np.array(np.shape(field)[3:]).prod() # data size n = np.array(np.shape(field)[3:]).prod() # data size
field_fourier = np.fft.rfftn(field,axes=(0,1,2),s=shapeFFT) field_fourier = np.fft.rfftn(field,axes=(0,1,2),s=grid)
grad_fourier = np.empty(field_fourier.shape+(3,),'c16')
k_sk = np.where(np.arange(grid[2])>grid[2]//2,np.arange(grid[2])-grid[2],np.arange(grid[2]))/size[0] k_sk = np.where(np.arange(grid[0])>grid[0]//2,np.arange(grid[0])-grid[0],np.arange(grid[0]))/size[0]
if grid[2]%2 == 0: k_sk[grid[2]//2] = 0 # Nyquist freq=0 for even grid (Johnson, MIT, 2011) if grid[0]%2 == 0: k_sk[grid[0]//2] = 0 # Nyquist freq=0 for even grid (Johnson, MIT, 2011)
k_sj = np.where(np.arange(grid[1])>grid[1]//2,np.arange(grid[1])-grid[1],np.arange(grid[1]))/size[1] k_sj = np.where(np.arange(grid[1])>grid[1]//2,np.arange(grid[1])-grid[1],np.arange(grid[1]))/size[1]
if grid[1]%2 == 0: k_sj[grid[1]//2] = 0 # Nyquist freq=0 for even grid (Johnson, MIT, 2011) if grid[1]%2 == 0: k_sj[grid[1]//2] = 0 # Nyquist freq=0 for even grid (Johnson, MIT, 2011)
k_si = np.arange(grid[0]//2+1)/size[2] k_si = np.arange(grid[2]//2+1)/size[2]
kk, kj, ki = np.meshgrid(k_sk,k_sj,k_si,indexing = 'ij') 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') k_s = np.concatenate((ki[:,:,:,None],kj[:,:,:,None],kk[:,:,:,None]),axis = 3).astype('c16')
@ -83,7 +74,7 @@ def gradient(size,field):
gradient = (np.einsum('ijkl,ijkm->ijkm', field_fourier,k_s)*2.0j*np.pi if n == 1 else # scalar, 1 -> 3 gradient = (np.einsum('ijkl,ijkm->ijkm', field_fourier,k_s)*2.0j*np.pi if n == 1 else # scalar, 1 -> 3
np.einsum('ijkl,ijkm->ijklm',field_fourier,k_s)*2.0j*np.pi) # vector, 3 -> 3x3 np.einsum('ijkl,ijkm->ijklm',field_fourier,k_s)*2.0j*np.pi) # vector, 3 -> 3x3
return np.fft.irfftn(grad_fourier,axes=(0,1,2),s=shapeFFT).reshape([N,3*n]) return np.fft.irfftn(grad_fourier,axes=(0,1,2),s=grid)
#-------------------------------------------------------------------------------------------------- #--------------------------------------------------------------------------------------------------