From d064be785aeb46cc3c34bc31b2905dc0457d7b0f Mon Sep 17 00:00:00 2001 From: Craig Warren Date: Tue, 23 Feb 2016 17:35:10 +0000 Subject: [PATCH] Performance improvements to Cython code. --- gprMax/fields_update.pyx | 222 +++++++++++++++-------------- gprMax/pml_1order_update.pyx | 264 +++++++++++++++++++---------------- gprMax/pml_2order_update.pyx | 264 +++++++++++++++++++---------------- 3 files changed, 408 insertions(+), 342 deletions(-) diff --git a/gprMax/fields_update.pyx b/gprMax/fields_update.pyx index 46898821..9e6df5ed 100644 --- a/gprMax/fields_update.pyx +++ b/gprMax/fields_update.pyx @@ -25,7 +25,7 @@ from gprMax.constants cimport floattype_t, complextype_t ######################################### # 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. 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 """ - cdef int i, j, k, listIndex + cdef Py_ssize_t i, j, k + cdef int material if ny == 1 or nz == 1: pass else: - for i in prange(0, nx, nogil=True, schedule='static', 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 k in range(1, nz): - listIndex = 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]) + material = ID[0, i, j, k] + Ex[i, j, k] = updatecoeffsE[material, 0] * Ex[i, j, k] + updatecoeffsE[material, 2] * (Hz[i, j, k] - Hz[i, j - 1, k]) - updatecoeffsE[material, 3] * (Hy[i, j, k] - Hy[i, j, k - 1]) -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. 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 """ - cdef int i, j, k, listIndex, p + cdef Py_ssize_t i, j, k, p + cdef int material cdef float phi = 0.0 if ny == 1 or nz == 1: pass 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 k in range(1, nz): - listIndex = ID[0, i, j, k] + material = ID[0, i, j, k] phi = 0.0 for p in range(0, maxpoles): - phi = phi + updatecoeffsdispersive[listIndex, 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] - 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 + phi = phi + updatecoeffsdispersive[material, p * 3].real * Tx[p, i, j, k].real + Tx[p, i, j, k] = updatecoeffsdispersive[material, 1 + (p * 3)] * Tx[p, i, j, k] + updatecoeffsdispersive[material, 2 + (p * 3)] * Ex[i, j, k] + 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. 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 """ - cdef int i, j, k, listIndex, p + cdef Py_ssize_t i, j, k, p + cdef int material if ny == 1 or nz == 1: pass 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 k in range(1, nz): - listIndex = ID[0, i, j, k] + material = ID[0, i, j, k] 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. 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 """ - cdef int i, j, k, listIndex + cdef Py_ssize_t i, j, k + cdef int material cdef float phi = 0.0 if ny == 1 or nz == 1: pass 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 k in range(1, nz): - listIndex = ID[0, i, j, k] - phi = updatecoeffsdispersive[listIndex, 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] - 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 + material = ID[0, i, j, k] + phi = updatecoeffsdispersive[material, 0].real * Tx[0, i, j, k].real + Tx[0, i, j, k] = updatecoeffsdispersive[material, 1] * Tx[0, i, j, k] + updatecoeffsdispersive[material, 2] * Ex[i, j, k] + Ex[i, j, k] = updatecoeffsE[material, 0] * Ex[i, j, k] + updatecoeffsE[material, 2] * (Hz[i, j, k] - Hz[i, j - 1, k]) - updatecoeffsE[material, 3] * (Hy[i, j, k] - Hy[i, j, k - 1]) - updatecoeffsE[material, 4] * phi -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. 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 """ - cdef int i, j, k, listIndex + cdef Py_ssize_t i, j, k + cdef int material if ny == 1 or nz == 1: pass else: - for i in prange(0, nx, nogil=True, schedule='static', 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 k in range(1, nz): - listIndex = ID[0, i, j, k] - Tx[0, i, j, k] = Tx[0, i, j, k] - updatecoeffsdispersive[listIndex, 2] * Ex[i, j, k] + material = ID[0, i, j, k] + Tx[0, i, j, k] = Tx[0, i, j, k] - updatecoeffsdispersive[material, 2] * Ex[i, j, k] ######################################### # Electric field updates - Ey component # ######################################### -cpdef 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. 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 """ - cdef int i, j, k, listIndex + cdef Py_ssize_t i, j, k + cdef int material if nx == 1 or nz == 1: pass else: - for i in prange(1, nx, nogil=True, schedule='static', 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 k in range(1, nz): - listIndex = 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]) + material = ID[1, i, j, k] + Ey[i, j, k] = updatecoeffsE[material, 0] * Ey[i, j, k] + updatecoeffsE[material, 3] * (Hx[i, j, k] - Hx[i, j, k - 1]) - updatecoeffsE[material, 1] * (Hz[i, j, k] - Hz[i - 1, j, k]) -cpdef 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. 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 """ - cdef int i, j, k, listIndex, p + cdef Py_ssize_t i, j, k, p + cdef int material cdef float phi = 0.0 if nx == 1 or nz == 1: pass else: - for i in prange(1, nx, nogil=True, schedule='static', 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 k in range(1, nz): - listIndex = ID[1, i, j, k] + material = ID[1, i, j, k] phi = 0.0 for p in range(0, maxpoles): - phi = phi + updatecoeffsdispersive[listIndex, 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] - 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 + phi = phi + updatecoeffsdispersive[material, p * 3].real * Ty[p, i, j, k].real + Ty[p, i, j, k] = updatecoeffsdispersive[material, 1 + (p * 3)] * Ty[p, i, j, k] + updatecoeffsdispersive[material, 2 + (p * 3)] * Ey[i, j, k] + Ey[i, j, k] = updatecoeffsE[material, 0] * Ey[i, j, k] + updatecoeffsE[material, 3] * (Hx[i, j, k] - Hx[i, j, k - 1]) - updatecoeffsE[material, 1] * (Hz[i, j, k] - Hz[i - 1, j, k]) - updatecoeffsE[material, 4] * phi -cpdef 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. 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 """ - cdef int i, j, k, listIndex, p + cdef Py_ssize_t i, j, k, p + cdef int material if nx == 1 or nz == 1: pass else: - for i in prange(1, nx, nogil=True, schedule='static', 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 k in range(1, nz): - listIndex = ID[1, i, j, k] + material = ID[1, i, j, k] 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. 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 """ - cdef int i, j, k, listIndex + cdef Py_ssize_t i, j, k + cdef int material cdef float phi = 0.0 if nx == 1 or nz == 1: pass else: - for i in prange(1, nx, nogil=True, schedule='static', 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 k in range(1, nz): - listIndex = ID[1, i, j, k] - phi = updatecoeffsdispersive[listIndex, 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] - 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 + material = ID[1, i, j, k] + phi = updatecoeffsdispersive[material, 0].real * Ty[0, i, j, k].real + Ty[0, i, j, k] = updatecoeffsdispersive[material, 1] * Ty[0, i, j, k] + updatecoeffsdispersive[material, 2] * Ey[i, j, k] + Ey[i, j, k] = updatecoeffsE[material, 0] * Ey[i, j, k] + updatecoeffsE[material, 3] * (Hx[i, j, k] - Hx[i, j, k - 1]) - updatecoeffsE[material, 1] * (Hz[i, j, k] - Hz[i - 1, j, k]) - updatecoeffsE[material, 4] * phi -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. 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 """ - cdef int i, j, k, listIndex + cdef Py_ssize_t i, j, k + cdef int material if nx == 1 or nz == 1: pass else: - for i in prange(1, nx, nogil=True, schedule='static', 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 k in range(1, nz): - listIndex = ID[1, i, j, k] - Ty[0, i, j, k] = Ty[0, i, j, k] - updatecoeffsdispersive[listIndex, 2] * Ey[i, j, k] + material = ID[1, i, j, k] + Ty[0, i, j, k] = Ty[0, i, j, k] - updatecoeffsdispersive[material, 2] * Ey[i, j, k] ######################################### # Electric field updates - Ez component # ######################################### -cpdef 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. 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 """ - cdef int i, j, k, listIndex + cdef Py_ssize_t i, j, k + cdef int material if nx == 1 or ny == 1: pass else: - for i in prange(1, nx, nogil=True, schedule='static', 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 k in range(0, nz): - listIndex = 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]) + material = ID[2, i, j, k] + Ez[i, j, k] = updatecoeffsE[material, 0] * Ez[i, j, k] + updatecoeffsE[material, 1] * (Hy[i, j, k] - Hy[i - 1, j, k]) - updatecoeffsE[material, 2] * (Hx[i, j, k] - Hx[i, j - 1, k]) -cpdef 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. 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 """ - cdef int i, j, k, listIndex, p + cdef Py_ssize_t i, j, k, p + cdef int material cdef float phi = 0.0 if nx == 1 or ny == 1: pass else: - for i in prange(1, nx, nogil=True, schedule='static', 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 k in range(0, nz): - listIndex = ID[2, i, j, k] + material = ID[2, i, j, k] phi = 0.0 for p in range(0, maxpoles): - phi = phi + updatecoeffsdispersive[listIndex, 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] - 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 + phi = phi + updatecoeffsdispersive[material, p * 3].real * Tz[p, i, j, k].real + Tz[p, i, j, k] = updatecoeffsdispersive[material, 1 + (p * 3)] * Tz[p, i, j, k] + updatecoeffsdispersive[material, 2 + (p * 3)] * Ez[i, j, k] + Ez[i, j, k] = updatecoeffsE[material, 0] * Ez[i, j, k] + updatecoeffsE[material, 1] * (Hy[i, j, k] - Hy[i - 1, j, k]) - updatecoeffsE[material, 2] * (Hx[i, j, k] - Hx[i, j - 1, k]) - updatecoeffsE[material, 4] * phi -cpdef 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. 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 """ - cdef int i, j, k, listIndex, p + cdef Py_ssize_t i, j, k, p + cdef int material if nx == 1 or ny == 1: pass else: - for i in prange(1, nx, nogil=True, schedule='static', 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 k in range(0, nz): - listIndex = ID[2, i, j, k] + material = ID[2, i, j, k] 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. 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 """ - cdef int i, j, k, listIndex + cdef Py_ssize_t i, j, k + cdef int material cdef float phi = 0.0 if nx == 1 or ny == 1: pass else: - for i in prange(1, nx, nogil=True, schedule='static', 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 k in range(0, nz): - listIndex = ID[2, i, j, k] - phi = updatecoeffsdispersive[listIndex, 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] - 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 + material = ID[2, i, j, k] + phi = updatecoeffsdispersive[material, 0].real * Tz[0, i, j, k].real + 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[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. 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 """ - cdef int i, j, k, listIndex + cdef Py_ssize_t i, j, k + cdef int material if nx == 1 or ny == 1: pass else: - for i in prange(1, nx, nogil=True, schedule='static', 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 k in range(0, nz): - listIndex = ID[2, i, j, k] - Tz[0, i, j, k] = Tz[0, i, j, k] - updatecoeffsdispersive[listIndex, 2] * Ez[i, j, k] + material = ID[2, 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 # ######################################### -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. 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 """ - cdef int i, j, k, listIndex + cdef Py_ssize_t i, j, k + cdef int material if nx == 1: pass 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 k in range(0, nz): - listIndex = 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]) + material = ID[3, i, j, k] + Hx[i, j, k] = updatecoeffsH[material, 0] * Hx[i, j, k] - updatecoeffsH[material, 2] * (Ez[i, j + 1, k] - Ez[i, j, k]) + updatecoeffsH[material, 3] * (Ey[i, j, k + 1] - Ey[i, j, k]) ######################################### # 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. 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 """ - cdef int i, j, k, listIndex + cdef Py_ssize_t i, j, k + cdef int material if ny == 1: pass 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 k in range(0, nz): - listIndex = 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]) + material = ID[4, i, j, k] + Hy[i, j, k] = updatecoeffsH[material, 0] * Hy[i, j, k] - updatecoeffsH[material, 3] * (Ex[i, j, k + 1] - Ex[i, j, k]) + updatecoeffsH[material, 1] * (Ez[i + 1, j, k] - Ez[i, j, k]) ######################################### # 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. 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 """ - cdef int i, j, k, listIndex + cdef Py_ssize_t i, j, k + cdef int material if nz == 1: pass 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 k in range(1, nz): - listIndex = 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]) + material = ID[5, 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]) diff --git a/gprMax/pml_1order_update.pyx b/gprMax/pml_1order_update.pyx index 85ff3de4..97924d5d 100644 --- a/gprMax/pml_1order_update.pyx +++ b/gprMax/pml_1order_update.pyx @@ -25,7 +25,7 @@ from gprMax.constants cimport floattype_t, complextype_t ############################################# # 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. 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 """ - 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 nx = xf - xs ny = yf - ys 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 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 - 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 -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. 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 """ - 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 nx = xf - xs ny = yf - ys 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 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 - 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 -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. 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 """ - 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 nx = xf - xs ny = yf - ys 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 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 - 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 -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. 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 """ - 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 nx = xf - xs ny = yf - ys 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 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 - 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 ############################################# # 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. 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 """ - 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 nx = xf - xs ny = yf - ys 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 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 - 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 -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. 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 """ - 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 nx = xf - xs ny = yf - ys 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 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 - 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 -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. 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 """ - 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 nx = xf - xs ny = yf - ys 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 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 - 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 -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. 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 """ - 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 nx = xf - xs ny = yf - ys 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 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 - 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 ############################################# # 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. 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 """ - 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 nx = xf - xs ny = yf - ys 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 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 - 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 -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. 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 """ - 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 nx = xf - xs ny = yf - ys 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 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 - 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 -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. 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 """ - 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 nx = xf - xs ny = yf - ys 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 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 - 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 -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. 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 """ - 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 nx = xf - xs ny = yf - ys 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 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 - 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 ############################################# # 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. 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 """ - 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 nx = xf - xs ny = yf - ys 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 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 - 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 -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. 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 """ - 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 nx = xf - xs ny = yf - ys 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 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 - 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 -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. 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 """ - 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 nx = xf - xs ny = yf - ys 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 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 - 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 -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. 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 """ - 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 nx = xf - xs ny = yf - ys 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 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 - 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 ############################################# # 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. 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 """ - 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 nx = xf - xs ny = yf - ys 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 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 - 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 -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. 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 """ - 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 nx = xf - xs ny = yf - ys 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 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 - 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 -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. 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 """ - 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 nx = xf - xs ny = yf - ys 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 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 - 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 -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. 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 """ - 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 nx = xf - xs ny = yf - ys 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 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 - 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 ############################################# # 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. 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 """ - 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 nx = xf - xs ny = yf - ys 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 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 - 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 -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. 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 """ - 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 nx = xf - xs ny = yf - ys 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 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 - 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 -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. 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 """ - 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 nx = xf - xs ny = yf - ys 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 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 - 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 -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. 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 """ - 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 nx = xf - xs ny = yf - ys 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 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 - 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 diff --git a/gprMax/pml_2order_update.pyx b/gprMax/pml_2order_update.pyx index 8d8f226f..4883c978 100644 --- a/gprMax/pml_2order_update.pyx +++ b/gprMax/pml_2order_update.pyx @@ -25,7 +25,7 @@ from gprMax.constants cimport floattype_t, complextype_t ############################################# # 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. 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 """ - 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 nx = xf - xs ny = yf - ys 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 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 - 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[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. 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 """ - 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 nx = xf - xs ny = yf - ys 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 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 - 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[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. 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 """ - 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 nx = xf - xs ny = yf - ys 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 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 - 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[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. 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 """ - 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 nx = xf - xs ny = yf - ys 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 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 - 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[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 # ############################################# -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. 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 """ - 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 nx = xf - xs ny = yf - ys 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 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 - 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[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. 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 """ - 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 nx = xf - xs ny = yf - ys 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 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 - 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[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. 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 """ - 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 nx = xf - xs ny = yf - ys 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 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 - 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[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. 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 """ - 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 nx = xf - xs ny = yf - ys 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 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 - 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[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 # ############################################# -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. 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 """ - 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 nx = xf - xs ny = yf - ys 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 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 - 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[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. 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 """ - 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 nx = xf - xs ny = yf - ys 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 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 - 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[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. 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 """ - 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 nx = xf - xs ny = yf - ys 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 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 - 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[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. 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 """ - 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 nx = xf - xs ny = yf - ys 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 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 - 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[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 # ############################################# -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. 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 """ - 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 nx = xf - xs ny = yf - ys 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 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 - 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[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. 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 """ - 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 nx = xf - xs ny = yf - ys 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 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 - 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[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. 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 """ - 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 nx = xf - xs ny = yf - ys 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 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 - 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[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. 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 """ - 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 nx = xf - xs ny = yf - ys 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 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 - 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[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 # ############################################# -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. 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 """ - 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 nx = xf - xs ny = yf - ys 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 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 - 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[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. 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 """ - 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 nx = xf - xs ny = yf - ys 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 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 - 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[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. 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 """ - 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 nx = xf - xs ny = yf - ys 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 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 - 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[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. 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 """ - 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 nx = xf - xs ny = yf - ys 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 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 - 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[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 # ############################################# -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. 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 """ - 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 nx = xf - xs ny = yf - ys 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 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 - 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[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. 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 """ - 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 nx = xf - xs ny = yf - ys 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 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 - 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[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. 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 """ - 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 nx = xf - xs ny = yf - ys 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 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 - 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[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. 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 """ - 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 nx = xf - xs ny = yf - ys 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 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 - 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[0, i, j, k] = RE[0, j] * HPhi[0, i, j, k] - RF[0, j] * dEx