Overhauled to put all component updates into single functions, for less function calls.

这个提交包含在:
Craig Warren
2016-02-25 18:02:11 +00:00
父节点 b40b908ff7
当前提交 3dbd912f17

查看文件

@@ -22,10 +22,10 @@ from cython.parallel import prange
from gprMax.constants cimport floattype_t, complextype_t from gprMax.constants cimport floattype_t, complextype_t
######################################### #####################################
# Electric field updates - Ex component # # Electric field updates - standard #
######################################### #####################################
cpdef void update_ex(int nx, int ny, int nz, int nthreads, floattype_t[:, ::1] updatecoeffsE, np.uint32_t[:, :, :, ::1] ID, floattype_t[:, :, ::1] Ex, floattype_t[:, :, ::1] Hy, floattype_t[:, :, ::1] Hz): cpdef void update_electric(int nx, int ny, int nz, int nthreads, floattype_t[:, ::1] updatecoeffsE, np.uint32_t[:, :, :, ::1] ID, floattype_t[:, :, ::1] Ex, floattype_t[:, :, ::1] Ey, floattype_t[:, :, ::1] Ez, floattype_t[:, :, ::1] Hx, floattype_t[:, :, ::1] Hy, floattype_t[:, :, ::1] Hz):
"""This function updates the Ex field components. """This function updates the Ex field components.
Args: Args:
@@ -37,6 +37,7 @@ cpdef void update_ex(int nx, int ny, int nz, int nthreads, floattype_t[:, ::1] u
cdef Py_ssize_t i, j, k cdef Py_ssize_t i, j, k
cdef int material cdef int material
# Ex component
if ny == 1 or nz == 1: if ny == 1 or nz == 1:
pass pass
else: else:
@@ -46,8 +47,31 @@ cpdef void update_ex(int nx, int ny, int nz, int nthreads, floattype_t[:, ::1] u
material = ID[0, i, j, k] material = ID[0, i, j, k]
Ex[i, j, k] = updatecoeffsE[material, 0] * Ex[i, j, k] + updatecoeffsE[material, 2] * (Hz[i, j, k] - Hz[i, j - 1, k]) - updatecoeffsE[material, 3] * (Hy[i, j, k] - Hy[i, j, k - 1]) Ex[i, j, k] = updatecoeffsE[material, 0] * Ex[i, j, k] + updatecoeffsE[material, 2] * (Hz[i, j, k] - Hz[i, j - 1, k]) - updatecoeffsE[material, 3] * (Hy[i, j, k] - Hy[i, j, k - 1])
# Ey component
if nx == 1 or nz == 1:
pass
else:
for i in prange(1, nx, nogil=True, schedule='static', num_threads=nthreads):
for j in range(0, ny):
for k in range(1, nz):
material = ID[1, i, j, k]
Ey[i, j, k] = updatecoeffsE[material, 0] * Ey[i, j, k] + updatecoeffsE[material, 3] * (Hx[i, j, k] - Hx[i, j, k - 1]) - updatecoeffsE[material, 1] * (Hz[i, j, k] - Hz[i - 1, j, k])
cpdef void update_ex_dispersive_multipole_A(int nx, int ny, int nz, int nthreads, int maxpoles, floattype_t[:, ::1] updatecoeffsE, complextype_t[:, ::1] updatecoeffsdispersive, np.uint32_t[:, :, :, ::1] ID, complextype_t[:, :, :, ::1] Tx, floattype_t[:, :, ::1] Ex, floattype_t[:, :, ::1] Hy, floattype_t[:, :, ::1] Hz): # Ez component
if nx == 1 or ny == 1:
pass
else:
for i in prange(1, nx, nogil=True, schedule='static', num_threads=nthreads):
for j in range(1, ny):
for k in range(0, nz):
material = ID[2, i, j, k]
Ez[i, j, k] = updatecoeffsE[material, 0] * Ez[i, j, k] + updatecoeffsE[material, 1] * (Hy[i, j, k] - Hy[i - 1, j, k]) - updatecoeffsE[material, 2] * (Hx[i, j, k] - Hx[i, j - 1, k])
#################################################
# Electric field updates - dispersive materials #
#################################################
cpdef void update_electric_dispersive_multipole_A(int nx, int ny, int nz, int nthreads, int maxpoles, floattype_t[:, ::1] updatecoeffsE, complextype_t[:, ::1] updatecoeffsdispersive, np.uint32_t[:, :, :, ::1] ID, complextype_t[:, :, :, ::1] Tx, complextype_t[:, :, :, ::1] Ty, complextype_t[:, :, :, ::1] Tz, floattype_t[:, :, ::1] Ex, floattype_t[:, :, ::1] Ey, floattype_t[:, :, ::1] Ez, floattype_t[:, :, ::1] Hx, floattype_t[:, :, ::1] Hy, floattype_t[:, :, ::1] Hz):
"""This function updates the Ex field components when dispersive materials (with multiple poles) are present. """This function updates the Ex field components when dispersive materials (with multiple poles) are present.
Args: Args:
@@ -57,10 +81,11 @@ cpdef void update_ex_dispersive_multipole_A(int nx, int ny, int nz, int nthreads
updatecoeffs, T, ID, E, H (memoryviews): Access to update coeffients, temporary, ID and field component arrays updatecoeffs, T, ID, E, H (memoryviews): Access to update coeffients, temporary, ID and field component arrays
""" """
cdef Py_ssize_t i, j, k, p cdef Py_ssize_t i, j, k, pole
cdef int material cdef int material
cdef float phi = 0.0 cdef float phi = 0
# Ex component
if ny == 1 or nz == 1: if ny == 1 or nz == 1:
pass pass
else: else:
@@ -68,13 +93,43 @@ cpdef void update_ex_dispersive_multipole_A(int nx, int ny, int nz, int nthreads
for j in range(1, ny): for j in range(1, ny):
for k in range(1, nz): for k in range(1, nz):
material = ID[0, i, j, k] material = ID[0, i, j, k]
phi = 0.0 phi = 0
for p in range(0, maxpoles): for pole in range(maxpoles):
phi = phi + updatecoeffsdispersive[material, p * 3].real * Tx[p, i, j, k].real phi = phi + updatecoeffsdispersive[material, pole * 3].real * Tx[pole, i, j, k].real
Tx[p, i, j, k] = updatecoeffsdispersive[material, 1 + (p * 3)] * Tx[p, i, j, k] + updatecoeffsdispersive[material, 2 + (p * 3)] * Ex[i, j, k] Tx[pole, i, j, k] = updatecoeffsdispersive[material, 1 + (pole * 3)] * Tx[pole, i, j, k] + updatecoeffsdispersive[material, 2 + (pole * 3)] * Ex[i, j, k]
Ex[i, j, k] = updatecoeffsE[material, 0] * Ex[i, j, k] + updatecoeffsE[material, 2] * (Hz[i, j, k] - Hz[i, j - 1, k]) - updatecoeffsE[material, 3] * (Hy[i, j, k] - Hy[i, j, k - 1]) - updatecoeffsE[material, 4] * phi Ex[i, j, k] = updatecoeffsE[material, 0] * Ex[i, j, k] + updatecoeffsE[material, 2] * (Hz[i, j, k] - Hz[i, j - 1, k]) - updatecoeffsE[material, 3] * (Hy[i, j, k] - Hy[i, j, k - 1]) - updatecoeffsE[material, 4] * phi
cpdef void update_ex_dispersive_multipole_B(int nx, int ny, int nz, int nthreads, int maxpoles, complextype_t[:, ::1] updatecoeffsdispersive, np.uint32_t[:, :, :, ::1] ID, complextype_t[:, :, :, ::1] Tx, floattype_t[:, :, ::1] Ex): # Ey component
if nx == 1 or nz == 1:
pass
else:
for i in prange(1, nx, nogil=True, schedule='static', num_threads=nthreads):
for j in range(0, ny):
for k in range(1, nz):
material = ID[1, i, j, k]
phi = 0
for pole in range(maxpoles):
phi = phi + updatecoeffsdispersive[material, pole * 3].real * Ty[pole, i, j, k].real
Ty[pole, i, j, k] = updatecoeffsdispersive[material, 1 + (pole * 3)] * Ty[pole, i, j, k] + updatecoeffsdispersive[material, 2 + (pole * 3)] * Ey[i, j, k]
Ey[i, j, k] = updatecoeffsE[material, 0] * Ey[i, j, k] + updatecoeffsE[material, 3] * (Hx[i, j, k] - Hx[i, j, k - 1]) - updatecoeffsE[material, 1] * (Hz[i, j, k] - Hz[i - 1, j, k]) - updatecoeffsE[material, 4] * phi
# Ez component
if nx == 1 or ny == 1:
pass
else:
for i in prange(1, nx, nogil=True, schedule='static', num_threads=nthreads):
for j in range(1, ny):
for k in range(0, nz):
material = ID[2, i, j, k]
phi = 0
for pole in range(maxpoles):
phi = phi + updatecoeffsdispersive[material, pole * 3].real * Tz[pole, i, j, k].real
Tz[pole, i, j, k] = updatecoeffsdispersive[material, 1 + (pole * 3)] * Tz[pole, i, j, k] + updatecoeffsdispersive[material, 2 + (pole * 3)] * Ez[i, j, k]
Ez[i, j, k] = updatecoeffsE[material, 0] * Ez[i, j, k] + updatecoeffsE[material, 1] * (Hy[i, j, k] - Hy[i - 1, j, k]) - updatecoeffsE[material, 2] * (Hx[i, j, k] - Hx[i, j - 1, k]) - updatecoeffsE[material, 4] * phi
cpdef void update_electric_dispersive_multipole_B(int nx, int ny, int nz, int nthreads, int maxpoles, complextype_t[:, ::1] updatecoeffsdispersive, np.uint32_t[:, :, :, ::1] ID, complextype_t[:, :, :, ::1] Tx, complextype_t[:, :, :, ::1] Ty, complextype_t[:, :, :, ::1] Tz, floattype_t[:, :, ::1] Ex, floattype_t[:, :, ::1] Ey, floattype_t[:, :, ::1] Ez):
"""This function updates the Ex field components when dispersive materials (with multiple poles) are present. """This function updates the Ex field components when dispersive materials (with multiple poles) are present.
Args: Args:
@@ -84,9 +139,10 @@ cpdef void update_ex_dispersive_multipole_B(int nx, int ny, int nz, int nthreads
updatecoeffs, T, ID, E (memoryviews): Access to update coeffients, temporary, ID and field component arrays updatecoeffs, T, ID, E (memoryviews): Access to update coeffients, temporary, ID and field component arrays
""" """
cdef Py_ssize_t i, j, k, p cdef Py_ssize_t i, j, k, pole
cdef int material cdef int material
# Ex component
if ny == 1 or nz == 1: if ny == 1 or nz == 1:
pass pass
else: else:
@@ -94,11 +150,33 @@ cpdef void update_ex_dispersive_multipole_B(int nx, int ny, int nz, int nthreads
for j in range(1, ny): for j in range(1, ny):
for k in range(1, nz): for k in range(1, nz):
material = ID[0, i, j, k] material = ID[0, i, j, k]
for p in range(0, maxpoles): for pole in range(maxpoles):
Tx[p, i, j, k] = Tx[p, i, j, k] - updatecoeffsdispersive[material, 2 + (p * 3)] * Ex[i, j, k] Tx[pole, i, j, k] = Tx[pole, i, j, k] - updatecoeffsdispersive[material, 2 + (pole * 3)] * Ex[i, j, k]
# Ey component
if nx == 1 or nz == 1:
pass
else:
for i in prange(1, nx, nogil=True, schedule='static', num_threads=nthreads):
for j in range(0, ny):
for k in range(1, nz):
material = ID[1, i, j, k]
for pole in range(maxpoles):
Ty[pole, i, j, k] = Ty[pole, i, j, k] - updatecoeffsdispersive[material, 2 + (pole * 3)] * Ey[i, j, k]
# Ez component
if nx == 1 or ny == 1:
pass
else:
for i in prange(1, nx, nogil=True, schedule='static', num_threads=nthreads):
for j in range(1, ny):
for k in range(0, nz):
material = ID[2, i, j, k]
for pole in range(maxpoles):
Tz[pole, i, j, k] = Tz[pole, i, j, k] - updatecoeffsdispersive[material, 2 + (pole * 3)] * Ez[i, j, k]
cpdef void update_ex_dispersive_1pole_A(int nx, int ny, int nz, int nthreads, floattype_t[:, ::1] updatecoeffsE, complextype_t[:, ::1] updatecoeffsdispersive, np.uint32_t[:, :, :, ::1] ID, complextype_t[:, :, :, ::1] Tx, floattype_t[:, :, ::1] Ex, floattype_t[:, :, ::1] Hy, floattype_t[:, :, ::1] Hz): cpdef void update_electric_dispersive_1pole_A(int nx, int ny, int nz, int nthreads, floattype_t[:, ::1] updatecoeffsE, complextype_t[:, ::1] updatecoeffsdispersive, np.uint32_t[:, :, :, ::1] ID, complextype_t[:, :, :, ::1] Tx, complextype_t[:, :, :, ::1] Ty, complextype_t[:, :, :, ::1] Tz, floattype_t[:, :, ::1] Ex, floattype_t[:, :, ::1] Ey, floattype_t[:, :, ::1] Ez, floattype_t[:, :, ::1] Hx, floattype_t[:, :, ::1] Hy, floattype_t[:, :, ::1] Hz):
"""This function updates the Ex field components when dispersive materials (with 1 pole) are present. """This function updates the Ex field components when dispersive materials (with 1 pole) are present.
Args: Args:
@@ -109,8 +187,9 @@ cpdef void update_ex_dispersive_1pole_A(int nx, int ny, int nz, int nthreads, fl
cdef Py_ssize_t i, j, k cdef Py_ssize_t i, j, k
cdef int material cdef int material
cdef float phi = 0.0 cdef float phi = 0
# Ex component
if ny == 1 or nz == 1: if ny == 1 or nz == 1:
pass pass
else: else:
@@ -122,119 +201,7 @@ cpdef void update_ex_dispersive_1pole_A(int nx, int ny, int nz, int nthreads, fl
Tx[0, i, j, k] = updatecoeffsdispersive[material, 1] * Tx[0, i, j, k] + updatecoeffsdispersive[material, 2] * Ex[i, j, k] Tx[0, i, j, k] = updatecoeffsdispersive[material, 1] * Tx[0, i, j, k] + updatecoeffsdispersive[material, 2] * Ex[i, j, k]
Ex[i, j, k] = updatecoeffsE[material, 0] * Ex[i, j, k] + updatecoeffsE[material, 2] * (Hz[i, j, k] - Hz[i, j - 1, k]) - updatecoeffsE[material, 3] * (Hy[i, j, k] - Hy[i, j, k - 1]) - updatecoeffsE[material, 4] * phi Ex[i, j, k] = updatecoeffsE[material, 0] * Ex[i, j, k] + updatecoeffsE[material, 2] * (Hz[i, j, k] - Hz[i, j - 1, k]) - updatecoeffsE[material, 3] * (Hy[i, j, k] - Hy[i, j, k - 1]) - updatecoeffsE[material, 4] * phi
# Ey component
cpdef void update_ex_dispersive_1pole_B(int nx, int ny, int nz, int nthreads, complextype_t[:, ::1] updatecoeffsdispersive, np.uint32_t[:, :, :, ::1] ID, complextype_t[:, :, :, ::1] Tx, floattype_t[:, :, ::1] Ex):
"""This function updates the Ex field components when dispersive materials (with 1 pole) are present.
Args:
nx, ny, nz (int): Grid size in cells
nthreads (int): Number of threads to use
updatecoeffs, T, ID, E (memoryviews): Access to update coeffients, temporary, ID and field component arrays
"""
cdef Py_ssize_t i, j, k
cdef int material
if ny == 1 or nz == 1:
pass
else:
for i in prange(0, nx, nogil=True, schedule='static', num_threads=nthreads):
for j in range(1, ny):
for k in range(1, nz):
material = ID[0, i, j, k]
Tx[0, i, j, k] = Tx[0, i, j, k] - updatecoeffsdispersive[material, 2] * Ex[i, j, k]
#########################################
# Electric field updates - Ey component #
#########################################
cpdef void update_ey(int nx, int ny, int nz, int nthreads, floattype_t[:, ::1] updatecoeffsE, np.uint32_t[:, :, :, ::1] ID, floattype_t[:, :, ::1] Ey, floattype_t[:, :, ::1] Hx, floattype_t[:, :, ::1] Hz):
"""This function updates the Ey field components.
Args:
nx, ny, nz (int): Grid size in cells
nthreads (int): Number of threads to use
updatecoeffs, ID, E, H (memoryviews): Access to update coeffients, ID and field component arrays
"""
cdef Py_ssize_t i, j, k
cdef int material
if nx == 1 or nz == 1:
pass
else:
for i in prange(1, nx, nogil=True, schedule='static', num_threads=nthreads):
for j in range(0, ny):
for k in range(1, nz):
material = ID[1, i, j, k]
Ey[i, j, k] = updatecoeffsE[material, 0] * Ey[i, j, k] + updatecoeffsE[material, 3] * (Hx[i, j, k] - Hx[i, j, k - 1]) - updatecoeffsE[material, 1] * (Hz[i, j, k] - Hz[i - 1, j, k])
cpdef void update_ey_dispersive_multipole_A(int nx, int ny, int nz, int nthreads, int maxpoles, floattype_t[:, ::1] updatecoeffsE, complextype_t[:, ::1] updatecoeffsdispersive, np.uint32_t[:, :, :, ::1] ID, complextype_t[:, :, :, ::1] Ty, floattype_t[:, :, ::1] Ey, floattype_t[:, :, ::1] Hx, floattype_t[:, :, ::1] Hz):
"""This function updates the Ey field components when dispersive materials (with multiple poles) are present.
Args:
nx, ny, nz (int): Grid size in cells
maxpoles (int): Maximum number of poles
nthreads (int): Number of threads to use
updatecoeffs, T, ID, E, H (memoryviews): Access to update coeffients, temporary, ID and field component arrays
"""
cdef Py_ssize_t i, j, k, p
cdef int material
cdef float phi = 0.0
if nx == 1 or nz == 1:
pass
else:
for i in prange(1, nx, nogil=True, schedule='static', num_threads=nthreads):
for j in range(0, ny):
for k in range(1, nz):
material = ID[1, i, j, k]
phi = 0.0
for p in range(0, maxpoles):
phi = phi + updatecoeffsdispersive[material, p * 3].real * Ty[p, i, j, k].real
Ty[p, i, j, k] = updatecoeffsdispersive[material, 1 + (p * 3)] * Ty[p, i, j, k] + updatecoeffsdispersive[material, 2 + (p * 3)] * Ey[i, j, k]
Ey[i, j, k] = updatecoeffsE[material, 0] * Ey[i, j, k] + updatecoeffsE[material, 3] * (Hx[i, j, k] - Hx[i, j, k - 1]) - updatecoeffsE[material, 1] * (Hz[i, j, k] - Hz[i - 1, j, k]) - updatecoeffsE[material, 4] * phi
cpdef void update_ey_dispersive_multipole_B(int nx, int ny, int nz, int nthreads, int maxpoles, complextype_t[:, ::1] updatecoeffsdispersive, np.uint32_t[:, :, :, ::1] ID, complextype_t[:, :, :, ::1] Ty, floattype_t[:, :, ::1] Ey):
"""This function updates the Ey field components when dispersive materials (with multiple poles) are present.
Args:
nx, ny, nz (int): Grid size in cells
maxpoles (int): Maximum number of poles
nthreads (int): Number of threads to use
updatecoeffs, T, ID, E, H (memoryviews): Access to update coeffients, temporary, ID and field component arrays
"""
cdef Py_ssize_t i, j, k, p
cdef int material
if nx == 1 or nz == 1:
pass
else:
for i in prange(1, nx, nogil=True, schedule='static', num_threads=nthreads):
for j in range(0, ny):
for k in range(1, nz):
material = ID[1, i, j, k]
for p in range(0, maxpoles):
Ty[p, i, j, k] = Ty[p, i, j, k] - updatecoeffsdispersive[material, 2 + (p * 3)] * Ey[i, j, k]
cpdef void update_ey_dispersive_1pole_A(int nx, int ny, int nz, int nthreads, floattype_t[:, ::1] updatecoeffsE, complextype_t[:, ::1] updatecoeffsdispersive, np.uint32_t[:, :, :, ::1] ID, complextype_t[:, :, :, ::1] Ty, floattype_t[:, :, ::1] Ey, floattype_t[:, :, ::1] Hx, floattype_t[:, :, ::1] Hz):
"""This function updates the Ey field components when dispersive materials (with 1 pole) are present.
Args:
nx, ny, nz (int): Grid size in cells
nthreads (int): Number of threads to use
updatecoeffs, T, ID, E, H (memoryviews): Access to update coeffients, temporary, ID and field component arrays
"""
cdef Py_ssize_t i, j, k
cdef int material
cdef float phi = 0.0
if nx == 1 or nz == 1: if nx == 1 or nz == 1:
pass pass
else: else:
@@ -246,119 +213,7 @@ cpdef void update_ey_dispersive_1pole_A(int nx, int ny, int nz, int nthreads, fl
Ty[0, i, j, k] = updatecoeffsdispersive[material, 1] * Ty[0, i, j, k] + updatecoeffsdispersive[material, 2] * Ey[i, j, k] Ty[0, i, j, k] = updatecoeffsdispersive[material, 1] * Ty[0, i, j, k] + updatecoeffsdispersive[material, 2] * Ey[i, j, k]
Ey[i, j, k] = updatecoeffsE[material, 0] * Ey[i, j, k] + updatecoeffsE[material, 3] * (Hx[i, j, k] - Hx[i, j, k - 1]) - updatecoeffsE[material, 1] * (Hz[i, j, k] - Hz[i - 1, j, k]) - updatecoeffsE[material, 4] * phi Ey[i, j, k] = updatecoeffsE[material, 0] * Ey[i, j, k] + updatecoeffsE[material, 3] * (Hx[i, j, k] - Hx[i, j, k - 1]) - updatecoeffsE[material, 1] * (Hz[i, j, k] - Hz[i - 1, j, k]) - updatecoeffsE[material, 4] * phi
# Ez component
cpdef void update_ey_dispersive_1pole_B(int nx, int ny, int nz, int nthreads, complextype_t[:, ::1] updatecoeffsdispersive, np.uint32_t[:, :, :, ::1] ID, complextype_t[:, :, :, ::1] Ty, floattype_t[:, :, ::1] Ey):
"""This function updates the Ey field components when dispersive materials (with 1 pole) are present.
Args:
nx, ny, nz (int): Grid size in cells
nthreads (int): Number of threads to use
updatecoeffs, T, ID, E (memoryviews): Access to update coeffients, temporary, ID and field component arrays
"""
cdef Py_ssize_t i, j, k
cdef int material
if nx == 1 or nz == 1:
pass
else:
for i in prange(1, nx, nogil=True, schedule='static', num_threads=nthreads):
for j in range(0, ny):
for k in range(1, nz):
material = ID[1, i, j, k]
Ty[0, i, j, k] = Ty[0, i, j, k] - updatecoeffsdispersive[material, 2] * Ey[i, j, k]
#########################################
# Electric field updates - Ez component #
#########################################
cpdef void update_ez(int nx, int ny, int nz, int nthreads, floattype_t[:, ::1] updatecoeffsE, np.uint32_t[:, :, :, ::1] ID, floattype_t[:, :, ::1] Ez, floattype_t[:, :, ::1] Hx, floattype_t[:, :, ::1] Hy):
"""This function updates the Ez field components.
Args:
nx, ny, nz (int): Grid size in cells
nthreads (int): Number of threads to use
updatecoeffs, ID, E, H (memoryviews): Access to update coeffients, ID and field component arrays
"""
cdef Py_ssize_t i, j, k
cdef int material
if nx == 1 or ny == 1:
pass
else:
for i in prange(1, nx, nogil=True, schedule='static', num_threads=nthreads):
for j in range(1, ny):
for k in range(0, nz):
material = ID[2, i, j, k]
Ez[i, j, k] = updatecoeffsE[material, 0] * Ez[i, j, k] + updatecoeffsE[material, 1] * (Hy[i, j, k] - Hy[i - 1, j, k]) - updatecoeffsE[material, 2] * (Hx[i, j, k] - Hx[i, j - 1, k])
cpdef void update_ez_dispersive_multipole_A(int nx, int ny, int nz, int nthreads, int maxpoles, floattype_t[:, ::1] updatecoeffsE, complextype_t[:, ::1] updatecoeffsdispersive, np.uint32_t[:, :, :, ::1] ID, complextype_t[:, :, :, ::1] Tz, floattype_t[:, :, ::1] Ez, floattype_t[:, :, ::1] Hx, floattype_t[:, :, ::1] Hy):
"""This function updates the Ez field components when dispersive materials (with multiple poles) are present.
Args:
nx, ny, nz (int): Grid size in cells
maxpoles (int): Maximum number of poles
nthreads (int): Number of threads to use
updatecoeffs, T, ID, E, H (memoryviews): Access to update coeffients, temporary, ID and field component arrays
"""
cdef Py_ssize_t i, j, k, p
cdef int material
cdef float phi = 0.0
if nx == 1 or ny == 1:
pass
else:
for i in prange(1, nx, nogil=True, schedule='static', num_threads=nthreads):
for j in range(1, ny):
for k in range(0, nz):
material = ID[2, i, j, k]
phi = 0.0
for p in range(0, maxpoles):
phi = phi + updatecoeffsdispersive[material, p * 3].real * Tz[p, i, j, k].real
Tz[p, i, j, k] = updatecoeffsdispersive[material, 1 + (p * 3)] * Tz[p, i, j, k] + updatecoeffsdispersive[material, 2 + (p * 3)] * Ez[i, j, k]
Ez[i, j, k] = updatecoeffsE[material, 0] * Ez[i, j, k] + updatecoeffsE[material, 1] * (Hy[i, j, k] - Hy[i - 1, j, k]) - updatecoeffsE[material, 2] * (Hx[i, j, k] - Hx[i, j - 1, k]) - updatecoeffsE[material, 4] * phi
cpdef void update_ez_dispersive_multipole_B(int nx, int ny, int nz, int nthreads, int maxpoles, complextype_t[:, ::1] updatecoeffsdispersive, np.uint32_t[:, :, :, ::1] ID, complextype_t[:, :, :, ::1] Tz, floattype_t[:, :, ::1] Ez):
"""This function updates the Ez field components when dispersive materials (with multiple poles) are present.
Args:
nx, ny, nz (int): Grid size in cells
maxpoles (int): Maximum number of poles
nthreads (int): Number of threads to use
updatecoeffs, T, ID, E (memoryviews): Access to update coeffients, temporary, ID and field component arrays
"""
cdef Py_ssize_t i, j, k, p
cdef int material
if nx == 1 or ny == 1:
pass
else:
for i in prange(1, nx, nogil=True, schedule='static', num_threads=nthreads):
for j in range(1, ny):
for k in range(0, nz):
material = ID[2, i, j, k]
for p in range(0, maxpoles):
Tz[p, i, j, k] = Tz[p, i, j, k] - updatecoeffsdispersive[material, 2 + (p * 3)] * Ez[i, j, k]
cpdef void update_ez_dispersive_1pole_A(int nx, int ny, int nz, int nthreads, floattype_t[:, ::1] updatecoeffsE, complextype_t[:, ::1] updatecoeffsdispersive, np.uint32_t[:, :, :, ::1] ID, complextype_t[:, :, :, ::1] Tz, floattype_t[:, :, ::1] Ez, floattype_t[:, :, ::1] Hx, floattype_t[:, :, ::1] Hy):
"""This function updates the Ez field components when dispersive materials (with 1 pole) are present.
Args:
nx, ny, nz (int): Grid size in cells
nthreads (int): Number of threads to use
updatecoeffs, T, ID, E, H (memoryviews): Access to update coeffients, temporary, ID and field component arrays
"""
cdef Py_ssize_t i, j, k
cdef int material
cdef float phi = 0.0
if nx == 1 or ny == 1: if nx == 1 or ny == 1:
pass pass
else: else:
@@ -371,8 +226,8 @@ cpdef void update_ez_dispersive_1pole_A(int nx, int ny, int nz, int nthreads, fl
Ez[i, j, k] = updatecoeffsE[material, 0] * Ez[i, j, k] + updatecoeffsE[material, 1] * (Hy[i, j, k] - Hy[i - 1, j, k]) - updatecoeffsE[material, 2] * (Hx[i, j, k] - Hx[i, j - 1, k]) - updatecoeffsE[material, 4] * phi Ez[i, j, k] = updatecoeffsE[material, 0] * Ez[i, j, k] + updatecoeffsE[material, 1] * (Hy[i, j, k] - Hy[i - 1, j, k]) - updatecoeffsE[material, 2] * (Hx[i, j, k] - Hx[i, j - 1, k]) - updatecoeffsE[material, 4] * phi
cpdef void update_ez_dispersive_1pole_B(int nx, int ny, int nz, int nthreads, complextype_t[:, ::1] updatecoeffsdispersive, np.uint32_t[:, :, :, ::1] ID, complextype_t[:, :, :, ::1] Tz, floattype_t[:, :, ::1] Ez): cpdef void update_electric_dispersive_1pole_B(int nx, int ny, int nz, int nthreads, complextype_t[:, ::1] updatecoeffsdispersive, np.uint32_t[:, :, :, ::1] ID, complextype_t[:, :, :, ::1] Tx, complextype_t[:, :, :, ::1] Ty, complextype_t[:, :, :, ::1] Tz, floattype_t[:, :, ::1] Ex, floattype_t[:, :, ::1] Ey, floattype_t[:, :, ::1] Ez):
"""This function updates the Ez field components when dispersive materials (with 1 pole) are present. """This function updates the Ex field components when dispersive materials (with 1 pole) are present.
Args: Args:
nx, ny, nz (int): Grid size in cells nx, ny, nz (int): Grid size in cells
@@ -383,6 +238,27 @@ cpdef void update_ez_dispersive_1pole_B(int nx, int ny, int nz, int nthreads, co
cdef Py_ssize_t i, j, k cdef Py_ssize_t i, j, k
cdef int material cdef int material
# Ex component
if ny == 1 or nz == 1:
pass
else:
for i in prange(0, nx, nogil=True, schedule='static', num_threads=nthreads):
for j in range(1, ny):
for k in range(1, nz):
material = ID[0, i, j, k]
Tx[0, i, j, k] = Tx[0, i, j, k] - updatecoeffsdispersive[material, 2] * Ex[i, j, k]
# Ey component
if nx == 1 or nz == 1:
pass
else:
for i in prange(1, nx, nogil=True, schedule='static', num_threads=nthreads):
for j in range(0, ny):
for k in range(1, nz):
material = ID[1, i, j, k]
Ty[0, i, j, k] = Ty[0, i, j, k] - updatecoeffsdispersive[material, 2] * Ey[i, j, k]
# Ez component
if nx == 1 or ny == 1: if nx == 1 or ny == 1:
pass pass
else: else:
@@ -393,10 +269,10 @@ cpdef void update_ez_dispersive_1pole_B(int nx, int ny, int nz, int nthreads, co
Tz[0, i, j, k] = Tz[0, i, j, k] - updatecoeffsdispersive[material, 2] * Ez[i, j, k] Tz[0, i, j, k] = Tz[0, i, j, k] - updatecoeffsdispersive[material, 2] * Ez[i, j, k]
######################################### ##########################
# Magnetic field updates - Hx component # # Magnetic field updates #
######################################### ##########################
cpdef void update_hx(int nx, int ny, int nz, int nthreads, floattype_t[:, ::1] updatecoeffsH, np.uint32_t[:, :, :, ::1] ID, floattype_t[:, :, ::1] Hx, floattype_t[:, :, ::1] Ey, floattype_t[:, :, ::1] Ez): cpdef void update_magnetic(int nx, int ny, int nz, int nthreads, floattype_t[:, ::1] updatecoeffsH, np.uint32_t[:, :, :, ::1] ID, floattype_t[:, :, ::1] Ex, floattype_t[:, :, ::1] Ey, floattype_t[:, :, ::1] Ez, floattype_t[:, :, ::1] Hx, floattype_t[:, :, ::1] Hy, floattype_t[:, :, ::1] Hz):
"""This function updates the Hx field components. """This function updates the Hx field components.
Args: Args:
@@ -408,6 +284,7 @@ cpdef void update_hx(int nx, int ny, int nz, int nthreads, floattype_t[:, ::1] u
cdef Py_ssize_t i, j, k cdef Py_ssize_t i, j, k
cdef int material cdef int material
# Hx component
if nx == 1: if nx == 1:
pass pass
else: else:
@@ -417,22 +294,7 @@ cpdef void update_hx(int nx, int ny, int nz, int nthreads, floattype_t[:, ::1] u
material = ID[3, i, j, k] material = ID[3, i, j, k]
Hx[i, j, k] = updatecoeffsH[material, 0] * Hx[i, j, k] - updatecoeffsH[material, 2] * (Ez[i, j + 1, k] - Ez[i, j, k]) + updatecoeffsH[material, 3] * (Ey[i, j, k + 1] - Ey[i, j, k]) Hx[i, j, k] = updatecoeffsH[material, 0] * Hx[i, j, k] - updatecoeffsH[material, 2] * (Ez[i, j + 1, k] - Ez[i, j, k]) + updatecoeffsH[material, 3] * (Ey[i, j, k + 1] - Ey[i, j, k])
# Hy component
#########################################
# Magnetic field updates - Hy component #
#########################################
cpdef void update_hy(int nx, int ny, int nz, int nthreads, floattype_t[:, ::1] updatecoeffsH, np.uint32_t[:, :, :, ::1] ID, floattype_t[:, :, ::1] Hy, floattype_t[:, :, ::1] Ex, floattype_t[:, :, ::1] Ez):
"""This function updates the Hy field components.
Args:
nx, ny, nz (int): Grid size in cells
nthreads (int): Number of threads to use
updatecoeffs, ID, E, H (memoryviews): Access to update coeffients, ID and field component arrays
"""
cdef Py_ssize_t i, j, k
cdef int material
if ny == 1: if ny == 1:
pass pass
else: else:
@@ -442,22 +304,7 @@ cpdef void update_hy(int nx, int ny, int nz, int nthreads, floattype_t[:, ::1] u
material = ID[4, i, j, k] material = ID[4, i, j, k]
Hy[i, j, k] = updatecoeffsH[material, 0] * Hy[i, j, k] - updatecoeffsH[material, 3] * (Ex[i, j, k + 1] - Ex[i, j, k]) + updatecoeffsH[material, 1] * (Ez[i + 1, j, k] - Ez[i, j, k]) Hy[i, j, k] = updatecoeffsH[material, 0] * Hy[i, j, k] - updatecoeffsH[material, 3] * (Ex[i, j, k + 1] - Ex[i, j, k]) + updatecoeffsH[material, 1] * (Ez[i + 1, j, k] - Ez[i, j, k])
# Hz component
#########################################
# Magnetic field updates - Hz component #
#########################################
cpdef void update_hz(int nx, int ny, int nz, int nthreads, floattype_t[:, ::1] updatecoeffsH, np.uint32_t[:, :, :, ::1] ID, floattype_t[:, :, ::1] Hz, floattype_t[:, :, ::1] Ex, floattype_t[:, :, ::1] Ey):
"""This function updates the Hz field components.
Args:
nx, ny, nz (int): Grid size in cells
nthreads (int): Number of threads to use
updatecoeffs, ID, E, H (memoryviews): Access to update coeffients, ID and field component arrays
"""
cdef Py_ssize_t i, j, k
cdef int material
if nz == 1: if nz == 1:
pass pass
else: else: