publicImage/gprMax
1
0
派生 0
镜像自地址 https://gitee.com/sunhf/gprMax.git 已同步 2025-08-07 04:56:51 +08:00
代码 工单 软件包 项目 版本发布 百科 动态

Performance improvements to Cython code.

浏览代码
这个提交包含在：
Craig Warren
2016-02-23 17:35:10 +00:00
父节点 927b637d33
当前提交 d064be785a
共有 3 个文件被更改，包括 408 次插入 和 342 次删除
下载 Patch 文件 下载 Diff 文件 展开所有文件 折叠所有文件

222
gprMax/fields_update.pyx
查看文件

@@ -25,7 +25,7 @@ from gprMax.constants cimport floattype_t, complextype_t
######################################### #########################################
# Electric field updates - Ex component # # Electric field updates - Ex component #
######################################### #########################################
cpdef update_ex(int nx, int ny, int nz, int nthreads, floattype_t[:, :] updatecoeffsE, np.uint32_t[:, :, :, :] ID, floattype_t[:, :, :] Ex, floattype_t[:, :, :] Hy, floattype_t[:, :, :] Hz): 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):
"""This function updates the Ex field components. """This function updates the Ex field components.
Args: Args:
@@ -34,19 +34,20 @@ cpdef update_ex(int nx, int ny, int nz, int nthreads, floattype_t[:, :] updateco
updatecoeffs, ID, E, H (memoryviews): Access to update coeffients, ID and field component arrays updatecoeffs, ID, E, H (memoryviews): Access to update coeffients, ID and field component arrays
""" """
cdef int i, j, k, listIndex cdef Py_ssize_t i, j, k
cdef int material
if ny == 1 or nz == 1: if ny == 1 or nz == 1:
pass pass
else: else:
for i in prange(0, nx, nogil=True, schedule='static', chunksize=1, num_threads=nthreads): for i in prange(0, nx, nogil=True, schedule='static', num_threads=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):
listIndex = ID[0, i, j, k] material = ID[0, i, j, k]
Ex[i, j, k] = updatecoeffsE[listIndex, 0] * Ex[i, j, k] + updatecoeffsE[listIndex, 2] * (Hz[i, j, k] - Hz[i, j - 1, k]) - updatecoeffsE[listIndex, 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])
cpdef update_ex_dispersive_multipole_A(int nx, int ny, int nz, int nthreads, int maxpoles, floattype_t[:, :] updatecoeffsE, complextype_t[:, :] updatecoeffsdispersive, np.uint32_t[:, :, :, :] ID, complextype_t[:, :, :, :] Tx, floattype_t[:, :, :] Ex, floattype_t[:, :, :] Hy, floattype_t[:, :, :] Hz): 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):
"""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:
@@ -56,23 +57,24 @@ cpdef update_ex_dispersive_multipole_A(int nx, int ny, int nz, int nthreads, int
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 int i, j, k, listIndex, p cdef Py_ssize_t i, j, k, p
cdef int material
cdef float phi = 0.0 cdef float phi = 0.0
if ny == 1 or nz == 1: if ny == 1 or nz == 1:
pass pass
else: else:
for i in prange(0, nx, nogil=True, schedule='static', chunksize=1, num_threads=nthreads): for i in prange(0, nx, nogil=True, schedule='static', num_threads=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):
listIndex = ID[0, i, j, k] material = ID[0, i, j, k]
phi = 0.0 phi = 0.0
for p in range(0, maxpoles): for p in range(0, maxpoles):
phi = phi + updatecoeffsdispersive[listIndex, p * 3].real * Tx[p, i, j, k].real phi = phi + updatecoeffsdispersive[material, p * 3].real * Tx[p, i, j, k].real
Tx[p, i, j, k] = updatecoeffsdispersive[listIndex, 1 + (p * 3)] * Tx[p, i, j, k] + updatecoeffsdispersive[listIndex, 2 + (p * 3)] * Ex[i, j, k] Tx[p, i, j, k] = updatecoeffsdispersive[material, 1 + (p * 3)] * Tx[p, i, j, k] + updatecoeffsdispersive[material, 2 + (p * 3)] * Ex[i, j, k]
Ex[i, j, k] = updatecoeffsE[listIndex, 0] * Ex[i, j, k] + updatecoeffsE[listIndex, 2] * (Hz[i, j, k] - Hz[i, j - 1, k]) - updatecoeffsE[listIndex, 3] * (Hy[i, j, k] - Hy[i, j, k - 1]) - updatecoeffsE[listIndex, 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 update_ex_dispersive_multipole_B(int nx, int ny, int nz, int nthreads, int maxpoles, complextype_t[:, :] updatecoeffsdispersive, np.uint32_t[:, :, :, :] ID, complextype_t[:, :, :, :] Tx, floattype_t[:, :, :] Ex): 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):
"""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:
@@ -82,20 +84,21 @@ cpdef update_ex_dispersive_multipole_B(int nx, int ny, int nz, int nthreads, int
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 int i, j, k, listIndex, p cdef Py_ssize_t i, j, k, p
cdef int material
if ny == 1 or nz == 1: if ny == 1 or nz == 1:
pass pass
else: else:
for i in prange(0, nx, nogil=True, schedule='static', chunksize=1, num_threads=nthreads): for i in prange(0, nx, nogil=True, schedule='static', num_threads=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):
listIndex = ID[0, i, j, k] material = ID[0, i, j, k]
for p in range(0, maxpoles): for p in range(0, maxpoles):
Tx[p, i, j, k] = Tx[p, i, j, k] - updatecoeffsdispersive[listIndex, 2 + (p * 3)] * Ex[i, j, k] Tx[p, i, j, k] = Tx[p, i, j, k] - updatecoeffsdispersive[material, 2 + (p * 3)] * Ex[i, j, k]
cpdef update_ex_dispersive_1pole_A(int nx, int ny, int nz, int nthreads, floattype_t[:, :] updatecoeffsE, complextype_t[:, :] updatecoeffsdispersive, np.uint32_t[:, :, :, :] ID, complextype_t[:, :, :, :] Tx, floattype_t[:, :, :] Ex, floattype_t[:, :, :] Hy, floattype_t[:, :, :] Hz): 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):
"""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:
@@ -104,22 +107,23 @@ cpdef update_ex_dispersive_1pole_A(int nx, int ny, int nz, int nthreads, floatty
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 int i, j, k, listIndex cdef Py_ssize_t i, j, k
cdef int material
cdef float phi = 0.0 cdef float phi = 0.0
if ny == 1 or nz == 1: if ny == 1 or nz == 1:
pass pass
else: else:
for i in prange(0, nx, nogil=True, schedule='static', chunksize=1, num_threads=nthreads): for i in prange(0, nx, nogil=True, schedule='static', num_threads=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):
listIndex = ID[0, i, j, k] material = ID[0, i, j, k]
phi = updatecoeffsdispersive[listIndex, 0].real * Tx[0, i, j, k].real phi = updatecoeffsdispersive[material, 0].real * Tx[0, i, j, k].real
Tx[0, i, j, k] = updatecoeffsdispersive[listIndex, 1] * Tx[0, i, j, k] + updatecoeffsdispersive[listIndex, 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[listIndex, 0] * Ex[i, j, k] + updatecoeffsE[listIndex, 2] * (Hz[i, j, k] - Hz[i, j - 1, k]) - updatecoeffsE[listIndex, 3] * (Hy[i, j, k] - Hy[i, j, k - 1]) - updatecoeffsE[listIndex, 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 update_ex_dispersive_1pole_B(int nx, int ny, int nz, int nthreads, complextype_t[:, :] updatecoeffsdispersive, np.uint32_t[:, :, :, :] ID, complextype_t[:, :, :, :] Tx, floattype_t[:, :, :] Ex): 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. """This function updates the Ex field components when dispersive materials (with 1 pole) are present.
Args: Args:
@@ -128,22 +132,23 @@ cpdef update_ex_dispersive_1pole_B(int nx, int ny, int nz, int nthreads, complex
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 int i, j, k, listIndex cdef Py_ssize_t i, j, k
cdef int material
if ny == 1 or nz == 1: if ny == 1 or nz == 1:
pass pass
else: else:
for i in prange(0, nx, nogil=True, schedule='static', chunksize=1, num_threads=nthreads): for i in prange(0, nx, nogil=True, schedule='static', num_threads=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):
listIndex = ID[0, i, j, k] material = ID[0, i, j, k]
Tx[0, i, j, k] = Tx[0, i, j, k] - updatecoeffsdispersive[listIndex, 2] * Ex[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 # # Electric field updates - Ey component #
######################################### #########################################
cpdef update_ey(int nx, int ny, int nz, int nthreads, floattype_t[:, :] updatecoeffsE, np.uint32_t[:, :, :, :] ID, floattype_t[:, :, :] Ey, floattype_t[:, :, :] Hx, floattype_t[:, :, :] Hz): 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. """This function updates the Ey field components.
Args: Args:
@@ -152,19 +157,20 @@ cpdef update_ey(int nx, int ny, int nz, int nthreads, floattype_t[:, :] updateco
updatecoeffs, ID, E, H (memoryviews): Access to update coeffients, ID and field component arrays updatecoeffs, ID, E, H (memoryviews): Access to update coeffients, ID and field component arrays
""" """
cdef int i, j, k, listIndex cdef Py_ssize_t i, j, k
cdef int material
if nx == 1 or nz == 1: if nx == 1 or nz == 1:
pass pass
else: else:
for i in prange(1, nx, nogil=True, schedule='static', chunksize=1, num_threads=nthreads): for i in prange(1, nx, nogil=True, schedule='static', num_threads=nthreads):
for j in range(0, ny): for j in range(0, ny):
for k in range(1, nz): for k in range(1, nz):
listIndex = ID[1, i, j, k] material = ID[1, i, j, k]
Ey[i, j, k] = updatecoeffsE[listIndex, 0] * Ey[i, j, k] + updatecoeffsE[listIndex, 3] * (Hx[i, j, k] - Hx[i, j, k - 1]) - updatecoeffsE[listIndex, 1] * (Hz[i, j, k] - Hz[i - 1, 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 update_ey_dispersive_multipole_A(int nx, int ny, int nz, int nthreads, int maxpoles, floattype_t[:, :] updatecoeffsE, complextype_t[:, :] updatecoeffsdispersive, np.uint32_t[:, :, :, :] ID, complextype_t[:, :, :, :] Ty, floattype_t[:, :, :] Ey, floattype_t[:, :, :] Hx, floattype_t[:, :, :] Hz): 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. """This function updates the Ey field components when dispersive materials (with multiple poles) are present.
Args: Args:
@@ -174,24 +180,25 @@ cpdef update_ey_dispersive_multipole_A(int nx, int ny, int nz, int nthreads, int
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 int i, j, k, listIndex, p cdef Py_ssize_t i, j, k, p
cdef int material
cdef float phi = 0.0 cdef float phi = 0.0
if nx == 1 or nz == 1: if nx == 1 or nz == 1:
pass pass
else: else:
for i in prange(1, nx, nogil=True, schedule='static', chunksize=1, num_threads=nthreads): for i in prange(1, nx, nogil=True, schedule='static', num_threads=nthreads):
for j in range(0, ny): for j in range(0, ny):
for k in range(1, nz): for k in range(1, nz):
listIndex = ID[1, i, j, k] material = ID[1, i, j, k]
phi = 0.0 phi = 0.0
for p in range(0, maxpoles): for p in range(0, maxpoles):
phi = phi + updatecoeffsdispersive[listIndex, p * 3].real * Ty[p, i, j, k].real phi = phi + updatecoeffsdispersive[material, p * 3].real * Ty[p, i, j, k].real
Ty[p, i, j, k] = updatecoeffsdispersive[listIndex, 1 + (p * 3)] * Ty[p, i, j, k] + updatecoeffsdispersive[listIndex, 2 + (p * 3)] * Ey[i, j, k] 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[listIndex, 0] * Ey[i, j, k] + updatecoeffsE[listIndex, 3] * (Hx[i, j, k] - Hx[i, j, k - 1]) - updatecoeffsE[listIndex, 1] * (Hz[i, j, k] - Hz[i - 1, j, k]) - updatecoeffsE[listIndex, 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
cpdef update_ey_dispersive_multipole_B(int nx, int ny, int nz, int nthreads, int maxpoles, complextype_t[:, :] updatecoeffsdispersive, np.uint32_t[:, :, :, :] ID, complextype_t[:, :, :, :] Ty, floattype_t[:, :, :] Ey): 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. """This function updates the Ey field components when dispersive materials (with multiple poles) are present.
Args: Args:
@@ -201,20 +208,21 @@ cpdef update_ey_dispersive_multipole_B(int nx, int ny, int nz, int nthreads, int
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 int i, j, k, listIndex, p cdef Py_ssize_t i, j, k, p
cdef int material
if nx == 1 or nz == 1: if nx == 1 or nz == 1:
pass pass
else: else:
for i in prange(1, nx, nogil=True, schedule='static', chunksize=1, num_threads=nthreads): for i in prange(1, nx, nogil=True, schedule='static', num_threads=nthreads):
for j in range(0, ny): for j in range(0, ny):
for k in range(1, nz): for k in range(1, nz):
listIndex = ID[1, i, j, k] material = ID[1, i, j, k]
for p in range(0, maxpoles): for p in range(0, maxpoles):
Ty[p, i, j, k] = Ty[p, i, j, k] - updatecoeffsdispersive[listIndex, 2 + (p * 3)] * Ey[i, j, k] Ty[p, i, j, k] = Ty[p, i, j, k] - updatecoeffsdispersive[material, 2 + (p * 3)] * Ey[i, j, k]
cpdef update_ey_dispersive_1pole_A(int nx, int ny, int nz, int nthreads, floattype_t[:, :] updatecoeffsE, complextype_t[:, :] updatecoeffsdispersive, np.uint32_t[:, :, :, :] ID, complextype_t[:, :, :, :] Ty, floattype_t[:, :, :] Ey, floattype_t[:, :, :] Hx, floattype_t[:, :, :] Hz): 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. """This function updates the Ey field components when dispersive materials (with 1 pole) are present.
Args: Args:
@@ -223,22 +231,23 @@ cpdef update_ey_dispersive_1pole_A(int nx, int ny, int nz, int nthreads, floatty
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 int i, j, k, listIndex cdef Py_ssize_t i, j, k
cdef int material
cdef float phi = 0.0 cdef float phi = 0.0
if nx == 1 or nz == 1: if nx == 1 or nz == 1:
pass pass
else: else:
for i in prange(1, nx, nogil=True, schedule='static', chunksize=1, num_threads=nthreads): for i in prange(1, nx, nogil=True, schedule='static', num_threads=nthreads):
for j in range(0, ny): for j in range(0, ny):
for k in range(1, nz): for k in range(1, nz):
listIndex = ID[1, i, j, k] material = ID[1, i, j, k]
phi = updatecoeffsdispersive[listIndex, 0].real * Ty[0, i, j, k].real phi = updatecoeffsdispersive[material, 0].real * Ty[0, i, j, k].real
Ty[0, i, j, k] = updatecoeffsdispersive[listIndex, 1] * Ty[0, i, j, k] + updatecoeffsdispersive[listIndex, 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[listIndex, 0] * Ey[i, j, k] + updatecoeffsE[listIndex, 3] * (Hx[i, j, k] - Hx[i, j, k - 1]) - updatecoeffsE[listIndex, 1] * (Hz[i, j, k] - Hz[i - 1, j, k]) - updatecoeffsE[listIndex, 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
cpdef update_ey_dispersive_1pole_B(int nx, int ny, int nz, int nthreads, complextype_t[:, :] updatecoeffsdispersive, np.uint32_t[:, :, :, :] ID, complextype_t[:, :, :, :] Ty, floattype_t[:, :, :] Ey): 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. """This function updates the Ey field components when dispersive materials (with 1 pole) are present.
Args: Args:
@@ -247,22 +256,23 @@ cpdef update_ey_dispersive_1pole_B(int nx, int ny, int nz, int nthreads, complex
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 int i, j, k, listIndex cdef Py_ssize_t i, j, k
cdef int material
if nx == 1 or nz == 1: if nx == 1 or nz == 1:
pass pass
else: else:
for i in prange(1, nx, nogil=True, schedule='static', chunksize=1, num_threads=nthreads): for i in prange(1, nx, nogil=True, schedule='static', num_threads=nthreads):
for j in range(0, ny): for j in range(0, ny):
for k in range(1, nz): for k in range(1, nz):
listIndex = ID[1, i, j, k] material = ID[1, i, j, k]
Ty[0, i, j, k] = Ty[0, i, j, k] - updatecoeffsdispersive[listIndex, 2] * Ey[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 # # Electric field updates - Ez component #
######################################### #########################################
cpdef update_ez(int nx, int ny, int nz, int nthreads, floattype_t[:, :] updatecoeffsE, np.uint32_t[:, :, :, :] ID, floattype_t[:, :, :] Ez, floattype_t[:, :, :] Hx, floattype_t[:, :, :] Hy): 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. """This function updates the Ez field components.
Args: Args:
@@ -271,19 +281,20 @@ cpdef update_ez(int nx, int ny, int nz, int nthreads, floattype_t[:, :] updateco
updatecoeffs, ID, E, H (memoryviews): Access to update coeffients, ID and field component arrays updatecoeffs, ID, E, H (memoryviews): Access to update coeffients, ID and field component arrays
""" """
cdef int i, j, k, listIndex cdef Py_ssize_t i, j, k
cdef int material
if nx == 1 or ny == 1: if nx == 1 or ny == 1:
pass pass
else: else:
for i in prange(1, nx, nogil=True, schedule='static', chunksize=1, num_threads=nthreads): for i in prange(1, nx, nogil=True, schedule='static', num_threads=nthreads):
for j in range(1, ny): for j in range(1, ny):
for k in range(0, nz): for k in range(0, nz):
listIndex = ID[2, i, j, k] material = ID[2, i, j, k]
Ez[i, j, k] = updatecoeffsE[listIndex, 0] * Ez[i, j, k] + updatecoeffsE[listIndex, 1] * (Hy[i, j, k] - Hy[i - 1, j, k]) - updatecoeffsE[listIndex, 2] * (Hx[i, j, k] - Hx[i, j - 1, 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 update_ez_dispersive_multipole_A(int nx, int ny, int nz, int nthreads, int maxpoles, floattype_t[:, :] updatecoeffsE, complextype_t[:, :] updatecoeffsdispersive, np.uint32_t[:, :, :, :] ID, complextype_t[:, :, :, :] Tz, floattype_t[:, :, :] Ez, floattype_t[:, :, :] Hx, floattype_t[:, :, :] Hy): 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. """This function updates the Ez field components when dispersive materials (with multiple poles) are present.
Args: Args:
@@ -293,24 +304,25 @@ cpdef update_ez_dispersive_multipole_A(int nx, int ny, int nz, int nthreads, int
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 int i, j, k, listIndex, p cdef Py_ssize_t i, j, k, p
cdef int material
cdef float phi = 0.0 cdef float phi = 0.0
if nx == 1 or ny == 1: if nx == 1 or ny == 1:
pass pass
else: else:
for i in prange(1, nx, nogil=True, schedule='static', chunksize=1, num_threads=nthreads): for i in prange(1, nx, nogil=True, schedule='static', num_threads=nthreads):
for j in range(1, ny): for j in range(1, ny):
for k in range(0, nz): for k in range(0, nz):
listIndex = ID[2, i, j, k] material = ID[2, i, j, k]
phi = 0.0 phi = 0.0
for p in range(0, maxpoles): for p in range(0, maxpoles):
phi = phi + updatecoeffsdispersive[listIndex, p * 3].real * Tz[p, i, j, k].real phi = phi + updatecoeffsdispersive[material, p * 3].real * Tz[p, i, j, k].real
Tz[p, i, j, k] = updatecoeffsdispersive[listIndex, 1 + (p * 3)] * Tz[p, i, j, k] + updatecoeffsdispersive[listIndex, 2 + (p * 3)] * Ez[i, j, k] 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[listIndex, 0] * Ez[i, j, k] + updatecoeffsE[listIndex, 1] * (Hy[i, j, k] - Hy[i - 1, j, k]) - updatecoeffsE[listIndex, 2] * (Hx[i, j, k] - Hx[i, j - 1, k]) - updatecoeffsE[listIndex, 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 update_ez_dispersive_multipole_B(int nx, int ny, int nz, int nthreads, int maxpoles, complextype_t[:, :] updatecoeffsdispersive, np.uint32_t[:, :, :, :] ID, complextype_t[:, :, :, :] Tz, floattype_t[:, :, :] Ez): 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. """This function updates the Ez field components when dispersive materials (with multiple poles) are present.
Args: Args:
@@ -320,20 +332,21 @@ cpdef update_ez_dispersive_multipole_B(int nx, int ny, int nz, int nthreads, int
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 int i, j, k, listIndex, p cdef Py_ssize_t i, j, k, p
cdef int material
if nx == 1 or ny == 1: if nx == 1 or ny == 1:
pass pass
else: else:
for i in prange(1, nx, nogil=True, schedule='static', chunksize=1, num_threads=nthreads): for i in prange(1, nx, nogil=True, schedule='static', num_threads=nthreads):
for j in range(1, ny): for j in range(1, ny):
for k in range(0, nz): for k in range(0, nz):
listIndex = ID[2, i, j, k] material = ID[2, i, j, k]
for p in range(0, maxpoles): for p in range(0, maxpoles):
Tz[p, i, j, k] = Tz[p, i, j, k] - updatecoeffsdispersive[listIndex, 2 + (p * 3)] * Ez[i, j, k] Tz[p, i, j, k] = Tz[p, i, j, k] - updatecoeffsdispersive[material, 2 + (p * 3)] * Ez[i, j, k]
cpdef update_ez_dispersive_1pole_A(int nx, int ny, int nz, int nthreads, floattype_t[:, :] updatecoeffsE, complextype_t[:, :] updatecoeffsdispersive, np.uint32_t[:, :, :, :] ID, complextype_t[:, :, :, :] Tz, floattype_t[:, :, :] Ez, floattype_t[:, :, :] Hx, floattype_t[:, :, :] Hy): 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. """This function updates the Ez field components when dispersive materials (with 1 pole) are present.
Args: Args:
@@ -342,22 +355,23 @@ cpdef update_ez_dispersive_1pole_A(int nx, int ny, int nz, int nthreads, floatty
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 int i, j, k, listIndex cdef Py_ssize_t i, j, k
cdef int material
cdef float phi = 0.0 cdef float phi = 0.0
if nx == 1 or ny == 1: if nx == 1 or ny == 1:
pass pass
else: else:
for i in prange(1, nx, nogil=True, schedule='static', chunksize=1, num_threads=nthreads): for i in prange(1, nx, nogil=True, schedule='static', num_threads=nthreads):
for j in range(1, ny): for j in range(1, ny):
for k in range(0, nz): for k in range(0, nz):
listIndex = ID[2, i, j, k] material = ID[2, i, j, k]
phi = updatecoeffsdispersive[listIndex, 0].real * Tz[0, i, j, k].real phi = updatecoeffsdispersive[material, 0].real * Tz[0, i, j, k].real
Tz[0, i, j, k] = updatecoeffsdispersive[listIndex, 1] * Tz[0, i, j, k] + updatecoeffsdispersive[listIndex, 2] * Ez[i, j, k] Tz[0, i, j, k] = updatecoeffsdispersive[material, 1] * Tz[0, i, j, k] + updatecoeffsdispersive[material, 2] * Ez[i, j, k]
Ez[i, j, k] = updatecoeffsE[listIndex, 0] * Ez[i, j, k] + updatecoeffsE[listIndex, 1] * (Hy[i, j, k] - Hy[i - 1, j, k]) - updatecoeffsE[listIndex, 2] * (Hx[i, j, k] - Hx[i, j - 1, k]) - updatecoeffsE[listIndex, 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 update_ez_dispersive_1pole_B(int nx, int ny, int nz, int nthreads, complextype_t[:, :] updatecoeffsdispersive, np.uint32_t[:, :, :, :] ID, complextype_t[:, :, :, :] Tz, floattype_t[:, :, :] Ez): 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):
"""This function updates the Ez field components when dispersive materials (with 1 pole) are present. """This function updates the Ez field components when dispersive materials (with 1 pole) are present.
Args: Args:
@@ -366,22 +380,23 @@ cpdef update_ez_dispersive_1pole_B(int nx, int ny, int nz, int nthreads, complex
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 int i, j, k, listIndex cdef Py_ssize_t i, j, k
cdef int material
if nx == 1 or ny == 1: if nx == 1 or ny == 1:
pass pass
else: else:
for i in prange(1, nx, nogil=True, schedule='static', chunksize=1, num_threads=nthreads): for i in prange(1, nx, nogil=True, schedule='static', num_threads=nthreads):
for j in range(1, ny): for j in range(1, ny):
for k in range(0, nz): for k in range(0, nz):
listIndex = ID[2, i, j, k] material = ID[2, i, j, k]
Tz[0, i, j, k] = Tz[0, i, j, k] - updatecoeffsdispersive[listIndex, 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 - Hx component #
######################################### #########################################
cpdef update_hx(int nx, int ny, int nz, int nthreads, floattype_t[:, :] updatecoeffsH, np.uint32_t[:, :, :, :] ID, floattype_t[:, :, :] Hx, floattype_t[:, :, :] Ey, floattype_t[:, :, :] Ez): 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):
"""This function updates the Hx field components. """This function updates the Hx field components.
Args: Args:
@@ -390,22 +405,23 @@ cpdef update_hx(int nx, int ny, int nz, int nthreads, floattype_t[:, :] updateco
updatecoeffs, ID, E, H (memoryviews): Access to update coeffients, ID and field component arrays updatecoeffs, ID, E, H (memoryviews): Access to update coeffients, ID and field component arrays
""" """
cdef int i, j, k, listIndex cdef Py_ssize_t i, j, k
cdef int material
if nx == 1: if nx == 1:
pass pass
else: else:
for i in prange(1, nx, nogil=True, schedule='static', chunksize=1, num_threads=nthreads): for i in prange(1, nx, nogil=True, schedule='static', num_threads=nthreads):
for j in range(0, ny): for j in range(0, ny):
for k in range(0, nz): for k in range(0, nz):
listIndex = ID[3, i, j, k] material = ID[3, i, j, k]
Hx[i, j, k] = updatecoeffsH[listIndex, 0] * Hx[i, j, k] - updatecoeffsH[listIndex, 2] * (Ez[i, j + 1, k] - Ez[i, j, k]) + updatecoeffsH[listIndex, 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])
######################################### #########################################
# Magnetic field updates - Hy component # # Magnetic field updates - Hy component #
######################################### #########################################
cpdef update_hy(int nx, int ny, int nz, int nthreads, floattype_t[:, :] updatecoeffsH, np.uint32_t[:, :, :, :] ID, floattype_t[:, :, :] Hy, floattype_t[:, :, :] Ex, floattype_t[:, :, :] Ez): 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. """This function updates the Hy field components.
Args: Args:
@@ -414,22 +430,23 @@ cpdef update_hy(int nx, int ny, int nz, int nthreads, floattype_t[:, :] updateco
updatecoeffs, ID, E, H (memoryviews): Access to update coeffients, ID and field component arrays updatecoeffs, ID, E, H (memoryviews): Access to update coeffients, ID and field component arrays
""" """
cdef int i, j, k, listIndex cdef Py_ssize_t i, j, k
cdef int material
if ny == 1: if ny == 1:
pass pass
else: else:
for i in prange(0, nx, nogil=True, schedule='static', chunksize=1, num_threads=nthreads): for i in prange(0, nx, nogil=True, schedule='static', num_threads=nthreads):
for j in range(1, ny): for j in range(1, ny):
for k in range(0, nz): for k in range(0, nz):
listIndex = ID[4, i, j, k] material = ID[4, i, j, k]
Hy[i, j, k] = updatecoeffsH[listIndex, 0] * Hy[i, j, k] - updatecoeffsH[listIndex, 3] * (Ex[i, j, k + 1] - Ex[i, j, k]) + updatecoeffsH[listIndex, 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])
######################################### #########################################
# Magnetic field updates - Hz component # # Magnetic field updates - Hz component #
######################################### #########################################
cpdef update_hz(int nx, int ny, int nz, int nthreads, floattype_t[:, :] updatecoeffsH, np.uint32_t[:, :, :, :] ID, floattype_t[:, :, :] Hz, floattype_t[:, :, :] Ex, floattype_t[:, :, :] Ey): 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. """This function updates the Hz field components.
Args: Args:
@@ -438,14 +455,15 @@ cpdef update_hz(int nx, int ny, int nz, int nthreads, floattype_t[:, :] updateco
updatecoeffs, ID, E, H (memoryviews): Access to update coeffients, ID and field component arrays updatecoeffs, ID, E, H (memoryviews): Access to update coeffients, ID and field component arrays
""" """
cdef int i, j, k, listIndex cdef Py_ssize_t i, j, k
cdef int material
if nz == 1: if nz == 1:
pass pass
else: else:
for i in prange(0, nx, nogil=True, schedule='static', chunksize=1, num_threads=nthreads): for i in prange(0, nx, nogil=True, schedule='static', num_threads=nthreads):
for j in range(0, ny): for j in range(0, ny):
for k in range(1, nz): for k in range(1, nz):
listIndex = ID[5, i, j, k] material = ID[5, i, j, k]
Hz[i, j, k] = updatecoeffsH[listIndex, 0] * Hz[i, j, k] - updatecoeffsH[listIndex, 1] * (Ey[i + 1, j, k] - Ey[i, j, k]) + updatecoeffsH[listIndex, 2] * (Ex[i, j + 1, k] - Ex[i, j, k]) Hz[i, j, k] = updatecoeffsH[material, 0] * Hz[i, j, k] - updatecoeffsH[material, 1] * (Ey[i + 1, j, k] - Ey[i, j, k]) + updatecoeffsH[material, 2] * (Ex[i, j + 1, k] - Ex[i, j, k])

264
gprMax/pml_1order_update.pyx
查看文件

@@ -25,7 +25,7 @@ from gprMax.constants cimport floattype_t, complextype_t
############################################# #############################################
# Electric field PML updates - Ex component # # Electric field PML updates - Ex component #
############################################# #############################################
cpdef update_pml_1order_ex_yplus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, :] updatecoeffsE, np.uint32_t[:, :, :, :] ID, floattype_t[:, :, :] Ex, floattype_t[:, :, :] Hz, floattype_t[:, :, :, :] EPhi, floattype_t[:, :] RA, floattype_t[:, :] RB, floattype_t[:, :] RE, floattype_t[:, :] RF, float dy): cpdef void update_pml_1order_ex_yplus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, ::1] updatecoeffsE, np.uint32_t[:, :, :, ::1] ID, floattype_t[:, :, ::1] Ex, floattype_t[:, :, ::1] Hz, floattype_t[:, :, :, ::1] EPhi, floattype_t[:, ::1] RA, floattype_t[:, ::1] RB, floattype_t[:, ::1] RE, floattype_t[:, ::1] RF, float dy):
"""This function updates the Ex field components in the y stretching direction. """This function updates the Ex field components in the y stretching direction.
Args: Args:
@@ -36,22 +36,23 @@ cpdef update_pml_1order_ex_yplus(int xs, int xf, int ys, int yf, int zs, int zf,
dx, dy, dz (float): Spatial discretisation dx, dy, dz (float): Spatial discretisation
""" """
cdef int i, j, k, nx, ny, nz, listIndex cdef Py_ssize_t i, j, k
cdef int nx, ny, nz, material
cdef float dHz cdef float dHz
nx = xf - xs nx = xf - xs
ny = yf - ys ny = yf - ys
nz = zf - zs nz = zf - zs
for i in prange(0, nx, nogil=True, schedule='static', chunksize=1, num_threads=nthreads): for i in prange(0, nx, nogil=True, schedule='static', num_threads=nthreads):
for j in range(0, ny): for j in range(0, ny):
for k in range(0, nz): for k in range(0, nz):
listIndex = ID[0, i + xs, j + ys, k + zs] material = ID[0, i + xs, j + ys, k + zs]
dHz = (Hz[i + xs, j + ys, k + zs] - Hz[i + xs, j - 1 + ys, k + zs]) / dy dHz = (Hz[i + xs, j + ys, k + zs] - Hz[i + xs, j - 1 + ys, k + zs]) / dy
Ex[i + xs, j + ys, k + zs] = Ex[i + xs, j + ys, k + zs] + updatecoeffsE[listIndex, 4] * ((RA[0, j] - 1) * dHz + RB[0, j] * EPhi[0, i, j, k]) Ex[i + xs, j + ys, k + zs] = Ex[i + xs, j + ys, k + zs] + updatecoeffsE[material, 4] * ((RA[0, j] - 1) * dHz + RB[0, j] * EPhi[0, i, j, k])
EPhi[0, i, j, k] = RE[0, j] * EPhi[0, i, j, k] - RF[0, j] * dHz EPhi[0, i, j, k] = RE[0, j] * EPhi[0, i, j, k] - RF[0, j] * dHz
cpdef update_pml_1order_ex_yminus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, :] updatecoeffsE, np.uint32_t[:, :, :, :] ID, floattype_t[:, :, :] Ex, floattype_t[:, :, :] Hz, floattype_t[:, :, :, :] EPhi, floattype_t[:, :] RA, floattype_t[:, :] RB, floattype_t[:, :] RE, floattype_t[:, :] RF, float dy): cpdef void update_pml_1order_ex_yminus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, ::1] updatecoeffsE, np.uint32_t[:, :, :, ::1] ID, floattype_t[:, :, ::1] Ex, floattype_t[:, :, ::1] Hz, floattype_t[:, :, :, ::1] EPhi, floattype_t[:, ::1] RA, floattype_t[:, ::1] RB, floattype_t[:, ::1] RE, floattype_t[:, ::1] RF, float dy):
"""This function updates the Ex field components in the y stretching direction. """This function updates the Ex field components in the y stretching direction.
Args: Args:
@@ -62,22 +63,23 @@ cpdef update_pml_1order_ex_yminus(int xs, int xf, int ys, int yf, int zs, int zf
dx, dy, dz (float): Spatial discretisation dx, dy, dz (float): Spatial discretisation
""" """
cdef int i, j, k, nx, ny, nz, listIndex cdef Py_ssize_t i, j, k
cdef int nx, ny, nz, material
cdef float dHz cdef float dHz
nx = xf - xs nx = xf - xs
ny = yf - ys ny = yf - ys
nz = zf - zs nz = zf - zs
for i in prange(0, nx, nogil=True, schedule='static', chunksize=1, num_threads=nthreads): for i in prange(0, nx, nogil=True, schedule='static', num_threads=nthreads):
for j in range(0, ny): for j in range(0, ny):
for k in range(0, nz): for k in range(0, nz):
listIndex = ID[0, i + xs, yf - j, k + zs] material = ID[0, i + xs, yf - j, k + zs]
dHz = (Hz[i + xs, yf - j, k + zs] - Hz[i + xs, yf - j - 1, k + zs]) / dy dHz = (Hz[i + xs, yf - j, k + zs] - Hz[i + xs, yf - j - 1, k + zs]) / dy
Ex[i + xs, yf - j, k + zs] = Ex[i + xs, yf - j, k + zs] + updatecoeffsE[listIndex, 4] * ((RA[0, j] - 1) * dHz + RB[0, j] * EPhi[0, i, j, k]) Ex[i + xs, yf - j, k + zs] = Ex[i + xs, yf - j, k + zs] + updatecoeffsE[material, 4] * ((RA[0, j] - 1) * dHz + RB[0, j] * EPhi[0, i, j, k])
EPhi[0, i, j, k] = RE[0, j] * EPhi[0, i, j, k] - RF[0, j] * dHz EPhi[0, i, j, k] = RE[0, j] * EPhi[0, i, j, k] - RF[0, j] * dHz
cpdef update_pml_1order_ex_zplus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, :] updatecoeffsE, np.uint32_t[:, :, :, :] ID, floattype_t[:, :, :] Ex, floattype_t[:, :, :] Hy, floattype_t[:, :, :, :] EPhi, floattype_t[:, :] RA, floattype_t[:, :] RB, floattype_t[:, :] RE, floattype_t[:, :] RF, float dz): cpdef void update_pml_1order_ex_zplus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, ::1] updatecoeffsE, np.uint32_t[:, :, :, ::1] ID, floattype_t[:, :, ::1] Ex, floattype_t[:, :, ::1] Hy, floattype_t[:, :, :, ::1] EPhi, floattype_t[:, ::1] RA, floattype_t[:, ::1] RB, floattype_t[:, ::1] RE, floattype_t[:, ::1] RF, float dz):
"""This function updates the Ex field components in the z stretching direction. """This function updates the Ex field components in the z stretching direction.
Args: Args:
@@ -88,22 +90,23 @@ cpdef update_pml_1order_ex_zplus(int xs, int xf, int ys, int yf, int zs, int zf,
dx, dy, dz (float): Spatial discretisation dx, dy, dz (float): Spatial discretisation
""" """
cdef int i, j, k, nx, ny, nz, listIndex cdef Py_ssize_t i, j, k
cdef int nx, ny, nz, material
cdef float dHy cdef float dHy
nx = xf - xs nx = xf - xs
ny = yf - ys ny = yf - ys
nz = zf - zs nz = zf - zs
for i in prange(0, nx, nogil=True, schedule='static', chunksize=1, num_threads=nthreads): for i in prange(0, nx, nogil=True, schedule='static', num_threads=nthreads):
for j in range(0, ny): for j in range(0, ny):
for k in range(0, nz): for k in range(0, nz):
listIndex = ID[0, i + xs, j + ys, k + zs] material = ID[0, i + xs, j + ys, k + zs]
dHy = (Hy[i + xs, j + ys, k + zs] - Hy[i + xs, j + ys, k - 1 + zs]) / dz dHy = (Hy[i + xs, j + ys, k + zs] - Hy[i + xs, j + ys, k - 1 + zs]) / dz
Ex[i + xs, j + ys, k + zs] = Ex[i + xs, j + ys, k + zs] - updatecoeffsE[listIndex, 4] * ((RA[0, k] - 1) * dHy + RB[0, k] * EPhi[0, i, j, k]) Ex[i + xs, j + ys, k + zs] = Ex[i + xs, j + ys, k + zs] - updatecoeffsE[material, 4] * ((RA[0, k] - 1) * dHy + RB[0, k] * EPhi[0, i, j, k])
EPhi[0, i, j, k] = RE[0, k] * EPhi[0, i, j, k] - RF[0, k] * dHy EPhi[0, i, j, k] = RE[0, k] * EPhi[0, i, j, k] - RF[0, k] * dHy
cpdef update_pml_1order_ex_zminus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, :] updatecoeffsE, np.uint32_t[:, :, :, :] ID, floattype_t[:, :, :] Ex, floattype_t[:, :, :] Hy, floattype_t[:, :, :, :] EPhi, floattype_t[:, :] RA, floattype_t[:, :] RB, floattype_t[:, :] RE, floattype_t[:, :] RF, float dz): cpdef void update_pml_1order_ex_zminus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, ::1] updatecoeffsE, np.uint32_t[:, :, :, ::1] ID, floattype_t[:, :, ::1] Ex, floattype_t[:, :, ::1] Hy, floattype_t[:, :, :, ::1] EPhi, floattype_t[:, ::1] RA, floattype_t[:, ::1] RB, floattype_t[:, ::1] RE, floattype_t[:, ::1] RF, float dz):
"""This function updates the Ex field components in the z stretching direction. """This function updates the Ex field components in the z stretching direction.
Args: Args:
@@ -114,25 +117,26 @@ cpdef update_pml_1order_ex_zminus(int xs, int xf, int ys, int yf, int zs, int zf
dx, dy, dz (float): Spatial discretisation dx, dy, dz (float): Spatial discretisation
""" """
cdef int i, j, k, nx, ny, nz, listIndex cdef Py_ssize_t i, j, k
cdef int nx, ny, nz, material
cdef float dHy cdef float dHy
nx = xf - xs nx = xf - xs
ny = yf - ys ny = yf - ys
nz = zf - zs nz = zf - zs
for i in prange(0, nx, nogil=True, schedule='static', chunksize=1, num_threads=nthreads): for i in prange(0, nx, nogil=True, schedule='static', num_threads=nthreads):
for j in range(0, ny): for j in range(0, ny):
for k in range(0, nz): for k in range(0, nz):
listIndex = ID[0, i + xs, j + ys, zf - k] material = ID[0, i + xs, j + ys, zf - k]
dHy = (Hy[i + xs, j + ys, zf - k] - Hy[i + xs, j + ys, zf - k - 1]) / dz dHy = (Hy[i + xs, j + ys, zf - k] - Hy[i + xs, j + ys, zf - k - 1]) / dz
Ex[i + xs, j + ys, zf - k] = Ex[i + xs, j + ys, zf - k] - updatecoeffsE[listIndex, 4] * ((RA[0, k] - 1) * dHy + RB[0, k] * EPhi[0, i, j, k]) Ex[i + xs, j + ys, zf - k] = Ex[i + xs, j + ys, zf - k] - updatecoeffsE[material, 4] * ((RA[0, k] - 1) * dHy + RB[0, k] * EPhi[0, i, j, k])
EPhi[0, i, j, k] = RE[0, k] * EPhi[0, i, j, k] - RF[0, k] * dHy EPhi[0, i, j, k] = RE[0, k] * EPhi[0, i, j, k] - RF[0, k] * dHy
############################################# #############################################
# Electric field PML updates - Ey component # # Electric field PML updates - Ey component #
############################################# #############################################
cpdef update_pml_1order_ey_xplus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, :] updatecoeffsE, np.uint32_t[:, :, :, :] ID, floattype_t[:, :, :] Ey, floattype_t[:, :, :] Hz, floattype_t[:, :, :, :] EPhi, floattype_t[:, :] RA, floattype_t[:, :] RB, floattype_t[:, :] RE, floattype_t[:, :] RF, float dx): cpdef void update_pml_1order_ey_xplus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, ::1] updatecoeffsE, np.uint32_t[:, :, :, ::1] ID, floattype_t[:, :, ::1] Ey, floattype_t[:, :, ::1] Hz, floattype_t[:, :, :, ::1] EPhi, floattype_t[:, ::1] RA, floattype_t[:, ::1] RB, floattype_t[:, ::1] RE, floattype_t[:, ::1] RF, float dx):
"""This function updates the Ey field components in the x stretching direction. """This function updates the Ey field components in the x stretching direction.
Args: Args:
@@ -143,22 +147,23 @@ cpdef update_pml_1order_ey_xplus(int xs, int xf, int ys, int yf, int zs, int zf,
dx, dy, dz (float): Spatial discretisation dx, dy, dz (float): Spatial discretisation
""" """
cdef int i, j, k, nx, ny, nz, listIndex cdef Py_ssize_t i, j, k
cdef int nx, ny, nz, material
cdef float dHz cdef float dHz
nx = xf - xs nx = xf - xs
ny = yf - ys ny = yf - ys
nz = zf - zs nz = zf - zs
for i in prange(0, nx, nogil=True, schedule='static', chunksize=1, num_threads=nthreads): for i in prange(0, nx, nogil=True, schedule='static', num_threads=nthreads):
for j in range(0, ny): for j in range(0, ny):
for k in range(0, nz): for k in range(0, nz):
listIndex = ID[1, i + xs, j + ys, k + zs] material = ID[1, i + xs, j + ys, k + zs]
dHz = (Hz[i + xs, j + ys, k + zs] - Hz[i - 1 + xs, j + ys, k + zs]) / dx dHz = (Hz[i + xs, j + ys, k + zs] - Hz[i - 1 + xs, j + ys, k + zs]) / dx
Ey[i + xs, j + ys, k + zs] = Ey[i + xs, j + ys, k + zs] - updatecoeffsE[listIndex, 4] * ((RA[0, i] - 1) * dHz + RB[0, i] * EPhi[0, i, j, k]) Ey[i + xs, j + ys, k + zs] = Ey[i + xs, j + ys, k + zs] - updatecoeffsE[material, 4] * ((RA[0, i] - 1) * dHz + RB[0, i] * EPhi[0, i, j, k])
EPhi[0, i, j, k] = RE[0, i] * EPhi[0, i, j, k] - RF[0, i] * dHz EPhi[0, i, j, k] = RE[0, i] * EPhi[0, i, j, k] - RF[0, i] * dHz
cpdef update_pml_1order_ey_xminus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, :] updatecoeffsE, np.uint32_t[:, :, :, :] ID, floattype_t[:, :, :] Ey, floattype_t[:, :, :] Hz, floattype_t[:, :, :, :] EPhi, floattype_t[:, :] RA, floattype_t[:, :] RB, floattype_t[:, :] RE, floattype_t[:, :] RF, float dx): cpdef void update_pml_1order_ey_xminus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, ::1] updatecoeffsE, np.uint32_t[:, :, :, ::1] ID, floattype_t[:, :, ::1] Ey, floattype_t[:, :, ::1] Hz, floattype_t[:, :, :, ::1] EPhi, floattype_t[:, ::1] RA, floattype_t[:, ::1] RB, floattype_t[:, ::1] RE, floattype_t[:, ::1] RF, float dx):
"""This function updates the Ey field components in the x stretching direction. """This function updates the Ey field components in the x stretching direction.
Args: Args:
@@ -169,22 +174,23 @@ cpdef update_pml_1order_ey_xminus(int xs, int xf, int ys, int yf, int zs, int zf
dx, dy, dz (float): Spatial discretisation dx, dy, dz (float): Spatial discretisation
""" """
cdef int i, j, k, nx, ny, nz, listIndex cdef Py_ssize_t i, j, k
cdef int nx, ny, nz, material
cdef float dHz cdef float dHz
nx = xf - xs nx = xf - xs
ny = yf - ys ny = yf - ys
nz = zf - zs nz = zf - zs
for i in prange(0, nx, nogil=True, schedule='static', chunksize=1, num_threads=nthreads): for i in prange(0, nx, nogil=True, schedule='static', num_threads=nthreads):
for j in range(0, ny): for j in range(0, ny):
for k in range(0, nz): for k in range(0, nz):
listIndex = ID[1, xf - i, j + ys, k + zs] material = ID[1, xf - i, j + ys, k + zs]
dHz = (Hz[xf - i, j + ys, k + zs] - Hz[xf - i - 1, j + ys, k + zs]) / dx dHz = (Hz[xf - i, j + ys, k + zs] - Hz[xf - i - 1, j + ys, k + zs]) / dx
Ey[xf - i, j + ys, k + zs] = Ey[xf - i, j + ys, k + zs] - updatecoeffsE[listIndex, 4] * ((RA[0, i] - 1) * dHz + RB[0, i] * EPhi[0, i, j, k]) Ey[xf - i, j + ys, k + zs] = Ey[xf - i, j + ys, k + zs] - updatecoeffsE[material, 4] * ((RA[0, i] - 1) * dHz + RB[0, i] * EPhi[0, i, j, k])
EPhi[0, i, j, k] = RE[0, i] * EPhi[0, i, j, k] - RF[0, i] * dHz EPhi[0, i, j, k] = RE[0, i] * EPhi[0, i, j, k] - RF[0, i] * dHz
cpdef update_pml_1order_ey_zplus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, :] updatecoeffsE, np.uint32_t[:, :, :, :] ID, floattype_t[:, :, :] Ey, floattype_t[:, :, :] Hx, floattype_t[:, :, :, :] EPhi, floattype_t[:, :] RA, floattype_t[:, :] RB, floattype_t[:, :] RE, floattype_t[:, :] RF, float dz): cpdef void update_pml_1order_ey_zplus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, ::1] updatecoeffsE, np.uint32_t[:, :, :, ::1] ID, floattype_t[:, :, ::1] Ey, floattype_t[:, :, ::1] Hx, floattype_t[:, :, :, ::1] EPhi, floattype_t[:, ::1] RA, floattype_t[:, ::1] RB, floattype_t[:, ::1] RE, floattype_t[:, ::1] RF, float dz):
"""This function updates the Ey field components in the z stretching direction. """This function updates the Ey field components in the z stretching direction.
Args: Args:
@@ -195,22 +201,23 @@ cpdef update_pml_1order_ey_zplus(int xs, int xf, int ys, int yf, int zs, int zf,
dx, dy, dz (float): Spatial discretisation dx, dy, dz (float): Spatial discretisation
""" """
cdef int i, j, k, nx, ny, nz, listIndex cdef Py_ssize_t i, j, k
cdef int nx, ny, nz, material
cdef float dHx cdef float dHx
nx = xf - xs nx = xf - xs
ny = yf - ys ny = yf - ys
nz = zf - zs nz = zf - zs
for i in prange(0, nx, nogil=True, schedule='static', chunksize=1, num_threads=nthreads): for i in prange(0, nx, nogil=True, schedule='static', num_threads=nthreads):
for j in range(0, ny): for j in range(0, ny):
for k in range(0, nz): for k in range(0, nz):
listIndex = ID[1, i + xs, j + ys, k + zs] material = ID[1, i + xs, j + ys, k + zs]
dHx = (Hx[i + xs, j + ys, k + zs] - Hx[i + xs, j + ys, k - 1 + zs]) / dz dHx = (Hx[i + xs, j + ys, k + zs] - Hx[i + xs, j + ys, k - 1 + zs]) / dz
Ey[i + xs, j + ys, k + zs] = Ey[i + xs, j + ys, k + zs] + updatecoeffsE[listIndex, 4] * ((RA[0, k] - 1) * dHx + RB[0, k] * EPhi[0, i, j, k]) Ey[i + xs, j + ys, k + zs] = Ey[i + xs, j + ys, k + zs] + updatecoeffsE[material, 4] * ((RA[0, k] - 1) * dHx + RB[0, k] * EPhi[0, i, j, k])
EPhi[0, i, j, k] = RE[0, k] * EPhi[0, i, j, k] - RF[0, k] * dHx EPhi[0, i, j, k] = RE[0, k] * EPhi[0, i, j, k] - RF[0, k] * dHx
cpdef update_pml_1order_ey_zminus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, :] updatecoeffsE, np.uint32_t[:, :, :, :] ID, floattype_t[:, :, :] Ey, floattype_t[:, :, :] Hx, floattype_t[:, :, :, :] EPhi, floattype_t[:, :] RA, floattype_t[:, :] RB, floattype_t[:, :] RE, floattype_t[:, :] RF, float dz): cpdef void update_pml_1order_ey_zminus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, ::1] updatecoeffsE, np.uint32_t[:, :, :, ::1] ID, floattype_t[:, :, ::1] Ey, floattype_t[:, :, ::1] Hx, floattype_t[:, :, :, ::1] EPhi, floattype_t[:, ::1] RA, floattype_t[:, ::1] RB, floattype_t[:, ::1] RE, floattype_t[:, ::1] RF, float dz):
"""This function updates the Ey field components in the z stretching direction. """This function updates the Ey field components in the z stretching direction.
Args: Args:
@@ -221,25 +228,26 @@ cpdef update_pml_1order_ey_zminus(int xs, int xf, int ys, int yf, int zs, int zf
dx, dy, dz (float): Spatial discretisation dx, dy, dz (float): Spatial discretisation
""" """
cdef int i, j, k, nx, ny, nz, listIndex cdef Py_ssize_t i, j, k
cdef int nx, ny, nz, material
cdef float dHx cdef float dHx
nx = xf - xs nx = xf - xs
ny = yf - ys ny = yf - ys
nz = zf - zs nz = zf - zs
for i in prange(0, nx, nogil=True, schedule='static', chunksize=1, num_threads=nthreads): for i in prange(0, nx, nogil=True, schedule='static', num_threads=nthreads):
for j in range(0, ny): for j in range(0, ny):
for k in range(0, nz): for k in range(0, nz):
listIndex = ID[1, i + xs, j + ys, zf - k] material = ID[1, i + xs, j + ys, zf - k]
dHx = (Hx[i + xs, j + ys, zf - k] - Hx[i + xs, j + ys, zf - k - 1]) / dz dHx = (Hx[i + xs, j + ys, zf - k] - Hx[i + xs, j + ys, zf - k - 1]) / dz
Ey[i + xs, j + ys, zf - k] = Ey[i + xs, j + ys, zf - k] + updatecoeffsE[listIndex, 4] * ((RA[0, k] - 1) * dHx + RB[0, k] * EPhi[0, i, j, k]) Ey[i + xs, j + ys, zf - k] = Ey[i + xs, j + ys, zf - k] + updatecoeffsE[material, 4] * ((RA[0, k] - 1) * dHx + RB[0, k] * EPhi[0, i, j, k])
EPhi[0, i, j, k] = RE[0, k] * EPhi[0, i, j, k] - RF[0, k] * dHx EPhi[0, i, j, k] = RE[0, k] * EPhi[0, i, j, k] - RF[0, k] * dHx
############################################# #############################################
# Electric field PML updates - Ez component # # Electric field PML updates - Ez component #
############################################# #############################################
cpdef update_pml_1order_ez_xplus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, :] updatecoeffsE, np.uint32_t[:, :, :, :] ID, floattype_t[:, :, :] Ez, floattype_t[:, :, :] Hy, floattype_t[:, :, :, :] EPhi, floattype_t[:, :] RA, floattype_t[:, :] RB, floattype_t[:, :] RE, floattype_t[:, :] RF, float dx): cpdef void update_pml_1order_ez_xplus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, ::1] updatecoeffsE, np.uint32_t[:, :, :, ::1] ID, floattype_t[:, :, ::1] Ez, floattype_t[:, :, ::1] Hy, floattype_t[:, :, :, ::1] EPhi, floattype_t[:, ::1] RA, floattype_t[:, ::1] RB, floattype_t[:, ::1] RE, floattype_t[:, ::1] RF, float dx):
"""This function updates the Ey field components in the z stretching direction. """This function updates the Ey field components in the z stretching direction.
Args: Args:
@@ -250,22 +258,23 @@ cpdef update_pml_1order_ez_xplus(int xs, int xf, int ys, int yf, int zs, int zf,
dx, dy, dz (float): Spatial discretisation dx, dy, dz (float): Spatial discretisation
""" """
cdef int i, j, k, nx, ny, nz, listIndex cdef Py_ssize_t i, j, k
cdef int nx, ny, nz, material
cdef float dHy cdef float dHy
nx = xf - xs nx = xf - xs
ny = yf - ys ny = yf - ys
nz = zf - zs nz = zf - zs
for i in prange(0, nx, nogil=True, schedule='static', chunksize=1, num_threads=nthreads): for i in prange(0, nx, nogil=True, schedule='static', num_threads=nthreads):
for j in range(0, ny): for j in range(0, ny):
for k in range(0, nz): for k in range(0, nz):
listIndex = ID[2, i + xs, j + ys, k + zs] material = ID[2, i + xs, j + ys, k + zs]
dHy = (Hy[i + xs, j + ys, k + zs] - Hy[i - 1 + xs, j + ys, k + zs]) / dx dHy = (Hy[i + xs, j + ys, k + zs] - Hy[i - 1 + xs, j + ys, k + zs]) / dx
Ez[i + xs, j + ys, k + zs] = Ez[i + xs, j + ys, k + zs] + updatecoeffsE[listIndex, 4] * ((RA[0, i] - 1) * dHy + RB[0, i] * EPhi[0, i, j, k]) Ez[i + xs, j + ys, k + zs] = Ez[i + xs, j + ys, k + zs] + updatecoeffsE[material, 4] * ((RA[0, i] - 1) * dHy + RB[0, i] * EPhi[0, i, j, k])
EPhi[0, i, j, k] = RE[0, i] * EPhi[0, i, j, k] - RF[0, i] * dHy EPhi[0, i, j, k] = RE[0, i] * EPhi[0, i, j, k] - RF[0, i] * dHy
cpdef update_pml_1order_ez_xminus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, :] updatecoeffsE, np.uint32_t[:, :, :, :] ID, floattype_t[:, :, :] Ez, floattype_t[:, :, :] Hy, floattype_t[:, :, :, :] EPhi, floattype_t[:, :] RA, floattype_t[:, :] RB, floattype_t[:, :] RE, floattype_t[:, :] RF, float dx): cpdef void update_pml_1order_ez_xminus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, ::1] updatecoeffsE, np.uint32_t[:, :, :, ::1] ID, floattype_t[:, :, ::1] Ez, floattype_t[:, :, ::1] Hy, floattype_t[:, :, :, ::1] EPhi, floattype_t[:, ::1] RA, floattype_t[:, ::1] RB, floattype_t[:, ::1] RE, floattype_t[:, ::1] RF, float dx):
"""This function updates the Ez field components in the x stretching direction. """This function updates the Ez field components in the x stretching direction.
Args: Args:
@@ -276,22 +285,23 @@ cpdef update_pml_1order_ez_xminus(int xs, int xf, int ys, int yf, int zs, int zf
dx, dy, dz (float): Spatial discretisation dx, dy, dz (float): Spatial discretisation
""" """
cdef int i, j, k, nx, ny, nz, listIndex cdef Py_ssize_t i, j, k
cdef int nx, ny, nz, material
cdef float dHy cdef float dHy
nx = xf - xs nx = xf - xs
ny = yf - ys ny = yf - ys
nz = zf - zs nz = zf - zs
for i in prange(0, nx, nogil=True, schedule='static', chunksize=1, num_threads=nthreads): for i in prange(0, nx, nogil=True, schedule='static', num_threads=nthreads):
for j in range(0, ny): for j in range(0, ny):
for k in range(0, nz): for k in range(0, nz):
listIndex = ID[2, xf - i, j + ys, k + zs] material = ID[2, xf - i, j + ys, k + zs]
dHy = (Hy[xf - i, j + ys, k + zs] - Hy[xf - i - 1, j + ys, k + zs]) / dx dHy = (Hy[xf - i, j + ys, k + zs] - Hy[xf - i - 1, j + ys, k + zs]) / dx
Ez[xf - i, j + ys, k + zs] = Ez[xf - i, j + ys, k + zs] + updatecoeffsE[listIndex, 4] * ((RA[0, i] - 1) * dHy + RB[0, i] * EPhi[0, i, j, k]) Ez[xf - i, j + ys, k + zs] = Ez[xf - i, j + ys, k + zs] + updatecoeffsE[material, 4] * ((RA[0, i] - 1) * dHy + RB[0, i] * EPhi[0, i, j, k])
EPhi[0, i, j, k] = RE[0, i] * EPhi[0, i, j, k] - RF[0, i] * dHy EPhi[0, i, j, k] = RE[0, i] * EPhi[0, i, j, k] - RF[0, i] * dHy
cpdef update_pml_1order_ez_yplus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, :] updatecoeffsE, np.uint32_t[:, :, :, :] ID, floattype_t[:, :, :] Ez, floattype_t[:, :, :] Hx, floattype_t[:, :, :, :] EPhi, floattype_t[:, :] RA, floattype_t[:, :] RB, floattype_t[:, :] RE, floattype_t[:, :] RF, float dy): cpdef void update_pml_1order_ez_yplus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, ::1] updatecoeffsE, np.uint32_t[:, :, :, ::1] ID, floattype_t[:, :, ::1] Ez, floattype_t[:, :, ::1] Hx, floattype_t[:, :, :, ::1] EPhi, floattype_t[:, ::1] RA, floattype_t[:, ::1] RB, floattype_t[:, ::1] RE, floattype_t[:, ::1] RF, float dy):
"""This function updates the Ez field components in the y stretching direction. """This function updates the Ez field components in the y stretching direction.
Args: Args:
@@ -302,22 +312,23 @@ cpdef update_pml_1order_ez_yplus(int xs, int xf, int ys, int yf, int zs, int zf,
dx, dy, dz (float): Spatial discretisation dx, dy, dz (float): Spatial discretisation
""" """
cdef int i, j, k, nx, ny, nz, listIndex cdef Py_ssize_t i, j, k
cdef int nx, ny, nz, material
cdef float dHx cdef float dHx
nx = xf - xs nx = xf - xs
ny = yf - ys ny = yf - ys
nz = zf - zs nz = zf - zs
for i in prange(0, nx, nogil=True, schedule='static', chunksize=1, num_threads=nthreads): for i in prange(0, nx, nogil=True, schedule='static', num_threads=nthreads):
for j in range(0, ny): for j in range(0, ny):
for k in range(0, nz): for k in range(0, nz):
listIndex = ID[2, i + xs, j + ys, k + zs] material = ID[2, i + xs, j + ys, k + zs]
dHx = (Hx[i + xs, j + ys, k + zs] - Hx[i + xs, j - 1 + ys, k + zs]) / dy dHx = (Hx[i + xs, j + ys, k + zs] - Hx[i + xs, j - 1 + ys, k + zs]) / dy
Ez[i + xs, j + ys, k + zs] = Ez[i + xs, j + ys, k + zs] - updatecoeffsE[listIndex, 4] * ((RA[0, j] - 1) * dHx + RB[0, j] * EPhi[0, i, j, k]) Ez[i + xs, j + ys, k + zs] = Ez[i + xs, j + ys, k + zs] - updatecoeffsE[material, 4] * ((RA[0, j] - 1) * dHx + RB[0, j] * EPhi[0, i, j, k])
EPhi[0, i, j, k] = RE[0, j] * EPhi[0, i, j, k] - RF[0, j] * dHx EPhi[0, i, j, k] = RE[0, j] * EPhi[0, i, j, k] - RF[0, j] * dHx
cpdef update_pml_1order_ez_yminus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, :] updatecoeffsE, np.uint32_t[:, :, :, :] ID, floattype_t[:, :, :] Ez, floattype_t[:, :, :] Hx, floattype_t[:, :, :, :] EPhi, floattype_t[:, :] RA, floattype_t[:, :] RB, floattype_t[:, :] RE, floattype_t[:, :] RF, float dy): cpdef void update_pml_1order_ez_yminus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, ::1] updatecoeffsE, np.uint32_t[:, :, :, ::1] ID, floattype_t[:, :, ::1] Ez, floattype_t[:, :, ::1] Hx, floattype_t[:, :, :, ::1] EPhi, floattype_t[:, ::1] RA, floattype_t[:, ::1] RB, floattype_t[:, ::1] RE, floattype_t[:, ::1] RF, float dy):
"""This function updates the Ez field components in the y stretching direction. """This function updates the Ez field components in the y stretching direction.
Args: Args:
@@ -328,25 +339,26 @@ cpdef update_pml_1order_ez_yminus(int xs, int xf, int ys, int yf, int zs, int zf
dx, dy, dz (float): Spatial discretisation dx, dy, dz (float): Spatial discretisation
""" """
cdef int i, j, k, nx, ny, nz, listIndex cdef Py_ssize_t i, j, k
cdef int nx, ny, nz, material
cdef float dHx cdef float dHx
nx = xf - xs nx = xf - xs
ny = yf - ys ny = yf - ys
nz = zf - zs nz = zf - zs
for i in prange(0, nx, nogil=True, schedule='static', chunksize=1, num_threads=nthreads): for i in prange(0, nx, nogil=True, schedule='static', num_threads=nthreads):
for j in range(0, ny): for j in range(0, ny):
for k in range(0, nz): for k in range(0, nz):
listIndex = ID[2, i + xs, yf - j, k + zs] material = ID[2, i + xs, yf - j, k + zs]
dHx = (Hx[i + xs, yf - j, k + zs] - Hx[i + xs, yf - j - 1, k + zs]) / dy dHx = (Hx[i + xs, yf - j, k + zs] - Hx[i + xs, yf - j - 1, k + zs]) / dy
Ez[i + xs, yf - j, k + zs] = Ez[i + xs, yf - j, k + zs] - updatecoeffsE[listIndex, 4] * ((RA[0, j] - 1) * dHx + RB[0, j] * EPhi[0, i, j, k]) Ez[i + xs, yf - j, k + zs] = Ez[i + xs, yf - j, k + zs] - updatecoeffsE[material, 4] * ((RA[0, j] - 1) * dHx + RB[0, j] * EPhi[0, i, j, k])
EPhi[0, i, j, k] = RE[0, j] * EPhi[0, i, j, k] - RF[0, j] * dHx EPhi[0, i, j, k] = RE[0, j] * EPhi[0, i, j, k] - RF[0, j] * dHx
############################################# #############################################
# Magnetic field PML updates - Hx component # # Magnetic field PML updates - Hx component #
############################################# #############################################
cpdef update_pml_1order_hx_yplus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, :] updatecoeffsH, np.uint32_t[:, :, :, :] ID, floattype_t[:, :, :] Hx, floattype_t[:, :, :] Ez, floattype_t[:, :, :, :] HPhi, floattype_t[:, :] RA, floattype_t[:, :] RB, floattype_t[:, :] RE, floattype_t[:, :] RF, float dy): cpdef void update_pml_1order_hx_yplus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, ::1] updatecoeffsH, np.uint32_t[:, :, :, ::1] ID, floattype_t[:, :, ::1] Hx, floattype_t[:, :, ::1] Ez, floattype_t[:, :, :, ::1] HPhi, floattype_t[:, ::1] RA, floattype_t[:, ::1] RB, floattype_t[:, ::1] RE, floattype_t[:, ::1] RF, float dy):
"""This function updates the Hx field components in the y stretching direction. """This function updates the Hx field components in the y stretching direction.
Args: Args:
@@ -357,22 +369,23 @@ cpdef update_pml_1order_hx_yplus(int xs, int xf, int ys, int yf, int zs, int zf,
dx, dy, dz (float): Spatial discretisation dx, dy, dz (float): Spatial discretisation
""" """
cdef int i, j, k, nx, ny, nz, listIndex cdef Py_ssize_t i, j, k
cdef int nx, ny, nz, material
cdef float dEz cdef float dEz
nx = xf - xs nx = xf - xs
ny = yf - ys ny = yf - ys
nz = zf - zs nz = zf - zs
for i in prange(0, nx, nogil=True, schedule='static', chunksize=1, num_threads=nthreads): for i in prange(0, nx, nogil=True, schedule='static', num_threads=nthreads):
for j in range(0, ny): for j in range(0, ny):
for k in range(0, nz): for k in range(0, nz):
listIndex = ID[3, i + xs, j + ys, k + zs] material = ID[3, i + xs, j + ys, k + zs]
dEz = (Ez[i + xs, j + 1 + ys, k + zs] - Ez[i + xs, j + ys, k + zs]) / dy dEz = (Ez[i + xs, j + 1 + ys, k + zs] - Ez[i + xs, j + ys, k + zs]) / dy
Hx[i + xs, j + ys, k + zs] = Hx[i + xs, j + ys, k + zs] - updatecoeffsH[listIndex, 4] * ((RA[0, j] - 1) * dEz + RB[0, j] * HPhi[0, i, j, k]) Hx[i + xs, j + ys, k + zs] = Hx[i + xs, j + ys, k + zs] - updatecoeffsH[material, 4] * ((RA[0, j] - 1) * dEz + RB[0, j] * HPhi[0, i, j, k])
HPhi[0, i, j, k] = RE[0, j] * HPhi[0, i, j, k] - RF[0, j] * dEz HPhi[0, i, j, k] = RE[0, j] * HPhi[0, i, j, k] - RF[0, j] * dEz
cpdef update_pml_1order_hx_yminus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, :] updatecoeffsH, np.uint32_t[:, :, :, :] ID, floattype_t[:, :, :] Hx, floattype_t[:, :, :] Ez, floattype_t[:, :, :, :] HPhi, floattype_t[:, :] RA, floattype_t[:, :] RB, floattype_t[:, :] RE, floattype_t[:, :] RF, float dy): cpdef void update_pml_1order_hx_yminus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, ::1] updatecoeffsH, np.uint32_t[:, :, :, ::1] ID, floattype_t[:, :, ::1] Hx, floattype_t[:, :, ::1] Ez, floattype_t[:, :, :, ::1] HPhi, floattype_t[:, ::1] RA, floattype_t[:, ::1] RB, floattype_t[:, ::1] RE, floattype_t[:, ::1] RF, float dy):
"""This function updates the Hx field components in the y stretching direction. """This function updates the Hx field components in the y stretching direction.
Args: Args:
@@ -383,22 +396,23 @@ cpdef update_pml_1order_hx_yminus(int xs, int xf, int ys, int yf, int zs, int zf
dx, dy, dz (float): Spatial discretisation dx, dy, dz (float): Spatial discretisation
""" """
cdef int i, j, k, nx, ny, nz, listIndex cdef Py_ssize_t i, j, k
cdef int nx, ny, nz, material
cdef float dEz cdef float dEz
nx = xf - xs nx = xf - xs
ny = yf - ys ny = yf - ys
nz = zf - zs nz = zf - zs
for i in prange(0, nx, nogil=True, schedule='static', chunksize=1, num_threads=nthreads): for i in prange(0, nx, nogil=True, schedule='static', num_threads=nthreads):
for j in range(0, ny): for j in range(0, ny):
for k in range(0, nz): for k in range(0, nz):
listIndex = ID[3, i + xs, yf - (j + 1), k + zs] material = ID[3, i + xs, yf - (j + 1), k + zs]
dEz = (Ez[i + xs, yf - j, k + zs] - Ez[i + xs, yf - (j + 1), k + zs]) / dy dEz = (Ez[i + xs, yf - j, k + zs] - Ez[i + xs, yf - (j + 1), k + zs]) / dy
Hx[i + xs, yf - (j + 1), k + zs] = Hx[i + xs, yf - (j + 1), k + zs] - updatecoeffsH[listIndex, 4] * ((RA[0, j] - 1) * dEz + RB[0, j] * HPhi[0, i, j, k]) Hx[i + xs, yf - (j + 1), k + zs] = Hx[i + xs, yf - (j + 1), k + zs] - updatecoeffsH[material, 4] * ((RA[0, j] - 1) * dEz + RB[0, j] * HPhi[0, i, j, k])
HPhi[0, i, j, k] = RE[0, j] * HPhi[0, i, j, k] - RF[0, j] * dEz HPhi[0, i, j, k] = RE[0, j] * HPhi[0, i, j, k] - RF[0, j] * dEz
cpdef update_pml_1order_hx_zplus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, :] updatecoeffsH, np.uint32_t[:, :, :, :] ID, floattype_t[:, :, :] Hx, floattype_t[:, :, :] Ey, floattype_t[:, :, :, :] HPhi, floattype_t[:, :] RA, floattype_t[:, :] RB, floattype_t[:, :] RE, floattype_t[:, :] RF, float dz): cpdef void update_pml_1order_hx_zplus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, ::1] updatecoeffsH, np.uint32_t[:, :, :, ::1] ID, floattype_t[:, :, ::1] Hx, floattype_t[:, :, ::1] Ey, floattype_t[:, :, :, ::1] HPhi, floattype_t[:, ::1] RA, floattype_t[:, ::1] RB, floattype_t[:, ::1] RE, floattype_t[:, ::1] RF, float dz):
"""This function updates the Hx field components in the z stretching direction. """This function updates the Hx field components in the z stretching direction.
Args: Args:
@@ -409,22 +423,23 @@ cpdef update_pml_1order_hx_zplus(int xs, int xf, int ys, int yf, int zs, int zf,
dx, dy, dz (float): Spatial discretisation dx, dy, dz (float): Spatial discretisation
""" """
cdef int i, j, k, nx, ny, nz, listIndex cdef Py_ssize_t i, j, k
cdef int nx, ny, nz, material
cdef float dEy cdef float dEy
nx = xf - xs nx = xf - xs
ny = yf - ys ny = yf - ys
nz = zf - zs nz = zf - zs
for i in prange(0, nx, nogil=True, schedule='static', chunksize=1, num_threads=nthreads): for i in prange(0, nx, nogil=True, schedule='static', num_threads=nthreads):
for j in range(0, ny): for j in range(0, ny):
for k in range(0, nz): for k in range(0, nz):
listIndex = ID[3, i + xs, j + ys, k + zs] material = ID[3, i + xs, j + ys, k + zs]
dEy = (Ey[i + xs, j + ys, k + 1 + zs] - Ey[i + xs, j + ys, k + zs]) / dz dEy = (Ey[i + xs, j + ys, k + 1 + zs] - Ey[i + xs, j + ys, k + zs]) / dz
Hx[i + xs, j + ys, k + zs] = Hx[i + xs, j + ys, k + zs] + updatecoeffsH[listIndex, 4] * ((RA[0, k] - 1) * dEy + RB[0, k] * HPhi[0, i, j, k]) Hx[i + xs, j + ys, k + zs] = Hx[i + xs, j + ys, k + zs] + updatecoeffsH[material, 4] * ((RA[0, k] - 1) * dEy + RB[0, k] * HPhi[0, i, j, k])
HPhi[0, i, j, k] = RE[0, k] * HPhi[0, i, j, k] - RF[0, k] * dEy HPhi[0, i, j, k] = RE[0, k] * HPhi[0, i, j, k] - RF[0, k] * dEy
cpdef update_pml_1order_hx_zminus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, :] updatecoeffsH, np.uint32_t[:, :, :, :] ID, floattype_t[:, :, :] Hx, floattype_t[:, :, :] Ey, floattype_t[:, :, :, :] HPhi, floattype_t[:, :] RA, floattype_t[:, :] RB, floattype_t[:, :] RE, floattype_t[:, :] RF, float dz): cpdef void update_pml_1order_hx_zminus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, ::1] updatecoeffsH, np.uint32_t[:, :, :, ::1] ID, floattype_t[:, :, ::1] Hx, floattype_t[:, :, ::1] Ey, floattype_t[:, :, :, ::1] HPhi, floattype_t[:, ::1] RA, floattype_t[:, ::1] RB, floattype_t[:, ::1] RE, floattype_t[:, ::1] RF, float dz):
"""This function updates the Hx field components in the z stretching direction. """This function updates the Hx field components in the z stretching direction.
Args: Args:
@@ -435,25 +450,26 @@ cpdef update_pml_1order_hx_zminus(int xs, int xf, int ys, int yf, int zs, int zf
dx, dy, dz (float): Spatial discretisation dx, dy, dz (float): Spatial discretisation
""" """
cdef int i, j, k, nx, ny, nz, listIndex cdef Py_ssize_t i, j, k
cdef int nx, ny, nz, material
cdef float dEy cdef float dEy
nx = xf - xs nx = xf - xs
ny = yf - ys ny = yf - ys
nz = zf - zs nz = zf - zs
for i in prange(0, nx, nogil=True, schedule='static', chunksize=1, num_threads=nthreads): for i in prange(0, nx, nogil=True, schedule='static', num_threads=nthreads):
for j in range(0, ny): for j in range(0, ny):
for k in range(0, nz): for k in range(0, nz):
listIndex = ID[3, i + xs, j + ys, zf - (k + 1)] material = ID[3, i + xs, j + ys, zf - (k + 1)]
dEy = (Ey[i + xs, j + ys, zf - k] - Ey[i + xs, j + ys, zf - (k + 1)]) / dz dEy = (Ey[i + xs, j + ys, zf - k] - Ey[i + xs, j + ys, zf - (k + 1)]) / dz
Hx[i + xs, j + ys, zf - (k + 1)] = Hx[i + xs, j + ys, zf - (k + 1)] + updatecoeffsH[listIndex, 4] * ((RA[0, k] - 1) * dEy + RB[0, k] * HPhi[0, i, j, k]) Hx[i + xs, j + ys, zf - (k + 1)] = Hx[i + xs, j + ys, zf - (k + 1)] + updatecoeffsH[material, 4] * ((RA[0, k] - 1) * dEy + RB[0, k] * HPhi[0, i, j, k])
HPhi[0, i, j, k] = RE[0, k] * HPhi[0, i, j, k] - RF[0, k] * dEy HPhi[0, i, j, k] = RE[0, k] * HPhi[0, i, j, k] - RF[0, k] * dEy
############################################# #############################################
# Magnetic field PML updates - Hy component # # Magnetic field PML updates - Hy component #
############################################# #############################################
cpdef update_pml_1order_hy_xplus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, :] updatecoeffsH, np.uint32_t[:, :, :, :] ID, floattype_t[:, :, :] Hy, floattype_t[:, :, :] Ez, floattype_t[:, :, :, :] HPhi, floattype_t[:, :] RA, floattype_t[:, :] RB, floattype_t[:, :] RE, floattype_t[:, :] RF, float dx): cpdef void update_pml_1order_hy_xplus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, ::1] updatecoeffsH, np.uint32_t[:, :, :, ::1] ID, floattype_t[:, :, ::1] Hy, floattype_t[:, :, ::1] Ez, floattype_t[:, :, :, ::1] HPhi, floattype_t[:, ::1] RA, floattype_t[:, ::1] RB, floattype_t[:, ::1] RE, floattype_t[:, ::1] RF, float dx):
"""This function updates the Hy field components in the x stretching direction. """This function updates the Hy field components in the x stretching direction.
Args: Args:
@@ -464,22 +480,23 @@ cpdef update_pml_1order_hy_xplus(int xs, int xf, int ys, int yf, int zs, int zf,
dx, dy, dz (float): Spatial discretisation dx, dy, dz (float): Spatial discretisation
""" """
cdef int i, j, k, nx, ny, nz, listIndex cdef Py_ssize_t i, j, k
cdef int nx, ny, nz, material
cdef float dEz cdef float dEz
nx = xf - xs nx = xf - xs
ny = yf - ys ny = yf - ys
nz = zf - zs nz = zf - zs
for i in prange(0, nx, nogil=True, schedule='static', chunksize=1, num_threads=nthreads): for i in prange(0, nx, nogil=True, schedule='static', num_threads=nthreads):
for j in range(0, ny): for j in range(0, ny):
for k in range(0, nz): for k in range(0, nz):
listIndex = ID[4, i + xs, j + ys, k + zs] material = ID[4, i + xs, j + ys, k + zs]
dEz = (Ez[i + 1 + xs, j + ys, k + zs] - Ez[i + xs, j + ys, k + zs]) / dx dEz = (Ez[i + 1 + xs, j + ys, k + zs] - Ez[i + xs, j + ys, k + zs]) / dx
Hy[i + xs, j + ys, k + zs] = Hy[i + xs, j + ys, k + zs] + updatecoeffsH[listIndex, 4] * ((RA[0, i] - 1) * dEz + RB[0, i] * HPhi[0, i, j, k]) Hy[i + xs, j + ys, k + zs] = Hy[i + xs, j + ys, k + zs] + updatecoeffsH[material, 4] * ((RA[0, i] - 1) * dEz + RB[0, i] * HPhi[0, i, j, k])
HPhi[0, i, j, k] = RE[0, i] * HPhi[0, i, j, k] - RF[0, i] * dEz HPhi[0, i, j, k] = RE[0, i] * HPhi[0, i, j, k] - RF[0, i] * dEz
cpdef update_pml_1order_hy_xminus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, :] updatecoeffsH, np.uint32_t[:, :, :, :] ID, floattype_t[:, :, :] Hy, floattype_t[:, :, :] Ez, floattype_t[:, :, :, :] HPhi, floattype_t[:, :] RA, floattype_t[:, :] RB, floattype_t[:, :] RE, floattype_t[:, :] RF, float dx): cpdef void update_pml_1order_hy_xminus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, ::1] updatecoeffsH, np.uint32_t[:, :, :, ::1] ID, floattype_t[:, :, ::1] Hy, floattype_t[:, :, ::1] Ez, floattype_t[:, :, :, ::1] HPhi, floattype_t[:, ::1] RA, floattype_t[:, ::1] RB, floattype_t[:, ::1] RE, floattype_t[:, ::1] RF, float dx):
"""This function updates the Hy field components in the x stretching direction. """This function updates the Hy field components in the x stretching direction.
Args: Args:
@@ -490,22 +507,23 @@ cpdef update_pml_1order_hy_xminus(int xs, int xf, int ys, int yf, int zs, int zf
dx, dy, dz (float): Spatial discretisation dx, dy, dz (float): Spatial discretisation
""" """
cdef int i, j, k, nx, ny, nz, listIndex cdef Py_ssize_t i, j, k
cdef int nx, ny, nz, material
cdef float dEz cdef float dEz
nx = xf - xs nx = xf - xs
ny = yf - ys ny = yf - ys
nz = zf - zs nz = zf - zs
for i in prange(0, nx, nogil=True, schedule='static', chunksize=1, num_threads=nthreads): for i in prange(0, nx, nogil=True, schedule='static', num_threads=nthreads):
for j in range(0, ny): for j in range(0, ny):
for k in range(0, nz): for k in range(0, nz):
listIndex = ID[4, xf - (i + 1), j + ys, k + zs] material = ID[4, xf - (i + 1), j + ys, k + zs]
dEz = (Ez[xf - i, j + ys, k + zs] - Ez[xf - (i + 1), j + ys, k + zs]) / dx dEz = (Ez[xf - i, j + ys, k + zs] - Ez[xf - (i + 1), j + ys, k + zs]) / dx
Hy[xf - (i + 1), j + ys, k + zs] = Hy[xf - (i + 1), j + ys, k + zs] + updatecoeffsH[listIndex, 4] * ((RA[0, i] - 1) * dEz + RB[0, i] * HPhi[0, i, j, k]) Hy[xf - (i + 1), j + ys, k + zs] = Hy[xf - (i + 1), j + ys, k + zs] + updatecoeffsH[material, 4] * ((RA[0, i] - 1) * dEz + RB[0, i] * HPhi[0, i, j, k])
HPhi[0, i, j, k] = RE[0, i] * HPhi[0, i, j, k] - RF[0, i] * dEz HPhi[0, i, j, k] = RE[0, i] * HPhi[0, i, j, k] - RF[0, i] * dEz
cpdef update_pml_1order_hy_zplus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, :] updatecoeffsH, np.uint32_t[:, :, :, :] ID, floattype_t[:, :, :] Hy, floattype_t[:, :, :] Ex, floattype_t[:, :, :, :] HPhi, floattype_t[:, :] RA, floattype_t[:, :] RB, floattype_t[:, :] RE, floattype_t[:, :] RF, float dz): cpdef void update_pml_1order_hy_zplus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, ::1] updatecoeffsH, np.uint32_t[:, :, :, ::1] ID, floattype_t[:, :, ::1] Hy, floattype_t[:, :, ::1] Ex, floattype_t[:, :, :, ::1] HPhi, floattype_t[:, ::1] RA, floattype_t[:, ::1] RB, floattype_t[:, ::1] RE, floattype_t[:, ::1] RF, float dz):
"""This function updates the Hy field components in the z stretching direction. """This function updates the Hy field components in the z stretching direction.
Args: Args:
@@ -516,22 +534,23 @@ cpdef update_pml_1order_hy_zplus(int xs, int xf, int ys, int yf, int zs, int zf,
dx, dy, dz (float): Spatial discretisation dx, dy, dz (float): Spatial discretisation
""" """
cdef int i, j, k, nx, ny, nz, listIndex cdef Py_ssize_t i, j, k
cdef int nx, ny, nz, material
cdef float dEx cdef float dEx
nx = xf - xs nx = xf - xs
ny = yf - ys ny = yf - ys
nz = zf - zs nz = zf - zs
for i in prange(0, nx, nogil=True, schedule='static', chunksize=1, num_threads=nthreads): for i in prange(0, nx, nogil=True, schedule='static', num_threads=nthreads):
for j in range(0, ny): for j in range(0, ny):
for k in range(0, nz): for k in range(0, nz):
listIndex = ID[4, i + xs, j + ys, k + zs] material = ID[4, i + xs, j + ys, k + zs]
dEx = (Ex[i + xs, j + ys, k + 1 + zs] - Ex[i + xs, j + ys, k + zs]) / dz dEx = (Ex[i + xs, j + ys, k + 1 + zs] - Ex[i + xs, j + ys, k + zs]) / dz
Hy[i + xs, j + ys, k + zs] = Hy[i + xs, j + ys, k + zs] - updatecoeffsH[listIndex, 4] * ((RA[0, k] - 1) * dEx + RB[0, k] * HPhi[0, i, j, k]) Hy[i + xs, j + ys, k + zs] = Hy[i + xs, j + ys, k + zs] - updatecoeffsH[material, 4] * ((RA[0, k] - 1) * dEx + RB[0, k] * HPhi[0, i, j, k])
HPhi[0, i, j, k] = RE[0, k] * HPhi[0, i, j, k] - RF[0, k] * dEx HPhi[0, i, j, k] = RE[0, k] * HPhi[0, i, j, k] - RF[0, k] * dEx
cpdef update_pml_1order_hy_zminus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, :] updatecoeffsH, np.uint32_t[:, :, :, :] ID, floattype_t[:, :, :] Hy, floattype_t[:, :, :] Ex, floattype_t[:, :, :, :] HPhi, floattype_t[:, :] RA, floattype_t[:, :] RB, floattype_t[:, :] RE, floattype_t[:, :] RF, float dz): cpdef void update_pml_1order_hy_zminus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, ::1] updatecoeffsH, np.uint32_t[:, :, :, ::1] ID, floattype_t[:, :, ::1] Hy, floattype_t[:, :, ::1] Ex, floattype_t[:, :, :, ::1] HPhi, floattype_t[:, ::1] RA, floattype_t[:, ::1] RB, floattype_t[:, ::1] RE, floattype_t[:, ::1] RF, float dz):
"""This function updates the Hy field components in the z stretching direction. """This function updates the Hy field components in the z stretching direction.
Args: Args:
@@ -542,25 +561,26 @@ cpdef update_pml_1order_hy_zminus(int xs, int xf, int ys, int yf, int zs, int zf
dx, dy, dz (float): Spatial discretisation dx, dy, dz (float): Spatial discretisation
""" """
cdef int i, j, k, nx, ny, nz, listIndex cdef Py_ssize_t i, j, k
cdef int nx, ny, nz, material
cdef float dEx cdef float dEx
nx = xf - xs nx = xf - xs
ny = yf - ys ny = yf - ys
nz = zf - zs nz = zf - zs
for i in prange(0, nx, nogil=True, schedule='static', chunksize=1, num_threads=nthreads): for i in prange(0, nx, nogil=True, schedule='static', num_threads=nthreads):
for j in range(0, ny): for j in range(0, ny):
for k in range(0, nz): for k in range(0, nz):
listIndex = ID[4, i + xs, j + ys, zf - (k + 1)] material = ID[4, i + xs, j + ys, zf - (k + 1)]
dEx = (Ex[i + xs, j + ys, zf - k] - Ex[i + xs, j + ys, zf - (k + 1)]) / dz dEx = (Ex[i + xs, j + ys, zf - k] - Ex[i + xs, j + ys, zf - (k + 1)]) / dz
Hy[i + xs, j + ys, zf - (k + 1)] = Hy[i + xs, j + ys, zf - (k + 1)] - updatecoeffsH[listIndex, 4] * ((RA[0, k] - 1) * dEx + RB[0, k] * HPhi[0, i, j, k]) Hy[i + xs, j + ys, zf - (k + 1)] = Hy[i + xs, j + ys, zf - (k + 1)] - updatecoeffsH[material, 4] * ((RA[0, k] - 1) * dEx + RB[0, k] * HPhi[0, i, j, k])
HPhi[0, i, j, k] = RE[0, k] * HPhi[0, i, j, k] - RF[0, k] * dEx HPhi[0, i, j, k] = RE[0, k] * HPhi[0, i, j, k] - RF[0, k] * dEx
############################################# #############################################
# Magnetic field PML updates - Hz component # # Magnetic field PML updates - Hz component #
############################################# #############################################
cpdef update_pml_1order_hz_xplus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, :] updatecoeffsH, np.uint32_t[:, :, :, :] ID, floattype_t[:, :, :] Hz, floattype_t[:, :, :] Ey, floattype_t[:, :, :, :] HPhi, floattype_t[:, :] RA, floattype_t[:, :] RB, floattype_t[:, :] RE, floattype_t[:, :] RF, float dx): cpdef void update_pml_1order_hz_xplus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, ::1] updatecoeffsH, np.uint32_t[:, :, :, ::1] ID, floattype_t[:, :, ::1] Hz, floattype_t[:, :, ::1] Ey, floattype_t[:, :, :, ::1] HPhi, floattype_t[:, ::1] RA, floattype_t[:, ::1] RB, floattype_t[:, ::1] RE, floattype_t[:, ::1] RF, float dx):
"""This function updates the Hz field components in the x stretching direction. """This function updates the Hz field components in the x stretching direction.
Args: Args:
@@ -571,22 +591,23 @@ cpdef update_pml_1order_hz_xplus(int xs, int xf, int ys, int yf, int zs, int zf,
dx, dy, dz (float): Spatial discretisation dx, dy, dz (float): Spatial discretisation
""" """
cdef int i, j, k, nx, ny, nz, listIndex cdef Py_ssize_t i, j, k
cdef int nx, ny, nz, material
cdef float dEy cdef float dEy
nx = xf - xs nx = xf - xs
ny = yf - ys ny = yf - ys
nz = zf - zs nz = zf - zs
for i in prange(0, nx, nogil=True, schedule='static', chunksize=1, num_threads=nthreads): for i in prange(0, nx, nogil=True, schedule='static', num_threads=nthreads):
for j in range(0, ny): for j in range(0, ny):
for k in range(0, nz): for k in range(0, nz):
listIndex = ID[5, i + xs, j + ys, k + zs] material = ID[5, i + xs, j + ys, k + zs]
dEy = (Ey[i + 1 + xs, j + ys, k + zs] - Ey[i + xs, j + ys, k + zs]) / dx dEy = (Ey[i + 1 + xs, j + ys, k + zs] - Ey[i + xs, j + ys, k + zs]) / dx
Hz[i + xs, j + ys, k + zs] = Hz[i + xs, j + ys, k + zs] - updatecoeffsH[listIndex, 4] * ((RA[0, i] - 1) * dEy + RB[0, i] * HPhi[0, i, j, k]) Hz[i + xs, j + ys, k + zs] = Hz[i + xs, j + ys, k + zs] - updatecoeffsH[material, 4] * ((RA[0, i] - 1) * dEy + RB[0, i] * HPhi[0, i, j, k])
HPhi[0, i, j, k] = RE[0, i] * HPhi[0, i, j, k] - RF[0, i] * dEy HPhi[0, i, j, k] = RE[0, i] * HPhi[0, i, j, k] - RF[0, i] * dEy
cpdef update_pml_1order_hz_xminus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, :] updatecoeffsH, np.uint32_t[:, :, :, :] ID, floattype_t[:, :, :] Hz, floattype_t[:, :, :] Ey, floattype_t[:, :, :, :] HPhi, floattype_t[:, :] RA, floattype_t[:, :] RB, floattype_t[:, :] RE, floattype_t[:, :] RF, float dx): cpdef void update_pml_1order_hz_xminus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, ::1] updatecoeffsH, np.uint32_t[:, :, :, ::1] ID, floattype_t[:, :, ::1] Hz, floattype_t[:, :, ::1] Ey, floattype_t[:, :, :, ::1] HPhi, floattype_t[:, ::1] RA, floattype_t[:, ::1] RB, floattype_t[:, ::1] RE, floattype_t[:, ::1] RF, float dx):
"""This function updates the Hz field components in the x stretching direction. """This function updates the Hz field components in the x stretching direction.
Args: Args:
@@ -597,22 +618,23 @@ cpdef update_pml_1order_hz_xminus(int xs, int xf, int ys, int yf, int zs, int zf
dx, dy, dz (float): Spatial discretisation dx, dy, dz (float): Spatial discretisation
""" """
cdef int i, j, k, nx, ny, nz, listIndex cdef Py_ssize_t i, j, k
cdef int nx, ny, nz, material
cdef float dEy cdef float dEy
nx = xf - xs nx = xf - xs
ny = yf - ys ny = yf - ys
nz = zf - zs nz = zf - zs
for i in prange(0, nx, nogil=True, schedule='static', chunksize=1, num_threads=nthreads): for i in prange(0, nx, nogil=True, schedule='static', num_threads=nthreads):
for j in range(0, ny): for j in range(0, ny):
for k in range(0, nz): for k in range(0, nz):
listIndex = ID[5, xf - (i + 1), j + ys, k + zs] material = ID[5, xf - (i + 1), j + ys, k + zs]
dEy = (Ey[xf - i, j + ys, k + zs] - Ey[xf - (i + 1), j + ys, k + zs]) / dx dEy = (Ey[xf - i, j + ys, k + zs] - Ey[xf - (i + 1), j + ys, k + zs]) / dx
Hz[xf - (i + 1), j + ys, k + zs] = Hz[xf - (i + 1), j + ys, k + zs] - updatecoeffsH[listIndex, 4] * ((RA[0, i] - 1) * dEy + RB[0, i] * HPhi[0, i, j, k]) Hz[xf - (i + 1), j + ys, k + zs] = Hz[xf - (i + 1), j + ys, k + zs] - updatecoeffsH[material, 4] * ((RA[0, i] - 1) * dEy + RB[0, i] * HPhi[0, i, j, k])
HPhi[0, i, j, k] = RE[0, i] * HPhi[0, i, j, k] - RF[0, i] * dEy HPhi[0, i, j, k] = RE[0, i] * HPhi[0, i, j, k] - RF[0, i] * dEy
cpdef update_pml_1order_hz_yplus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, :] updatecoeffsH, np.uint32_t[:, :, :, :] ID, floattype_t[:, :, :] Hz, floattype_t[:, :, :] Ex, floattype_t[:, :, :, :] HPhi, floattype_t[:, :] RA, floattype_t[:, :] RB, floattype_t[:, :] RE, floattype_t[:, :] RF, float dy): cpdef void update_pml_1order_hz_yplus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, ::1] updatecoeffsH, np.uint32_t[:, :, :, ::1] ID, floattype_t[:, :, ::1] Hz, floattype_t[:, :, ::1] Ex, floattype_t[:, :, :, ::1] HPhi, floattype_t[:, ::1] RA, floattype_t[:, ::1] RB, floattype_t[:, ::1] RE, floattype_t[:, ::1] RF, float dy):
"""This function updates the Hz field components in the y stretching direction. """This function updates the Hz field components in the y stretching direction.
Args: Args:
@@ -623,22 +645,23 @@ cpdef update_pml_1order_hz_yplus(int xs, int xf, int ys, int yf, int zs, int zf,
dx, dy, dz (float): Spatial discretisation dx, dy, dz (float): Spatial discretisation
""" """
cdef int i, j, k, nx, ny, nz, listIndex cdef Py_ssize_t i, j, k
cdef int nx, ny, nz, material
cdef float dEx cdef float dEx
nx = xf - xs nx = xf - xs
ny = yf - ys ny = yf - ys
nz = zf - zs nz = zf - zs
for i in prange(0, nx, nogil=True, schedule='static', chunksize=1, num_threads=nthreads): for i in prange(0, nx, nogil=True, schedule='static', num_threads=nthreads):
for j in range(0, ny): for j in range(0, ny):
for k in range(0, nz): for k in range(0, nz):
listIndex = ID[5, i + xs, j + ys, k + zs] material = ID[5, i + xs, j + ys, k + zs]
dEx = (Ex[i + xs, j + 1 + ys, k + zs] - Ex[i + xs, j + ys, k + zs]) / dy dEx = (Ex[i + xs, j + 1 + ys, k + zs] - Ex[i + xs, j + ys, k + zs]) / dy
Hz[i + xs, j + ys, k + zs] = Hz[i + xs, j + ys, k + zs] + updatecoeffsH[listIndex, 4] * ((RA[0, j] - 1) * dEx + RB[0, j] * HPhi[0, i, j, k]) Hz[i + xs, j + ys, k + zs] = Hz[i + xs, j + ys, k + zs] + updatecoeffsH[material, 4] * ((RA[0, j] - 1) * dEx + RB[0, j] * HPhi[0, i, j, k])
HPhi[0, i, j, k] = RE[0, j] * HPhi[0, i, j, k] - RF[0, j] * dEx HPhi[0, i, j, k] = RE[0, j] * HPhi[0, i, j, k] - RF[0, j] * dEx
cpdef update_pml_1order_hz_yminus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, :] updatecoeffsH, np.uint32_t[:, :, :, :] ID, floattype_t[:, :, :] Hz, floattype_t[:, :, :] Ex, floattype_t[:, :, :, :] HPhi, floattype_t[:, :] RA, floattype_t[:, :] RB, floattype_t[:, :] RE, floattype_t[:, :] RF, float dy): cpdef void update_pml_1order_hz_yminus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, ::1] updatecoeffsH, np.uint32_t[:, :, :, ::1] ID, floattype_t[:, :, ::1] Hz, floattype_t[:, :, ::1] Ex, floattype_t[:, :, :, ::1] HPhi, floattype_t[:, ::1] RA, floattype_t[:, ::1] RB, floattype_t[:, ::1] RE, floattype_t[:, ::1] RF, float dy):
"""This function updates the Hz field components in the y stretching direction. """This function updates the Hz field components in the y stretching direction.
Args: Args:
@@ -649,17 +672,18 @@ cpdef update_pml_1order_hz_yminus(int xs, int xf, int ys, int yf, int zs, int zf
dx, dy, dz (float): Spatial discretisation dx, dy, dz (float): Spatial discretisation
""" """
cdef int i, j, k, nx, ny, nz, listIndex cdef Py_ssize_t i, j, k
cdef int nx, ny, nz, material
cdef float dEx cdef float dEx
nx = xf - xs nx = xf - xs
ny = yf - ys ny = yf - ys
nz = zf - zs nz = zf - zs
for i in prange(0, nx, nogil=True, schedule='static', chunksize=1, num_threads=nthreads): for i in prange(0, nx, nogil=True, schedule='static', num_threads=nthreads):
for j in range(0, ny): for j in range(0, ny):
for k in range(0, nz): for k in range(0, nz):
listIndex = ID[5, i + xs, yf - (j + 1), k + zs] material = ID[5, i + xs, yf - (j + 1), k + zs]
dEx = (Ex[i + xs, yf - j, k + zs] - Ex[i + xs, yf - (j + 1), k + zs]) / dy dEx = (Ex[i + xs, yf - j, k + zs] - Ex[i + xs, yf - (j + 1), k + zs]) / dy
Hz[i + xs, yf - (j + 1), k + zs] = Hz[i + xs, yf - (j + 1), k + zs] + updatecoeffsH[listIndex, 4] * ((RA[0, j] - 1) * dEx + RB[0, j] * HPhi[0, i, j, k]) Hz[i + xs, yf - (j + 1), k + zs] = Hz[i + xs, yf - (j + 1), k + zs] + updatecoeffsH[material, 4] * ((RA[0, j] - 1) * dEx + RB[0, j] * HPhi[0, i, j, k])
HPhi[0, i, j, k] = RE[0, j] * HPhi[0, i, j, k] - RF[0, j] * dEx HPhi[0, i, j, k] = RE[0, j] * HPhi[0, i, j, k] - RF[0, j] * dEx

264
gprMax/pml_2order_update.pyx
查看文件

@@ -25,7 +25,7 @@ from gprMax.constants cimport floattype_t, complextype_t
############################################# #############################################
# Electric field PML updates - Ex component # # Electric field PML updates - Ex component #
############################################# #############################################
cpdef update_pml_2order_ex_yplus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, :] updatecoeffsE, np.uint32_t[:, :, :, :] ID, floattype_t[:, :, :] Ex, floattype_t[:, :, :] Hz, floattype_t[:, :, :, :] EPhi, floattype_t[:, :] RA, floattype_t[:, :] RB, floattype_t[:, :] RE, floattype_t[:, :] RF, float dy): cpdef update_pml_2order_ex_yplus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, ::1] updatecoeffsE, np.uint32_t[:, :, :, ::1] ID, floattype_t[:, :, ::1] Ex, floattype_t[:, :, ::1] Hz, floattype_t[:, :, :, ::1] EPhi, floattype_t[:, ::1] RA, floattype_t[:, ::1] RB, floattype_t[:, ::1] RE, floattype_t[:, ::1] RF, float dy):
"""This function updates the Ex field components in the y stretching direction. """This function updates the Ex field components in the y stretching direction.
Args: Args:
@@ -36,23 +36,24 @@ cpdef update_pml_2order_ex_yplus(int xs, int xf, int ys, int yf, int zs, int zf,
dx, dy, dz (float): Spatial discretisation dx, dy, dz (float): Spatial discretisation
""" """
cdef int i, j, k, nx, ny, nz, listIndex cdef Py_ssize_t i, j, k
cdef int nx, ny, nz, material
cdef float dHz cdef float dHz
nx = xf - xs nx = xf - xs
ny = yf - ys ny = yf - ys
nz = zf - zs nz = zf - zs
for i in prange(0, nx, nogil=True, schedule='static', chunksize=1, num_threads=nthreads): for i in prange(0, nx, nogil=True, schedule='static', num_threads=nthreads):
for j in range(0, ny): for j in range(0, ny):
for k in range(0, nz): for k in range(0, nz):
listIndex = ID[0, i + xs, j + ys, k + zs] material = ID[0, i + xs, j + ys, k + zs]
dHz = (Hz[i + xs, j + ys, k + zs] - Hz[i + xs, j - 1 + ys, k + zs]) / dy dHz = (Hz[i + xs, j + ys, k + zs] - Hz[i + xs, j - 1 + ys, k + zs]) / dy
Ex[i + xs, j + ys, k + zs] = Ex[i + xs, j + ys, k + zs] + updatecoeffsE[listIndex, 4] * ((RA[0, j] * RA[1, j] - 1) * dHz + RA[1, j] * RB[0, j] * EPhi[0, i, j, k] + RB[1, j] * EPhi[1, i, j, k]) Ex[i + xs, j + ys, k + zs] = Ex[i + xs, j + ys, k + zs] + updatecoeffsE[material, 4] * ((RA[0, j] * RA[1, j] - 1) * dHz + RA[1, j] * RB[0, j] * EPhi[0, i, j, k] + RB[1, j] * EPhi[1, i, j, k])
EPhi[1, i, j, k] = RE[1, j] * EPhi[1, i, j, k] - RF[1, j] * (RA[0, j] * dHz + RB[0, j] * EPhi[0, i, j, k]) EPhi[1, i, j, k] = RE[1, j] * EPhi[1, i, j, k] - RF[1, j] * (RA[0, j] * dHz + RB[0, j] * EPhi[0, i, j, k])
EPhi[0, i, j, k] = RE[0, j] * EPhi[0, i, j, k] - RF[0, j] * dHz EPhi[0, i, j, k] = RE[0, j] * EPhi[0, i, j, k] - RF[0, j] * dHz
cpdef update_pml_2order_ex_yminus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, :] updatecoeffsE, np.uint32_t[:, :, :, :] ID, floattype_t[:, :, :] Ex, floattype_t[:, :, :] Hz, floattype_t[:, :, :, :] EPhi, floattype_t[:, :] RA, floattype_t[:, :] RB, floattype_t[:, :] RE, floattype_t[:, :] RF, float dy): cpdef update_pml_2order_ex_yminus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, ::1] updatecoeffsE, np.uint32_t[:, :, :, ::1] ID, floattype_t[:, :, ::1] Ex, floattype_t[:, :, ::1] Hz, floattype_t[:, :, :, ::1] EPhi, floattype_t[:, ::1] RA, floattype_t[:, ::1] RB, floattype_t[:, ::1] RE, floattype_t[:, ::1] RF, float dy):
"""This function updates the Ex field components in the y stretching direction. """This function updates the Ex field components in the y stretching direction.
Args: Args:
@@ -63,23 +64,24 @@ cpdef update_pml_2order_ex_yminus(int xs, int xf, int ys, int yf, int zs, int zf
dx, dy, dz (float): Spatial discretisation dx, dy, dz (float): Spatial discretisation
""" """
cdef int i, j, k, nx, ny, nz, listIndex cdef Py_ssize_t i, j, k
cdef int nx, ny, nz, material
cdef float dHz cdef float dHz
nx = xf - xs nx = xf - xs
ny = yf - ys ny = yf - ys
nz = zf - zs nz = zf - zs
for i in prange(0, nx, nogil=True, schedule='static', chunksize=1, num_threads=nthreads): for i in prange(0, nx, nogil=True, schedule='static', num_threads=nthreads):
for j in range(0, ny): for j in range(0, ny):
for k in range(0, nz): for k in range(0, nz):
listIndex = ID[0, i + xs, yf - j, k + zs] material = ID[0, i + xs, yf - j, k + zs]
dHz = (Hz[i + xs, yf - j, k + zs] - Hz[i + xs, yf - j - 1, k + zs]) / dy dHz = (Hz[i + xs, yf - j, k + zs] - Hz[i + xs, yf - j - 1, k + zs]) / dy
Ex[i + xs, yf - j, k + zs] = Ex[i + xs, yf - j, k + zs] + updatecoeffsE[listIndex, 4] * ((RA[0, j] * RA[1, j] - 1) * dHz + RA[1, j] * RB[0, j] * EPhi[0, i, j, k] + RB[1, j] * EPhi[1, i, j, k]) Ex[i + xs, yf - j, k + zs] = Ex[i + xs, yf - j, k + zs] + updatecoeffsE[material, 4] * ((RA[0, j] * RA[1, j] - 1) * dHz + RA[1, j] * RB[0, j] * EPhi[0, i, j, k] + RB[1, j] * EPhi[1, i, j, k])
EPhi[1, i, j, k] = RE[1, j] * EPhi[1, i, j, k] - RF[1, j] * (RA[0, j] * dHz + RB[0, j] * EPhi[0, i, j, k]) EPhi[1, i, j, k] = RE[1, j] * EPhi[1, i, j, k] - RF[1, j] * (RA[0, j] * dHz + RB[0, j] * EPhi[0, i, j, k])
EPhi[0, i, j, k] = RE[0, j] * EPhi[0, i, j, k] - RF[0, j] * dHz EPhi[0, i, j, k] = RE[0, j] * EPhi[0, i, j, k] - RF[0, j] * dHz
cpdef update_pml_2order_ex_zplus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, :] updatecoeffsE, np.uint32_t[:, :, :, :] ID, floattype_t[:, :, :] Ex, floattype_t[:, :, :] Hy, floattype_t[:, :, :, :] EPhi, floattype_t[:, :] RA, floattype_t[:, :] RB, floattype_t[:, :] RE, floattype_t[:, :] RF, float dz): cpdef update_pml_2order_ex_zplus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, ::1] updatecoeffsE, np.uint32_t[:, :, :, ::1] ID, floattype_t[:, :, ::1] Ex, floattype_t[:, :, ::1] Hy, floattype_t[:, :, :, ::1] EPhi, floattype_t[:, ::1] RA, floattype_t[:, ::1] RB, floattype_t[:, ::1] RE, floattype_t[:, ::1] RF, float dz):
"""This function updates the Ex field components in the z stretching direction. """This function updates the Ex field components in the z stretching direction.
Args: Args:
@@ -90,23 +92,24 @@ cpdef update_pml_2order_ex_zplus(int xs, int xf, int ys, int yf, int zs, int zf,
dx, dy, dz (float): Spatial discretisation dx, dy, dz (float): Spatial discretisation
""" """
cdef int i, j, k, nx, ny, nz, listIndex cdef Py_ssize_t i, j, k
cdef int nx, ny, nz, material
cdef float dHy cdef float dHy
nx = xf - xs nx = xf - xs
ny = yf - ys ny = yf - ys
nz = zf - zs nz = zf - zs
for i in prange(0, nx, nogil=True, schedule='static', chunksize=1, num_threads=nthreads): for i in prange(0, nx, nogil=True, schedule='static', num_threads=nthreads):
for j in range(0, ny): for j in range(0, ny):
for k in range(0, nz): for k in range(0, nz):
listIndex = ID[0, i + xs, j + ys, k + zs] material = ID[0, i + xs, j + ys, k + zs]
dHy = (Hy[i + xs, j + ys, k + zs] - Hy[i + xs, j + ys, k - 1 + zs]) / dz dHy = (Hy[i + xs, j + ys, k + zs] - Hy[i + xs, j + ys, k - 1 + zs]) / dz
Ex[i + xs, j + ys, k + zs] = Ex[i + xs, j + ys, k + zs] - updatecoeffsE[listIndex, 4] * ((RA[0, k] * RA[1, k] - 1) * dHy + RA[1, k] * RB[0, k] * EPhi[0, i, j, k] + RB[1, k] * EPhi[1, i, j, k]) Ex[i + xs, j + ys, k + zs] = Ex[i + xs, j + ys, k + zs] - updatecoeffsE[material, 4] * ((RA[0, k] * RA[1, k] - 1) * dHy + RA[1, k] * RB[0, k] * EPhi[0, i, j, k] + RB[1, k] * EPhi[1, i, j, k])
EPhi[1, i, j, k] = RE[1, k] * EPhi[1, i, j, k] - RF[1, k] * (RA[0, k] * dHy + RB[0, k] * EPhi[0, i, j, k]) EPhi[1, i, j, k] = RE[1, k] * EPhi[1, i, j, k] - RF[1, k] * (RA[0, k] * dHy + RB[0, k] * EPhi[0, i, j, k])
EPhi[0, i, j, k] = RE[0, k] * EPhi[0, i, j, k] - RF[0, k] * dHy EPhi[0, i, j, k] = RE[0, k] * EPhi[0, i, j, k] - RF[0, k] * dHy
cpdef update_pml_2order_ex_zminus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, :] updatecoeffsE, np.uint32_t[:, :, :, :] ID, floattype_t[:, :, :] Ex, floattype_t[:, :, :] Hy, floattype_t[:, :, :, :] EPhi, floattype_t[:, :] RA, floattype_t[:, :] RB, floattype_t[:, :] RE, floattype_t[:, :] RF, float dz): cpdef update_pml_2order_ex_zminus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, ::1] updatecoeffsE, np.uint32_t[:, :, :, ::1] ID, floattype_t[:, :, ::1] Ex, floattype_t[:, :, ::1] Hy, floattype_t[:, :, :, ::1] EPhi, floattype_t[:, ::1] RA, floattype_t[:, ::1] RB, floattype_t[:, ::1] RE, floattype_t[:, ::1] RF, float dz):
"""This function updates the Ex field components in the z stretching direction. """This function updates the Ex field components in the z stretching direction.
Args: Args:
@@ -117,18 +120,19 @@ cpdef update_pml_2order_ex_zminus(int xs, int xf, int ys, int yf, int zs, int zf
dx, dy, dz (float): Spatial discretisation dx, dy, dz (float): Spatial discretisation
""" """
cdef int i, j, k, nx, ny, nz, listIndex cdef Py_ssize_t i, j, k
cdef int nx, ny, nz, material
cdef float dHy cdef float dHy
nx = xf - xs nx = xf - xs
ny = yf - ys ny = yf - ys
nz = zf - zs nz = zf - zs
for i in prange(0, nx, nogil=True, schedule='static', chunksize=1, num_threads=nthreads): for i in prange(0, nx, nogil=True, schedule='static', num_threads=nthreads):
for j in range(0, ny): for j in range(0, ny):
for k in range(0, nz): for k in range(0, nz):
listIndex = ID[0, i + xs, j + ys, zf - k] material = ID[0, i + xs, j + ys, zf - k]
dHy = (Hy[i + xs, j + ys, zf - k] - Hy[i + xs, j + ys, zf - k - 1]) / dz dHy = (Hy[i + xs, j + ys, zf - k] - Hy[i + xs, j + ys, zf - k - 1]) / dz
Ex[i + xs, j + ys, zf - k] = Ex[i + xs, j + ys, zf - k] - updatecoeffsE[listIndex, 4] * ((RA[0, k] * RA[1, k] - 1) * dHy + RA[1, k] * RB[0, k] * EPhi[0, i, j, k] + RB[1, k] * EPhi[1, i, j, k]) Ex[i + xs, j + ys, zf - k] = Ex[i + xs, j + ys, zf - k] - updatecoeffsE[material, 4] * ((RA[0, k] * RA[1, k] - 1) * dHy + RA[1, k] * RB[0, k] * EPhi[0, i, j, k] + RB[1, k] * EPhi[1, i, j, k])
EPhi[1, i, j, k] = RE[1, k] * EPhi[1, i, j, k] - RF[1, k] * (RA[0, k] * dHy + RB[0, k] * EPhi[0, i, j, k]) EPhi[1, i, j, k] = RE[1, k] * EPhi[1, i, j, k] - RF[1, k] * (RA[0, k] * dHy + RB[0, k] * EPhi[0, i, j, k])
EPhi[0, i, j, k] = RE[0, k] * EPhi[0, i, j, k] - RF[0, k] * dHy EPhi[0, i, j, k] = RE[0, k] * EPhi[0, i, j, k] - RF[0, k] * dHy
@@ -136,7 +140,7 @@ cpdef update_pml_2order_ex_zminus(int xs, int xf, int ys, int yf, int zs, int zf
############################################# #############################################
# Electric field PML updates - Ey component # # Electric field PML updates - Ey component #
############################################# #############################################
cpdef update_pml_2order_ey_xplus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, :] updatecoeffsE, np.uint32_t[:, :, :, :] ID, floattype_t[:, :, :] Ey, floattype_t[:, :, :] Hz, floattype_t[:, :, :, :] EPhi, floattype_t[:, :] RA, floattype_t[:, :] RB, floattype_t[:, :] RE, floattype_t[:, :] RF, float dx): cpdef update_pml_2order_ey_xplus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, ::1] updatecoeffsE, np.uint32_t[:, :, :, ::1] ID, floattype_t[:, :, ::1] Ey, floattype_t[:, :, ::1] Hz, floattype_t[:, :, :, ::1] EPhi, floattype_t[:, ::1] RA, floattype_t[:, ::1] RB, floattype_t[:, ::1] RE, floattype_t[:, ::1] RF, float dx):
"""This function updates the Ey field components in the x stretching direction. """This function updates the Ey field components in the x stretching direction.
Args: Args:
@@ -147,23 +151,24 @@ cpdef update_pml_2order_ey_xplus(int xs, int xf, int ys, int yf, int zs, int zf,
dx, dy, dz (float): Spatial discretisation dx, dy, dz (float): Spatial discretisation
""" """
cdef int i, j, k, nx, ny, nz, listIndex cdef Py_ssize_t i, j, k
cdef int nx, ny, nz, material
cdef float dHz cdef float dHz
nx = xf - xs nx = xf - xs
ny = yf - ys ny = yf - ys
nz = zf - zs nz = zf - zs
for i in prange(0, nx, nogil=True, schedule='static', chunksize=1, num_threads=nthreads): for i in prange(0, nx, nogil=True, schedule='static', num_threads=nthreads):
for j in range(0, ny): for j in range(0, ny):
for k in range(0, nz): for k in range(0, nz):
listIndex = ID[1, i + xs, j + ys, k + zs] material = ID[1, i + xs, j + ys, k + zs]
dHz = (Hz[i + xs, j + ys, k + zs] - Hz[i - 1 + xs, j + ys, k + zs]) / dx dHz = (Hz[i + xs, j + ys, k + zs] - Hz[i - 1 + xs, j + ys, k + zs]) / dx
Ey[i + xs, j + ys, k + zs] = Ey[i + xs, j + ys, k + zs] - updatecoeffsE[listIndex, 4] * ((RA[0, i] * RA[1, i] - 1) * dHz + RA[1, i] * RB[0, i] * EPhi[0, i, j, k] + RB[1, i] * EPhi[1, i, j, k]) Ey[i + xs, j + ys, k + zs] = Ey[i + xs, j + ys, k + zs] - updatecoeffsE[material, 4] * ((RA[0, i] * RA[1, i] - 1) * dHz + RA[1, i] * RB[0, i] * EPhi[0, i, j, k] + RB[1, i] * EPhi[1, i, j, k])
EPhi[1, i, j, k] = RE[1, i] * EPhi[1, i, j, k] - RF[1, i] * (RA[0, i] * dHz + RB[0, i] * EPhi[0, i, j, k]) EPhi[1, i, j, k] = RE[1, i] * EPhi[1, i, j, k] - RF[1, i] * (RA[0, i] * dHz + RB[0, i] * EPhi[0, i, j, k])
EPhi[0, i, j, k] = RE[0, i] * EPhi[0, i, j, k] - RF[0, i] * dHz EPhi[0, i, j, k] = RE[0, i] * EPhi[0, i, j, k] - RF[0, i] * dHz
cpdef update_pml_2order_ey_xminus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, :] updatecoeffsE, np.uint32_t[:, :, :, :] ID, floattype_t[:, :, :] Ey, floattype_t[:, :, :] Hz, floattype_t[:, :, :, :] EPhi, floattype_t[:, :] RA, floattype_t[:, :] RB, floattype_t[:, :] RE, floattype_t[:, :] RF, float dx): cpdef update_pml_2order_ey_xminus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, ::1] updatecoeffsE, np.uint32_t[:, :, :, ::1] ID, floattype_t[:, :, ::1] Ey, floattype_t[:, :, ::1] Hz, floattype_t[:, :, :, ::1] EPhi, floattype_t[:, ::1] RA, floattype_t[:, ::1] RB, floattype_t[:, ::1] RE, floattype_t[:, ::1] RF, float dx):
"""This function updates the Ey field components in the x stretching direction. """This function updates the Ey field components in the x stretching direction.
Args: Args:
@@ -174,23 +179,24 @@ cpdef update_pml_2order_ey_xminus(int xs, int xf, int ys, int yf, int zs, int zf
dx, dy, dz (float): Spatial discretisation dx, dy, dz (float): Spatial discretisation
""" """
cdef int i, j, k, nx, ny, nz, listIndex cdef Py_ssize_t i, j, k
cdef int nx, ny, nz, material
cdef float dHz cdef float dHz
nx = xf - xs nx = xf - xs
ny = yf - ys ny = yf - ys
nz = zf - zs nz = zf - zs
for i in prange(0, nx, nogil=True, schedule='static', chunksize=1, num_threads=nthreads): for i in prange(0, nx, nogil=True, schedule='static', num_threads=nthreads):
for j in range(0, ny): for j in range(0, ny):
for k in range(0, nz): for k in range(0, nz):
listIndex = ID[1, xf - i, j + ys, k + zs] material = ID[1, xf - i, j + ys, k + zs]
dHz = (Hz[xf - i, j + ys, k + zs] - Hz[xf - i - 1, j + ys, k + zs]) / dx dHz = (Hz[xf - i, j + ys, k + zs] - Hz[xf - i - 1, j + ys, k + zs]) / dx
Ey[xf - i, j + ys, k + zs] = Ey[xf - i, j + ys, k + zs] - updatecoeffsE[listIndex, 4] * ((RA[0, i] * RA[1, i] - 1) * dHz + RA[1, i] * RB[0, i] * EPhi[0, i, j, k] + RB[1, i] * EPhi[1, i, j, k]) Ey[xf - i, j + ys, k + zs] = Ey[xf - i, j + ys, k + zs] - updatecoeffsE[material, 4] * ((RA[0, i] * RA[1, i] - 1) * dHz + RA[1, i] * RB[0, i] * EPhi[0, i, j, k] + RB[1, i] * EPhi[1, i, j, k])
EPhi[1, i, j, k] = RE[1, i] * EPhi[1, i, j, k] - RF[1, i] * (RA[0, i] * dHz + RB[0, i] * EPhi[0, i, j, k]) EPhi[1, i, j, k] = RE[1, i] * EPhi[1, i, j, k] - RF[1, i] * (RA[0, i] * dHz + RB[0, i] * EPhi[0, i, j, k])
EPhi[0, i, j, k] = RE[0, i] * EPhi[0, i, j, k] - RF[0, i] * dHz EPhi[0, i, j, k] = RE[0, i] * EPhi[0, i, j, k] - RF[0, i] * dHz
cpdef update_pml_2order_ey_zplus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, :] updatecoeffsE, np.uint32_t[:, :, :, :] ID, floattype_t[:, :, :] Ey, floattype_t[:, :, :] Hx, floattype_t[:, :, :, :] EPhi, floattype_t[:, :] RA, floattype_t[:, :] RB, floattype_t[:, :] RE, floattype_t[:, :] RF, float dz): cpdef update_pml_2order_ey_zplus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, ::1] updatecoeffsE, np.uint32_t[:, :, :, ::1] ID, floattype_t[:, :, ::1] Ey, floattype_t[:, :, ::1] Hx, floattype_t[:, :, :, ::1] EPhi, floattype_t[:, ::1] RA, floattype_t[:, ::1] RB, floattype_t[:, ::1] RE, floattype_t[:, ::1] RF, float dz):
"""This function updates the Ey field components in the z stretching direction. """This function updates the Ey field components in the z stretching direction.
Args: Args:
@@ -201,23 +207,24 @@ cpdef update_pml_2order_ey_zplus(int xs, int xf, int ys, int yf, int zs, int zf,
dx, dy, dz (float): Spatial discretisation dx, dy, dz (float): Spatial discretisation
""" """
cdef int i, j, k, nx, ny, nz, listIndex cdef Py_ssize_t i, j, k
cdef int nx, ny, nz, material
cdef float dHx cdef float dHx
nx = xf - xs nx = xf - xs
ny = yf - ys ny = yf - ys
nz = zf - zs nz = zf - zs
for i in prange(0, nx, nogil=True, schedule='static', chunksize=1, num_threads=nthreads): for i in prange(0, nx, nogil=True, schedule='static', num_threads=nthreads):
for j in range(0, ny): for j in range(0, ny):
for k in range(0, nz): for k in range(0, nz):
listIndex = ID[1, i + xs, j + ys, k + zs] material = ID[1, i + xs, j + ys, k + zs]
dHx = (Hx[i + xs, j + ys, k + zs] - Hx[i + xs, j + ys, k - 1 + zs]) / dz dHx = (Hx[i + xs, j + ys, k + zs] - Hx[i + xs, j + ys, k - 1 + zs]) / dz
Ey[i + xs, j + ys, k + zs] = Ey[i + xs, j + ys, k + zs] + updatecoeffsE[listIndex, 4] * ((RA[0, k] * RA[1, k] - 1) * dHx + RA[1, k] * RB[0, k] * EPhi[0, i, j, k] + RB[1, k] * EPhi[1, i, j, k]) Ey[i + xs, j + ys, k + zs] = Ey[i + xs, j + ys, k + zs] + updatecoeffsE[material, 4] * ((RA[0, k] * RA[1, k] - 1) * dHx + RA[1, k] * RB[0, k] * EPhi[0, i, j, k] + RB[1, k] * EPhi[1, i, j, k])
EPhi[1, i, j, k] = RE[1, k] * EPhi[1, i, j, k] - RF[1, k] * (RA[0, k] * dHx + RB[0, k] * EPhi[0, i, j, k]) EPhi[1, i, j, k] = RE[1, k] * EPhi[1, i, j, k] - RF[1, k] * (RA[0, k] * dHx + RB[0, k] * EPhi[0, i, j, k])
EPhi[0, i, j, k] = RE[0, k] * EPhi[0, i, j, k] - RF[0, k] * dHx EPhi[0, i, j, k] = RE[0, k] * EPhi[0, i, j, k] - RF[0, k] * dHx
cpdef update_pml_2order_ey_zminus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, :] updatecoeffsE, np.uint32_t[:, :, :, :] ID, floattype_t[:, :, :] Ey, floattype_t[:, :, :] Hx, floattype_t[:, :, :, :] EPhi, floattype_t[:, :] RA, floattype_t[:, :] RB, floattype_t[:, :] RE, floattype_t[:, :] RF, float dz): cpdef update_pml_2order_ey_zminus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, ::1] updatecoeffsE, np.uint32_t[:, :, :, ::1] ID, floattype_t[:, :, ::1] Ey, floattype_t[:, :, ::1] Hx, floattype_t[:, :, :, ::1] EPhi, floattype_t[:, ::1] RA, floattype_t[:, ::1] RB, floattype_t[:, ::1] RE, floattype_t[:, ::1] RF, float dz):
"""This function updates the Ey field components in the z stretching direction. """This function updates the Ey field components in the z stretching direction.
Args: Args:
@@ -228,18 +235,19 @@ cpdef update_pml_2order_ey_zminus(int xs, int xf, int ys, int yf, int zs, int zf
dx, dy, dz (float): Spatial discretisation dx, dy, dz (float): Spatial discretisation
""" """
cdef int i, j, k, nx, ny, nz, listIndex cdef Py_ssize_t i, j, k
cdef int nx, ny, nz, material
cdef float dHx cdef float dHx
nx = xf - xs nx = xf - xs
ny = yf - ys ny = yf - ys
nz = zf - zs nz = zf - zs
for i in prange(0, nx, nogil=True, schedule='static', chunksize=1, num_threads=nthreads): for i in prange(0, nx, nogil=True, schedule='static', num_threads=nthreads):
for j in range(0, ny): for j in range(0, ny):
for k in range(0, nz): for k in range(0, nz):
listIndex = ID[1, i + xs, j + ys, zf - k] material = ID[1, i + xs, j + ys, zf - k]
dHx = (Hx[i + xs, j + ys, zf - k] - Hx[i + xs, j + ys, zf - k - 1]) / dz dHx = (Hx[i + xs, j + ys, zf - k] - Hx[i + xs, j + ys, zf - k - 1]) / dz
Ey[i + xs, j + ys, zf - k] = Ey[i + xs, j + ys, zf - k] + updatecoeffsE[listIndex, 4] * ((RA[0, k] * RA[1, k] - 1) * dHx + RA[1, k] * RB[0, k] * EPhi[0, i, j, k] + RB[1, k] * EPhi[1, i, j, k]) Ey[i + xs, j + ys, zf - k] = Ey[i + xs, j + ys, zf - k] + updatecoeffsE[material, 4] * ((RA[0, k] * RA[1, k] - 1) * dHx + RA[1, k] * RB[0, k] * EPhi[0, i, j, k] + RB[1, k] * EPhi[1, i, j, k])
EPhi[1, i, j, k] = RE[1, k] * EPhi[1, i, j, k] - RF[1, k] * (RA[0, k] * dHx + RB[0, k] * EPhi[0, i, j, k]) EPhi[1, i, j, k] = RE[1, k] * EPhi[1, i, j, k] - RF[1, k] * (RA[0, k] * dHx + RB[0, k] * EPhi[0, i, j, k])
EPhi[0, i, j, k] = RE[0, k] * EPhi[0, i, j, k] - RF[0, k] * dHx EPhi[0, i, j, k] = RE[0, k] * EPhi[0, i, j, k] - RF[0, k] * dHx
@@ -247,7 +255,7 @@ cpdef update_pml_2order_ey_zminus(int xs, int xf, int ys, int yf, int zs, int zf
############################################# #############################################
# Electric field PML updates - Ez component # # Electric field PML updates - Ez component #
############################################# #############################################
cpdef update_pml_2order_ez_xplus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, :] updatecoeffsE, np.uint32_t[:, :, :, :] ID, floattype_t[:, :, :] Ez, floattype_t[:, :, :] Hy, floattype_t[:, :, :, :] EPhi, floattype_t[:, :] RA, floattype_t[:, :] RB, floattype_t[:, :] RE, floattype_t[:, :] RF, float dx): cpdef update_pml_2order_ez_xplus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, ::1] updatecoeffsE, np.uint32_t[:, :, :, ::1] ID, floattype_t[:, :, ::1] Ez, floattype_t[:, :, ::1] Hy, floattype_t[:, :, :, ::1] EPhi, floattype_t[:, ::1] RA, floattype_t[:, ::1] RB, floattype_t[:, ::1] RE, floattype_t[:, ::1] RF, float dx):
"""This function updates the Ey field components in the z stretching direction. """This function updates the Ey field components in the z stretching direction.
Args: Args:
@@ -258,23 +266,24 @@ cpdef update_pml_2order_ez_xplus(int xs, int xf, int ys, int yf, int zs, int zf,
dx, dy, dz (float): Spatial discretisation dx, dy, dz (float): Spatial discretisation
""" """
cdef int i, j, k, nx, ny, nz, listIndex cdef Py_ssize_t i, j, k
cdef int nx, ny, nz, material
cdef float dHy cdef float dHy
nx = xf - xs nx = xf - xs
ny = yf - ys ny = yf - ys
nz = zf - zs nz = zf - zs
for i in prange(0, nx, nogil=True, schedule='static', chunksize=1, num_threads=nthreads): for i in prange(0, nx, nogil=True, schedule='static', num_threads=nthreads):
for j in range(0, ny): for j in range(0, ny):
for k in range(0, nz): for k in range(0, nz):
listIndex = ID[2, i + xs, j + ys, k + zs] material = ID[2, i + xs, j + ys, k + zs]
dHy = (Hy[i + xs, j + ys, k + zs] - Hy[i - 1 + xs, j + ys, k + zs]) / dx dHy = (Hy[i + xs, j + ys, k + zs] - Hy[i - 1 + xs, j + ys, k + zs]) / dx
Ez[i + xs, j + ys, k + zs] = Ez[i + xs, j + ys, k + zs] + updatecoeffsE[listIndex, 4] * ((RA[0, i] * RA[1, i] - 1) * dHy + RA[1, i] * RB[0, i] * EPhi[0, i, j, k] + RB[1, i] * EPhi[1, i, j, k]) Ez[i + xs, j + ys, k + zs] = Ez[i + xs, j + ys, k + zs] + updatecoeffsE[material, 4] * ((RA[0, i] * RA[1, i] - 1) * dHy + RA[1, i] * RB[0, i] * EPhi[0, i, j, k] + RB[1, i] * EPhi[1, i, j, k])
EPhi[1, i, j, k] = RE[1, i] * EPhi[1, i, j, k] - RF[1, i] * (RA[0, i] * dHy + RB[0, i] * EPhi[0, i, j, k]) EPhi[1, i, j, k] = RE[1, i] * EPhi[1, i, j, k] - RF[1, i] * (RA[0, i] * dHy + RB[0, i] * EPhi[0, i, j, k])
EPhi[0, i, j, k] = RE[0, i] * EPhi[0, i, j, k] - RF[0, i] * dHy EPhi[0, i, j, k] = RE[0, i] * EPhi[0, i, j, k] - RF[0, i] * dHy
cpdef update_pml_2order_ez_xminus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, :] updatecoeffsE, np.uint32_t[:, :, :, :] ID, floattype_t[:, :, :] Ez, floattype_t[:, :, :] Hy, floattype_t[:, :, :, :] EPhi, floattype_t[:, :] RA, floattype_t[:, :] RB, floattype_t[:, :] RE, floattype_t[:, :] RF, float dx): cpdef update_pml_2order_ez_xminus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, ::1] updatecoeffsE, np.uint32_t[:, :, :, ::1] ID, floattype_t[:, :, ::1] Ez, floattype_t[:, :, ::1] Hy, floattype_t[:, :, :, ::1] EPhi, floattype_t[:, ::1] RA, floattype_t[:, ::1] RB, floattype_t[:, ::1] RE, floattype_t[:, ::1] RF, float dx):
"""This function updates the Ez field components in the x stretching direction. """This function updates the Ez field components in the x stretching direction.
Args: Args:
@@ -285,23 +294,24 @@ cpdef update_pml_2order_ez_xminus(int xs, int xf, int ys, int yf, int zs, int zf
dx, dy, dz (float): Spatial discretisation dx, dy, dz (float): Spatial discretisation
""" """
cdef int i, j, k, nx, ny, nz, listIndex cdef Py_ssize_t i, j, k
cdef int nx, ny, nz, material
cdef float dHy cdef float dHy
nx = xf - xs nx = xf - xs
ny = yf - ys ny = yf - ys
nz = zf - zs nz = zf - zs
for i in prange(0, nx, nogil=True, schedule='static', chunksize=1, num_threads=nthreads): for i in prange(0, nx, nogil=True, schedule='static', num_threads=nthreads):
for j in range(0, ny): for j in range(0, ny):
for k in range(0, nz): for k in range(0, nz):
listIndex = ID[2, xf - i, j + ys, k + zs] material = ID[2, xf - i, j + ys, k + zs]
dHy = (Hy[xf - i, j + ys, k + zs] - Hy[xf - i - 1, j + ys, k + zs]) / dx dHy = (Hy[xf - i, j + ys, k + zs] - Hy[xf - i - 1, j + ys, k + zs]) / dx
Ez[xf - i, j + ys, k + zs] = Ez[xf - i, j + ys, k + zs] + updatecoeffsE[listIndex, 4] * ((RA[0, i] * RA[1, i] - 1) * dHy + RA[1, i] * RB[0, i] * EPhi[0, i, j, k] + RB[1, i] * EPhi[1, i, j, k]) Ez[xf - i, j + ys, k + zs] = Ez[xf - i, j + ys, k + zs] + updatecoeffsE[material, 4] * ((RA[0, i] * RA[1, i] - 1) * dHy + RA[1, i] * RB[0, i] * EPhi[0, i, j, k] + RB[1, i] * EPhi[1, i, j, k])
EPhi[1, i, j, k] = RE[1, i] * EPhi[1, i, j, k] - RF[1, i] * (RA[0, i] * dHy + RB[0, i] * EPhi[0, i, j, k]) EPhi[1, i, j, k] = RE[1, i] * EPhi[1, i, j, k] - RF[1, i] * (RA[0, i] * dHy + RB[0, i] * EPhi[0, i, j, k])
EPhi[0, i, j, k] = RE[0, i] * EPhi[0, i, j, k] - RF[0, i] * dHy EPhi[0, i, j, k] = RE[0, i] * EPhi[0, i, j, k] - RF[0, i] * dHy
cpdef update_pml_2order_ez_yplus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, :] updatecoeffsE, np.uint32_t[:, :, :, :] ID, floattype_t[:, :, :] Ez, floattype_t[:, :, :] Hx, floattype_t[:, :, :, :] EPhi, floattype_t[:, :] RA, floattype_t[:, :] RB, floattype_t[:, :] RE, floattype_t[:, :] RF, float dy): cpdef update_pml_2order_ez_yplus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, ::1] updatecoeffsE, np.uint32_t[:, :, :, ::1] ID, floattype_t[:, :, ::1] Ez, floattype_t[:, :, ::1] Hx, floattype_t[:, :, :, ::1] EPhi, floattype_t[:, ::1] RA, floattype_t[:, ::1] RB, floattype_t[:, ::1] RE, floattype_t[:, ::1] RF, float dy):
"""This function updates the Ez field components in the y stretching direction. """This function updates the Ez field components in the y stretching direction.
Args: Args:
@@ -312,23 +322,24 @@ cpdef update_pml_2order_ez_yplus(int xs, int xf, int ys, int yf, int zs, int zf,
dx, dy, dz (float): Spatial discretisation dx, dy, dz (float): Spatial discretisation
""" """
cdef int i, j, k, nx, ny, nz, listIndex cdef Py_ssize_t i, j, k
cdef int nx, ny, nz, material
cdef float dHx cdef float dHx
nx = xf - xs nx = xf - xs
ny = yf - ys ny = yf - ys
nz = zf - zs nz = zf - zs
for i in prange(0, nx, nogil=True, schedule='static', chunksize=1, num_threads=nthreads): for i in prange(0, nx, nogil=True, schedule='static', num_threads=nthreads):
for j in range(0, ny): for j in range(0, ny):
for k in range(0, nz): for k in range(0, nz):
listIndex = ID[2, i + xs, j + ys, k + zs] material = ID[2, i + xs, j + ys, k + zs]
dHx = (Hx[i + xs, j + ys, k + zs] - Hx[i + xs, j - 1 + ys, k + zs]) / dy dHx = (Hx[i + xs, j + ys, k + zs] - Hx[i + xs, j - 1 + ys, k + zs]) / dy
Ez[i + xs, j + ys, k + zs] = Ez[i + xs, j + ys, k + zs] - updatecoeffsE[listIndex, 4] * ((RA[0, j] * RA[1, j] - 1) * dHx + RA[1, j] * RB[0, j] * EPhi[0, i, j, k] + RB[1, j] * EPhi[1, i, j, k]) Ez[i + xs, j + ys, k + zs] = Ez[i + xs, j + ys, k + zs] - updatecoeffsE[material, 4] * ((RA[0, j] * RA[1, j] - 1) * dHx + RA[1, j] * RB[0, j] * EPhi[0, i, j, k] + RB[1, j] * EPhi[1, i, j, k])
EPhi[1, i, j, k] = RE[1, j] * EPhi[1, i, j, k] - RF[1, j] * (RA[0, j] * dHx + RB[0, j] * EPhi[0, i, j, k]) EPhi[1, i, j, k] = RE[1, j] * EPhi[1, i, j, k] - RF[1, j] * (RA[0, j] * dHx + RB[0, j] * EPhi[0, i, j, k])
EPhi[0, i, j, k] = RE[0, j] * EPhi[0, i, j, k] - RF[0, j] * dHx EPhi[0, i, j, k] = RE[0, j] * EPhi[0, i, j, k] - RF[0, j] * dHx
cpdef update_pml_2order_ez_yminus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, :] updatecoeffsE, np.uint32_t[:, :, :, :] ID, floattype_t[:, :, :] Ez, floattype_t[:, :, :] Hx, floattype_t[:, :, :, :] EPhi, floattype_t[:, :] RA, floattype_t[:, :] RB, floattype_t[:, :] RE, floattype_t[:, :] RF, float dy): cpdef update_pml_2order_ez_yminus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, ::1] updatecoeffsE, np.uint32_t[:, :, :, ::1] ID, floattype_t[:, :, ::1] Ez, floattype_t[:, :, ::1] Hx, floattype_t[:, :, :, ::1] EPhi, floattype_t[:, ::1] RA, floattype_t[:, ::1] RB, floattype_t[:, ::1] RE, floattype_t[:, ::1] RF, float dy):
"""This function updates the Ez field components in the y stretching direction. """This function updates the Ez field components in the y stretching direction.
Args: Args:
@@ -339,18 +350,19 @@ cpdef update_pml_2order_ez_yminus(int xs, int xf, int ys, int yf, int zs, int zf
dx, dy, dz (float): Spatial discretisation dx, dy, dz (float): Spatial discretisation
""" """
cdef int i, j, k, nx, ny, nz, listIndex cdef Py_ssize_t i, j, k
cdef int nx, ny, nz, material
cdef float dHx cdef float dHx
nx = xf - xs nx = xf - xs
ny = yf - ys ny = yf - ys
nz = zf - zs nz = zf - zs
for i in prange(0, nx, nogil=True, schedule='static', chunksize=1, num_threads=nthreads): for i in prange(0, nx, nogil=True, schedule='static', num_threads=nthreads):
for j in range(0, ny): for j in range(0, ny):
for k in range(0, nz): for k in range(0, nz):
listIndex = ID[2, i + xs, yf - j, k + zs] material = ID[2, i + xs, yf - j, k + zs]
dHx = (Hx[i + xs, yf - j, k + zs] - Hx[i + xs, yf - j - 1, k + zs]) / dy dHx = (Hx[i + xs, yf - j, k + zs] - Hx[i + xs, yf - j - 1, k + zs]) / dy
Ez[i + xs, yf - j, k + zs] = Ez[i + xs, yf - j, k + zs] - updatecoeffsE[listIndex, 4] * ((RA[0, j] * RA[1, j] - 1) * dHx + RA[1, j] * RB[0, j] * EPhi[0, i, j, k] + RB[1, j] * EPhi[1, i, j, k]) Ez[i + xs, yf - j, k + zs] = Ez[i + xs, yf - j, k + zs] - updatecoeffsE[material, 4] * ((RA[0, j] * RA[1, j] - 1) * dHx + RA[1, j] * RB[0, j] * EPhi[0, i, j, k] + RB[1, j] * EPhi[1, i, j, k])
EPhi[1, i, j, k] = RE[1, j] * EPhi[1, i, j, k] - RF[1, j] * (RA[0, j] * dHx + RB[0, j] * EPhi[0, i, j, k]) EPhi[1, i, j, k] = RE[1, j] * EPhi[1, i, j, k] - RF[1, j] * (RA[0, j] * dHx + RB[0, j] * EPhi[0, i, j, k])
EPhi[0, i, j, k] = RE[0, j] * EPhi[0, i, j, k] - RF[0, j] * dHx EPhi[0, i, j, k] = RE[0, j] * EPhi[0, i, j, k] - RF[0, j] * dHx
@@ -358,7 +370,7 @@ cpdef update_pml_2order_ez_yminus(int xs, int xf, int ys, int yf, int zs, int zf
############################################# #############################################
# Magnetic field PML updates - Hx component # # Magnetic field PML updates - Hx component #
############################################# #############################################
cpdef update_pml_2order_hx_yplus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, :] updatecoeffsH, np.uint32_t[:, :, :, :] ID, floattype_t[:, :, :] Hx, floattype_t[:, :, :] Ez, floattype_t[:, :, :, :] HPhi, floattype_t[:, :] RA, floattype_t[:, :] RB, floattype_t[:, :] RE, floattype_t[:, :] RF, float dy): cpdef update_pml_2order_hx_yplus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, ::1] updatecoeffsH, np.uint32_t[:, :, :, ::1] ID, floattype_t[:, :, ::1] Hx, floattype_t[:, :, ::1] Ez, floattype_t[:, :, :, ::1] HPhi, floattype_t[:, ::1] RA, floattype_t[:, ::1] RB, floattype_t[:, ::1] RE, floattype_t[:, ::1] RF, float dy):
"""This function updates the Hx field components in the y stretching direction. """This function updates the Hx field components in the y stretching direction.
Args: Args:
@@ -369,23 +381,24 @@ cpdef update_pml_2order_hx_yplus(int xs, int xf, int ys, int yf, int zs, int zf,
dx, dy, dz (float): Spatial discretisation dx, dy, dz (float): Spatial discretisation
""" """
cdef int i, j, k, nx, ny, nz, listIndex cdef Py_ssize_t i, j, k
cdef int nx, ny, nz, material
cdef float dEz cdef float dEz
nx = xf - xs nx = xf - xs
ny = yf - ys ny = yf - ys
nz = zf - zs nz = zf - zs
for i in prange(0, nx, nogil=True, schedule='static', chunksize=1, num_threads=nthreads): for i in prange(0, nx, nogil=True, schedule='static', num_threads=nthreads):
for j in range(0, ny): for j in range(0, ny):
for k in range(0, nz): for k in range(0, nz):
listIndex = ID[3, i + xs, j + ys, k + zs] material = ID[3, i + xs, j + ys, k + zs]
dEz = (Ez[i + xs, j + 1 + ys, k + zs] - Ez[i + xs, j + ys, k + zs]) / dy dEz = (Ez[i + xs, j + 1 + ys, k + zs] - Ez[i + xs, j + ys, k + zs]) / dy
Hx[i + xs, j + ys, k + zs] = Hx[i + xs, j + ys, k + zs] - updatecoeffsH[listIndex, 4] * ((RA[0, j] * RA[1, j] - 1) * dEz + RA[1, j] * RB[0, j] * HPhi[0, i, j, k] + RB[1, j] * HPhi[1, i, j, k]) Hx[i + xs, j + ys, k + zs] = Hx[i + xs, j + ys, k + zs] - updatecoeffsH[material, 4] * ((RA[0, j] * RA[1, j] - 1) * dEz + RA[1, j] * RB[0, j] * HPhi[0, i, j, k] + RB[1, j] * HPhi[1, i, j, k])
HPhi[1, i, j, k] = RE[1, j] * HPhi[1, i, j, k] - RF[1, j] * (RA[0, j] * dEz + RB[0, j] * HPhi[0, i, j, k]) HPhi[1, i, j, k] = RE[1, j] * HPhi[1, i, j, k] - RF[1, j] * (RA[0, j] * dEz + RB[0, j] * HPhi[0, i, j, k])
HPhi[0, i, j, k] = RE[0, j] * HPhi[0, i, j, k] - RF[0, j] * dEz HPhi[0, i, j, k] = RE[0, j] * HPhi[0, i, j, k] - RF[0, j] * dEz
cpdef update_pml_2order_hx_yminus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, :] updatecoeffsH, np.uint32_t[:, :, :, :] ID, floattype_t[:, :, :] Hx, floattype_t[:, :, :] Ez, floattype_t[:, :, :, :] HPhi, floattype_t[:, :] RA, floattype_t[:, :] RB, floattype_t[:, :] RE, floattype_t[:, :] RF, float dy): cpdef update_pml_2order_hx_yminus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, ::1] updatecoeffsH, np.uint32_t[:, :, :, ::1] ID, floattype_t[:, :, ::1] Hx, floattype_t[:, :, ::1] Ez, floattype_t[:, :, :, ::1] HPhi, floattype_t[:, ::1] RA, floattype_t[:, ::1] RB, floattype_t[:, ::1] RE, floattype_t[:, ::1] RF, float dy):
"""This function updates the Hx field components in the y stretching direction. """This function updates the Hx field components in the y stretching direction.
Args: Args:
@@ -396,23 +409,24 @@ cpdef update_pml_2order_hx_yminus(int xs, int xf, int ys, int yf, int zs, int zf
dx, dy, dz (float): Spatial discretisation dx, dy, dz (float): Spatial discretisation
""" """
cdef int i, j, k, nx, ny, nz, listIndex cdef Py_ssize_t i, j, k
cdef int nx, ny, nz, material
cdef float dEz cdef float dEz
nx = xf - xs nx = xf - xs
ny = yf - ys ny = yf - ys
nz = zf - zs nz = zf - zs
for i in prange(0, nx, nogil=True, schedule='static', chunksize=1, num_threads=nthreads): for i in prange(0, nx, nogil=True, schedule='static', num_threads=nthreads):
for j in range(0, ny): for j in range(0, ny):
for k in range(0, nz): for k in range(0, nz):
listIndex = ID[3, i + xs, yf - (j + 1), k + zs] material = ID[3, i + xs, yf - (j + 1), k + zs]
dEz = (Ez[i + xs, yf - j, k + zs] - Ez[i + xs, yf - (j + 1), k + zs]) / dy dEz = (Ez[i + xs, yf - j, k + zs] - Ez[i + xs, yf - (j + 1), k + zs]) / dy
Hx[i + xs, yf - (j + 1), k + zs] = Hx[i + xs, yf - (j + 1), k + zs] - updatecoeffsH[listIndex, 4] * ((RA[0, j] * RA[1, j] - 1) * dEz + RA[1, j] * RB[0, j] * HPhi[0, i, j, k] + RB[1, j] * HPhi[1, i, j, k]) Hx[i + xs, yf - (j + 1), k + zs] = Hx[i + xs, yf - (j + 1), k + zs] - updatecoeffsH[material, 4] * ((RA[0, j] * RA[1, j] - 1) * dEz + RA[1, j] * RB[0, j] * HPhi[0, i, j, k] + RB[1, j] * HPhi[1, i, j, k])
HPhi[1, i, j, k] = RE[1, j] * HPhi[1, i, j, k] - RF[1, j] * (RA[0, j] * dEz + RB[0, j] * HPhi[0, i, j, k]) HPhi[1, i, j, k] = RE[1, j] * HPhi[1, i, j, k] - RF[1, j] * (RA[0, j] * dEz + RB[0, j] * HPhi[0, i, j, k])
HPhi[0, i, j, k] = RE[0, j] * HPhi[0, i, j, k] - RF[0, j] * dEz HPhi[0, i, j, k] = RE[0, j] * HPhi[0, i, j, k] - RF[0, j] * dEz
cpdef update_pml_2order_hx_zplus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, :] updatecoeffsH, np.uint32_t[:, :, :, :] ID, floattype_t[:, :, :] Hx, floattype_t[:, :, :] Ey, floattype_t[:, :, :, :] HPhi, floattype_t[:, :] RA, floattype_t[:, :] RB, floattype_t[:, :] RE, floattype_t[:, :] RF, float dz): cpdef update_pml_2order_hx_zplus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, ::1] updatecoeffsH, np.uint32_t[:, :, :, ::1] ID, floattype_t[:, :, ::1] Hx, floattype_t[:, :, ::1] Ey, floattype_t[:, :, :, ::1] HPhi, floattype_t[:, ::1] RA, floattype_t[:, ::1] RB, floattype_t[:, ::1] RE, floattype_t[:, ::1] RF, float dz):
"""This function updates the Hx field components in the z stretching direction. """This function updates the Hx field components in the z stretching direction.
Args: Args:
@@ -423,23 +437,24 @@ cpdef update_pml_2order_hx_zplus(int xs, int xf, int ys, int yf, int zs, int zf,
dx, dy, dz (float): Spatial discretisation dx, dy, dz (float): Spatial discretisation
""" """
cdef int i, j, k, nx, ny, nz, listIndex cdef Py_ssize_t i, j, k
cdef int nx, ny, nz, material
cdef float dEy cdef float dEy
nx = xf - xs nx = xf - xs
ny = yf - ys ny = yf - ys
nz = zf - zs nz = zf - zs
for i in prange(0, nx, nogil=True, schedule='static', chunksize=1, num_threads=nthreads): for i in prange(0, nx, nogil=True, schedule='static', num_threads=nthreads):
for j in range(0, ny): for j in range(0, ny):
for k in range(0, nz): for k in range(0, nz):
listIndex = ID[3, i + xs, j + ys, k + zs] material = ID[3, i + xs, j + ys, k + zs]
dEy = (Ey[i + xs, j + ys, k + 1 + zs] - Ey[i + xs, j + ys, k + zs]) / dz dEy = (Ey[i + xs, j + ys, k + 1 + zs] - Ey[i + xs, j + ys, k + zs]) / dz
Hx[i + xs, j + ys, k + zs] = Hx[i + xs, j + ys, k + zs] + updatecoeffsH[listIndex, 4] * ((RA[0, k] * RA[1, k] - 1) * dEy + RA[1, k] * RB[0, k] * HPhi[0, i, j, k] + RB[1, k] * HPhi[1, i, j, k]) Hx[i + xs, j + ys, k + zs] = Hx[i + xs, j + ys, k + zs] + updatecoeffsH[material, 4] * ((RA[0, k] * RA[1, k] - 1) * dEy + RA[1, k] * RB[0, k] * HPhi[0, i, j, k] + RB[1, k] * HPhi[1, i, j, k])
HPhi[1, i, j, k] = RE[1, k] * HPhi[1, i, j, k] - RF[1, k] * (RA[0, k] * dEy + RB[0, k] * HPhi[0, i, j, k]) HPhi[1, i, j, k] = RE[1, k] * HPhi[1, i, j, k] - RF[1, k] * (RA[0, k] * dEy + RB[0, k] * HPhi[0, i, j, k])
HPhi[0, i, j, k] = RE[0, k] * HPhi[0, i, j, k] - RF[0, k] * dEy HPhi[0, i, j, k] = RE[0, k] * HPhi[0, i, j, k] - RF[0, k] * dEy
cpdef update_pml_2order_hx_zminus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, :] updatecoeffsH, np.uint32_t[:, :, :, :] ID, floattype_t[:, :, :] Hx, floattype_t[:, :, :] Ey, floattype_t[:, :, :, :] HPhi, floattype_t[:, :] RA, floattype_t[:, :] RB, floattype_t[:, :] RE, floattype_t[:, :] RF, float dz): cpdef update_pml_2order_hx_zminus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, ::1] updatecoeffsH, np.uint32_t[:, :, :, ::1] ID, floattype_t[:, :, ::1] Hx, floattype_t[:, :, ::1] Ey, floattype_t[:, :, :, ::1] HPhi, floattype_t[:, ::1] RA, floattype_t[:, ::1] RB, floattype_t[:, ::1] RE, floattype_t[:, ::1] RF, float dz):
"""This function updates the Hx field components in the z stretching direction. """This function updates the Hx field components in the z stretching direction.
Args: Args:
@@ -450,18 +465,19 @@ cpdef update_pml_2order_hx_zminus(int xs, int xf, int ys, int yf, int zs, int zf
dx, dy, dz (float): Spatial discretisation dx, dy, dz (float): Spatial discretisation
""" """
cdef int i, j, k, nx, ny, nz, listIndex cdef Py_ssize_t i, j, k
cdef int nx, ny, nz, material
cdef float dEy cdef float dEy
nx = xf - xs nx = xf - xs
ny = yf - ys ny = yf - ys
nz = zf - zs nz = zf - zs
for i in prange(0, nx, nogil=True, schedule='static', chunksize=1, num_threads=nthreads): for i in prange(0, nx, nogil=True, schedule='static', num_threads=nthreads):
for j in range(0, ny): for j in range(0, ny):
for k in range(0, nz): for k in range(0, nz):
listIndex = ID[3, i + xs, j + ys, zf - (k + 1)] material = ID[3, i + xs, j + ys, zf - (k + 1)]
dEy = (Ey[i + xs, j + ys, zf - k] - Ey[i + xs, j + ys, zf - (k + 1)]) / dz dEy = (Ey[i + xs, j + ys, zf - k] - Ey[i + xs, j + ys, zf - (k + 1)]) / dz
Hx[i + xs, j + ys, zf - (k + 1)] = Hx[i + xs, j + ys, zf - (k + 1)] + updatecoeffsH[listIndex, 4] * ((RA[0, k] * RA[1, k] - 1) * dEy + RA[1, k] * RB[0, k] * HPhi[0, i, j, k] + RB[1, k] * HPhi[1, i, j, k]) Hx[i + xs, j + ys, zf - (k + 1)] = Hx[i + xs, j + ys, zf - (k + 1)] + updatecoeffsH[material, 4] * ((RA[0, k] * RA[1, k] - 1) * dEy + RA[1, k] * RB[0, k] * HPhi[0, i, j, k] + RB[1, k] * HPhi[1, i, j, k])
HPhi[1, i, j, k] = RE[1, k] * HPhi[1, i, j, k] - RF[1, k] * (RA[0, k] * dEy + RB[0, k] * HPhi[0, i, j, k]) HPhi[1, i, j, k] = RE[1, k] * HPhi[1, i, j, k] - RF[1, k] * (RA[0, k] * dEy + RB[0, k] * HPhi[0, i, j, k])
HPhi[0, i, j, k] = RE[0, k] * HPhi[0, i, j, k] - RF[0, k] * dEy HPhi[0, i, j, k] = RE[0, k] * HPhi[0, i, j, k] - RF[0, k] * dEy
@@ -469,7 +485,7 @@ cpdef update_pml_2order_hx_zminus(int xs, int xf, int ys, int yf, int zs, int zf
############################################# #############################################
# Magnetic field PML updates - Hy component # # Magnetic field PML updates - Hy component #
############################################# #############################################
cpdef update_pml_2order_hy_xplus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, :] updatecoeffsH, np.uint32_t[:, :, :, :] ID, floattype_t[:, :, :] Hy, floattype_t[:, :, :] Ez, floattype_t[:, :, :, :] HPhi, floattype_t[:, :] RA, floattype_t[:, :] RB, floattype_t[:, :] RE, floattype_t[:, :] RF, float dx): cpdef update_pml_2order_hy_xplus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, ::1] updatecoeffsH, np.uint32_t[:, :, :, ::1] ID, floattype_t[:, :, ::1] Hy, floattype_t[:, :, ::1] Ez, floattype_t[:, :, :, ::1] HPhi, floattype_t[:, ::1] RA, floattype_t[:, ::1] RB, floattype_t[:, ::1] RE, floattype_t[:, ::1] RF, float dx):
"""This function updates the Hy field components in the x stretching direction. """This function updates the Hy field components in the x stretching direction.
Args: Args:
@@ -480,23 +496,24 @@ cpdef update_pml_2order_hy_xplus(int xs, int xf, int ys, int yf, int zs, int zf,
dx, dy, dz (float): Spatial discretisation dx, dy, dz (float): Spatial discretisation
""" """
cdef int i, j, k, nx, ny, nz, listIndex cdef Py_ssize_t i, j, k
cdef int nx, ny, nz, material
cdef float dEz cdef float dEz
nx = xf - xs nx = xf - xs
ny = yf - ys ny = yf - ys
nz = zf - zs nz = zf - zs
for i in prange(0, nx, nogil=True, schedule='static', chunksize=1, num_threads=nthreads): for i in prange(0, nx, nogil=True, schedule='static', num_threads=nthreads):
for j in range(0, ny): for j in range(0, ny):
for k in range(0, nz): for k in range(0, nz):
listIndex = ID[4, i + xs, j + ys, k + zs] material = ID[4, i + xs, j + ys, k + zs]
dEz = (Ez[i + 1 + xs, j + ys, k + zs] - Ez[i + xs, j + ys, k + zs]) / dx dEz = (Ez[i + 1 + xs, j + ys, k + zs] - Ez[i + xs, j + ys, k + zs]) / dx
Hy[i + xs, j + ys, k + zs] = Hy[i + xs, j + ys, k + zs] + updatecoeffsH[listIndex, 4] * ((RA[0, i] * RA[1, i] - 1) * dEz + RA[1, i] * RB[0, i] * HPhi[0, i, j, k] + RB[1, i] * HPhi[1, i, j, k]) Hy[i + xs, j + ys, k + zs] = Hy[i + xs, j + ys, k + zs] + updatecoeffsH[material, 4] * ((RA[0, i] * RA[1, i] - 1) * dEz + RA[1, i] * RB[0, i] * HPhi[0, i, j, k] + RB[1, i] * HPhi[1, i, j, k])
HPhi[1, i, j, k] = RE[1, i] * HPhi[1, i, j, k] - RF[1, i] * (RA[0, i] * dEz + RB[0, i] * HPhi[0, i, j, k]) HPhi[1, i, j, k] = RE[1, i] * HPhi[1, i, j, k] - RF[1, i] * (RA[0, i] * dEz + RB[0, i] * HPhi[0, i, j, k])
HPhi[0, i, j, k] = RE[0, i] * HPhi[0, i, j, k] - RF[0, i] * dEz HPhi[0, i, j, k] = RE[0, i] * HPhi[0, i, j, k] - RF[0, i] * dEz
cpdef update_pml_2order_hy_xminus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, :] updatecoeffsH, np.uint32_t[:, :, :, :] ID, floattype_t[:, :, :] Hy, floattype_t[:, :, :] Ez, floattype_t[:, :, :, :] HPhi, floattype_t[:, :] RA, floattype_t[:, :] RB, floattype_t[:, :] RE, floattype_t[:, :] RF, float dx): cpdef update_pml_2order_hy_xminus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, ::1] updatecoeffsH, np.uint32_t[:, :, :, ::1] ID, floattype_t[:, :, ::1] Hy, floattype_t[:, :, ::1] Ez, floattype_t[:, :, :, ::1] HPhi, floattype_t[:, ::1] RA, floattype_t[:, ::1] RB, floattype_t[:, ::1] RE, floattype_t[:, ::1] RF, float dx):
"""This function updates the Hy field components in the x stretching direction. """This function updates the Hy field components in the x stretching direction.
Args: Args:
@@ -507,23 +524,24 @@ cpdef update_pml_2order_hy_xminus(int xs, int xf, int ys, int yf, int zs, int zf
dx, dy, dz (float): Spatial discretisation dx, dy, dz (float): Spatial discretisation
""" """
cdef int i, j, k, nx, ny, nz, listIndex cdef Py_ssize_t i, j, k
cdef int nx, ny, nz, material
cdef float dEz cdef float dEz
nx = xf - xs nx = xf - xs
ny = yf - ys ny = yf - ys
nz = zf - zs nz = zf - zs
for i in prange(0, nx, nogil=True, schedule='static', chunksize=1, num_threads=nthreads): for i in prange(0, nx, nogil=True, schedule='static', num_threads=nthreads):
for j in range(0, ny): for j in range(0, ny):
for k in range(0, nz): for k in range(0, nz):
listIndex = ID[4, xf - (i + 1), j + ys, k + zs] material = ID[4, xf - (i + 1), j + ys, k + zs]
dEz = (Ez[xf - i, j + ys, k + zs] - Ez[xf - (i + 1), j + ys, k + zs]) / dx dEz = (Ez[xf - i, j + ys, k + zs] - Ez[xf - (i + 1), j + ys, k + zs]) / dx
Hy[xf - (i + 1), j + ys, k + zs] = Hy[xf - (i + 1), j + ys, k + zs] + updatecoeffsH[listIndex, 4] * ((RA[0, i] * RA[1, i] - 1) * dEz + RA[1, i] * RB[0, i] * HPhi[0, i, j, k] + RB[1, i] * HPhi[1, i, j, k]) Hy[xf - (i + 1), j + ys, k + zs] = Hy[xf - (i + 1), j + ys, k + zs] + updatecoeffsH[material, 4] * ((RA[0, i] * RA[1, i] - 1) * dEz + RA[1, i] * RB[0, i] * HPhi[0, i, j, k] + RB[1, i] * HPhi[1, i, j, k])
HPhi[1, i, j, k] = RE[1, i] * HPhi[1, i, j, k] - RF[1, i] * (RA[0, i] * dEz + RB[0, i] * HPhi[0, i, j, k]) HPhi[1, i, j, k] = RE[1, i] * HPhi[1, i, j, k] - RF[1, i] * (RA[0, i] * dEz + RB[0, i] * HPhi[0, i, j, k])
HPhi[0, i, j, k] = RE[0, i] * HPhi[0, i, j, k] - RF[0, i] * dEz HPhi[0, i, j, k] = RE[0, i] * HPhi[0, i, j, k] - RF[0, i] * dEz
cpdef update_pml_2order_hy_zplus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, :] updatecoeffsH, np.uint32_t[:, :, :, :] ID, floattype_t[:, :, :] Hy, floattype_t[:, :, :] Ex, floattype_t[:, :, :, :] HPhi, floattype_t[:, :] RA, floattype_t[:, :] RB, floattype_t[:, :] RE, floattype_t[:, :] RF, float dz): cpdef update_pml_2order_hy_zplus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, ::1] updatecoeffsH, np.uint32_t[:, :, :, ::1] ID, floattype_t[:, :, ::1] Hy, floattype_t[:, :, ::1] Ex, floattype_t[:, :, :, ::1] HPhi, floattype_t[:, ::1] RA, floattype_t[:, ::1] RB, floattype_t[:, ::1] RE, floattype_t[:, ::1] RF, float dz):
"""This function updates the Hy field components in the z stretching direction. """This function updates the Hy field components in the z stretching direction.
Args: Args:
@@ -534,23 +552,24 @@ cpdef update_pml_2order_hy_zplus(int xs, int xf, int ys, int yf, int zs, int zf,
dx, dy, dz (float): Spatial discretisation dx, dy, dz (float): Spatial discretisation
""" """
cdef int i, j, k, nx, ny, nz, listIndex cdef Py_ssize_t i, j, k
cdef int nx, ny, nz, material
cdef float dEx cdef float dEx
nx = xf - xs nx = xf - xs
ny = yf - ys ny = yf - ys
nz = zf - zs nz = zf - zs
for i in prange(0, nx, nogil=True, schedule='static', chunksize=1, num_threads=nthreads): for i in prange(0, nx, nogil=True, schedule='static', num_threads=nthreads):
for j in range(0, ny): for j in range(0, ny):
for k in range(0, nz): for k in range(0, nz):
listIndex = ID[4, i + xs, j + ys, k + zs] material = ID[4, i + xs, j + ys, k + zs]
dEx = (Ex[i + xs, j + ys, k + 1 + zs] - Ex[i + xs, j + ys, k + zs]) / dz dEx = (Ex[i + xs, j + ys, k + 1 + zs] - Ex[i + xs, j + ys, k + zs]) / dz
Hy[i + xs, j + ys, k + zs] = Hy[i + xs, j + ys, k + zs] - updatecoeffsH[listIndex, 4] * ((RA[0, k] * RA[1, k] - 1) * dEx + RA[1, k] * RB[0, k] * HPhi[0, i, j, k] + RB[1, k] * HPhi[1, i, j, k]) Hy[i + xs, j + ys, k + zs] = Hy[i + xs, j + ys, k + zs] - updatecoeffsH[material, 4] * ((RA[0, k] * RA[1, k] - 1) * dEx + RA[1, k] * RB[0, k] * HPhi[0, i, j, k] + RB[1, k] * HPhi[1, i, j, k])
HPhi[1, i, j, k] = RE[1, k] * HPhi[1, i, j, k] - RF[1, k] * (RA[0, k] * dEx + RB[0, k] * HPhi[0, i, j, k]) HPhi[1, i, j, k] = RE[1, k] * HPhi[1, i, j, k] - RF[1, k] * (RA[0, k] * dEx + RB[0, k] * HPhi[0, i, j, k])
HPhi[0, i, j, k] = RE[0, k] * HPhi[0, i, j, k] - RF[0, k] * dEx HPhi[0, i, j, k] = RE[0, k] * HPhi[0, i, j, k] - RF[0, k] * dEx
cpdef update_pml_2order_hy_zminus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, :] updatecoeffsH, np.uint32_t[:, :, :, :] ID, floattype_t[:, :, :] Hy, floattype_t[:, :, :] Ex, floattype_t[:, :, :, :] HPhi, floattype_t[:, :] RA, floattype_t[:, :] RB, floattype_t[:, :] RE, floattype_t[:, :] RF, float dz): cpdef update_pml_2order_hy_zminus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, ::1] updatecoeffsH, np.uint32_t[:, :, :, ::1] ID, floattype_t[:, :, ::1] Hy, floattype_t[:, :, ::1] Ex, floattype_t[:, :, :, ::1] HPhi, floattype_t[:, ::1] RA, floattype_t[:, ::1] RB, floattype_t[:, ::1] RE, floattype_t[:, ::1] RF, float dz):
"""This function updates the Hy field components in the z stretching direction. """This function updates the Hy field components in the z stretching direction.
Args: Args:
@@ -561,18 +580,19 @@ cpdef update_pml_2order_hy_zminus(int xs, int xf, int ys, int yf, int zs, int zf
dx, dy, dz (float): Spatial discretisation dx, dy, dz (float): Spatial discretisation
""" """
cdef int i, j, k, nx, ny, nz, listIndex cdef Py_ssize_t i, j, k
cdef int nx, ny, nz, material
cdef float dEx cdef float dEx
nx = xf - xs nx = xf - xs
ny = yf - ys ny = yf - ys
nz = zf - zs nz = zf - zs
for i in prange(0, nx, nogil=True, schedule='static', chunksize=1, num_threads=nthreads): for i in prange(0, nx, nogil=True, schedule='static', num_threads=nthreads):
for j in range(0, ny): for j in range(0, ny):
for k in range(0, nz): for k in range(0, nz):
listIndex = ID[4, i + xs, j + ys, zf - (k + 1)] material = ID[4, i + xs, j + ys, zf - (k + 1)]
dEx = (Ex[i + xs, j + ys, zf - k] - Ex[i + xs, j + ys, zf - (k + 1)]) / dz dEx = (Ex[i + xs, j + ys, zf - k] - Ex[i + xs, j + ys, zf - (k + 1)]) / dz
Hy[i + xs, j + ys, zf - (k + 1)] = Hy[i + xs, j + ys, zf - (k + 1)] - updatecoeffsH[listIndex, 4] * ((RA[0, k] * RA[1, k] - 1) * dEx + RA[1, k] * RB[0, k] * HPhi[0, i, j, k] + RB[1, k] * HPhi[1, i, j, k]) Hy[i + xs, j + ys, zf - (k + 1)] = Hy[i + xs, j + ys, zf - (k + 1)] - updatecoeffsH[material, 4] * ((RA[0, k] * RA[1, k] - 1) * dEx + RA[1, k] * RB[0, k] * HPhi[0, i, j, k] + RB[1, k] * HPhi[1, i, j, k])
HPhi[1, i, j, k] = RE[1, k] * HPhi[1, i, j, k] - RF[1, k] * (RA[0, k] * dEx + RB[0, k] * HPhi[0, i, j, k]) HPhi[1, i, j, k] = RE[1, k] * HPhi[1, i, j, k] - RF[1, k] * (RA[0, k] * dEx + RB[0, k] * HPhi[0, i, j, k])
HPhi[0, i, j, k] = RE[0, k] * HPhi[0, i, j, k] - RF[0, k] * dEx HPhi[0, i, j, k] = RE[0, k] * HPhi[0, i, j, k] - RF[0, k] * dEx
@@ -580,7 +600,7 @@ cpdef update_pml_2order_hy_zminus(int xs, int xf, int ys, int yf, int zs, int zf
############################################# #############################################
# Magnetic field PML updates - Hz component # # Magnetic field PML updates - Hz component #
############################################# #############################################
cpdef update_pml_2order_hz_xplus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, :] updatecoeffsH, np.uint32_t[:, :, :, :] ID, floattype_t[:, :, :] Hz, floattype_t[:, :, :] Ey, floattype_t[:, :, :, :] HPhi, floattype_t[:, :] RA, floattype_t[:, :] RB, floattype_t[:, :] RE, floattype_t[:, :] RF, float dx): cpdef update_pml_2order_hz_xplus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, ::1] updatecoeffsH, np.uint32_t[:, :, :, ::1] ID, floattype_t[:, :, ::1] Hz, floattype_t[:, :, ::1] Ey, floattype_t[:, :, :, ::1] HPhi, floattype_t[:, ::1] RA, floattype_t[:, ::1] RB, floattype_t[:, ::1] RE, floattype_t[:, ::1] RF, float dx):
"""This function updates the Hz field components in the x stretching direction. """This function updates the Hz field components in the x stretching direction.
Args: Args:
@@ -591,23 +611,24 @@ cpdef update_pml_2order_hz_xplus(int xs, int xf, int ys, int yf, int zs, int zf,
dx, dy, dz (float): Spatial discretisation dx, dy, dz (float): Spatial discretisation
""" """
cdef int i, j, k, nx, ny, nz, listIndex cdef Py_ssize_t i, j, k
cdef int nx, ny, nz, material
cdef float dEy cdef float dEy
nx = xf - xs nx = xf - xs
ny = yf - ys ny = yf - ys
nz = zf - zs nz = zf - zs
for i in prange(0, nx, nogil=True, schedule='static', chunksize=1, num_threads=nthreads): for i in prange(0, nx, nogil=True, schedule='static', num_threads=nthreads):
for j in range(0, ny): for j in range(0, ny):
for k in range(0, nz): for k in range(0, nz):
listIndex = ID[5, i + xs, j + ys, k + zs] material = ID[5, i + xs, j + ys, k + zs]
dEy = (Ey[i + 1 + xs, j + ys, k + zs] - Ey[i + xs, j + ys, k + zs]) / dx dEy = (Ey[i + 1 + xs, j + ys, k + zs] - Ey[i + xs, j + ys, k + zs]) / dx
Hz[i + xs, j + ys, k + zs] = Hz[i + xs, j + ys, k + zs] - updatecoeffsH[listIndex, 4] * ((RA[0, i] * RA[1, i] - 1) * dEy + RA[1, i] * RB[0, i] * HPhi[0, i, j, k] + RB[1, i] * HPhi[1, i, j, k]) Hz[i + xs, j + ys, k + zs] = Hz[i + xs, j + ys, k + zs] - updatecoeffsH[material, 4] * ((RA[0, i] * RA[1, i] - 1) * dEy + RA[1, i] * RB[0, i] * HPhi[0, i, j, k] + RB[1, i] * HPhi[1, i, j, k])
HPhi[1, i, j, k] = RE[1, i] * HPhi[1, i, j, k] - RF[1, i] * (RA[0, i] * dEy + RB[0, i] * HPhi[0, i, j, k]) HPhi[1, i, j, k] = RE[1, i] * HPhi[1, i, j, k] - RF[1, i] * (RA[0, i] * dEy + RB[0, i] * HPhi[0, i, j, k])
HPhi[0, i, j, k] = RE[0, i] * HPhi[0, i, j, k] - RF[0, i] * dEy HPhi[0, i, j, k] = RE[0, i] * HPhi[0, i, j, k] - RF[0, i] * dEy
cpdef update_pml_2order_hz_xminus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, :] updatecoeffsH, np.uint32_t[:, :, :, :] ID, floattype_t[:, :, :] Hz, floattype_t[:, :, :] Ey, floattype_t[:, :, :, :] HPhi, floattype_t[:, :] RA, floattype_t[:, :] RB, floattype_t[:, :] RE, floattype_t[:, :] RF, float dx): cpdef update_pml_2order_hz_xminus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, ::1] updatecoeffsH, np.uint32_t[:, :, :, ::1] ID, floattype_t[:, :, ::1] Hz, floattype_t[:, :, ::1] Ey, floattype_t[:, :, :, ::1] HPhi, floattype_t[:, ::1] RA, floattype_t[:, ::1] RB, floattype_t[:, ::1] RE, floattype_t[:, ::1] RF, float dx):
"""This function updates the Hz field components in the x stretching direction. """This function updates the Hz field components in the x stretching direction.
Args: Args:
@@ -618,23 +639,24 @@ cpdef update_pml_2order_hz_xminus(int xs, int xf, int ys, int yf, int zs, int zf
dx, dy, dz (float): Spatial discretisation dx, dy, dz (float): Spatial discretisation
""" """
cdef int i, j, k, nx, ny, nz, listIndex cdef Py_ssize_t i, j, k
cdef int nx, ny, nz, material
cdef float dEy cdef float dEy
nx = xf - xs nx = xf - xs
ny = yf - ys ny = yf - ys
nz = zf - zs nz = zf - zs
for i in prange(0, nx, nogil=True, schedule='static', chunksize=1, num_threads=nthreads): for i in prange(0, nx, nogil=True, schedule='static', num_threads=nthreads):
for j in range(0, ny): for j in range(0, ny):
for k in range(0, nz): for k in range(0, nz):
listIndex = ID[5, xf - (i + 1), j + ys, k + zs] material = ID[5, xf - (i + 1), j + ys, k + zs]
dEy = (Ey[xf - i, j + ys, k + zs] - Ey[xf - (i + 1), j + ys, k + zs]) / dx dEy = (Ey[xf - i, j + ys, k + zs] - Ey[xf - (i + 1), j + ys, k + zs]) / dx
Hz[xf - (i + 1), j + ys, k + zs] = Hz[xf - (i + 1), j + ys, k + zs] - updatecoeffsH[listIndex, 4] * ((RA[0, i] * RA[1, i] - 1) * dEy + RA[1, i] * RB[0, i] * HPhi[0, i, j, k] + RB[1, i] * HPhi[1, i, j, k]) Hz[xf - (i + 1), j + ys, k + zs] = Hz[xf - (i + 1), j + ys, k + zs] - updatecoeffsH[material, 4] * ((RA[0, i] * RA[1, i] - 1) * dEy + RA[1, i] * RB[0, i] * HPhi[0, i, j, k] + RB[1, i] * HPhi[1, i, j, k])
HPhi[1, i, j, k] = RE[1, i] * HPhi[1, i, j, k] - RF[1, i] * (RA[0, i] * dEy + RB[0, i] * HPhi[0, i, j, k]) HPhi[1, i, j, k] = RE[1, i] * HPhi[1, i, j, k] - RF[1, i] * (RA[0, i] * dEy + RB[0, i] * HPhi[0, i, j, k])
HPhi[0, i, j, k] = RE[0, i] * HPhi[0, i, j, k] - RF[0, i] * dEy HPhi[0, i, j, k] = RE[0, i] * HPhi[0, i, j, k] - RF[0, i] * dEy
cpdef update_pml_2order_hz_yplus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, :] updatecoeffsH, np.uint32_t[:, :, :, :] ID, floattype_t[:, :, :] Hz, floattype_t[:, :, :] Ex, floattype_t[:, :, :, :] HPhi, floattype_t[:, :] RA, floattype_t[:, :] RB, floattype_t[:, :] RE, floattype_t[:, :] RF, float dy): cpdef update_pml_2order_hz_yplus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, ::1] updatecoeffsH, np.uint32_t[:, :, :, ::1] ID, floattype_t[:, :, ::1] Hz, floattype_t[:, :, ::1] Ex, floattype_t[:, :, :, ::1] HPhi, floattype_t[:, ::1] RA, floattype_t[:, ::1] RB, floattype_t[:, ::1] RE, floattype_t[:, ::1] RF, float dy):
"""This function updates the Hz field components in the y stretching direction. """This function updates the Hz field components in the y stretching direction.
Args: Args:
@@ -645,23 +667,24 @@ cpdef update_pml_2order_hz_yplus(int xs, int xf, int ys, int yf, int zs, int zf,
dx, dy, dz (float): Spatial discretisation dx, dy, dz (float): Spatial discretisation
""" """
cdef int i, j, k, nx, ny, nz, listIndex cdef Py_ssize_t i, j, k
cdef int nx, ny, nz, material
cdef float dEx cdef float dEx
nx = xf - xs nx = xf - xs
ny = yf - ys ny = yf - ys
nz = zf - zs nz = zf - zs
for i in prange(0, nx, nogil=True, schedule='static', chunksize=1, num_threads=nthreads): for i in prange(0, nx, nogil=True, schedule='static', num_threads=nthreads):
for j in range(0, ny): for j in range(0, ny):
for k in range(0, nz): for k in range(0, nz):
listIndex = ID[5, i + xs, j + ys, k + zs] material = ID[5, i + xs, j + ys, k + zs]
dEx = (Ex[i + xs, j + 1 + ys, k + zs] - Ex[i + xs, j + ys, k + zs]) / dy dEx = (Ex[i + xs, j + 1 + ys, k + zs] - Ex[i + xs, j + ys, k + zs]) / dy
Hz[i + xs, j + ys, k + zs] = Hz[i + xs, j + ys, k + zs] + updatecoeffsH[listIndex, 4] * ((RA[0, j] * RA[1, j] - 1) * dEx + RA[1, j] * RB[0, j] * HPhi[0, i, j, k] + RB[1, j] * HPhi[1, i, j, k]) Hz[i + xs, j + ys, k + zs] = Hz[i + xs, j + ys, k + zs] + updatecoeffsH[material, 4] * ((RA[0, j] * RA[1, j] - 1) * dEx + RA[1, j] * RB[0, j] * HPhi[0, i, j, k] + RB[1, j] * HPhi[1, i, j, k])
HPhi[1, i, j, k] = RE[1, j] * HPhi[1, i, j, k] - RF[1, j] * (RA[0, j] * dEx + RB[0, j] * HPhi[0, i, j, k]) HPhi[1, i, j, k] = RE[1, j] * HPhi[1, i, j, k] - RF[1, j] * (RA[0, j] * dEx + RB[0, j] * HPhi[0, i, j, k])
HPhi[0, i, j, k] = RE[0, j] * HPhi[0, i, j, k] - RF[0, j] * dEx HPhi[0, i, j, k] = RE[0, j] * HPhi[0, i, j, k] - RF[0, j] * dEx
cpdef update_pml_2order_hz_yminus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, :] updatecoeffsH, np.uint32_t[:, :, :, :] ID, floattype_t[:, :, :] Hz, floattype_t[:, :, :] Ex, floattype_t[:, :, :, :] HPhi, floattype_t[:, :] RA, floattype_t[:, :] RB, floattype_t[:, :] RE, floattype_t[:, :] RF, float dy): cpdef update_pml_2order_hz_yminus(int xs, int xf, int ys, int yf, int zs, int zf, int nthreads, floattype_t[:, ::1] updatecoeffsH, np.uint32_t[:, :, :, ::1] ID, floattype_t[:, :, ::1] Hz, floattype_t[:, :, ::1] Ex, floattype_t[:, :, :, ::1] HPhi, floattype_t[:, ::1] RA, floattype_t[:, ::1] RB, floattype_t[:, ::1] RE, floattype_t[:, ::1] RF, float dy):
"""This function updates the Hz field components in the y stretching direction. """This function updates the Hz field components in the y stretching direction.
Args: Args:
@@ -672,18 +695,19 @@ cpdef update_pml_2order_hz_yminus(int xs, int xf, int ys, int yf, int zs, int zf
dx, dy, dz (float): Spatial discretisation dx, dy, dz (float): Spatial discretisation
""" """
cdef int i, j, k, nx, ny, nz, listIndex cdef Py_ssize_t i, j, k
cdef int nx, ny, nz, material
cdef float dEx cdef float dEx
nx = xf - xs nx = xf - xs
ny = yf - ys ny = yf - ys
nz = zf - zs nz = zf - zs
for i in prange(0, nx, nogil=True, schedule='static', chunksize=1, num_threads=nthreads): for i in prange(0, nx, nogil=True, schedule='static', num_threads=nthreads):
for j in range(0, ny): for j in range(0, ny):
for k in range(0, nz): for k in range(0, nz):
listIndex = ID[5, i + xs, yf - (j + 1), k + zs] material = ID[5, i + xs, yf - (j + 1), k + zs]
dEx = (Ex[i + xs, yf - j, k + zs] - Ex[i + xs, yf - (j + 1), k + zs]) / dy dEx = (Ex[i + xs, yf - j, k + zs] - Ex[i + xs, yf - (j + 1), k + zs]) / dy
Hz[i + xs, yf - (j + 1), k + zs] = Hz[i + xs, yf - (j + 1), k + zs] + updatecoeffsH[listIndex, 4] * ((RA[0, j] * RA[1, j] - 1) * dEx + RA[1, j] * RB[0, j] * HPhi[0, i, j, k] + RB[1, j] * HPhi[1, i, j, k]) Hz[i + xs, yf - (j + 1), k + zs] = Hz[i + xs, yf - (j + 1), k + zs] + updatecoeffsH[material, 4] * ((RA[0, j] * RA[1, j] - 1) * dEx + RA[1, j] * RB[0, j] * HPhi[0, i, j, k] + RB[1, j] * HPhi[1, i, j, k])
HPhi[1, i, j, k] = RE[1, j] * HPhi[1, i, j, k] - RF[1, j] * (RA[0, j] * dEx + RB[0, j] * HPhi[0, i, j, k]) HPhi[1, i, j, k] = RE[1, j] * HPhi[1, i, j, k] - RF[1, j] * (RA[0, j] * dEx + RB[0, j] * HPhi[0, i, j, k])
HPhi[0, i, j, k] = RE[0, j] * HPhi[0, i, j, k] - RF[0, j] * dEx HPhi[0, i, j, k] = RE[0, j] * HPhi[0, i, j, k] - RF[0, j] * dEx