diff --git a/gprMax/pml_1order_update.pyx b/gprMax/pml_1order_update.pyx index 4f016418..52c68c94 100644 --- a/gprMax/pml_1order_update.pyx +++ b/gprMax/pml_1order_update.pyx @@ -23,668 +23,529 @@ from cython.parallel import prange from gprMax.constants cimport floattype_t, complextype_t -############################################# -# Electric field PML updates - Ex component # -############################################# -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. +############################################ +# Electric field PML updates - xminus slab # +############################################ +cpdef void update_pml_1order_electric_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] Ex, floattype_t[:, :, ::1] Ey, floattype_t[:, :, ::1] Ez, floattype_t[:, :, ::1] Hx, floattype_t[:, :, ::1] Hy, floattype_t[:, :, ::1] Hz, floattype_t[:, :, :, ::1] EPhi1, floattype_t[:, :, :, ::1] EPhi2, floattype_t[:, ::1] ERA, floattype_t[:, ::1] ERB, floattype_t[:, ::1] ERE, floattype_t[:, ::1] ERF, float d): + """This function updates the Ey and Ez field components for the xminus slab. - Args: + Args: xs, xf, ys, yf, zs, zf (int): Cell coordinates of entire box nthreads (int): Number of threads to use updatecoeffs, ID, E, H (memoryviews): Access to update coeffients, ID and field component arrays - EPhi, HPhi, RA, RB, RE, RF (memoryviews): Access to PML coefficient arrays - dx, dy, dz (float): Spatial discretisation - """ + EPhi, HPhi, ERA, ERB, ERE, ERF (memoryviews): Access to PML coefficient arrays + d (float): Spatial discretisation, e.g. dx, dy or dz + """ - 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', num_threads=nthreads): - for j in range(0, ny): - for k in range(0, nz): - 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[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 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: - xs, xf, ys, yf, zs, zf (int): Cell coordinates of entire box - nthreads (int): Number of threads to use - updatecoeffs, ID, E, H (memoryviews): Access to update coeffients, ID and field component arrays - EPhi, HPhi, RA, RB, RE, RF (memoryviews): Access to PML coefficient arrays - dx, dy, dz (float): Spatial discretisation - """ - - 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', num_threads=nthreads): - for j in range(0, ny): - for k in range(0, nz): - 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[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 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: - xs, xf, ys, yf, zs, zf (int): Cell coordinates of entire box - nthreads (int): Number of threads to use - updatecoeffs, ID, E, H (memoryviews): Access to update coeffients, ID and field component arrays - EPhi, HPhi, RA, RB, RE, RF (memoryviews): Access to PML coefficient arrays - dx, dy, dz (float): Spatial discretisation - """ - - 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', num_threads=nthreads): - for j in range(0, ny): - for k in range(0, nz): - 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[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 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: - xs, xf, ys, yf, zs, zf (int): Cell coordinates of entire box - nthreads (int): Number of threads to use - updatecoeffs, ID, E, H (memoryviews): Access to update coeffients, ID and field component arrays - EPhi, HPhi, RA, RB, RE, RF (memoryviews): Access to PML coefficient arrays - dx, dy, dz (float): Spatial discretisation - """ - - 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', num_threads=nthreads): - for j in range(0, ny): - for k in range(0, nz): - 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[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 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: - xs, xf, ys, yf, zs, zf (int): Cell coordinates of entire box - nthreads (int): Number of threads to use - updatecoeffs, ID, E, H (memoryviews): Access to update coeffients, ID and field component arrays - EPhi, HPhi, RA, RB, RE, RF (memoryviews): Access to PML coefficient arrays - dx, dy, dz (float): Spatial discretisation - """ - - 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', num_threads=nthreads): - for j in range(0, ny): - for k in range(0, nz): - 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[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 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: - xs, xf, ys, yf, zs, zf (int): Cell coordinates of entire box - nthreads (int): Number of threads to use - updatecoeffs, ID, E, H (memoryviews): Access to update coeffients, ID and field component arrays - EPhi, HPhi, RA, RB, RE, RF (memoryviews): Access to PML coefficient arrays - dx, dy, dz (float): Spatial discretisation - """ - - 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', num_threads=nthreads): - for j in range(0, ny): - for k in range(0, nz): - 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[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 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: - xs, xf, ys, yf, zs, zf (int): Cell coordinates of entire box - nthreads (int): Number of threads to use - updatecoeffs, ID, E, H (memoryviews): Access to update coeffients, ID and field component arrays - EPhi, HPhi, RA, RB, RE, RF (memoryviews): Access to PML coefficient arrays - dx, dy, dz (float): Spatial discretisation - """ - - 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', num_threads=nthreads): - for j in range(0, ny): - for k in range(0, nz): - 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[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 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: - xs, xf, ys, yf, zs, zf (int): Cell coordinates of entire box - nthreads (int): Number of threads to use - updatecoeffs, ID, E, H (memoryviews): Access to update coeffients, ID and field component arrays - EPhi, HPhi, RA, RB, RE, RF (memoryviews): Access to PML coefficient arrays - dx, dy, dz (float): Spatial discretisation - """ - - 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', num_threads=nthreads): - for j in range(0, ny): - for k in range(0, nz): - 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[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 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: - xs, xf, ys, yf, zs, zf (int): Cell coordinates of entire box - nthreads (int): Number of threads to use - updatecoeffs, ID, E, H (memoryviews): Access to update coeffients, ID and field component arrays - EPhi, HPhi, RA, RB, RE, RF (memoryviews): Access to PML coefficient arrays - dx, dy, dz (float): Spatial discretisation - """ - - 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', num_threads=nthreads): - for j in range(0, ny): - for k in range(0, nz): - 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[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 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: - xs, xf, ys, yf, zs, zf (int): Cell coordinates of entire box - nthreads (int): Number of threads to use - updatecoeffs, ID, E, H (memoryviews): Access to update coeffients, ID and field component arrays - EPhi, HPhi, RA, RB, RE, RF (memoryviews): Access to PML coefficient arrays - dx, dy, dz (float): Spatial discretisation - """ - - 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', num_threads=nthreads): - for j in range(0, ny): - for k in range(0, nz): - 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[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 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: - xs, xf, ys, yf, zs, zf (int): Cell coordinates of entire box - nthreads (int): Number of threads to use - updatecoeffs, ID, E, H (memoryviews): Access to update coeffients, ID and field component arrays - EPhi, HPhi, RA, RB, RE, RF (memoryviews): Access to PML coefficient arrays - dx, dy, dz (float): Spatial discretisation - """ - - 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', num_threads=nthreads): - for j in range(0, ny): - for k in range(0, nz): - 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[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 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: - xs, xf, ys, yf, zs, zf (int): Cell coordinates of entire box - nthreads (int): Number of threads to use - updatecoeffs, ID, E, H (memoryviews): Access to update coeffients, ID and field component arrays - EPhi, HPhi, RA, RB, RE, RF (memoryviews): Access to PML coefficient arrays - dx, dy, dz (float): Spatial discretisation - """ - - 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', num_threads=nthreads): - for j in range(0, ny): - for k in range(0, nz): - 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[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 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: - xs, xf, ys, yf, zs, zf (int): Cell coordinates of entire box - nthreads (int): Number of threads to use - updatecoeffs, ID, E, H (memoryviews): Access to update coeffients, ID and field component arrays - EPhi, HPhi, RA, RB, RE, RF (memoryviews): Access to PML coefficient arrays - dx, dy, dz (float): Spatial discretisation - """ - - 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', num_threads=nthreads): - for j in range(0, ny): - for k in range(0, nz): - 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[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 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: - xs, xf, ys, yf, zs, zf (int): Cell coordinates of entire box - nthreads (int): Number of threads to use - updatecoeffs, ID, E, H (memoryviews): Access to update coeffients, ID and field component arrays - EPhi, HPhi, RA, RB, RE, RF (memoryviews): Access to PML coefficient arrays - dx, dy, dz (float): Spatial discretisation - """ - - 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', num_threads=nthreads): - for j in range(0, ny): - for k in range(0, nz): - 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[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 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: - xs, xf, ys, yf, zs, zf (int): Cell coordinates of entire box - nthreads (int): Number of threads to use - updatecoeffs, ID, E, H (memoryviews): Access to update coeffients, ID and field component arrays - EPhi, HPhi, RA, RB, RE, RF (memoryviews): Access to PML coefficient arrays - dx, dy, dz (float): Spatial discretisation - """ - - 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', num_threads=nthreads): - for j in range(0, ny): - for k in range(0, nz): - 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[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 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: - xs, xf, ys, yf, zs, zf (int): Cell coordinates of entire box - nthreads (int): Number of threads to use - updatecoeffs, ID, E, H (memoryviews): Access to update coeffients, ID and field component arrays - EPhi, HPhi, RA, RB, RE, RF (memoryviews): Access to PML coefficient arrays - dx, dy, dz (float): Spatial discretisation - """ - - 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', num_threads=nthreads): - for j in range(0, ny): - for k in range(0, nz): - 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[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 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: - xs, xf, ys, yf, zs, zf (int): Cell coordinates of entire box - nthreads (int): Number of threads to use - updatecoeffs, ID, E, H (memoryviews): Access to update coeffients, ID and field component arrays - EPhi, HPhi, RA, RB, RE, RF (memoryviews): Access to PML coefficient arrays - dx, dy, dz (float): Spatial discretisation - """ - - 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', num_threads=nthreads): - for j in range(0, ny): - for k in range(0, nz): - 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[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 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: - xs, xf, ys, yf, zs, zf (int): Cell coordinates of entire box - nthreads (int): Number of threads to use - updatecoeffs, ID, E, H (memoryviews): Access to update coeffients, ID and field component arrays - EPhi, HPhi, RA, RB, RE, RF (memoryviews): Access to PML coefficient arrays - dx, dy, dz (float): Spatial discretisation - """ - - 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', num_threads=nthreads): - for j in range(0, ny): - for k in range(0, nz): - 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[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 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: - xs, xf, ys, yf, zs, zf (int): Cell coordinates of entire box - nthreads (int): Number of threads to use - updatecoeffs, ID, E, H (memoryviews): Access to update coeffients, ID and field component arrays - EPhi, HPhi, RA, RB, RE, RF (memoryviews): Access to PML coefficient arrays - dx, dy, dz (float): Spatial discretisation - """ - - 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', num_threads=nthreads): - for j in range(0, ny): - for k in range(0, nz): - 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[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 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: - xs, xf, ys, yf, zs, zf (int): Cell coordinates of entire box - nthreads (int): Number of threads to use - updatecoeffs, ID, E, H (memoryviews): Access to update coeffients, ID and field component arrays - EPhi, HPhi, RA, RB, RE, RF (memoryviews): Access to PML coefficient arrays - dx, dy, dz (float): Spatial discretisation - """ - - cdef Py_ssize_t i, j, k - cdef int nx, ny, nz, material - cdef float dEx + cdef Py_ssize_t i, j, k, ii, jj, kk + cdef int nx, ny, nz, materialEy, materialEz + cdef float dx, dHy, dHz, RA0, RB0, RE0, RF0 + dx = d nx = xf - xs ny = yf - ys nz = zf - zs for i in prange(0, nx, nogil=True, schedule='static', num_threads=nthreads): + RA0 = (ERA[0, i] - 1) + RB0 = ERB[0, i] + RE0 = ERE[0, i] + RF0 = ERF[0, i] + ii = xf - i for j in range(0, ny): + jj = j + ys for k in range(0, nz): - 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[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 + kk = k + zs + # Ey + materialEy = ID[1, ii, jj, kk] + dHz = (Hz[ii, jj, kk] - Hz[ii - 1, jj, kk]) / dx + Ey[ii, jj, kk] = Ey[ii, jj, kk] - updatecoeffsE[materialEy, 4] * (RA0 * dHz + RB0 * EPhi1[0, i, j, k]) + EPhi1[0, i, j, k] = RE0 * EPhi1[0, i, j, k] - RF0 * dHz + # Ez + materialEz = ID[2, ii, jj, kk] + dHy = (Hy[ii, jj, kk] - Hy[ii - 1, jj, kk]) / dx + Ez[ii, jj, kk] = Ez[ii, jj, kk] + updatecoeffsE[materialEz, 4] * (RA0 * dHy + RB0 * EPhi2[0, i, j, k]) + EPhi2[0, i, j, k] = RE0 * EPhi2[0, i, j, k] - RF0 * dHy -############################################# -# Magnetic field PML updates - Hz component # -############################################# -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. +########################################### +# Electric field PML updates - xplus slab # +########################################### +cpdef void update_pml_1order_electric_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] Ex, floattype_t[:, :, ::1] Ey, floattype_t[:, :, ::1] Ez, floattype_t[:, :, ::1] Hx, floattype_t[:, :, ::1] Hy, floattype_t[:, :, ::1] Hz, floattype_t[:, :, :, ::1] EPhi1, floattype_t[:, :, :, ::1] EPhi2, floattype_t[:, ::1] ERA, floattype_t[:, ::1] ERB, floattype_t[:, ::1] ERE, floattype_t[:, ::1] ERF, float d): + """This function updates the Ey and Ez field components for the xplus slab. Args: xs, xf, ys, yf, zs, zf (int): Cell coordinates of entire box nthreads (int): Number of threads to use updatecoeffs, ID, E, H (memoryviews): Access to update coeffients, ID and field component arrays - EPhi, HPhi, RA, RB, RE, RF (memoryviews): Access to PML coefficient arrays - dx, dy, dz (float): Spatial discretisation + EPhi, HPhi, ERA, ERB, ERE, ERF (memoryviews): Access to PML coefficient arrays + d (float): Spatial discretisation, e.g. dx, dy or dz """ - cdef Py_ssize_t i, j, k - cdef int nx, ny, nz, material - cdef float dEy + cdef Py_ssize_t i, j, k, ii, jj, kk + cdef int nx, ny, nz, materialEy, materialEz + cdef float dx, dHy, dHz, RA0, RB0, RE0, RF0 + dx = d nx = xf - xs ny = yf - ys nz = zf - zs for i in prange(0, nx, nogil=True, schedule='static', num_threads=nthreads): + RA0 = (ERA[0, i] - 1) + RB0 = ERB[0, i] + RE0 = ERE[0, i] + RF0 = ERF[0, i] + ii = i + xs for j in range(0, ny): + jj = j + ys for k in range(0, nz): - 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[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 + kk = k + zs + # Ey + materialEy = ID[1, ii, jj, kk] + dHz = (Hz[ii, jj, kk] - Hz[ii - 1, jj, kk]) / dx + Ey[ii, jj, kk] = Ey[ii, jj, kk] - updatecoeffsE[materialEy, 4] * (RA0 * dHz + RB0 * EPhi1[0, i, j, k]) + EPhi1[0, i, j, k] = RE0 * EPhi1[0, i, j, k] - RF0 * dHz + # Ez + materialEz = ID[2, ii, jj, kk] + dHy = (Hy[ii, jj, kk] - Hy[ii - 1, jj, kk]) / dx + Ez[ii, jj, kk] = Ez[ii, jj, kk] + updatecoeffsE[materialEz, 4] * (RA0 * dHy + RB0 * EPhi2[0, i, j, k]) + EPhi2[0, i, j, k] = RE0 * EPhi2[0, i, j, k] - RF0 * dHy -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. +############################################ +# Electric field PML updates - yminus slab # +############################################ +cpdef void update_pml_1order_electric_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] Ey, floattype_t[:, :, ::1] Ez, floattype_t[:, :, ::1] Hx, floattype_t[:, :, ::1] Hy, floattype_t[:, :, ::1] Hz, floattype_t[:, :, :, ::1] EPhi1, floattype_t[:, :, :, ::1] EPhi2, floattype_t[:, ::1] ERA, floattype_t[:, ::1] ERB, floattype_t[:, ::1] ERE, floattype_t[:, ::1] ERF, float d): + """This function updates the Ex and Ez field components for the yminus slab. Args: xs, xf, ys, yf, zs, zf (int): Cell coordinates of entire box nthreads (int): Number of threads to use updatecoeffs, ID, E, H (memoryviews): Access to update coeffients, ID and field component arrays - EPhi, HPhi, RA, RB, RE, RF (memoryviews): Access to PML coefficient arrays - dx, dy, dz (float): Spatial discretisation + EPhi, HPhi, ERA, ERB, ERE, ERF (memoryviews): Access to PML coefficient arrays + d (float): Spatial discretisation, e.g. dx, dy or dz """ - cdef Py_ssize_t i, j, k - cdef int nx, ny, nz, material - cdef float dEy + cdef Py_ssize_t i, j, k, ii, jj, kk + cdef int nx, ny, nz, materialEx, materialEz + cdef float dy, dHx, dHz, RA0, RB0, RE0, RF0 + dy = d nx = xf - xs ny = yf - ys nz = zf - zs for i in prange(0, nx, nogil=True, schedule='static', num_threads=nthreads): + ii = i + xs for j in range(0, ny): + jj = yf - j + RA0 = (ERA[0, j] - 1) + RB0 = ERB[0, j] + RE0 = ERE[0, j] + RF0 = ERF[0, j] for k in range(0, nz): - 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[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 + kk = k + zs + # Ex + materialEx = ID[0, ii, jj, kk] + dHz = (Hz[ii, jj, kk] - Hz[ii, jj - 1, kk]) / dy + Ex[ii, jj, kk] = Ex[ii, jj, kk] + updatecoeffsE[materialEx, 4] * (RA0 * dHz + RB0 * EPhi1[0, i, j, k]) + EPhi1[0, i, j, k] = RE0 * EPhi1[0, i, j, k] - RF0 * dHz + # Ez + materialEz = ID[2, ii, jj, kk] + dHx = (Hx[ii, jj, kk] - Hx[ii, jj - 1, kk]) / dy + Ez[ii, jj, kk] = Ez[ii, jj, kk] - updatecoeffsE[materialEz, 4] * (RA0 * dHx + RB0 * EPhi2[0, i, j, k]) + EPhi2[0, i, j, k] = RE0 * EPhi2[0, i, j, k] - RF0 * dHx -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. +########################################### +# Electric field PML updates - yplus slab # +########################################### +cpdef void update_pml_1order_electric_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] Ey, floattype_t[:, :, ::1] Ez, floattype_t[:, :, ::1] Hx, floattype_t[:, :, ::1] Hy, floattype_t[:, :, ::1] Hz, floattype_t[:, :, :, ::1] EPhi1, floattype_t[:, :, :, ::1] EPhi2, floattype_t[:, ::1] ERA, floattype_t[:, ::1] ERB, floattype_t[:, ::1] ERE, floattype_t[:, ::1] ERF, float d): + """This function updates the Ex and Ez field components for the yplus slab. Args: xs, xf, ys, yf, zs, zf (int): Cell coordinates of entire box nthreads (int): Number of threads to use updatecoeffs, ID, E, H (memoryviews): Access to update coeffients, ID and field component arrays - EPhi, HPhi, RA, RB, RE, RF (memoryviews): Access to PML coefficient arrays - dx, dy, dz (float): Spatial discretisation + EPhi, HPhi, ERA, ERB, ERE, ERF (memoryviews): Access to PML coefficient arrays + d (float): Spatial discretisation, e.g. dx, dy or dz """ - cdef Py_ssize_t i, j, k - cdef int nx, ny, nz, material - cdef float dEx + cdef Py_ssize_t i, j, k, ii, jj, kk + cdef int nx, ny, nz, materialEx, materialEz + cdef float dy, dHx, dHz, RA0, RB0, RE0, RF0 + dy = d nx = xf - xs ny = yf - ys nz = zf - zs for i in prange(0, nx, nogil=True, schedule='static', num_threads=nthreads): + ii = i + xs for j in range(0, ny): + jj = j + ys + RA0 = (ERA[0, j] - 1) + RB0 = ERB[0, j] + RE0 = ERE[0, j] + RF0 = ERF[0, j] for k in range(0, nz): - 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[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 + kk = k + zs + # Ex + materialEx = ID[0, ii, jj, kk] + dHz = (Hz[ii, jj, kk] - Hz[ii, jj - 1, kk]) / dy + Ex[ii, jj, kk] = Ex[ii, jj, kk] + updatecoeffsE[materialEx, 4] * (RA0 * dHz + RB0 * EPhi1[0, i, j, k]) + EPhi1[0, i, j, k] = RE0 * EPhi1[0, i, j, k] - RF0 * dHz + # Ez + materialEz = ID[2, ii, jj, kk] + dHx = (Hx[ii, jj, kk] - Hx[ii, jj - 1, kk]) / dy + Ez[ii, jj, kk] = Ez[ii, jj, kk] - updatecoeffsE[materialEz, 4] * (RA0 * dHx + RB0 * EPhi2[0, i, j, k]) + EPhi2[0, i, j, k] = RE0 * EPhi2[0, i, j, k] - RF0 * dHx -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. +############################################ +# Electric field PML updates - zminus slab # +############################################ +cpdef void update_pml_1order_electric_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] Ey, floattype_t[:, :, ::1] Ez, floattype_t[:, :, ::1] Hx, floattype_t[:, :, ::1] Hy, floattype_t[:, :, ::1] Hz, floattype_t[:, :, :, ::1] EPhi1, floattype_t[:, :, :, ::1] EPhi2, floattype_t[:, ::1] ERA, floattype_t[:, ::1] ERB, floattype_t[:, ::1] ERE, floattype_t[:, ::1] ERF, float d): + """This function updates the Ex and Ey field components for the zminus slab. + + Args: + xs, xf, ys, yf, zs, zf (int): Cell coordinates of entire box + nthreads (int): Number of threads to use + updatecoeffs, ID, E, H (memoryviews): Access to update coeffients, ID and field component arrays + EPhi, HPhi, ERA, ERB, ERE, ERF (memoryviews): Access to PML coefficient arrays + d (float): Spatial discretisation, e.g. dx, dy or dz + """ + + cdef Py_ssize_t i, j, k, ii, jj, kk + cdef int nx, ny, nz, materialEx, materialEy + cdef float dz, dHx, dHy, RA0, RB0, RE0, RF0 + dz = d + nx = xf - xs + ny = yf - ys + nz = zf - zs + + for i in prange(0, nx, nogil=True, schedule='static', num_threads=nthreads): + ii = i + xs + for j in range(0, ny): + jj = j + ys + for k in range(0, nz): + kk = zf - k + RA0 = (ERA[0, k] - 1) + RB0 = ERB[0, k] + RE0 = ERE[0, k] + RF0 = ERF[0, k] + # Ex + materialEx = ID[0, ii, jj, kk] + dHy = (Hy[ii, jj, kk] - Hy[ii, jj, kk - 1]) / dz + Ex[ii, jj, kk] = Ex[ii, jj, kk] - updatecoeffsE[materialEx, 4] * (RA0 * dHy + RB0 * EPhi1[0, i, j, k]) + EPhi1[0, i, j, k] = RE0 * EPhi1[0, i, j, k] - RF0 * dHy + # Ey + materialEy = ID[1, ii, jj, kk] + dHx = (Hx[ii, jj, kk] - Hx[ii, jj, kk - 1]) / dz + Ey[ii, jj, kk] = Ey[ii, jj, kk] + updatecoeffsE[materialEy, 4] * (RA0 * dHx + RB0 * EPhi2[0, i, j, k]) + EPhi2[0, i, j, k] = RE0 * EPhi2[0, i, j, k] - RF0 * dHx + + +########################################### +# Electric field PML updates - zplus slab # +########################################### +cpdef void update_pml_1order_electric_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] Ey, floattype_t[:, :, ::1] Ez, floattype_t[:, :, ::1] Hx, floattype_t[:, :, ::1] Hy, floattype_t[:, :, ::1] Hz, floattype_t[:, :, :, ::1] EPhi1, floattype_t[:, :, :, ::1] EPhi2, floattype_t[:, ::1] ERA, floattype_t[:, ::1] ERB, floattype_t[:, ::1] ERE, floattype_t[:, ::1] ERF, float d): + """This function updates the Ex and Ey field components for the zplus slab. Args: xs, xf, ys, yf, zs, zf (int): Cell coordinates of entire box nthreads (int): Number of threads to use updatecoeffs, ID, E, H (memoryviews): Access to update coeffients, ID and field component arrays - EPhi, HPhi, RA, RB, RE, RF (memoryviews): Access to PML coefficient arrays - dx, dy, dz (float): Spatial discretisation + EPhi, HPhi, ERA, ERB, ERE, ERF (memoryviews): Access to PML coefficient arrays + d (float): Spatial discretisation, e.g. dx, dy or dz """ - cdef Py_ssize_t i, j, k - cdef int nx, ny, nz, material - cdef float dEx + cdef Py_ssize_t i, j, k, ii, jj, kk + cdef int nx, ny, nz, materialEx, materialEy + cdef float dz, dHx, dHy, RA0, RB0, RE0, RF0 + dz = d nx = xf - xs ny = yf - ys nz = zf - zs for i in prange(0, nx, nogil=True, schedule='static', num_threads=nthreads): + ii = i + xs for j in range(0, ny): + jj = j + ys for k in range(0, nz): - 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[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 + kk = k + zs + RA0 = (ERA[0, k] - 1) + RB0 = ERB[0, k] + RE0 = ERE[0, k] + RF0 = ERF[0, k] + # Ex + materialEx = ID[0, ii, jj, kk] + dHy = (Hy[ii, jj, kk] - Hy[ii, jj, kk - 1]) / dz + Ex[ii, jj, kk] = Ex[ii, jj, kk] - updatecoeffsE[materialEx, 4] * (RA0 * dHy + RB0 * EPhi1[0, i, j, k]) + EPhi1[0, i, j, k] = RE0 * EPhi1[0, i, j, k] - RF0 * dHy + # Ey + materialEy = ID[1, ii, jj, kk] + dHx = (Hx[ii, jj, kk] - Hx[ii, jj, kk - 1]) / dz + Ey[ii, jj, kk] = Ey[ii, jj, kk] + updatecoeffsE[materialEy, 4] * (RA0 * dHx + RB0 * EPhi2[0, i, j, k]) + EPhi2[0, i, j, k] = RE0 * EPhi2[0, i, j, k] - RF0 * dHx + +############################################ +# Magnetic field PML updates - xminus slab # +############################################ +cpdef void update_pml_1order_magnetic_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] Ex, floattype_t[:, :, ::1] Ey, floattype_t[:, :, ::1] Ez, floattype_t[:, :, ::1] Hx, floattype_t[:, :, ::1] Hy, floattype_t[:, :, ::1] Hz, floattype_t[:, :, :, ::1] HPhi1, floattype_t[:, :, :, ::1] HPhi2, floattype_t[:, ::1] HRA, floattype_t[:, ::1] HRB, floattype_t[:, ::1] HRE, floattype_t[:, ::1] HRF, float d): + """This function updates the Hy and Hz field components for the xminus slab. + + Args: + xs, xf, ys, yf, zs, zf (int): Cell coordinates of entire box + nthreads (int): Number of threads to use + updatecoeffs, ID, E, H (memoryviews): Access to update coeffients, ID and field component arrays + EPhi, HPhi, HRA, HRB, ERE, HRF (memoryviews): Access to PML coefficient arrays + d (float): Spatial discretisation, e.g. dx, dy or dz + """ + + cdef Py_ssize_t i, j, k, ii, jj, kk + cdef int nx, ny, nz, materialHy, materialHz + cdef float dx, dEy, dEz, RA0, RB0, RE0, RF0 + dx = d + nx = xf - xs + ny = yf - ys + nz = zf - zs + + for i in prange(0, nx, nogil=True, schedule='static', num_threads=nthreads): + ii = xf - (i + 1) + RA0 = (HRA[0, i] - 1) + RB0 = HRB[0, i] + RE0 = HRE[0, i] + RF0 = HRF[0, i] + for j in range(0, ny): + jj = j + ys + for k in range(0, nz): + kk = k + zs + # Hy + materialHy = ID[4, ii, jj, kk] + dEz = (Ez[ii + 1, jj, kk] - Ez[ii, jj, kk]) / dx + Hy[ii, jj, kk] = Hy[ii, jj, kk] + updatecoeffsH[materialHy, 4] * (RA0 * dEz + RB0 * HPhi1[0, i, j, k]) + HPhi1[0, i, j, k] = RE0 * HPhi1[0, i, j, k] - RF0 * dEz + # Hz + materialHz = ID[5, ii, jj, kk] + dEy = (Ey[ii + 1, jj, kk] - Ey[ii, jj, kk]) / dx + Hz[ii, jj, kk] = Hz[ii, jj, kk] - updatecoeffsH[materialHz, 4] * (RA0 * dEy + RB0 * HPhi2[0, i, j, k]) + HPhi2[0, i, j, k] = RE0 * HPhi2[0, i, j, k] - RF0 * dEy + + +########################################### +# Magnetic field PML updates - xplus slab # +########################################### +cpdef void update_pml_1order_magnetic_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] Ex, floattype_t[:, :, ::1] Ey, floattype_t[:, :, ::1] Ez, floattype_t[:, :, ::1] Hx, floattype_t[:, :, ::1] Hy, floattype_t[:, :, ::1] Hz, floattype_t[:, :, :, ::1] HPhi1, floattype_t[:, :, :, ::1] HPhi2, floattype_t[:, ::1] HRA, floattype_t[:, ::1] HRB, floattype_t[:, ::1] HRE, floattype_t[:, ::1] HRF, float d): + """This function updates the Hy and Hz field components for the xplus slab. + + Args: + xs, xf, ys, yf, zs, zf (int): Cell coordinates of entire box + nthreads (int): Number of threads to use + updatecoeffs, ID, E, H (memoryviews): Access to update coeffients, ID and field component arrays + EPhi, HPhi, HRA, HRB, HRE, HRF (memoryviews): Access to PML coefficient arrays + d (float): Spatial discretisation, e.g. dx, dy or dz + """ + + cdef Py_ssize_t i, j, k, ii, jj, kk + cdef int nx, ny, nz, materialHy, materialHz + cdef float dx, dEy, dEz, RA0, RB0, RE0, RF0 + dx = d + nx = xf - xs + ny = yf - ys + nz = zf - zs + + for i in prange(0, nx, nogil=True, schedule='static', num_threads=nthreads): + ii = i + xs + RA0 = (HRA[0, i] - 1) + RB0 = HRB[0, i] + RE0 = HRE[0, i] + RF0 = HRF[0, i] + for j in range(0, ny): + jj = j + ys + for k in range(0, nz): + kk = k + zs + # Hy + materialHy = ID[4, ii, jj, kk] + dEz = (Ez[ii + 1, jj, kk] - Ez[ii, jj, kk]) / dx + Hy[ii, jj, kk] = Hy[ii, jj, kk] + updatecoeffsH[materialHy, 4] * (RA0 * dEz + RB0 * HPhi1[0, i, j, k]) + HPhi1[0, i, j, k] = RE0 * HPhi1[0, i, j, k] - RF0 * dEz + # Hz + materialHz = ID[5, ii, jj, kk] + dEy = (Ey[ii + 1, jj, kk] - Ey[ii, jj, kk]) / dx + Hz[ii, jj, kk] = Hz[ii, jj, kk] - updatecoeffsH[materialHz, 4] * (RA0 * dEy + RB0 * HPhi2[0, i, j, k]) + HPhi2[0, i, j, k] = RE0 * HPhi2[0, i, j, k] - RF0 * dEy + + +############################################ +# Magnetic field PML updates - yminus slab # +############################################ +cpdef void update_pml_1order_magnetic_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] Ex, floattype_t[:, :, ::1] Ey, floattype_t[:, :, ::1] Ez, floattype_t[:, :, ::1] Hx, floattype_t[:, :, ::1] Hy, floattype_t[:, :, ::1] Hz, floattype_t[:, :, :, ::1] HPhi1, floattype_t[:, :, :, ::1] HPhi2, floattype_t[:, ::1] HRA, floattype_t[:, ::1] HRB, floattype_t[:, ::1] HRE, floattype_t[:, ::1] HRF, float d): + """This function updates the Hx and Hz field components for the yminus slab. + + Args: + xs, xf, ys, yf, zs, zf (int): Cell coordinates of entire box + nthreads (int): Number of threads to use + updatecoeffs, ID, E, H (memoryviews): Access to update coeffients, ID and field component arrays + EPhi, HPhi, HRA, HRB, HRE, HRF (memoryviews): Access to PML coefficient arrays + d (float): Spatial discretisation, e.g. dx, dy or dz + """ + + cdef Py_ssize_t i, j, k, ii, jj, kk + cdef int nx, ny, nz, materialHx, materialHz + cdef float dy, dEx, dEz, RA0, RB0, RE0, RF0 + dy = d + nx = xf - xs + ny = yf - ys + nz = zf - zs + + for i in prange(0, nx, nogil=True, schedule='static', num_threads=nthreads): + ii = i + xs + for j in range(0, ny): + jj = yf - (j + 1) + RA0 = (HRA[0, j] - 1) + RB0 = HRB[0, j] + RE0 = HRE[0, j] + RF0 = HRF[0, j] + for k in range(0, nz): + kk = k + zs + # Hx + materialHx = ID[3, ii, jj, kk] + dEz = (Ez[ii, jj + 1, kk] - Ez[ii, jj, kk]) / dy + Hx[ii, jj, kk] = Hx[ii, jj, kk] - updatecoeffsH[materialHx, 4] * (RA0 * dEz + RB0 * HPhi1[0, i, j, k]) + HPhi1[0, i, j, k] = RE0 * HPhi1[0, i, j, k] - RF0 * dEz + # Hz + materialHz = ID[5, ii, jj, kk] + dEx = (Ex[ii, jj + 1, kk] - Ex[ii, jj, kk]) / dy + Hz[ii, jj, kk] = Hz[ii, jj, kk] + updatecoeffsH[materialHz, 4] * (RA0 * dEx + RB0 * HPhi2[0, i, j, k]) + HPhi2[0, i, j, k] = RE0 * HPhi2[0, i, j, k] - RF0 * dEx + + +########################################### +# Magnetic field PML updates - yplus slab # +########################################### +cpdef void update_pml_1order_magnetic_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] Ex, floattype_t[:, :, ::1] Ey, floattype_t[:, :, ::1] Ez, floattype_t[:, :, ::1] Hx, floattype_t[:, :, ::1] Hy, floattype_t[:, :, ::1] Hz, floattype_t[:, :, :, ::1] HPhi1, floattype_t[:, :, :, ::1] HPhi2, floattype_t[:, ::1] HRA, floattype_t[:, ::1] HRB, floattype_t[:, ::1] HRE, floattype_t[:, ::1] HRF, float d): + """This function updates the Hx and Hz field components for the yplus slab. + + Args: + xs, xf, ys, yf, zs, zf (int): Cell coordinates of entire box + nthreads (int): Number of threads to use + updatecoeffs, ID, E, H (memoryviews): Access to update coeffients, ID and field component arrays + EPhi, HPhi, HRA, HRB, HRE, HRF (memoryviews): Access to PML coefficient arrays + d (float): Spatial discretisation, e.g. dx, dy or dz + """ + + cdef Py_ssize_t i, j, k, ii, jj, kk + cdef int nx, ny, nz, materialHx, materialHz + cdef float dy, dEx, dEz, RA0, RB0, RE0, RF0 + dy = d + nx = xf - xs + ny = yf - ys + nz = zf - zs + + for i in prange(0, nx, nogil=True, schedule='static', num_threads=nthreads): + ii = i + xs + for j in range(0, ny): + jj = j + ys + RA0 = (HRA[0, j] - 1) + RB0 = HRB[0, j] + RE0 = HRE[0, j] + RF0 = HRF[0, j] + for k in range(0, nz): + kk = k + zs + # Hx + materialHx = ID[3, ii, jj, kk] + dEz = (Ez[ii, jj + 1, kk] - Ez[ii, jj, kk]) / dy + Hx[ii, jj, kk] = Hx[ii, jj, kk] - updatecoeffsH[materialHx, 4] * (RA0 * dEz + RB0 * HPhi1[0, i, j, k]) + HPhi1[0, i, j, k] = RE0 * HPhi1[0, i, j, k] - RF0 * dEz + # Hz + materialHz = ID[5, ii, jj, kk] + dEx = (Ex[ii, jj + 1, kk] - Ex[ii, jj, kk]) / dy + Hz[ii, jj, kk] = Hz[ii, jj, kk] + updatecoeffsH[materialHz, 4] * (RA0 * dEx + RB0 * HPhi2[0, i, j, k]) + HPhi2[0, i, j, k] = RE0 * HPhi2[0, i, j, k] - RF0 * dEx + + +############################################ +# Magnetic field PML updates - zminus slab # +############################################ +cpdef void update_pml_1order_magnetic_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] Ex, floattype_t[:, :, ::1] Ey, floattype_t[:, :, ::1] Ez, floattype_t[:, :, ::1] Hx, floattype_t[:, :, ::1] Hy, floattype_t[:, :, ::1] Hz, floattype_t[:, :, :, ::1] HPhi1, floattype_t[:, :, :, ::1] HPhi2, floattype_t[:, ::1] HRA, floattype_t[:, ::1] HRB, floattype_t[:, ::1] HRE, floattype_t[:, ::1] HRF, float d): + """This function updates the Hx and Hy field components for the zminus slab. + + Args: + xs, xf, ys, yf, zs, zf (int): Cell coordinates of entire box + nthreads (int): Number of threads to use + updatecoeffs, ID, E, H (memoryviews): Access to update coeffients, ID and field component arrays + EPhi, HPhi, HRA, HRB, HRE, HRF (memoryviews): Access to PML coefficient arrays + d (float): Spatial discretisation, e.g. dx, dy or dz + """ + + cdef Py_ssize_t i, j, k, ii, jj, kk + cdef int nx, ny, nz, materialHx, materialHy + cdef float dz, dEx, dEy, RA0, RB0, RE0, RF0 + dz = d + nx = xf - xs + ny = yf - ys + nz = zf - zs + + for i in prange(0, nx, nogil=True, schedule='static', num_threads=nthreads): + ii = i + xs + for j in range(0, ny): + jj = j + ys + for k in range(0, nz): + kk = zf - (k + 1) + RA0 = (HRA[0, k] - 1) + RB0 = HRB[0, k] + RE0 = HRE[0, k] + RF0 = HRF[0, k] + # Hx + materialHx = ID[3, ii, jj, kk] + dEy = (Ey[ii, jj, kk + 1] - Ey[ii, jj, kk]) / dz + Hx[ii, jj, kk] = Hx[ii, jj, kk] + updatecoeffsH[materialHx, 4] * (RA0 * dEy + RB0 * HPhi1[0, i, j, k]) + HPhi1[0, i, j, k] = RE0 * HPhi1[0, i, j, k] - RF0 * dEy + # Hy + materialHy = ID[4, ii, jj, kk] + dEx = (Ex[ii, jj, kk + 1] - Ex[ii, jj, kk]) / dz + Hy[ii, jj, kk] = Hy[ii, jj, kk] - updatecoeffsH[materialHy, 4] * (RA0 * dEx + RB0 * HPhi2[0, i, j, k]) + HPhi2[0, i, j, k] = RE0 * HPhi2[0, i, j, k] - RF0 * dEx + + +########################################### +# Magnetic field PML updates - zplus slab # +########################################### +cpdef void update_pml_1order_magnetic_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] Ex, floattype_t[:, :, ::1] Ey, floattype_t[:, :, ::1] Ez, floattype_t[:, :, ::1] Hx, floattype_t[:, :, ::1] Hy, floattype_t[:, :, ::1] Hz, floattype_t[:, :, :, ::1] HPhi1, floattype_t[:, :, :, ::1] HPhi2, floattype_t[:, ::1] HRA, floattype_t[:, ::1] HRB, floattype_t[:, ::1] HRE, floattype_t[:, ::1] HRF, float d): + """This function updates the Hx and Hy field components for the zplus slab. + + Args: + xs, xf, ys, yf, zs, zf (int): Cell coordinates of entire box + nthreads (int): Number of threads to use + updatecoeffs, ID, E, H (memoryviews): Access to update coeffients, ID and field component arrays + EPhi, HPhi, HRA, HRB, HRE, HRF (memoryviews): Access to PML coefficient arrays + d (float): Spatial discretisation, e.g. dx, dy or dz + """ + + cdef Py_ssize_t i, j, k, ii, jj, kk + cdef int nx, ny, nz, materialHx, materialHy + cdef float dz, dEx, dEy, RA0, RB0, RE0, RF0 + dz = d + nx = xf - xs + ny = yf - ys + nz = zf - zs + + for i in prange(0, nx, nogil=True, schedule='static', num_threads=nthreads): + ii = i + xs + for j in range(0, ny): + jj = j + ys + for k in range(0, nz): + kk = k + zs + RA0 = (HRA[0, k] - 1) + RB0 = HRB[0, k] + RE0 = HRE[0, k] + RF0 = HRF[0, k] + # Hx + materialHx = ID[3, ii, jj, kk] + dEy = (Ey[ii, jj, kk + 1] - Ey[ii, jj, kk]) / dz + Hx[ii, jj, kk] = Hx[ii, jj, kk] + updatecoeffsH[materialHx, 4] * (RA0 * dEy + RB0 * HPhi1[0, i, j, k]) + HPhi1[0, i, j, k] = RE0 * HPhi1[0, i, j, k] - RF0 * dEy + # Hy + materialHy = ID[4, ii, jj, kk] + dEx = (Ex[ii, jj, kk + 1] - Ex[ii, jj, kk]) / dz + Hy[ii, jj, kk] = Hy[ii, jj, kk] - updatecoeffsH[materialHy, 4] * (RA0 * dEx + RB0 * HPhi2[0, i, j, k]) + HPhi2[0, i, j, k] = RE0 * HPhi2[0, i, j, k] - RF0 * dEx diff --git a/gprMax/pml_2order_update.pyx b/gprMax/pml_2order_update.pyx index e7001abb..2f64048f 100644 --- a/gprMax/pml_2order_update.pyx +++ b/gprMax/pml_2order_update.pyx @@ -23,692 +23,613 @@ from cython.parallel import prange 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[:, ::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. +############################################ +# Electric field PML updates - xminus slab # +############################################ +cpdef void update_pml_2order_electric_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] Ex, floattype_t[:, :, ::1] Ey, floattype_t[:, :, ::1] Ez, floattype_t[:, :, ::1] Hx, floattype_t[:, :, ::1] Hy, floattype_t[:, :, ::1] Hz, floattype_t[:, :, :, ::1] EPhi1, floattype_t[:, :, :, ::1] EPhi2, floattype_t[:, ::1] ERA, floattype_t[:, ::1] ERB, floattype_t[:, ::1] ERE, floattype_t[:, ::1] ERF, float d): + """This function updates the Ey and Ez field components for the xminus slab. - Args: + Args: xs, xf, ys, yf, zs, zf (int): Cell coordinates of entire box nthreads (int): Number of threads to use updatecoeffs, ID, E, H (memoryviews): Access to update coeffients, ID and field component arrays - EPhi, HPhi, RA, RB, RE, RF (memoryviews): Access to PML coefficient arrays - dx, dy, dz (float): Spatial discretisation - """ + EPhi, HPhi, ERA, ERB, ERE, ERF (memoryviews): Access to PML coefficient arrays + d (float): Spatial discretisation, e.g. dx, dy or dz + """ - 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', num_threads=nthreads): - for j in range(0, ny): - for k in range(0, nz): - 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[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[:, ::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: - xs, xf, ys, yf, zs, zf (int): Cell coordinates of entire box - nthreads (int): Number of threads to use - updatecoeffs, ID, E, H (memoryviews): Access to update coeffients, ID and field component arrays - EPhi, HPhi, RA, RB, RE, RF (memoryviews): Access to PML coefficient arrays - dx, dy, dz (float): Spatial discretisation - """ - - 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', num_threads=nthreads): - for j in range(0, ny): - for k in range(0, nz): - 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[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[:, ::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: - xs, xf, ys, yf, zs, zf (int): Cell coordinates of entire box - nthreads (int): Number of threads to use - updatecoeffs, ID, E, H (memoryviews): Access to update coeffients, ID and field component arrays - EPhi, HPhi, RA, RB, RE, RF (memoryviews): Access to PML coefficient arrays - dx, dy, dz (float): Spatial discretisation - """ - - 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', num_threads=nthreads): - for j in range(0, ny): - for k in range(0, nz): - 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[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[:, ::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: - xs, xf, ys, yf, zs, zf (int): Cell coordinates of entire box - nthreads (int): Number of threads to use - updatecoeffs, ID, E, H (memoryviews): Access to update coeffients, ID and field component arrays - EPhi, HPhi, RA, RB, RE, RF (memoryviews): Access to PML coefficient arrays - dx, dy, dz (float): Spatial discretisation - """ - - 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', num_threads=nthreads): - for j in range(0, ny): - for k in range(0, nz): - 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[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 - - -############################################# -# 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[:, ::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: - xs, xf, ys, yf, zs, zf (int): Cell coordinates of entire box - nthreads (int): Number of threads to use - updatecoeffs, ID, E, H (memoryviews): Access to update coeffients, ID and field component arrays - EPhi, HPhi, RA, RB, RE, RF (memoryviews): Access to PML coefficient arrays - dx, dy, dz (float): Spatial discretisation - """ - - 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', num_threads=nthreads): - for j in range(0, ny): - for k in range(0, nz): - 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[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[:, ::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: - xs, xf, ys, yf, zs, zf (int): Cell coordinates of entire box - nthreads (int): Number of threads to use - updatecoeffs, ID, E, H (memoryviews): Access to update coeffients, ID and field component arrays - EPhi, HPhi, RA, RB, RE, RF (memoryviews): Access to PML coefficient arrays - dx, dy, dz (float): Spatial discretisation - """ - - 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', num_threads=nthreads): - for j in range(0, ny): - for k in range(0, nz): - 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[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[:, ::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: - xs, xf, ys, yf, zs, zf (int): Cell coordinates of entire box - nthreads (int): Number of threads to use - updatecoeffs, ID, E, H (memoryviews): Access to update coeffients, ID and field component arrays - EPhi, HPhi, RA, RB, RE, RF (memoryviews): Access to PML coefficient arrays - dx, dy, dz (float): Spatial discretisation - """ - - 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', num_threads=nthreads): - for j in range(0, ny): - for k in range(0, nz): - 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[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[:, ::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: - xs, xf, ys, yf, zs, zf (int): Cell coordinates of entire box - nthreads (int): Number of threads to use - updatecoeffs, ID, E, H (memoryviews): Access to update coeffients, ID and field component arrays - EPhi, HPhi, RA, RB, RE, RF (memoryviews): Access to PML coefficient arrays - dx, dy, dz (float): Spatial discretisation - """ - - 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', num_threads=nthreads): - for j in range(0, ny): - for k in range(0, nz): - 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[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 - - -############################################# -# 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[:, ::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: - xs, xf, ys, yf, zs, zf (int): Cell coordinates of entire box - nthreads (int): Number of threads to use - updatecoeffs, ID, E, H (memoryviews): Access to update coeffients, ID and field component arrays - EPhi, HPhi, RA, RB, RE, RF (memoryviews): Access to PML coefficient arrays - dx, dy, dz (float): Spatial discretisation - """ - - 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', num_threads=nthreads): - for j in range(0, ny): - for k in range(0, nz): - 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[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[:, ::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: - xs, xf, ys, yf, zs, zf (int): Cell coordinates of entire box - nthreads (int): Number of threads to use - updatecoeffs, ID, E, H (memoryviews): Access to update coeffients, ID and field component arrays - EPhi, HPhi, RA, RB, RE, RF (memoryviews): Access to PML coefficient arrays - dx, dy, dz (float): Spatial discretisation - """ - - 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', num_threads=nthreads): - for j in range(0, ny): - for k in range(0, nz): - 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[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[:, ::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: - xs, xf, ys, yf, zs, zf (int): Cell coordinates of entire box - nthreads (int): Number of threads to use - updatecoeffs, ID, E, H (memoryviews): Access to update coeffients, ID and field component arrays - EPhi, HPhi, RA, RB, RE, RF (memoryviews): Access to PML coefficient arrays - dx, dy, dz (float): Spatial discretisation - """ - - 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', num_threads=nthreads): - for j in range(0, ny): - for k in range(0, nz): - 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[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[:, ::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: - xs, xf, ys, yf, zs, zf (int): Cell coordinates of entire box - nthreads (int): Number of threads to use - updatecoeffs, ID, E, H (memoryviews): Access to update coeffients, ID and field component arrays - EPhi, HPhi, RA, RB, RE, RF (memoryviews): Access to PML coefficient arrays - dx, dy, dz (float): Spatial discretisation - """ - - 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', num_threads=nthreads): - for j in range(0, ny): - for k in range(0, nz): - 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[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 - - -############################################# -# 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[:, ::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: - xs, xf, ys, yf, zs, zf (int): Cell coordinates of entire box - nthreads (int): Number of threads to use - updatecoeffs, ID, E, H (memoryviews): Access to update coeffients, ID and field component arrays - EPhi, HPhi, RA, RB, RE, RF (memoryviews): Access to PML coefficient arrays - dx, dy, dz (float): Spatial discretisation - """ - - 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', num_threads=nthreads): - for j in range(0, ny): - for k in range(0, nz): - 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[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[:, ::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: - xs, xf, ys, yf, zs, zf (int): Cell coordinates of entire box - nthreads (int): Number of threads to use - updatecoeffs, ID, E, H (memoryviews): Access to update coeffients, ID and field component arrays - EPhi, HPhi, RA, RB, RE, RF (memoryviews): Access to PML coefficient arrays - dx, dy, dz (float): Spatial discretisation - """ - - 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', num_threads=nthreads): - for j in range(0, ny): - for k in range(0, nz): - 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[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[:, ::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: - xs, xf, ys, yf, zs, zf (int): Cell coordinates of entire box - nthreads (int): Number of threads to use - updatecoeffs, ID, E, H (memoryviews): Access to update coeffients, ID and field component arrays - EPhi, HPhi, RA, RB, RE, RF (memoryviews): Access to PML coefficient arrays - dx, dy, dz (float): Spatial discretisation - """ - - 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', num_threads=nthreads): - for j in range(0, ny): - for k in range(0, nz): - 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[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[:, ::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: - xs, xf, ys, yf, zs, zf (int): Cell coordinates of entire box - nthreads (int): Number of threads to use - updatecoeffs, ID, E, H (memoryviews): Access to update coeffients, ID and field component arrays - EPhi, HPhi, RA, RB, RE, RF (memoryviews): Access to PML coefficient arrays - dx, dy, dz (float): Spatial discretisation - """ - - 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', num_threads=nthreads): - for j in range(0, ny): - for k in range(0, nz): - 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[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 - - -############################################# -# 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[:, ::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: - xs, xf, ys, yf, zs, zf (int): Cell coordinates of entire box - nthreads (int): Number of threads to use - updatecoeffs, ID, E, H (memoryviews): Access to update coeffients, ID and field component arrays - EPhi, HPhi, RA, RB, RE, RF (memoryviews): Access to PML coefficient arrays - dx, dy, dz (float): Spatial discretisation - """ - - 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', num_threads=nthreads): - for j in range(0, ny): - for k in range(0, nz): - 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[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[:, ::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: - xs, xf, ys, yf, zs, zf (int): Cell coordinates of entire box - nthreads (int): Number of threads to use - updatecoeffs, ID, E, H (memoryviews): Access to update coeffients, ID and field component arrays - EPhi, HPhi, RA, RB, RE, RF (memoryviews): Access to PML coefficient arrays - dx, dy, dz (float): Spatial discretisation - """ - - 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', num_threads=nthreads): - for j in range(0, ny): - for k in range(0, nz): - 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[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[:, ::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: - xs, xf, ys, yf, zs, zf (int): Cell coordinates of entire box - nthreads (int): Number of threads to use - updatecoeffs, ID, E, H (memoryviews): Access to update coeffients, ID and field component arrays - EPhi, HPhi, RA, RB, RE, RF (memoryviews): Access to PML coefficient arrays - dx, dy, dz (float): Spatial discretisation - """ - - 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', num_threads=nthreads): - for j in range(0, ny): - for k in range(0, nz): - 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[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[:, ::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: - xs, xf, ys, yf, zs, zf (int): Cell coordinates of entire box - nthreads (int): Number of threads to use - updatecoeffs, ID, E, H (memoryviews): Access to update coeffients, ID and field component arrays - EPhi, HPhi, RA, RB, RE, RF (memoryviews): Access to PML coefficient arrays - dx, dy, dz (float): Spatial discretisation - """ - - cdef Py_ssize_t i, j, k - cdef int nx, ny, nz, material - cdef float dEx + cdef Py_ssize_t i, j, k, ii, jj, kk + cdef int nx, ny, nz, materialEy, materialEz + cdef float dx, dHy, dHz, RA0, RB0, RE0, RF0, RA1, RB1, RE1, RF1, RA01 + dx = d nx = xf - xs ny = yf - ys nz = zf - zs for i in prange(0, nx, nogil=True, schedule='static', num_threads=nthreads): + RA0 = ERA[0, i] + RB0 = ERB[0, i] + RE0 = ERE[0, i] + RF0 = ERF[0, i] + RA1 = ERA[1, i] + RB1 = ERB[1, i] + RE1 = ERE[1, i] + RF1 = ERF[1, i] + RA01 = ERA[0, i] * ERA[1, i] - 1 + ii = xf - i for j in range(0, ny): + jj = j + ys for k in range(0, nz): - 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[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 + kk = k + zs + # Ey + materialEy = ID[1, ii, jj, kk] + dHz = (Hz[ii, jj, kk] - Hz[ii - 1, jj, kk]) / dx + Ey[ii, jj, kk] = Ey[ii, jj, kk] - updatecoeffsE[materialEy, 4] * (RA01 * dHz + RA1 * RB0 * EPhi1[0, i, j, k] + RB1 * EPhi1[1, i, j, k]) + EPhi1[1, i, j, k] = RE1 * EPhi1[1, i, j, k] - RF1 * (RA0 * dHz + RB0 * EPhi1[0, i, j, k]) + EPhi1[0, i, j, k] = RE0 * EPhi1[0, i, j, k] - RF0 * dHz + # Ez + materialEz = ID[2, ii, jj, kk] + dHy = (Hy[ii, jj, kk] - Hy[ii - 1, jj, kk]) / dx + Ez[ii, jj, kk] = Ez[ii, jj, kk] + updatecoeffsE[materialEz, 4] * (RA01 * dHy + RA1 * RB0 * EPhi2[0, i, j, k] + RB1 * EPhi2[1, i, j, k]) + EPhi2[1, i, j, k] = RE1 * EPhi2[1, i, j, k] - RF1 * (RA0 * dHy + RB0 * EPhi2[0, i, j, k]) + EPhi2[0, i, j, k] = RE0 * EPhi2[0, i, j, k] - RF0 * dHy -############################################# -# 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[:, ::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. +########################################### +# Electric field PML updates - xplus slab # +########################################### +cpdef void update_pml_2order_electric_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] Ex, floattype_t[:, :, ::1] Ey, floattype_t[:, :, ::1] Ez, floattype_t[:, :, ::1] Hx, floattype_t[:, :, ::1] Hy, floattype_t[:, :, ::1] Hz, floattype_t[:, :, :, ::1] EPhi1, floattype_t[:, :, :, ::1] EPhi2, floattype_t[:, ::1] ERA, floattype_t[:, ::1] ERB, floattype_t[:, ::1] ERE, floattype_t[:, ::1] ERF, float d): + """This function updates the Ey and Ez field components for the xplus slab. Args: xs, xf, ys, yf, zs, zf (int): Cell coordinates of entire box nthreads (int): Number of threads to use updatecoeffs, ID, E, H (memoryviews): Access to update coeffients, ID and field component arrays - EPhi, HPhi, RA, RB, RE, RF (memoryviews): Access to PML coefficient arrays - dx, dy, dz (float): Spatial discretisation + EPhi, HPhi, ERA, ERB, ERE, ERF (memoryviews): Access to PML coefficient arrays + d (float): Spatial discretisation, e.g. dx, dy or dz """ - cdef Py_ssize_t i, j, k - cdef int nx, ny, nz, material - cdef float dEy + cdef Py_ssize_t i, j, k, ii, jj, kk + cdef int nx, ny, nz, materialEy, materialEz + cdef float dx, dHy, dHz, RA0, RB0, RE0, RF0, RA1, RB1, RE1, RF1, RA01 + dx = d nx = xf - xs ny = yf - ys nz = zf - zs for i in prange(0, nx, nogil=True, schedule='static', num_threads=nthreads): + RA0 = ERA[0, i] + RB0 = ERB[0, i] + RE0 = ERE[0, i] + RF0 = ERF[0, i] + RA1 = ERA[1, i] + RB1 = ERB[1, i] + RE1 = ERE[1, i] + RF1 = ERF[1, i] + RA01 = ERA[0, i] * ERA[1, i] - 1 + ii = i + xs for j in range(0, ny): + jj = j + ys for k in range(0, nz): - 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[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 + kk = k + zs + # Ey + materialEy = ID[1, ii, jj, kk] + dHz = (Hz[ii, jj, kk] - Hz[ii - 1, jj, kk]) / dx + Ey[ii, jj, kk] = Ey[ii, jj, kk] - updatecoeffsE[materialEy, 4] * (RA01 * dHz + RA1 * RB0 * EPhi1[0, i, j, k] + RB1 * EPhi1[1, i, j, k]) + EPhi1[1, i, j, k] = RE1 * EPhi1[1, i, j, k] - RF1 * (RA0 * dHz + RB0 * EPhi1[0, i, j, k]) + EPhi1[0, i, j, k] = RE0 * EPhi1[0, i, j, k] - RF0 * dHz + # Ez + materialEz = ID[2, ii, jj, kk] + dHy = (Hy[ii, jj, kk] - Hy[ii - 1, jj, kk]) / dx + Ez[ii, jj, kk] = Ez[ii, jj, kk] + updatecoeffsE[materialEz, 4] * (RA01 * dHy + RA1 * RB0 * EPhi2[0, i, j, k] + RB1 * EPhi2[1, i, j, k]) + EPhi2[1, i, j, k] = RE1 * EPhi2[1, i, j, k] - RF1 * (RA0 * dHy + RB0 * EPhi2[0, i, j, k]) + EPhi2[0, i, j, k] = RE0 * EPhi2[0, i, j, k] - RF0 * dHy -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. +############################################ +# Electric field PML updates - yminus slab # +############################################ +cpdef void update_pml_2order_electric_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] Ey, floattype_t[:, :, ::1] Ez, floattype_t[:, :, ::1] Hx, floattype_t[:, :, ::1] Hy, floattype_t[:, :, ::1] Hz, floattype_t[:, :, :, ::1] EPhi1, floattype_t[:, :, :, ::1] EPhi2, floattype_t[:, ::1] ERA, floattype_t[:, ::1] ERB, floattype_t[:, ::1] ERE, floattype_t[:, ::1] ERF, float d): + """This function updates the Ex and Ez field components for the yminus slab. Args: xs, xf, ys, yf, zs, zf (int): Cell coordinates of entire box nthreads (int): Number of threads to use updatecoeffs, ID, E, H (memoryviews): Access to update coeffients, ID and field component arrays - EPhi, HPhi, RA, RB, RE, RF (memoryviews): Access to PML coefficient arrays - dx, dy, dz (float): Spatial discretisation + EPhi, HPhi, ERA, ERB, ERE, ERF (memoryviews): Access to PML coefficient arrays + d (float): Spatial discretisation, e.g. dx, dy or dz """ - cdef Py_ssize_t i, j, k - cdef int nx, ny, nz, material - cdef float dEy + cdef Py_ssize_t i, j, k, ii, jj, kk + cdef int nx, ny, nz, materialEx, materialEz + cdef float dy, dHx, dHz, RA0, RB0, RE0, RF0, RA1, RB1, RE1, RF1, RA01 + dy = d nx = xf - xs ny = yf - ys nz = zf - zs for i in prange(0, nx, nogil=True, schedule='static', num_threads=nthreads): + ii = i + xs for j in range(0, ny): + jj = yf - j + RA0 = ERA[0, j] + RB0 = ERB[0, j] + RE0 = ERE[0, j] + RF0 = ERF[0, j] + RA1 = ERA[1, j] + RB1 = ERB[1, j] + RE1 = ERE[1, j] + RF1 = ERF[1, j] + RA01 = ERA[0, j] * ERA[1, j] - 1 for k in range(0, nz): - 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[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 + kk = k + zs + # Ex + materialEx = ID[0, ii, jj, kk] + dHz = (Hz[ii, jj, kk] - Hz[ii, jj - 1, kk]) / dy + Ex[ii, jj, kk] = Ex[ii, jj, kk] + updatecoeffsE[materialEx, 4] * (RA01 * dHz + RA1 * RB0 * EPhi1[0, i, j, k] + RB1 * EPhi1[1, i, j, k]) + EPhi1[1, i, j, k] = RE1 * EPhi1[1, i, j, k] - RF1 * (RA0 * dHz + RB0 * EPhi1[0, i, j, k]) + EPhi1[0, i, j, k] = RE0 * EPhi1[0, i, j, k] - RF0 * dHz + # Ez + materialEz = ID[2, ii, jj, kk] + dHx = (Hx[ii, jj, kk] - Hx[ii, jj - 1, kk]) / dy + Ez[ii, jj, kk] = Ez[ii, jj, kk] - updatecoeffsE[materialEz, 4] * (RA01 * dHx + RA1 * RB0 * EPhi2[0, i, j, k] + RB1 * EPhi2[1, i, j, k]) + EPhi2[1, i, j, k] = RE1 * EPhi2[1, i, j, k] - RF1 * (RA0 * dHx + RB0 * EPhi2[0, i, j, k]) + EPhi2[0, i, j, k] = RE0 * EPhi2[0, i, j, k] - RF0 * dHx -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. +########################################### +# Electric field PML updates - yplus slab # +########################################### +cpdef void update_pml_2order_electric_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] Ey, floattype_t[:, :, ::1] Ez, floattype_t[:, :, ::1] Hx, floattype_t[:, :, ::1] Hy, floattype_t[:, :, ::1] Hz, floattype_t[:, :, :, ::1] EPhi1, floattype_t[:, :, :, ::1] EPhi2, floattype_t[:, ::1] ERA, floattype_t[:, ::1] ERB, floattype_t[:, ::1] ERE, floattype_t[:, ::1] ERF, float d): + """This function updates the Ex and Ez field components for the yplus slab. Args: xs, xf, ys, yf, zs, zf (int): Cell coordinates of entire box nthreads (int): Number of threads to use updatecoeffs, ID, E, H (memoryviews): Access to update coeffients, ID and field component arrays - EPhi, HPhi, RA, RB, RE, RF (memoryviews): Access to PML coefficient arrays - dx, dy, dz (float): Spatial discretisation + EPhi, HPhi, ERA, ERB, ERE, ERF (memoryviews): Access to PML coefficient arrays + d (float): Spatial discretisation, e.g. dx, dy or dz """ - cdef Py_ssize_t i, j, k - cdef int nx, ny, nz, material - cdef float dEx + cdef Py_ssize_t i, j, k, ii, jj, kk + cdef int nx, ny, nz, materialEx, materialEz + cdef float dy, dHx, dHz, RA0, RB0, RE0, RF0, RA1, RB1, RE1, RF1, RA01 + dy = d nx = xf - xs ny = yf - ys nz = zf - zs for i in prange(0, nx, nogil=True, schedule='static', num_threads=nthreads): + ii = i + xs for j in range(0, ny): + jj = j + ys + RA0 = ERA[0, j] + RB0 = ERB[0, j] + RE0 = ERE[0, j] + RF0 = ERF[0, j] + RA1 = ERA[1, j] + RB1 = ERB[1, j] + RE1 = ERE[1, j] + RF1 = ERF[1, j] + RA01 = ERA[0, j] * ERA[1, j] - 1 for k in range(0, nz): - 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[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 + kk = k + zs + # Ex + materialEx = ID[0, ii, jj, kk] + dHz = (Hz[ii, jj, kk] - Hz[ii, jj - 1, kk]) / dy + Ex[ii, jj, kk] = Ex[ii, jj, kk] + updatecoeffsE[materialEx, 4] * (RA01 * dHz + RA1 * RB0 * EPhi1[0, i, j, k] + RB1 * EPhi1[1, i, j, k]) + EPhi1[1, i, j, k] = RE1 * EPhi1[1, i, j, k] - RF1 * (RA0 * dHz + RB0 * EPhi1[0, i, j, k]) + EPhi1[0, i, j, k] = RE0 * EPhi1[0, i, j, k] - RF0 * dHz + # Ez + materialEz = ID[2, ii, jj, kk] + dHx = (Hx[ii, jj, kk] - Hx[ii, jj - 1, kk]) / dy + Ez[ii, jj, kk] = Ez[ii, jj, kk] - updatecoeffsE[materialEz, 4] * (RA01 * dHx + RA1 * RB0 * EPhi2[0, i, j, k] + RB1 * EPhi2[1, i, j, k]) + EPhi2[1, i, j, k] = RE1 * EPhi2[1, i, j, k] - RF1 * (RA0 * dHx + RB0 * EPhi2[0, i, j, k]) + EPhi2[0, i, j, k] = RE0 * EPhi2[0, i, j, k] - RF0 * dHx -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. +############################################ +# Electric field PML updates - zminus slab # +############################################ +cpdef void update_pml_2order_electric_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] Ey, floattype_t[:, :, ::1] Ez, floattype_t[:, :, ::1] Hx, floattype_t[:, :, ::1] Hy, floattype_t[:, :, ::1] Hz, floattype_t[:, :, :, ::1] EPhi1, floattype_t[:, :, :, ::1] EPhi2, floattype_t[:, ::1] ERA, floattype_t[:, ::1] ERB, floattype_t[:, ::1] ERE, floattype_t[:, ::1] ERF, float d): + """This function updates the Ex and Ey field components for the zminus slab. + + Args: + xs, xf, ys, yf, zs, zf (int): Cell coordinates of entire box + nthreads (int): Number of threads to use + updatecoeffs, ID, E, H (memoryviews): Access to update coeffients, ID and field component arrays + EPhi, HPhi, ERA, ERB, ERE, ERF (memoryviews): Access to PML coefficient arrays + d (float): Spatial discretisation, e.g. dx, dy or dz + """ + + cdef Py_ssize_t i, j, k, ii, jj, kk + cdef int nx, ny, nz, materialEx, materialEy + cdef float dz, dHx, dHy, RA0, RB0, RE0, RF0, RA1, RB1, RE1, RF1, RA01 + dz = d + nx = xf - xs + ny = yf - ys + nz = zf - zs + + for i in prange(0, nx, nogil=True, schedule='static', num_threads=nthreads): + ii = i + xs + for j in range(0, ny): + jj = j + ys + for k in range(0, nz): + kk = zf - k + RA0 = ERA[0, k] + RB0 = ERB[0, k] + RE0 = ERE[0, k] + RF0 = ERF[0, k] + RA1 = ERA[1, k] + RB1 = ERB[1, k] + RE1 = ERE[1, k] + RF1 = ERF[1, k] + RA01 = ERA[0, k] * ERA[1, k] - 1 + # Ex + materialEx = ID[0, ii, jj, kk] + dHy = (Hy[ii, jj, kk] - Hy[ii, jj, kk - 1]) / dz + Ex[ii, jj, kk] = Ex[ii, jj, kk] - updatecoeffsE[materialEx, 4] * (RA01 * dHy + RA1 * RB0 * EPhi1[0, i, j, k] + RB1 * EPhi1[1, i, j, k]) + EPhi1[1, i, j, k] = RE1 * EPhi1[1, i, j, k] - RF1 * (RA0 * dHy + RB0 * EPhi1[0, i, j, k]) + EPhi1[0, i, j, k] = RE0 * EPhi1[0, i, j, k] - RF0 * dHy + # Ey + materialEy = ID[1, ii, jj, kk] + dHx = (Hx[ii, jj, kk] - Hx[ii, jj, kk - 1]) / dz + Ey[ii, jj, kk] = Ey[ii, jj, kk] + updatecoeffsE[materialEy, 4] * (RA01 * dHx + RA1 * RB0 * EPhi2[0, i, j, k] + RB1 * EPhi2[1, i, j, k]) + EPhi2[1, i, j, k] = RE1 * EPhi2[1, i, j, k] - RF1 * (RA0 * dHx + RB0 * EPhi2[0, i, j, k]) + EPhi2[0, i, j, k] = RE0 * EPhi2[0, i, j, k] - RF0 * dHx + + +########################################### +# Electric field PML updates - zplus slab # +########################################### +cpdef void update_pml_2order_electric_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] Ey, floattype_t[:, :, ::1] Ez, floattype_t[:, :, ::1] Hx, floattype_t[:, :, ::1] Hy, floattype_t[:, :, ::1] Hz, floattype_t[:, :, :, ::1] EPhi1, floattype_t[:, :, :, ::1] EPhi2, floattype_t[:, ::1] ERA, floattype_t[:, ::1] ERB, floattype_t[:, ::1] ERE, floattype_t[:, ::1] ERF, float d): + """This function updates the Ex and Ey field components for the zplus slab. Args: xs, xf, ys, yf, zs, zf (int): Cell coordinates of entire box nthreads (int): Number of threads to use updatecoeffs, ID, E, H (memoryviews): Access to update coeffients, ID and field component arrays - EPhi, HPhi, RA, RB, RE, RF (memoryviews): Access to PML coefficient arrays - dx, dy, dz (float): Spatial discretisation + EPhi, HPhi, ERA, ERB, ERE, ERF (memoryviews): Access to PML coefficient arrays + d (float): Spatial discretisation, e.g. dx, dy or dz """ - cdef Py_ssize_t i, j, k - cdef int nx, ny, nz, material - cdef float dEx + cdef Py_ssize_t i, j, k, ii, jj, kk + cdef int nx, ny, nz, materialEx, materialEy + cdef float dz, dHx, dHy, RA0, RB0, RE0, RF0, RA1, RB1, RE1, RF1, RA01 + dz = d nx = xf - xs ny = yf - ys nz = zf - zs for i in prange(0, nx, nogil=True, schedule='static', num_threads=nthreads): + ii = i + xs for j in range(0, ny): + jj = j + ys for k in range(0, nz): - 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[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 + kk = k + zs + RA0 = ERA[0, k] + RB0 = ERB[0, k] + RE0 = ERE[0, k] + RF0 = ERF[0, k] + RA1 = ERA[1, k] + RB1 = ERB[1, k] + RE1 = ERE[1, k] + RF1 = ERF[1, k] + RA01 = ERA[0, k] * ERA[1, k] - 1 + # Ex + materialEx = ID[0, ii, jj, kk] + dHy = (Hy[ii, jj, kk] - Hy[ii, jj, kk - 1]) / dz + Ex[ii, jj, kk] = Ex[ii, jj, kk] - updatecoeffsE[materialEx, 4] * (RA01 * dHy + RA1 * RB0 * EPhi1[0, i, j, k] + RB1 * EPhi1[1, i, j, k]) + EPhi1[1, i, j, k] = RE1 * EPhi1[1, i, j, k] - RF1 * (RA0 * dHy + RB0 * EPhi1[0, i, j, k]) + EPhi1[0, i, j, k] = RE0 * EPhi1[0, i, j, k] - RF0 * dHy + # Ey + materialEy = ID[1, ii, jj, kk] + dHx = (Hx[ii, jj, kk] - Hx[ii, jj, kk - 1]) / dz + Ey[ii, jj, kk] = Ey[ii, jj, kk] + updatecoeffsE[materialEy, 4] * (RA01 * dHx + RA1 * RB0 * EPhi2[0, i, j, k] + RB1 * EPhi2[1, i, j, k]) + EPhi2[1, i, j, k] = RE1 * EPhi2[1, i, j, k] - RF1 * (RA0 * dHx + RB0 * EPhi2[0, i, j, k]) + EPhi2[0, i, j, k] = RE0 * EPhi2[0, i, j, k] - RF0 * dHx + +############################################ +# Magnetic field PML updates - xminus slab # +############################################ +cpdef void update_pml_2order_magnetic_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] Ex, floattype_t[:, :, ::1] Ey, floattype_t[:, :, ::1] Ez, floattype_t[:, :, ::1] Hx, floattype_t[:, :, ::1] Hy, floattype_t[:, :, ::1] Hz, floattype_t[:, :, :, ::1] HPhi1, floattype_t[:, :, :, ::1] HPhi2, floattype_t[:, ::1] HRA, floattype_t[:, ::1] HRB, floattype_t[:, ::1] HRE, floattype_t[:, ::1] HRF, float d): + """This function updates the Hy and Hz field components for the xminus slab. + + Args: + xs, xf, ys, yf, zs, zf (int): Cell coordinates of entire box + nthreads (int): Number of threads to use + updatecoeffs, ID, E, H (memoryviews): Access to update coeffients, ID and field component arrays + EPhi, HPhi, HRA, HRB, ERE, HRF (memoryviews): Access to PML coefficient arrays + d (float): Spatial discretisation, e.g. dx, dy or dz + """ + + cdef Py_ssize_t i, j, k, ii, jj, kk + cdef int nx, ny, nz, materialHy, materialHz + cdef float dx, dEy, dEz, RA0, RB0, RE0, RF0, RA1, RB1, RE1, RF1, RA01 + dx = d + nx = xf - xs + ny = yf - ys + nz = zf - zs + + for i in prange(0, nx, nogil=True, schedule='static', num_threads=nthreads): + ii = xf - (i + 1) + RA0 = HRA[0, i] + RB0 = HRB[0, i] + RE0 = HRE[0, i] + RF0 = HRF[0, i] + RA1 = HRA[1, i] + RB1 = HRB[1, i] + RE1 = HRE[1, i] + RF1 = HRF[1, i] + RA01 = HRA[0, i] * HRA[1, i] - 1 + for j in range(0, ny): + jj = j + ys + for k in range(0, nz): + kk = k + zs + # Hy + materialHy = ID[4, ii, jj, kk] + dEz = (Ez[ii + 1, jj, kk] - Ez[ii, jj, kk]) / dx + Hy[ii, jj, kk] = Hy[ii, jj, kk] + updatecoeffsH[materialHy, 4] * (RA01 * dEz + RA1 * RB0 * HPhi1[0, i, j, k] + RB1 * HPhi1[1, i, j, k]) + HPhi1[1, i, j, k] = RE1 * HPhi1[1, i, j, k] - RF1 * (RA0 * dEz + RB0 * HPhi1[0, i, j, k]) + HPhi1[0, i, j, k] = RE0 * HPhi1[0, i, j, k] - RF0 * dEz + # Hz + materialHz = ID[5, ii, jj, kk] + dEy = (Ey[ii + 1, jj, kk] - Ey[ii, jj, kk]) / dx + Hz[ii, jj, kk] = Hz[ii, jj, kk] - updatecoeffsH[materialHz, 4] * (RA01 * dEy + RA1 * RB0 * HPhi2[0, i, j, k] + RB1 * HPhi2[1, i, j, k]) + HPhi2[1, i, j, k] = RE1 * HPhi2[1, i, j, k] - RF1 * (RA0 * dEy + RB0 * HPhi2[0, i, j, k]) + HPhi2[0, i, j, k] = RE0 * HPhi2[0, i, j, k] - RF0 * dEy + + +########################################### +# Magnetic field PML updates - xplus slab # +########################################### +cpdef void update_pml_2order_magnetic_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] Ex, floattype_t[:, :, ::1] Ey, floattype_t[:, :, ::1] Ez, floattype_t[:, :, ::1] Hx, floattype_t[:, :, ::1] Hy, floattype_t[:, :, ::1] Hz, floattype_t[:, :, :, ::1] HPhi1, floattype_t[:, :, :, ::1] HPhi2, floattype_t[:, ::1] HRA, floattype_t[:, ::1] HRB, floattype_t[:, ::1] HRE, floattype_t[:, ::1] HRF, float d): + """This function updates the Hy and Hz field components for the xplus slab. + + Args: + xs, xf, ys, yf, zs, zf (int): Cell coordinates of entire box + nthreads (int): Number of threads to use + updatecoeffs, ID, E, H (memoryviews): Access to update coeffients, ID and field component arrays + EPhi, HPhi, HRA, HRB, HRE, HRF (memoryviews): Access to PML coefficient arrays + d (float): Spatial discretisation, e.g. dx, dy or dz + """ + + cdef Py_ssize_t i, j, k, ii, jj, kk + cdef int nx, ny, nz, materialHy, materialHz + cdef float dx, dEy, dEz, RA0, RB0, RE0, RF0, RA1, RB1, RE1, RF1, RA01 + dx = d + nx = xf - xs + ny = yf - ys + nz = zf - zs + + for i in prange(0, nx, nogil=True, schedule='static', num_threads=nthreads): + ii = i + xs + RA0 = HRA[0, i] + RB0 = HRB[0, i] + RE0 = HRE[0, i] + RF0 = HRF[0, i] + RA1 = HRA[1, i] + RB1 = HRB[1, i] + RE1 = HRE[1, i] + RF1 = HRF[1, i] + RA01 = HRA[0, i] * HRA[1, i] - 1 + for j in range(0, ny): + jj = j + ys + for k in range(0, nz): + kk = k + zs + # Hy + materialHy = ID[4, ii, jj, kk] + dEz = (Ez[ii + 1, jj, kk] - Ez[ii, jj, kk]) / dx + Hy[ii, jj, kk] = Hy[ii, jj, kk] + updatecoeffsH[materialHy, 4] * (RA01 * dEz + RA1 * RB0 * HPhi1[0, i, j, k] + RB1 * HPhi1[1, i, j, k]) + HPhi1[1, i, j, k] = RE1 * HPhi1[1, i, j, k] - RF1 * (RA0 * dEz + RB0 * HPhi1[0, i, j, k]) + HPhi1[0, i, j, k] = RE0 * HPhi1[0, i, j, k] - RF0 * dEz + # Hz + materialHz = ID[5, ii, jj, kk] + dEy = (Ey[ii + 1, jj, kk] - Ey[ii, jj, kk]) / dx + Hz[ii, jj, kk] = Hz[ii, jj, kk] - updatecoeffsH[materialHz, 4] * (RA01 * dEy + RA1 * RB0 * HPhi2[0, i, j, k] + RB1 * HPhi2[1, i, j, k]) + HPhi2[1, i, j, k] = RE1 * HPhi2[1, i, j, k] - RF1 * (RA0 * dEy + RB0 * HPhi2[0, i, j, k]) + HPhi2[0, i, j, k] = RE0 * HPhi2[0, i, j, k] - RF0 * dEy + + +############################################ +# Magnetic field PML updates - yminus slab # +############################################ +cpdef void update_pml_2order_magnetic_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] Ex, floattype_t[:, :, ::1] Ey, floattype_t[:, :, ::1] Ez, floattype_t[:, :, ::1] Hx, floattype_t[:, :, ::1] Hy, floattype_t[:, :, ::1] Hz, floattype_t[:, :, :, ::1] HPhi1, floattype_t[:, :, :, ::1] HPhi2, floattype_t[:, ::1] HRA, floattype_t[:, ::1] HRB, floattype_t[:, ::1] HRE, floattype_t[:, ::1] HRF, float d): + """This function updates the Hx and Hz field components for the yminus slab. + + Args: + xs, xf, ys, yf, zs, zf (int): Cell coordinates of entire box + nthreads (int): Number of threads to use + updatecoeffs, ID, E, H (memoryviews): Access to update coeffients, ID and field component arrays + EPhi, HPhi, HRA, HRB, HRE, HRF (memoryviews): Access to PML coefficient arrays + d (float): Spatial discretisation, e.g. dx, dy or dz + """ + + cdef Py_ssize_t i, j, k, ii, jj, kk + cdef int nx, ny, nz, materialHx, materialHz + cdef float dy, dEx, dEz, RA0, RB0, RE0, RF0, RA1, RB1, RE1, RF1, RA01 + dy = d + nx = xf - xs + ny = yf - ys + nz = zf - zs + + for i in prange(0, nx, nogil=True, schedule='static', num_threads=nthreads): + ii = i + xs + for j in range(0, ny): + jj = yf - (j + 1) + RA0 = HRA[0, j] + RB0 = HRB[0, j] + RE0 = HRE[0, j] + RF0 = HRF[0, j] + RA1 = HRA[1, j] + RB1 = HRB[1, j] + RE1 = HRE[1, j] + RF1 = HRF[1, j] + RA01 = HRA[0, j] * HRA[1, j] - 1 + for k in range(0, nz): + kk = k + zs + # Hx + materialHx = ID[3, ii, jj, kk] + dEz = (Ez[ii, jj + 1, kk] - Ez[ii, jj, kk]) / dy + Hx[ii, jj, kk] = Hx[ii, jj, kk] - updatecoeffsH[materialHx, 4] * (RA01 * dEz + RA1 * RB0 * HPhi1[0, i, j, k] + RB1 * HPhi1[1, i, j, k]) + HPhi1[1, i, j, k] = RE1 * HPhi1[1, i, j, k] - RF1 * (RA0 * dEz + RB0 * HPhi1[0, i, j, k]) + HPhi1[0, i, j, k] = RE0 * HPhi1[0, i, j, k] - RF0 * dEz + # Hz + materialHz = ID[5, ii, jj, kk] + dEx = (Ex[ii, jj + 1, kk] - Ex[ii, jj, kk]) / dy + Hz[ii, jj, kk] = Hz[ii, jj, kk] + updatecoeffsH[materialHz, 4] * (RA01 * dEx + RA1 * RB0 * HPhi2[0, i, j, k] + RB1 * HPhi2[1, i, j, k]) + HPhi2[1, i, j, k] = RE1 * HPhi2[1, i, j, k] - RF1 * (RA0 * dEx + RB0 * HPhi2[0, i, j, k]) + HPhi2[0, i, j, k] = RE0 * HPhi2[0, i, j, k] - RF0 * dEx + + +########################################### +# Magnetic field PML updates - yplus slab # +########################################### +cpdef void update_pml_2order_magnetic_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] Ex, floattype_t[:, :, ::1] Ey, floattype_t[:, :, ::1] Ez, floattype_t[:, :, ::1] Hx, floattype_t[:, :, ::1] Hy, floattype_t[:, :, ::1] Hz, floattype_t[:, :, :, ::1] HPhi1, floattype_t[:, :, :, ::1] HPhi2, floattype_t[:, ::1] HRA, floattype_t[:, ::1] HRB, floattype_t[:, ::1] HRE, floattype_t[:, ::1] HRF, float d): + """This function updates the Hx and Hz field components for the yplus slab. + + Args: + xs, xf, ys, yf, zs, zf (int): Cell coordinates of entire box + nthreads (int): Number of threads to use + updatecoeffs, ID, E, H (memoryviews): Access to update coeffients, ID and field component arrays + EPhi, HPhi, HRA, HRB, HRE, HRF (memoryviews): Access to PML coefficient arrays + d (float): Spatial discretisation, e.g. dx, dy or dz + """ + + cdef Py_ssize_t i, j, k, ii, jj, kk + cdef int nx, ny, nz, materialHx, materialHz + cdef float dy, dEx, dEz, RA0, RB0, RE0, RF0, RA1, RB1, RE1, RF1, RA01 + dy = d + nx = xf - xs + ny = yf - ys + nz = zf - zs + + for i in prange(0, nx, nogil=True, schedule='static', num_threads=nthreads): + ii = i + xs + for j in range(0, ny): + jj = j + ys + RA0 = HRA[0, j] + RB0 = HRB[0, j] + RE0 = HRE[0, j] + RF0 = HRF[0, j] + RA1 = HRA[1, j] + RB1 = HRB[1, j] + RE1 = HRE[1, j] + RF1 = HRF[1, j] + RA01 = HRA[0, j] * HRA[1, j] - 1 + for k in range(0, nz): + kk = k + zs + # Hx + materialHx = ID[3, ii, jj, kk] + dEz = (Ez[ii, jj + 1, kk] - Ez[ii, jj, kk]) / dy + Hx[ii, jj, kk] = Hx[ii, jj, kk] - updatecoeffsH[materialHx, 4] * (RA01 * dEz + RA1 * RB0 * HPhi1[0, i, j, k] + RB1 * HPhi1[1, i, j, k]) + HPhi1[1, i, j, k] = RE1 * HPhi1[1, i, j, k] - RF1 * (RA0 * dEz + RB0 * HPhi1[0, i, j, k]) + HPhi1[0, i, j, k] = RE0 * HPhi1[0, i, j, k] - RF0 * dEz + # Hz + materialHz = ID[5, ii, jj, kk] + dEx = (Ex[ii, jj + 1, kk] - Ex[ii, jj, kk]) / dy + Hz[ii, jj, kk] = Hz[ii, jj, kk] + updatecoeffsH[materialHz, 4] * (RA01 * dEx + RA1 * RB0 * HPhi2[0, i, j, k] + RB1 * HPhi2[1, i, j, k]) + HPhi2[1, i, j, k] = RE1 * HPhi2[1, i, j, k] - RF1 * (RA0 * dEx + RB0 * HPhi2[0, i, j, k]) + HPhi2[0, i, j, k] = RE0 * HPhi2[0, i, j, k] - RF0 * dEx + + +############################################ +# Magnetic field PML updates - zminus slab # +############################################ +cpdef void update_pml_2order_magnetic_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] Ex, floattype_t[:, :, ::1] Ey, floattype_t[:, :, ::1] Ez, floattype_t[:, :, ::1] Hx, floattype_t[:, :, ::1] Hy, floattype_t[:, :, ::1] Hz, floattype_t[:, :, :, ::1] HPhi1, floattype_t[:, :, :, ::1] HPhi2, floattype_t[:, ::1] HRA, floattype_t[:, ::1] HRB, floattype_t[:, ::1] HRE, floattype_t[:, ::1] HRF, float d): + """This function updates the Hx and Hy field components for the zminus slab. + + Args: + xs, xf, ys, yf, zs, zf (int): Cell coordinates of entire box + nthreads (int): Number of threads to use + updatecoeffs, ID, E, H (memoryviews): Access to update coeffients, ID and field component arrays + EPhi, HPhi, HRA, HRB, HRE, HRF (memoryviews): Access to PML coefficient arrays + d (float): Spatial discretisation, e.g. dx, dy or dz + """ + + cdef Py_ssize_t i, j, k, ii, jj, kk + cdef int nx, ny, nz, materialHx, materialHy + cdef float dz, dEx, dEy, RA0, RB0, RE0, RF0, RA1, RB1, RE1, RF1, RA01 + dz = d + nx = xf - xs + ny = yf - ys + nz = zf - zs + + for i in prange(0, nx, nogil=True, schedule='static', num_threads=nthreads): + ii = i + xs + for j in range(0, ny): + jj = j + ys + for k in range(0, nz): + kk = zf - (k + 1) + RA0 = HRA[0, k] + RB0 = HRB[0, k] + RE0 = HRE[0, k] + RF0 = HRF[0, k] + RA1 = HRA[1, k] + RB1 = HRB[1, k] + RE1 = HRE[1, k] + RF1 = HRF[1, k] + RA01 = HRA[0, k] * HRA[1, k] - 1 + # Hx + materialHx = ID[3, ii, jj, kk] + dEy = (Ey[ii, jj, kk + 1] - Ey[ii, jj, kk]) / dz + Hx[ii, jj, kk] = Hx[ii, jj, kk] + updatecoeffsH[materialHx, 4] * (RA01 * dEy + RA1 * RB0 * HPhi1[0, i, j, k] + RB1 * HPhi1[1, i, j, k]) + HPhi1[1, i, j, k] = RE1 * HPhi1[1, i, j, k] - RF1 * (RA0 * dEy + RB0 * HPhi1[0, i, j, k]) + HPhi1[0, i, j, k] = RE0 * HPhi1[0, i, j, k] - RF0 * dEy + # Hy + materialHy = ID[4, ii, jj, kk] + dEx = (Ex[ii, jj, kk + 1] - Ex[ii, jj, kk]) / dz + Hy[ii, jj, kk] = Hy[ii, jj, kk] - updatecoeffsH[materialHy, 4] * (RA01 * dEx + RA1 * RB0 * HPhi2[0, i, j, k] + RB1 * HPhi2[1, i, j, k]) + HPhi2[1, i, j, k] = RE1 * HPhi2[1, i, j, k] - RF1 * (RA0 * dEx + RB0 * HPhi2[0, i, j, k]) + HPhi2[0, i, j, k] = RE0 * HPhi2[0, i, j, k] - RF0 * dEx + + +########################################### +# Magnetic field PML updates - zplus slab # +########################################### +cpdef void update_pml_2order_magnetic_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] Ex, floattype_t[:, :, ::1] Ey, floattype_t[:, :, ::1] Ez, floattype_t[:, :, ::1] Hx, floattype_t[:, :, ::1] Hy, floattype_t[:, :, ::1] Hz, floattype_t[:, :, :, ::1] HPhi1, floattype_t[:, :, :, ::1] HPhi2, floattype_t[:, ::1] HRA, floattype_t[:, ::1] HRB, floattype_t[:, ::1] HRE, floattype_t[:, ::1] HRF, float d): + """This function updates the Hx and Hy field components for the zplus slab. + + Args: + xs, xf, ys, yf, zs, zf (int): Cell coordinates of entire box + nthreads (int): Number of threads to use + updatecoeffs, ID, E, H (memoryviews): Access to update coeffients, ID and field component arrays + EPhi, HPhi, HRA, HRB, HRE, HRF (memoryviews): Access to PML coefficient arrays + d (float): Spatial discretisation, e.g. dx, dy or dz + """ + + cdef Py_ssize_t i, j, k, ii, jj, kk + cdef int nx, ny, nz, materialHx, materialHy + cdef float dz, dEx, dEy, RA0, RB0, RE0, RF0, RA1, RB1, RE1, RF1, RA01 + dz = d + nx = xf - xs + ny = yf - ys + nz = zf - zs + + for i in prange(0, nx, nogil=True, schedule='static', num_threads=nthreads): + ii = i + xs + for j in range(0, ny): + jj = j + ys + for k in range(0, nz): + kk = k + zs + RA0 = HRA[0, k] + RB0 = HRB[0, k] + RE0 = HRE[0, k] + RF0 = HRF[0, k] + RA1 = HRA[1, k] + RB1 = HRB[1, k] + RE1 = HRE[1, k] + RF1 = HRF[1, k] + RA01 = HRA[0, k] * HRA[1, k] - 1 + # Hx + materialHx = ID[3, ii, jj, kk] + dEy = (Ey[ii, jj, kk + 1] - Ey[ii, jj, kk]) / dz + Hx[ii, jj, kk] = Hx[ii, jj, kk] + updatecoeffsH[materialHx, 4] * (RA01 * dEy + RA1 * RB0 * HPhi1[0, i, j, k] + RB1 * HPhi1[1, i, j, k]) + HPhi1[1, i, j, k] = RE1 * HPhi1[1, i, j, k] - RF1 * (RA0 * dEy + RB0 * HPhi1[0, i, j, k]) + HPhi1[0, i, j, k] = RE0 * HPhi1[0, i, j, k] - RF0 * dEy + # Hy + materialHy = ID[4, ii, jj, kk] + dEx = (Ex[ii, jj, kk + 1] - Ex[ii, jj, kk]) / dz + Hy[ii, jj, kk] = Hy[ii, jj, kk] - updatecoeffsH[materialHy, 4] * (RA01 * dEx + RA1 * RB0 * HPhi2[0, i, j, k] + RB1 * HPhi2[1, i, j, k]) + HPhi2[1, i, j, k] = RE1 * HPhi2[1, i, j, k] - RF1 * (RA0 * dEx + RB0 * HPhi2[0, i, j, k]) + HPhi2[0, i, j, k] = RE0 * HPhi2[0, i, j, k] - RF0 * dEx