diff --git a/gprMax/grid.py b/gprMax/grid.py index ee93fd18..bb565662 100644 --- a/gprMax/grid.py +++ b/gprMax/grid.py @@ -127,6 +127,7 @@ class FDTDGrid(Grid): self.pmlthickness = OrderedDict((key, 10) for key in PML.boundaryIDs) self.cfs = [] self.pmls = [] + self.pmlformulation = 'HORIPML' self.materials = [] self.mixingmodels = [] diff --git a/gprMax/input_cmds_file.py b/gprMax/input_cmds_file.py index 7f36668a..e11165f3 100644 --- a/gprMax/input_cmds_file.py +++ b/gprMax/input_cmds_file.py @@ -192,7 +192,7 @@ def check_cmd_names(processedlines, checkessential=True): essentialcmds = ['#domain', '#dx_dy_dz', '#time_window'] # Commands that there should only be one instance of in a model - singlecmds = dict.fromkeys(['#domain', '#dx_dy_dz', '#time_window', '#title', '#messages', '#num_threads', '#time_step_stability_factor', '#pml_cells', '#excitation_file', '#src_steps', '#rx_steps', '#taguchi', '#end_taguchi'], None) + singlecmds = dict.fromkeys(['#domain', '#dx_dy_dz', '#time_window', '#title', '#messages', '#num_threads', '#time_step_stability_factor', '#pml_cells', '#excitation_file', '#src_steps', '#rx_steps', '#taguchi', '#end_taguchi', '#pml_formulation'], None) # Commands that there can be multiple instances of in a model - these will be lists within the dictionary multiplecmds = {key: [] for key in ['#geometry_view', '#geometry_objects_write', '#material', '#soil_peplinski', '#add_dispersion_debye', '#add_dispersion_lorentz', '#add_dispersion_drude', '#waveform', '#voltage_source', '#hertzian_dipole', '#magnetic_dipole', '#transmission_line', '#rx', '#rx_array', '#snapshot', '#pml_cfs', '#include_file']} diff --git a/gprMax/input_cmds_multiuse.py b/gprMax/input_cmds_multiuse.py index 3ba653cc..63f5d915 100644 --- a/gprMax/input_cmds_multiuse.py +++ b/gprMax/input_cmds_multiuse.py @@ -795,8 +795,8 @@ def process_multicmds(multicmds, G): if tmp[1] not in CFSParameter.scalingdirections or tmp[5] not in CFSParameter.scalingdirections or tmp[9] not in CFSParameter.scalingdirections: raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' must have scaling type {}'.format(','.join(CFSParameter.scalingdirections))) if float(tmp[2]) < 0 or float(tmp[3]) < 0 or float(tmp[6]) < 0 or float(tmp[7]) < 0 or float(tmp[10]) < 0: - raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' minimum and maximum scaling values must be greater than zero') - if float(tmp[6]) < 1: + raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' minimum and maximum scaling values must be positive') + if float(tmp[6]) < 1 and G.pmlformulation == 'HORIPML': raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' minimum scaling value for kappa must be greater than or equal to one') cfsalpha = CFSParameter() diff --git a/gprMax/input_cmds_singleuse.py b/gprMax/input_cmds_singleuse.py index 693f1dc5..7f7abbdb 100644 --- a/gprMax/input_cmds_singleuse.py +++ b/gprMax/input_cmds_singleuse.py @@ -232,6 +232,20 @@ def process_singlecmds(singlecmds, G): if 2 * G.pmlthickness['x0'] >= G.nx or 2 * G.pmlthickness['y0'] >= G.ny or 2 * G.pmlthickness['z0'] >= G.nz or 2 * G.pmlthickness['xmax'] >= G.nx or 2 * G.pmlthickness['ymax'] >= G.ny or 2 * G.pmlthickness['zmax'] >= G.nz: raise CmdInputError(cmd + ' has too many cells for the domain size') + cmd = '#pml_formulation' + if singlecmds[cmd] is not None: + pmlformulations = ['MRIPML', 'HORIPML'] + tmp = singlecmds[cmd].split() + if len(tmp) != 1: + raise CmdInputError(cmd + ' requires exactly one parameter') + else: + G.pmlformulation = singlecmds[cmd].upper() + if singlecmds[cmd].upper() in pmlformulations: + if G.messages: + print('PML formulation: {}'.format(G.pmlformulation)) + else: + raise CmdInputError(cmd + ' PML formulation is not found') + # src_steps cmd = '#src_steps' if singlecmds[cmd] is not None: diff --git a/gprMax/pml.py b/gprMax/pml.py index 1074d8c6..4fedeec0 100644 --- a/gprMax/pml.py +++ b/gprMax/pml.py @@ -24,13 +24,14 @@ from tqdm import tqdm from gprMax.constants import e0 from gprMax.constants import z0 from gprMax.constants import floattype +from gprMax.exceptions import GeneralError class CFSParameter(object): """Individual CFS parameter (e.g. alpha, kappa, or sigma).""" # Allowable scaling profiles and directions - scalingprofiles = {'constant': 0, 'linear': 1, 'quadratic': 2, 'cubic': 3, 'quartic': 4, 'quintic': 5, 'sextic': 6} + scalingprofiles = {'constant': 0, 'linear': 1, 'quadratic': 2, 'cubic': 3, 'quartic': 4, 'quintic': 5, 'sextic': 6, 'dodeka': 12} scalingdirections = ['forward', 'reverse'] def __init__(self, ID=None, scaling='polynomial', scalingprofile=None, scalingdirection='forward', min=0, max=0): @@ -232,19 +233,41 @@ class PML(object): Ekappa, Hkappa = cfs.calculate_values(self.thickness, cfs.kappa) Esigma, Hsigma = cfs.calculate_values(self.thickness, cfs.sigma) - # Electric PML update coefficients - tmp = (2 * e0 * Ekappa) + G.dt * (Ealpha * Ekappa + Esigma) - self.ERA[x, :] = (2 * e0 + G.dt * Ealpha) / tmp - self.ERB[x, :] = (2 * e0 * Ekappa) / tmp - self.ERE[x, :] = ((2 * e0 * Ekappa) - G.dt * (Ealpha * Ekappa + Esigma)) / tmp - self.ERF[x, :] = (2 * Esigma * G.dt) / (Ekappa * tmp) + # Define different parameters depending on PML formulation + if G.pmlformulation == 'HORIPML': + # HORIPML Electric PML update coefficients + tmp = (2 * e0 * Ekappa) + G.dt * (Ealpha * Ekappa + Esigma) + self.ERA[x, :] = (2 * e0 + G.dt * Ealpha) / tmp + self.ERB[x, :] = (2 * e0 * Ekappa) / tmp + self.ERE[x, :] = ((2 * e0 * Ekappa) - G.dt * (Ealpha * Ekappa + Esigma)) / tmp + self.ERF[x, :] = (2 * Esigma * G.dt) / (Ekappa * tmp) + + # HORIPML Magnetic PML update coefficients + tmp = (2 * e0 * Hkappa) + G.dt * (Halpha * Hkappa + Hsigma) + self.HRA[x, :] = (2 * e0 + G.dt * Halpha) / tmp + self.HRB[x, :] = (2 * e0 * Hkappa) / tmp + self.HRE[x, :] = ((2 * e0 * Hkappa) - G.dt * (Halpha * Hkappa + Hsigma)) / tmp + self.HRF[x, :] = (2 * Hsigma * G.dt) / (Hkappa * tmp) + + elif G.pmlformulation == 'MRIPML': + # MRIPML Electric PML update coefficients + tmp = 2 * e0 + G.dt * Ealpha + self.ERA[x, :] = Ekappa + (G.dt * Esigma) / tmp + self.ERB[x, :] = (2 * e0) / tmp + self.ERE[x, :] = ((2 * e0) - G.dt * Ealpha) / tmp + self.ERF[x, :] = (2 * Esigma * G.dt) / tmp + + # MRIPML Magnetic PML update coefficients + tmp = 2 * e0 + G.dt * Halpha + self.HRA[x, :] = Hkappa + (G.dt * Hsigma) / tmp + self.HRB[x, :] = (2 * e0) / tmp + self.HRE[x, :] = ((2 * e0) - G.dt * Halpha) / tmp + self.HRF[x, :] = (2 * Hsigma * G.dt) / tmp + + else: + GeneralError('Do not know about ' + G.pmlformulation + ' No PML update coefficients calculated') + - # Magnetic PML update coefficients - tmp = (2 * e0 * Hkappa) + G.dt * (Halpha * Hkappa + Hsigma) - self.HRA[x, :] = (2 * e0 + G.dt * Halpha) / tmp - self.HRB[x, :] = (2 * e0 * Hkappa) / tmp - self.HRE[x, :] = ((2 * e0 * Hkappa) - G.dt * (Halpha * Hkappa + Hsigma)) / tmp - self.HRF[x, :] = (2 * Hsigma * G.dt) / (Hkappa * tmp) def update_electric(self, G): """This functions updates electric field components with the PML correction. @@ -253,7 +276,7 @@ class PML(object): G (class): Grid class instance - holds essential parameters describing the model. """ - func = getattr(import_module('gprMax.pml_updates_ext'), 'update_pml_' + str(len(self.CFS)) + 'order_electric_' + self.direction) + func = getattr(import_module('gprMax.pml_updates_ext'), 'update_pml_' + str(len(self.CFS)) + 'order_electric_' + G.pmlformulation + '_' + self.direction) func(self.xs, self.xf, self.ys, self.yf, self.zs, self.zf, G.nthreads, G.updatecoeffsE, G.ID, G.Ex, G.Ey, G.Ez, G.Hx, G.Hy, G.Hz, self.EPhi1, self.EPhi2, self.ERA, self.ERB, self.ERE, self.ERF, self.d) def update_magnetic(self, G): @@ -263,7 +286,7 @@ class PML(object): G (class): Grid class instance - holds essential parameters describing the model. """ - func = getattr(import_module('gprMax.pml_updates_ext'), 'update_pml_' + str(len(self.CFS)) + 'order_magnetic_' + self.direction) + func = getattr(import_module('gprMax.pml_updates_ext'), 'update_pml_' + str(len(self.CFS)) + 'order_magnetic_' + G.pmlformulation + '_' + self.direction) func(self.xs, self.xf, self.ys, self.yf, self.zs, self.zf, G.nthreads, G.updatecoeffsH, G.ID, G.Ex, G.Ey, G.Ez, G.Hx, G.Hy, G.Hz, self.HPhi1, self.HPhi2, self.HRA, self.HRB, self.HRE, self.HRF, self.d) def gpu_set_blocks_per_grid(self, G): diff --git a/gprMax/pml_updates_ext.pyx b/gprMax/pml_updates_ext.pyx index d6b3ffd9..53897065 100644 --- a/gprMax/pml_updates_ext.pyx +++ b/gprMax/pml_updates_ext.pyx @@ -26,7 +26,7 @@ from gprMax.constants cimport floattype_t, complextype_t ######################################################## # Electric field PML updates - 1st order - xminus slab # ######################################################## -cpdef void update_pml_1order_electric_xminus( +cpdef void update_pml_1order_electric_HORIPML_xminus( int xs, int xf, int ys, @@ -90,10 +90,76 @@ cpdef void update_pml_1order_electric_xminus( EPhi2[0, i, j, k] = RE0 * EPhi2[0, i, j, k] - RF0 * dHy +cpdef void update_pml_1order_electric_MRIPML_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: + 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, materialEy, materialEz + cdef float dx, dHy, dHz, RA0, RA1, RB0, RC0, 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): + RA1 = ((1 / ERA[0, i]) - 1) + RA0 = (1 / ERA[0, i]) + RB0 = ERB[0, i] + RE0 = ERE[0, i] + RF0 = ERF[0, i] + RC0 = RF0*RB0*RA0 + ii = xf - i + for j in range(0, ny): + jj = j + ys + for k in range(0, nz): + 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] * (RA1 * dHz - RA0 * EPhi1[0, i, j, k]) + EPhi1[0, i, j, k] = RE0 * EPhi1[0, i, j, k] + RC0 * dHz - RC0 * EPhi1[0, i, j, k] + # 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] * (RA1 * dHy - RA0 * EPhi2[0, i, j, k]) + EPhi2[0, i, j, k] = RE0 * EPhi2[0, i, j, k] + RC0 * dHy - RC0 * EPhi2[0, i, j, k] + + ####################################################### # Electric field PML updates - 1st order - xplus slab # ####################################################### -cpdef void update_pml_1order_electric_xplus( +cpdef void update_pml_1order_electric_HORIPML_xplus( int xs, int xf, int ys, @@ -156,11 +222,76 @@ cpdef void update_pml_1order_electric_xplus( 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_electric_MRIPML_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, 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, materialEy, materialEz + cdef float dx, dHy, dHz, RA0, RA1, RB0, RC0, 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): + RA1 = (1 / ERA[0, i]) - 1 + RA0 = 1 / ERA[0, i] + RB0 = ERB[0, i] + RE0 = ERE[0, i] + RF0 = ERF[0, i] + RC0 = RF0*RB0*RA0 + ii = i + xs + for j in range(0, ny): + jj = j + ys + for k in range(0, nz): + 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] * (RA1 * dHz - RA0 * EPhi1[0, i, j, k]) + EPhi1[0, i, j, k] = RE0 * EPhi1[0, i, j, k] + RC0 * dHz - RC0 * EPhi1[0, i, j, k] + # 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] * (RA1 * dHy - RA0 * EPhi2[0, i, j, k]) + EPhi2[0, i, j, k] = RE0 * EPhi2[0, i, j, k] + RC0 * dHy - RC0 * EPhi2[0, i, j, k] + ######################################################## # Electric field PML updates - 1st order - yminus slab # ######################################################## -cpdef void update_pml_1order_electric_yminus( +cpdef void update_pml_1order_electric_HORIPML_yminus( int xs, int xf, int ys, @@ -223,11 +354,76 @@ cpdef void update_pml_1order_electric_yminus( 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_electric_MRIPML_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, 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, materialEz + cdef float dy, dHx, dHz, RA0, RA1, RB0, RC0, 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 + RA1 = (1 / ERA[0, j]) - 1 + RA0 = 1 / ERA[0, j] + RB0 = ERB[0, j] + RE0 = ERE[0, j] + RF0 = ERF[0, j] + RC0 = RF0*RB0*RA0 + for k in range(0, nz): + 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] * (RA1 * dHz - RA0 * EPhi1[0, i, j, k]) + EPhi1[0, i, j, k] = RE0 * EPhi1[0, i, j, k] + RC0 * dHz - RC0 * EPhi1[0, i, j, k] + # 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] * (RA1 * dHx - RA0 * EPhi2[0, i, j, k]) + EPhi2[0, i, j, k] = RE0 * EPhi2[0, i, j, k] + RC0 * dHx - RC0 * EPhi2[0, i, j, k] + ####################################################### # Electric field PML updates - 1st order - yplus slab # ####################################################### -cpdef void update_pml_1order_electric_yplus( +cpdef void update_pml_1order_electric_HORIPML_yplus( int xs, int xf, int ys, @@ -290,11 +486,77 @@ cpdef void update_pml_1order_electric_yplus( 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_electric_MRIPML_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, 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, materialEz + cdef float dy, dHx, dHz, RA0, RA1, RB0, RC0, 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 + RA1 = (1 / ERA[0, j]) - 1 + RA0 = 1 / ERA[0, j] + RB0 = ERB[0, j] + RE0 = ERE[0, j] + RF0 = ERF[0, j] + RC0 = RF0*RB0*RA0 + for k in range(0, nz): + 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] * (RA1 * dHz - RA0 * EPhi1[0, i, j, k]) + EPhi1[0, i, j, k] = RE0 * EPhi1[0, i, j, k] + RC0 * dHz - RC0 * EPhi1[0, i, j, k] + # 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] * (RA1 * dHx - RA0 * EPhi2[0, i, j, k]) + EPhi2[0, i, j, k] = RE0 * EPhi2[0, i, j, k] + RC0 * dHx - RC0 * EPhi2[0, i, j, k] + + ######################################################## # Electric field PML updates - 1st order - zminus slab # ######################################################## -cpdef void update_pml_1order_electric_zminus( +cpdef void update_pml_1order_electric_HORIPML_zminus( int xs, int xf, int ys, @@ -357,11 +619,76 @@ cpdef void update_pml_1order_electric_zminus( 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 +cpdef void update_pml_1order_electric_MRIPML_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, RA1, RB0, RC0, 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 + RA1 = (1 / ERA[0, k]) - 1 + RA0 = 1 / ERA[0, k] + RB0 = ERB[0, k] + RE0 = ERE[0, k] + RF0 = ERF[0, k] + RC0 = RF0*RB0*RA0 + # 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] * (RA1 * dHy - RA0 * EPhi1[0, i, j, k]) + EPhi1[0, i, j, k] = RE0 * EPhi1[0, i, j, k] + RC0 * dHy - RC0 * EPhi1[0, i, j, k] + # 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] * (RA1 * dHx - RA0 * EPhi2[0, i, j, k]) + EPhi2[0, i, j, k] = RE0 * EPhi2[0, i, j, k] + RC0 * dHx - RC0 * EPhi2[0, i, j, k] + ####################################################### # Electric field PML updates - 1st order - zplus slab # ####################################################### -cpdef void update_pml_1order_electric_zplus( +cpdef void update_pml_1order_electric_HORIPML_zplus( int xs, int xf, int ys, @@ -424,11 +751,76 @@ cpdef void update_pml_1order_electric_zplus( 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 +cpdef void update_pml_1order_electric_MRIPML_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, 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, RA1, RB0, RC0, 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 + RA1 = (1 / ERA[0, k]) - 1 + RA0 = 1 / ERA[0, k] + RB0 = ERB[0, k] + RE0 = ERE[0, k] + RF0 = ERF[0, k] + RC0 = RF0*RB0*RA0 + # 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] * (RA1 * dHy - RA0 * EPhi1[0, i, j, k]) + EPhi1[0, i, j, k] = RE0 * EPhi1[0, i, j, k] + RC0 * dHy - RC0 * EPhi1[0, i, j, k] + # 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] * (RA1 * dHx - RA0 * EPhi2[0, i, j, k]) + EPhi2[0, i, j, k] = RE0 * EPhi2[0, i, j, k] + RC0 * dHx - RC0 * EPhi2[0, i, j, k] + ######################################################## # Magnetic field PML updates - 1st order - xminus slab # ######################################################## -cpdef void update_pml_1order_magnetic_xminus( +cpdef void update_pml_1order_magnetic_HORIPML_xminus( int xs, int xf, int ys, @@ -491,11 +883,76 @@ cpdef void update_pml_1order_magnetic_xminus( 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 +cpdef void update_pml_1order_magnetic_MRIPML_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, RA1, RB0, RC0, 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) + RA1 = (1 / HRA[0, i]) - 1 + RA0 = 1 / HRA[0, i] + RB0 = HRB[0, i] + RE0 = HRE[0, i] + RF0 = HRF[0, i] + RC0 = RF0*RB0*RA0 + 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] * (RA1 * dEz - RA0 * HPhi1[0, i, j, k]) + HPhi1[0, i, j, k] = RE0 * HPhi1[0, i, j, k] + RC0 * dEz - RC0 * HPhi1[0, i, j, k] + # 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] * (RA1 * dEy - RA0 * HPhi2[0, i, j, k]) + HPhi2[0, i, j, k] = RE0 * HPhi2[0, i, j, k] + RC0 * dEy - RC0 * HPhi2[0, i, j, k] + ####################################################### # Magnetic field PML updates - 1st order - xplus slab # ####################################################### -cpdef void update_pml_1order_magnetic_xplus( +cpdef void update_pml_1order_magnetic_HORIPML_xplus( int xs, int xf, int ys, @@ -558,11 +1015,76 @@ cpdef void update_pml_1order_magnetic_xplus( 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 +cpdef void update_pml_1order_magnetic_MRIPML_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, RA1, RB0, RC0, 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 + RA1 = (1 / HRA[0, i]) - 1 + RA0 = 1 / HRA[0, i] + RB0 = HRB[0, i] + RE0 = HRE[0, i] + RF0 = HRF[0, i] + RC0 = RF0*RB0*RA0 + 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] * (RA1 * dEz - RA0 * HPhi1[0, i, j, k]) + HPhi1[0, i, j, k] = RE0 * HPhi1[0, i, j, k] + RC0 * dEz - RC0 * HPhi1[0, i, j, k] + # 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] * (RA1 * dEy - RA0 * HPhi2[0, i, j, k]) + HPhi2[0, i, j, k] = RE0 * HPhi2[0, i, j, k] + RC0 * dEy - RC0 * HPhi2[0, i, j, k] + ######################################################## # Magnetic field PML updates - 1st order - yminus slab # ######################################################## -cpdef void update_pml_1order_magnetic_yminus( +cpdef void update_pml_1order_magnetic_HORIPML_yminus( int xs, int xf, int ys, @@ -625,11 +1147,76 @@ cpdef void update_pml_1order_magnetic_yminus( 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 +cpdef void update_pml_1order_magnetic_MRIPML_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, RA1, RB0, RC0, 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) + RA1 = (1 / HRA[0, j]) - 1 + RA0 = 1 / HRA[0, j] + RB0 = HRB[0, j] + RE0 = HRE[0, j] + RF0 = HRF[0, j] + RC0 = RF0*RB0*RA0 + 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] * (RA1 * dEz - RA0 * HPhi1[0, i, j, k]) + HPhi1[0, i, j, k] = RE0 * HPhi1[0, i, j, k] + RC0 * dEz - RC0 * HPhi1[0, i, j, k] + # 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] * (RA1 * dEx - RA0 * HPhi2[0, i, j, k]) + HPhi2[0, i, j, k] = RE0 * HPhi2[0, i, j, k] + RC0 * dEx - RC0 * HPhi2[0, i, j, k] + ####################################################### # Magnetic field PML updates - 1st order - yplus slab # ####################################################### -cpdef void update_pml_1order_magnetic_yplus( +cpdef void update_pml_1order_magnetic_HORIPML_yplus( int xs, int xf, int ys, @@ -692,11 +1279,76 @@ cpdef void update_pml_1order_magnetic_yplus( 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 +cpdef void update_pml_1order_magnetic_MRIPML_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, RA1, RB0, RC0, 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 + RA1 = (1 / HRA[0, j]) - 1 + RA0 = 1 / HRA[0, j] + RB0 = HRB[0, j] + RE0 = HRE[0, j] + RF0 = HRF[0, j] + RC0 = RF0*RB0*RA0 + 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] * (RA1 * dEz - RA0 * HPhi1[0, i, j, k]) + HPhi1[0, i, j, k] = RE0 * HPhi1[0, i, j, k] + RC0 * dEz - RC0 * HPhi1[0, i, j, k] + # 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] * (RA1 * dEx - RA0 * HPhi2[0, i, j, k]) + HPhi2[0, i, j, k] = RE0 * HPhi2[0, i, j, k] + RC0 * dEx - RC0 * HPhi2[0, i, j, k] + ######################################################## # Magnetic field PML updates - 1st order - zminus slab # ######################################################## -cpdef void update_pml_1order_magnetic_zminus( +cpdef void update_pml_1order_magnetic_HORIPML_zminus( int xs, int xf, int ys, @@ -759,11 +1411,77 @@ cpdef void update_pml_1order_magnetic_zminus( 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 +cpdef void update_pml_1order_magnetic_MRIPML_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, RA1, RB0, RC0, 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) + RA1 = (1 / HRA[0, k]) - 1 + RA0 = 1 / HRA[0, k] + RB0 = HRB[0, k] + RE0 = HRE[0, k] + RF0 = HRF[0, k] + RC0 = RF0*RB0*RA0 + # 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] * (RA1 * dEy - RA0 * HPhi1[0, i, j, k]) + HPhi1[0, i, j, k] = RE0 * HPhi1[0, i, j, k] + RC0 * dEy - RC0 * HPhi1[0, i, j, k] + # 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] * (RA1 * dEx - RA0 * HPhi2[0, i, j, k]) + HPhi2[0, i, j, k] = RE0 * HPhi2[0, i, j, k] + RC0 * dEx - RC0 * HPhi2[0, i, j, k] + + ####################################################### # Magnetic field PML updates - 1st order - zplus slab # ####################################################### -cpdef void update_pml_1order_magnetic_zplus( +cpdef void update_pml_1order_magnetic_HORIPML_zplus( int xs, int xf, int ys, @@ -826,11 +1544,76 @@ cpdef void update_pml_1order_magnetic_zplus( 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 +cpdef void update_pml_1order_magnetic_MRIPML_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, RA1, RB0, RC0, 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 + RA1 = (1 / HRA[0, k]) - 1 + RA0 = 1 / HRA[0, k] + RB0 = HRB[0, k] + RE0 = HRE[0, k] + RF0 = HRF[0, k] + RC0 = RF0*RB0*RA0 + # 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] * (RA1 * dEy - RA0 * HPhi1[0, i, j, k]) + HPhi1[0, i, j, k] = RE0 * HPhi1[0, i, j, k] + RC0 * dEy - RC0 * HPhi1[0, i, j, k] + # 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] * (RA1 * dEx - RA0 * HPhi2[0, i, j, k]) + HPhi2[0, i, j, k] = RE0 * HPhi2[0, i, j, k] + RC0 * dEx - RC0 * HPhi2[0, i, j, k] + ######################################################## # Electric field PML updates - 2nd order - xminus slab # ######################################################## -cpdef void update_pml_2order_electric_xminus( +cpdef void update_pml_2order_electric_HORIPML_xminus( int xs, int xf, int ys, @@ -900,11 +1683,86 @@ cpdef void update_pml_2order_electric_xminus( 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 void update_pml_2order_electric_MRIPML_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: + 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, materialEy, materialEz + cdef float dx, dHy, dHz, RA, IRA1, IRA, RB0, RC0, RE0, RF0, RB1, RE1, RC1, RF1, EPsi1, EPsi2 + dx = d + nx = xf - xs + ny = yf - ys + nz = zf - zs + + for i in prange(0, nx, nogil=True, schedule='static', num_threads=nthreads): + RA = ERA[0, i] + ERA[1, i] + IRA1 = (1 / RA) - 1 + IRA = 1 / RA + RB0 = ERB[0, i] + RE0 = ERE[0, i] + RF0 = ERF[0, i] + RC0 = RF0*RB0*IRA + RB1 = ERB[1, i] + RE1 = ERE[1, i] + RF1 = ERF[1, i] + RC1 = RF1*RB1*IRA + ii = xf - i + for j in range(0, ny): + jj = j + ys + for k in range(0, nz): + kk = k + zs + # Ey + materialEy = ID[1, ii, jj, kk] + dHz = (Hz[ii, jj, kk] - Hz[ii - 1, jj, kk]) / dx + EPsi1 = EPhi1[0, i, j, k] + EPhi1[1, i, j, k] + Ey[ii, jj, kk] = Ey[ii, jj, kk] - updatecoeffsE[materialEy, 4] * (IRA1 * dHz - IRA * EPsi1) + EPhi1[1, i, j, k] = RE1 * EPhi1[1, i, j, k] + RC1 * dHz - RC1 * EPsi1 + EPhi1[0, i, j, k] = RE0 * EPhi1[0, i, j, k] + RC0 * dHz - RC0 * EPsi1 + # Ez + materialEz = ID[2, ii, jj, kk] + dHy = (Hy[ii, jj, kk] - Hy[ii - 1, jj, kk]) / dx + EPsi2 = EPhi2[0, i, j, k] + EPhi2[1, i, j, k] + Ez[ii, jj, kk] = Ez[ii, jj, kk] + updatecoeffsE[materialEz, 4] * (IRA1 * dHy - IRA * EPsi2) + EPhi2[1, i, j, k] = RE1 * EPhi2[1, i, j, k] + RC1 * dHy - RC1 * EPsi2 + EPhi2[0, i, j, k] = RE0 * EPhi2[0, i, j, k] + RC0 * dHy - RC0 * EPsi2 + + ####################################################### # Electric field PML updates - 2nd order - xplus slab # ####################################################### -cpdef void update_pml_2order_electric_xplus( +cpdef void update_pml_2order_electric_HORIPML_xplus( int xs, int xf, int ys, @@ -974,11 +1832,85 @@ cpdef void update_pml_2order_electric_xplus( 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 void update_pml_2order_electric_MRIPML_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, 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, materialEy, materialEz + cdef float dx, dHy, dHz, RA, IRA1, IRA, RB0, RC0, RE0, RF0, RB1, RC1, RE1, RF1, EPsi1, EPsi2 + dx = d + nx = xf - xs + ny = yf - ys + nz = zf - zs + + for i in prange(0, nx, nogil=True, schedule='static', num_threads=nthreads): + RA = ERA[0, i] + ERA[1, i] + IRA1 = (1 / RA) - 1 + IRA = 1 / RA + RB0 = ERB[0, i] + RE0 = ERE[0, i] + RF0 = ERF[0, i] + RC0 = RF0*RB0*IRA + RB1 = ERB[1, i] + RE1 = ERE[1, i] + RF1 = ERF[1, i] + RC1 = RF1*RB1*IRA + ii = i + xs + for j in range(0, ny): + jj = j + ys + for k in range(0, nz): + kk = k + zs + # Ey + materialEy = ID[1, ii, jj, kk] + dHz = (Hz[ii, jj, kk] - Hz[ii - 1, jj, kk]) / dx + EPsi1 = EPhi1[0, i, j, k] + EPhi1[1, i, j, k] + Ey[ii, jj, kk] = Ey[ii, jj, kk] - updatecoeffsE[materialEy, 4] * (IRA1 * dHz - IRA * EPsi1) + EPhi1[1, i, j, k] = RE1 * EPhi1[1, i, j, k] + RC1 * dHz - RC1 * EPsi1 + EPhi1[0, i, j, k] = RE0 * EPhi1[0, i, j, k] + RC0 * dHz - RC0 * EPsi1 + # Ez + materialEz = ID[2, ii, jj, kk] + dHy = (Hy[ii, jj, kk] - Hy[ii - 1, jj, kk]) / dx + EPsi2 = EPhi2[0, i, j, k] + EPhi2[1, i, j, k] + Ez[ii, jj, kk] = Ez[ii, jj, kk] + updatecoeffsE[materialEz, 4] * (IRA1 * dHy - IRA * EPsi2) + EPhi2[1, i, j, k] = RE1 * EPhi2[1, i, j, k] + RC1 * dHy - RC1 * EPsi2 + EPhi2[0, i, j, k] = RE0 * EPhi2[0, i, j, k] + RC0 * dHy - RC0 * EPsi2 + ######################################################## # Electric field PML updates - 2nd order - yminus slab # ######################################################## -cpdef void update_pml_2order_electric_yminus( +cpdef void update_pml_2order_electric_HORIPML_yminus( int xs, int xf, int ys, @@ -1048,11 +1980,85 @@ cpdef void update_pml_2order_electric_yminus( 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 void update_pml_2order_electric_MRIPML_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, 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, materialEz + cdef float dy, dHx, dHz, RA, IRA1, IRA, RB0, RC0, RE0, RF0, RB1, RC1, RE1, RF1, EPsi1, EPsi2 + 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 + RA = ERA[0, j] + ERA[1, j] + IRA1 = (1 / RA) - 1 + IRA = 1 / RA + RB0 = ERB[0, j] + RE0 = ERE[0, j] + RF0 = ERF[0, j] + RC0 = RF0*RB0*IRA + RB1 = ERB[1, j] + RE1 = ERE[1, j] + RF1 = ERF[1, j] + RC1 = RF1*RB1*IRA + for k in range(0, nz): + kk = k + zs + # Ex + materialEx = ID[0, ii, jj, kk] + dHz = (Hz[ii, jj, kk] - Hz[ii, jj - 1, kk]) / dy + EPsi1 = EPhi1[0, i, j, k] + EPhi1[1, i, j, k] + Ex[ii, jj, kk] = Ex[ii, jj, kk] + updatecoeffsE[materialEx, 4] * (IRA1 * dHz - IRA * EPsi1) + EPhi1[1, i, j, k] = RE1 * EPhi1[1, i, j, k] + RC1 * dHz - RC1 * EPsi1 + EPhi1[0, i, j, k] = RE0 * EPhi1[0, i, j, k] + RC0 * dHz - RC0 * EPsi1 + # Ez + materialEz = ID[2, ii, jj, kk] + dHx = (Hx[ii, jj, kk] - Hx[ii, jj - 1, kk]) / dy + EPsi2 = EPhi2[0, i, j, k] + EPhi2[1, i, j, k] + Ez[ii, jj, kk] = Ez[ii, jj, kk] - updatecoeffsE[materialEz, 4] * (IRA1 * dHx - IRA * EPsi2) + EPhi2[1, i, j, k] = RE1 * EPhi2[1, i, j, k] + RC1 * dHx - RC1 * EPsi2 + EPhi2[0, i, j, k] = RE0 * EPhi2[0, i, j, k] + RC0 * dHx - RC0 * EPsi2 + ####################################################### # Electric field PML updates - 2nd order - yplus slab # ####################################################### -cpdef void update_pml_2order_electric_yplus( +cpdef void update_pml_2order_electric_HORIPML_yplus( int xs, int xf, int ys, @@ -1122,11 +2128,85 @@ cpdef void update_pml_2order_electric_yplus( 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 void update_pml_2order_electric_MRIPML_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, 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, materialEz + cdef float dy, dHx, dHz, RA, IRA1, IRA, RB0, RC0, RE0, RF0, RB1, RC1, RE1, RF1, EPsi1, EPsi2 + 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 + RA = ERA[0, j] + ERA[1, j] + IRA1 = (1 / RA) - 1 + IRA = 1 / RA + RB0 = ERB[0, j] + RE0 = ERE[0, j] + RF0 = ERF[0, j] + RC0 = RF0*RB0*IRA + RB1 = ERB[1, j] + RE1 = ERE[1, j] + RF1 = ERF[1, j] + RC1 = RF1*RB1*IRA + for k in range(0, nz): + kk = k + zs + # Ex + materialEx = ID[0, ii, jj, kk] + dHz = (Hz[ii, jj, kk] - Hz[ii, jj - 1, kk]) / dy + EPsi1 = EPhi1[0, i, j, k] + EPhi1[1, i, j, k] + Ex[ii, jj, kk] = Ex[ii, jj, kk] + updatecoeffsE[materialEx, 4] * (IRA1 * dHz - IRA * EPsi1) + EPhi1[1, i, j, k] = RE1 * EPhi1[1, i, j, k] + RC1 * dHz - RC1 * EPsi1 + EPhi1[0, i, j, k] = RE0 * EPhi1[0, i, j, k] + RC0 * dHz - RC0 * EPsi1 + # Ez + materialEz = ID[2, ii, jj, kk] + dHx = (Hx[ii, jj, kk] - Hx[ii, jj - 1, kk]) / dy + EPsi2 = EPhi2[0, i, j, k] + EPhi2[1, i, j, k] + Ez[ii, jj, kk] = Ez[ii, jj, kk] - updatecoeffsE[materialEz, 4] * (IRA1 * dHx - IRA * EPsi2) + EPhi2[1, i, j, k] = RE1 * EPhi2[1, i, j, k] + RC1 * dHx - RC1 * EPsi2 + EPhi2[0, i, j, k] = RE0 * EPhi2[0, i, j, k] + RC0 * dHx - RC0 * EPsi2 + ######################################################## # Electric field PML updates - 2nd order - zminus slab # ######################################################## -cpdef void update_pml_2order_electric_zminus( +cpdef void update_pml_2order_electric_HORIPML_zminus( int xs, int xf, int ys, @@ -1196,11 +2276,85 @@ cpdef void update_pml_2order_electric_zminus( 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 void update_pml_2order_electric_MRIPML_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, RA, IRA1, IRA, RB0, RC0, RE0, RF0, RB1, RC1, RE1, RF1, EPsi1, EPsi2 + 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 + RA = ERA[0, k] + ERA[1, k] + IRA1 = (1 /RA) - 1 + IRA = 1 / RA + RB0 = ERB[0, k] + RE0 = ERE[0, k] + RF0 = ERF[0, k] + RC0 = RF0*RB0*IRA + RB1 = ERB[1, k] + RE1 = ERE[1, k] + RF1 = ERF[1, k] + RC1 = RF1*RB1*IRA + # Ex + materialEx = ID[0, ii, jj, kk] + dHy = (Hy[ii, jj, kk] - Hy[ii, jj, kk - 1]) / dz + EPsi1 = EPhi1[0, i, j, k] + EPhi1[1, i, j, k] + Ex[ii, jj, kk] = Ex[ii, jj, kk] - updatecoeffsE[materialEx, 4] * (IRA1 * dHy - IRA * EPsi1) + EPhi1[1, i, j, k] = RE1 * EPhi1[1, i, j, k] + RC1 * dHy - RC1 * EPsi1 + EPhi1[0, i, j, k] = RE0 * EPhi1[0, i, j, k] + RC0 * dHy - RC0 * EPsi1 + # Ey + materialEy = ID[1, ii, jj, kk] + dHx = (Hx[ii, jj, kk] - Hx[ii, jj, kk - 1]) / dz + EPsi2 = EPhi2[0, i, j, k] + EPhi2[1, i, j, k] + Ey[ii, jj, kk] = Ey[ii, jj, kk] + updatecoeffsE[materialEy, 4] * (IRA1 * dHx - IRA * EPsi2) + EPhi2[1, i, j, k] = RE1 * EPhi2[1, i, j, k] + RC1 * dHx - RC1 * EPsi2 + EPhi2[0, i, j, k] = RE0 * EPhi2[0, i, j, k] + RC0 * dHx - RC0 * EPsi2 + ####################################################### # Electric field PML updates - 2nd order - zplus slab # ####################################################### -cpdef void update_pml_2order_electric_zplus( +cpdef void update_pml_2order_electric_HORIPML_zplus( int xs, int xf, int ys, @@ -1270,11 +2424,85 @@ cpdef void update_pml_2order_electric_zplus( 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 void update_pml_2order_electric_MRIPML_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, 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, RA, IRA1, IRA, RB0, RC0, RE0, RF0, RB1, RC1, RE1, RF1, EPsi1, EPsi2 + 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 + RA = ERA[0, k] + ERA[1, k] + IRA1 = (1 / RA) - 1 + IRA = 1 / RA + RB0 = ERB[0, k] + RE0 = ERE[0, k] + RF0 = ERF[0, k] + RC0 = RF0*RB0*IRA + RB1 = ERB[1, k] + RE1 = ERE[1, k] + RF1 = ERF[1, k] + RC1 = RF1*RB1*IRA + # Ex + materialEx = ID[0, ii, jj, kk] + dHy = (Hy[ii, jj, kk] - Hy[ii, jj, kk - 1]) / dz + EPsi1 = EPhi1[0, i, j, k] + EPhi1[1, i, j, k] + Ex[ii, jj, kk] = Ex[ii, jj, kk] - updatecoeffsE[materialEx, 4] * (IRA1 * dHy - IRA * EPsi1) + EPhi1[1, i, j, k] = RE1 * EPhi1[1, i, j, k] + RC1 * dHy - RC1 * EPsi1 + EPhi1[0, i, j, k] = RE0 * EPhi1[0, i, j, k] + RC0 * dHy - RC0 * EPsi1 + # Ey + materialEy = ID[1, ii, jj, kk] + dHx = (Hx[ii, jj, kk] - Hx[ii, jj, kk - 1]) / dz + EPsi2 = EPhi2[0, i, j, k] + EPhi2[1, i, j, k] + Ey[ii, jj, kk] = Ey[ii, jj, kk] + updatecoeffsE[materialEy, 4] * (IRA1 * dHx - IRA * EPsi2) + EPhi2[1, i, j, k] = RE1 * EPhi2[1, i, j, k] + RC1 * dHx - RC1 * EPsi2 + EPhi2[0, i, j, k] = RE0 * EPhi2[0, i, j, k] + RC0 * dHx - RC0 * EPsi2 + ######################################################## # Magnetic field PML updates - 2nd order - xminus slab # ######################################################## -cpdef void update_pml_2order_magnetic_xminus( +cpdef void update_pml_2order_magnetic_HORIPML_xminus( int xs, int xf, int ys, @@ -1344,11 +2572,85 @@ cpdef void update_pml_2order_magnetic_xminus( 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 +cpdef void update_pml_2order_magnetic_MRIPML_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, RA, IRA1, IRA, RB0, RC0, RE0, RF0, RB1, RC1, RE1, RF1, HPsi1, HPsi2 + 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) + RA = HRA[0, i] + HRA[1, i] + IRA1 = (1 / RA) - 1 + IRA = 1 / RA + RB0 = HRB[0, i] + RE0 = HRE[0, i] + RF0 = HRF[0, i] + RC0 = RF0*RB0*IRA + RB1 = HRB[1, i] + RE1 = HRE[1, i] + RF1 = HRF[1, i] + RC1 = RF1*RB1*IRA + 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 + HPsi1 = HPhi1[0, i, j, k] + HPhi1[1, i, j, k] + Hy[ii, jj, kk] = Hy[ii, jj, kk] + updatecoeffsH[materialHy, 4] * (IRA1 * dEz - IRA * HPsi1) + HPhi1[1, i, j, k] = RE1 * HPhi1[1, i, j, k] + RC1 * dEz - RC1 * HPsi1 + HPhi1[0, i, j, k] = RE0 * HPhi1[0, i, j, k] + RC0 * dEz - RC0 * HPsi1 + # Hz + materialHz = ID[5, ii, jj, kk] + dEy = (Ey[ii + 1, jj, kk] - Ey[ii, jj, kk]) / dx + HPsi2 = HPhi2[0, i, j, k] + HPhi2[1, i, j, k] + Hz[ii, jj, kk] = Hz[ii, jj, kk] - updatecoeffsH[materialHz, 4] * (IRA1 * dEy - IRA * HPsi2) + HPhi2[1, i, j, k] = RE1 * HPhi2[1, i, j, k] + RC1 * dEy - RC1 * HPsi2 + HPhi2[0, i, j, k] = RE0 * HPhi2[0, i, j, k] + RC0 * dEy - RC0 * HPsi2 + ####################################################### # Magnetic field PML updates - 2nd order - xplus slab # ####################################################### -cpdef void update_pml_2order_magnetic_xplus( +cpdef void update_pml_2order_magnetic_HORIPML_xplus( int xs, int xf, int ys, @@ -1418,11 +2720,85 @@ cpdef void update_pml_2order_magnetic_xplus( 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 +cpdef void update_pml_2order_magnetic_MRIPML_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, RA, IRA1, IRA, RB0, RC0, RE0, RF0, RB1, RC1, RE1, RF1, HPsi1, HPsi2 + 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 + RA = HRA[0, i] + HRA[1, i] + IRA1 = (1 / RA) - 1 + IRA = 1 / RA + RB0 = HRB[0, i] + RE0 = HRE[0, i] + RF0 = HRF[0, i] + RC0 = RF0*RB0*IRA + RB1 = HRB[1, i] + RE1 = HRE[1, i] + RF1 = HRF[1, i] + RC1 = RF1*RB1*IRA + 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 + HPsi1 = HPhi1[0, i, j, k] + HPhi1[1, i, j, k] + Hy[ii, jj, kk] = Hy[ii, jj, kk] + updatecoeffsH[materialHy, 4] * (IRA1 * dEz - IRA * HPsi1) + HPhi1[1, i, j, k] = RE1 * HPhi1[1, i, j, k] + RC1 * dEz - RC1 * HPsi1 + HPhi1[0, i, j, k] = RE0 * HPhi1[0, i, j, k] + RC0 * dEz - RC0 * HPsi1 + # Hz + materialHz = ID[5, ii, jj, kk] + dEy = (Ey[ii + 1, jj, kk] - Ey[ii, jj, kk]) / dx + HPsi2 = HPhi2[0, i, j, k] + HPhi2[1, i, j, k] + Hz[ii, jj, kk] = Hz[ii, jj, kk] - updatecoeffsH[materialHz, 4] * (IRA1 * dEy - IRA * HPsi2) + HPhi2[1, i, j, k] = RE1 * HPhi2[1, i, j, k] + RC1 * dEy - RC1 * HPsi2 + HPhi2[0, i, j, k] = RE0 * HPhi2[0, i, j, k] + RC0 * dEy - RC0 * HPsi2 + ######################################################## # Magnetic field PML updates - 2nd order - yminus slab # ######################################################## -cpdef void update_pml_2order_magnetic_yminus( +cpdef void update_pml_2order_magnetic_HORIPML_yminus( int xs, int xf, int ys, @@ -1492,11 +2868,85 @@ cpdef void update_pml_2order_magnetic_yminus( 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 +cpdef void update_pml_2order_magnetic_MRIPML_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, RA, IRA1, IRA, RB0, RC0, RE0, RF0, RB1, RC1, RE1, RF1, HPsi1, HPsi2 + 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) + RA = HRA[0, j] + HRA[1, j] + IRA1 = (1 / RA) - 1 + IRA = 1 / RA + RB0 = HRB[0, j] + RE0 = HRE[0, j] + RF0 = HRF[0, j] + RC0 = RF0*RB0*IRA + RB1 = HRB[1, j] + RE1 = HRE[1, j] + RF1 = HRF[1, j] + RC1 = RF1*RB1*IRA + 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 + HPsi1 = HPhi1[0, i, j, k] + HPhi1[1, i, j, k] + Hx[ii, jj, kk] = Hx[ii, jj, kk] - updatecoeffsH[materialHx, 4] * (IRA1 * dEz - IRA * HPsi1) + HPhi1[1, i, j, k] = RE1 * HPhi1[1, i, j, k] + RC1 * dEz - RC1 * HPsi1 + HPhi1[0, i, j, k] = RE0 * HPhi1[0, i, j, k] + RC0 * dEz - RC0 * HPsi1 + # Hz + materialHz = ID[5, ii, jj, kk] + dEx = (Ex[ii, jj + 1, kk] - Ex[ii, jj, kk]) / dy + HPsi2 = HPhi2[0, i, j, k] + HPhi2[1, i, j, k] + Hz[ii, jj, kk] = Hz[ii, jj, kk] + updatecoeffsH[materialHz, 4] * (IRA1 * dEx - IRA * HPsi2) + HPhi2[1, i, j, k] = RE1 * HPhi2[1, i, j, k] + RC1 * dEx - RC1 * HPsi2 + HPhi2[0, i, j, k] = RE0 * HPhi2[0, i, j, k] + RC0 * dEx - RC0 * HPsi2 + ####################################################### # Magnetic field PML updates - 2nd order - yplus slab # ####################################################### -cpdef void update_pml_2order_magnetic_yplus( +cpdef void update_pml_2order_magnetic_HORIPML_yplus( int xs, int xf, int ys, @@ -1566,11 +3016,85 @@ cpdef void update_pml_2order_magnetic_yplus( 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 +cpdef void update_pml_2order_magnetic_MRIPML_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, RA, IRA1, IRA, RB0, RC0, RE0, RF0, RB1, RC1, RE1, RF1, HPsi1, HPsi2 + 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 + RA = HRA[0, j] + HRA[1, j] + IRA1 = (1 / RA) - 1 + IRA = 1 / RA + RB0 = HRB[0, j] + RE0 = HRE[0, j] + RF0 = HRF[0, j] + RC0 = RF0*RB0*IRA + RB1 = HRB[1, j] + RE1 = HRE[1, j] + RF1 = HRF[1, j] + RC1 = RF1*RB1*IRA + 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 + HPsi1 = HPhi1[0, i, j, k] + HPhi1[1, i, j, k] + Hx[ii, jj, kk] = Hx[ii, jj, kk] - updatecoeffsH[materialHx, 4] * (IRA1 * dEz - IRA * HPsi1) + HPhi1[1, i, j, k] = RE1 * HPhi1[1, i, j, k] + RC1 * dEz - RC1 * HPsi1 + HPhi1[0, i, j, k] = RE0 * HPhi1[0, i, j, k] + RC0 * dEz - RC0 * HPsi1 + # Hz + materialHz = ID[5, ii, jj, kk] + dEx = (Ex[ii, jj + 1, kk] - Ex[ii, jj, kk]) / dy + HPsi2 = HPhi2[0, i, j, k] + HPhi2[1, i, j, k] + Hz[ii, jj, kk] = Hz[ii, jj, kk] + updatecoeffsH[materialHz, 4] * (IRA1 * dEx - IRA * HPsi2) + HPhi2[1, i, j, k] = RE1 * HPhi2[1, i, j, k] + RC1 * dEx - RC1 * HPsi2 + HPhi2[0, i, j, k] = RE0 * HPhi2[0, i, j, k] + RC0 * dEx - RC0 * HPsi2 + ######################################################## # Magnetic field PML updates - 2nd order - zminus slab # ######################################################## -cpdef void update_pml_2order_magnetic_zminus( +cpdef void update_pml_2order_magnetic_HORIPML_zminus( int xs, int xf, int ys, @@ -1640,11 +3164,85 @@ cpdef void update_pml_2order_magnetic_zminus( 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 +cpdef void update_pml_2order_magnetic_MRIPML_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, RA, IRA1, IRA, RB0, RC0, RE0, RF0, RB1, RC1, RE1, RF1, HPsi1, HPsi2 + 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) + RA = HRA[0, k] + HRA[1, k] + IRA1 = (1 / RA) - 1 + IRA = 1 / RA + RB0 = HRB[0, k] + RE0 = HRE[0, k] + RF0 = HRF[0, k] + RC0 = RF0*RB0*IRA + RB1 = HRB[1, k] + RE1 = HRE[1, k] + RF1 = HRF[1, k] + RC1 = RF1*RB1*IRA + # Hx + materialHx = ID[3, ii, jj, kk] + dEy = (Ey[ii, jj, kk + 1] - Ey[ii, jj, kk]) / dz + HPsi1 = HPhi1[0, i, j, k] + HPhi1[1, i, j, k] + Hx[ii, jj, kk] = Hx[ii, jj, kk] + updatecoeffsH[materialHx, 4] * (IRA1 * dEy - IRA * HPsi1) + HPhi1[1, i, j, k] = RE1 * HPhi1[1, i, j, k] + RC1 * dEy - RC1 * HPsi1 + HPhi1[0, i, j, k] = RE0 * HPhi1[0, i, j, k] + RC0 * dEy - RC0 * HPsi1 + # Hy + materialHy = ID[4, ii, jj, kk] + dEx = (Ex[ii, jj, kk + 1] - Ex[ii, jj, kk]) / dz + HPsi2 = HPhi2[0, i, j, k] + HPhi2[1, i, j, k] + Hy[ii, jj, kk] = Hy[ii, jj, kk] - updatecoeffsH[materialHy, 4] * (IRA1 * dEx - IRA * HPsi2) + HPhi2[1, i, j, k] = RE1 * HPhi2[1, i, j, k] + RC1 * dEx - RC1 * HPsi2 + HPhi2[0, i, j, k] = RE0 * HPhi2[0, i, j, k] + RC0 * dEx - RC0 * HPsi2 + ####################################################### # Magnetic field PML updates - 2nd order - zplus slab # ####################################################### -cpdef void update_pml_2order_magnetic_zplus( +cpdef void update_pml_2order_magnetic_HORIPML_zplus( int xs, int xf, int ys, @@ -1713,3 +3311,77 @@ cpdef void update_pml_2order_magnetic_zplus( 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 + +cpdef void update_pml_2order_magnetic_MRIPML_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, RA, IRA1, IRA, RB0, RC0, RE0, RF0, RB1, RC1, RE1, RF1, HPsi1, HPsi2 + 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 + RA = HRA[0, k] + HRA[1, k] + IRA1 = (1 / RA) - 1 + IRA = 1 / RA + RB0 = HRB[0, k] + RE0 = HRE[0, k] + RF0 = HRF[0, k] + RC0 = RF0*RB0*IRA + RB1 = HRB[1, k] + RE1 = HRE[1, k] + RF1 = HRF[1, k] + RC1 = RF1*RB1*IRA + # Hx + materialHx = ID[3, ii, jj, kk] + dEy = (Ey[ii, jj, kk + 1] - Ey[ii, jj, kk]) / dz + HPsi1 = HPhi1[0, i, j, k] + HPhi1[1, i, j, k] + Hx[ii, jj, kk] = Hx[ii, jj, kk] + updatecoeffsH[materialHx, 4] * (IRA1 * dEy - IRA * HPsi1) + HPhi1[1, i, j, k] = RE1 * HPhi1[1, i, j, k] + RC1 * dEy - RC1 * HPsi1 + HPhi1[0, i, j, k] = RE0 * HPhi1[0, i, j, k] + RC0 * dEy - RC0 * HPsi1 + # Hy + materialHy = ID[4, ii, jj, kk] + dEx = (Ex[ii, jj, kk + 1] - Ex[ii, jj, kk]) / dz + HPsi2 = HPhi2[0, i, j, k] + HPhi2[1, i, j, k] + Hy[ii, jj, kk] = Hy[ii, jj, kk] - updatecoeffsH[materialHy, 4] * (IRA1 * dEx - IRA * HPsi2) + HPhi2[1, i, j, k] = RE1 * HPhi2[1, i, j, k] + RC1 * dEx - RC1 * HPsi2 + HPhi2[0, i, j, k] = RE0 * HPhi2[0, i, j, k] + RC0 * dEx - RC0 * HPsi2 diff --git a/user_models/cylinder_Ascan_2D.in b/user_models/cylinder_Ascan_2D.in index edd88986..5da2702e 100644 --- a/user_models/cylinder_Ascan_2D.in +++ b/user_models/cylinder_Ascan_2D.in @@ -1,15 +1,20 @@ #title: A-scan from a metal cylinder buried in a dielectric half-space #domain: 0.240 0.210 0.002 #dx_dy_dz: 0.002 0.002 0.002 -#time_window: 3e-9 +#time_window: 300e-9 + +#pml_cfs: constant forward 0 0 cubic forward 1 4 dodeka forward 0 13.7934 +#pml_cfs: constant forward 0.1 0.1 constant forward 0 0 cubic forward 0 8.4883 +#material: 8 0.01 1 0 wet #material: 6 0 1 0 half_space #waveform: ricker 1 1.5e9 my_ricker #hertzian_dipole: z 0.100 0.170 0 my_ricker #rx: 0.140 0.170 0 - +#pml_formulation: MRIPML #box: 0 0 0 0.240 0.170 0.002 half_space #cylinder: 0.120 0.080 0 0.120 0.080 0.002 0.010 pec -#geometry_view: 0 0 0 0.240 0.210 0.002 0.002 0.002 0.002 cylinder_half_space n \ No newline at end of file +#geometry_view: 0 0 0 0.240 0.210 0.002 0.002 0.002 0.002 cylinder_half_space n +#messages: y \ No newline at end of file diff --git a/user_models/cylinder_Ascan_2D.out b/user_models/cylinder_Ascan_2D.out index 1b0d394c..de9e3368 100644 Binary files a/user_models/cylinder_Ascan_2D.out and b/user_models/cylinder_Ascan_2D.out differ