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