From 6b36b597bde68978e6325c9c1a3f3845f44fc4a0 Mon Sep 17 00:00:00 2001 From: fake-name Date: Tue, 25 Apr 2017 16:40:55 -0700 Subject: [PATCH] Break the comments to be mostly 80 chars wide, with a 100 char ceiling (mostly, probably missed some). No actual code modification yet, only whitespace. --- gprMax/constants.py | 3 +- gprMax/fields_outputs.py | 3 +- gprMax/fields_updates.pyx | 100 ++++++++++- gprMax/fractals.py | 10 +- gprMax/geometry_outputs.py | 17 +- gprMax/geometry_primitives.pyx | 262 +++++++++++++++++++++++++--- gprMax/gprMax.py | 50 ++++-- gprMax/grid.py | 28 ++- gprMax/input_cmd_funcs.py | 72 +++++--- gprMax/input_cmds_file.py | 44 +++-- gprMax/input_cmds_geometry.py | 14 +- gprMax/input_cmds_multiuse.py | 4 +- gprMax/materials.py | 19 ++- gprMax/model_build_run.py | 40 +++-- gprMax/optimisation_taguchi.py | 52 ++++-- gprMax/pml.py | 7 +- gprMax/pml_1order_update.pyx | 300 +++++++++++++++++++++++++++++++-- gprMax/pml_2order_update.pyx | 300 +++++++++++++++++++++++++++++++-- gprMax/snapshots.py | 3 +- gprMax/sources.py | 24 ++- gprMax/utilities.py | 10 +- 21 files changed, 1197 insertions(+), 165 deletions(-) diff --git a/gprMax/constants.py b/gprMax/constants.py index cf5f2e46..a98f4498 100644 --- a/gprMax/constants.py +++ b/gprMax/constants.py @@ -27,7 +27,8 @@ z0 = np.sqrt(m0 / e0) # Data types: # Solid and ID arrays use 32-bit integers (0 to 4294967295) # Rigid arrays use 8-bit integers (the smallest available type to store true/false) -# Fractal and dispersive coefficient arrays use complex numbers (complextype) which are represented as two floats +# Fractal and dispersive coefficient arrays use complex numbers (complextype) +# which are represented as two floats # Main field arrays use floats (floattype) and complex numbers (complextype) # Single precision diff --git a/gprMax/fields_outputs.py b/gprMax/fields_outputs.py index 0db3ad57..03bf8ab6 100644 --- a/gprMax/fields_outputs.py +++ b/gprMax/fields_outputs.py @@ -75,7 +75,8 @@ def write_hdf5_outputfile(outputfile, Ex, Ey, Ez, Hx, Hy, Hz, G): grp.attrs['Type'] = type(src).__name__ grp.attrs['Position'] = (src.xcoord * G.dx, src.ycoord * G.dy, src.zcoord * G.dz) - # Create group for transmission lines; add positional data, line resistance and line discretisation attributes; write arrays for line voltages and currents + # Create group for transmission lines; add positional data, line resistance and + # line discretisation attributes; write arrays for line voltages and currents for tlindex, tl in enumerate(G.transmissionlines): grp = f.create_group('/tls/tl' + str(tlindex + 1)) grp.attrs['Position'] = (tl.xcoord * G.dx, tl.ycoord * G.dy, tl.zcoord * G.dz) diff --git a/gprMax/fields_updates.pyx b/gprMax/fields_updates.pyx index 5ff58514..b742e54a 100644 --- a/gprMax/fields_updates.pyx +++ b/gprMax/fields_updates.pyx @@ -27,7 +27,20 @@ from gprMax.constants cimport complextype_t ############################################### # Electric field updates - standard materials # ############################################### -cpdef void update_electric(int nx, int ny, int nz, int nthreads, floattype_t[:, ::1] updatecoeffsE, np.uint32_t[:, :, :, ::1] ID, floattype_t[:, :, ::1] Ex, floattype_t[:, :, ::1] Ey, floattype_t[:, :, ::1] Ez, floattype_t[:, :, ::1] Hx, floattype_t[:, :, ::1] Hy, floattype_t[:, :, ::1] Hz): +cpdef void update_electric( + int nx, + int ny, + int nz, + int nthreads, + floattype_t[:, ::1] updatecoeffsE, + np.uint32_t[:, :, :, ::1] ID, + floattype_t[:, :, ::1] Ex, + floattype_t[:, :, ::1] Ey, + floattype_t[:, :, ::1] Ez, + floattype_t[:, :, ::1] Hx, + floattype_t[:, :, ::1] Hy, + floattype_t[:, :, ::1] Hz + ): """This function updates the electric field components. Args: @@ -97,7 +110,25 @@ cpdef void update_electric(int nx, int ny, int nz, int nthreads, floattype_t[:, ################################################# # Electric field updates - dispersive materials # ################################################# -cpdef void update_electric_dispersive_multipole_A(int nx, int ny, int nz, int nthreads, int maxpoles, floattype_t[:, ::1] updatecoeffsE, complextype_t[:, ::1] updatecoeffsdispersive, np.uint32_t[:, :, :, ::1] ID, complextype_t[:, :, :, ::1] Tx, complextype_t[:, :, :, ::1] Ty, complextype_t[:, :, :, ::1] Tz, floattype_t[:, :, ::1] Ex, floattype_t[:, :, ::1] Ey, floattype_t[:, :, ::1] Ez, floattype_t[:, :, ::1] Hx, floattype_t[:, :, ::1] Hy, floattype_t[:, :, ::1] Hz): +cpdef void update_electric_dispersive_multipole_A( + int nx, + int ny, + int nz, + int nthreads, + int maxpoles, + floattype_t[:, ::1] updatecoeffsE, + complextype_t[:, ::1] updatecoeffsdispersive, + np.uint32_t[:, :, :, ::1] ID, + complextype_t[:, :, :, ::1] Tx, + complextype_t[:, :, :, ::1] Ty, + complextype_t[:, :, :, ::1] Tz, + floattype_t[:, :, ::1] Ex, + floattype_t[:, :, ::1] Ey, + floattype_t[:, :, ::1] Ez, + floattype_t[:, :, ::1] Hx, + floattype_t[:, :, ::1] Hy, + floattype_t[:, :, ::1] Hz + ): """This function updates the electric field components when dispersive materials (with multiple poles) are present. Args: @@ -149,7 +180,21 @@ cpdef void update_electric_dispersive_multipole_A(int nx, int ny, int nz, int nt -cpdef void update_electric_dispersive_multipole_B(int nx, int ny, int nz, int nthreads, int maxpoles, complextype_t[:, ::1] updatecoeffsdispersive, np.uint32_t[:, :, :, ::1] ID, complextype_t[:, :, :, ::1] Tx, complextype_t[:, :, :, ::1] Ty, complextype_t[:, :, :, ::1] Tz, floattype_t[:, :, ::1] Ex, floattype_t[:, :, ::1] Ey, floattype_t[:, :, ::1] Ez): +cpdef void update_electric_dispersive_multipole_B( + int nx, + int ny, + int nz, + int nthreads, + int maxpoles, + complextype_t[:, ::1] updatecoeffsdispersive, + np.uint32_t[:, :, :, ::1] ID, + complextype_t[:, :, :, ::1] Tx, + complextype_t[:, :, :, ::1] Ty, + complextype_t[:, :, :, ::1] Tz, + floattype_t[:, :, ::1] Ex, + floattype_t[:, :, ::1] Ey, + floattype_t[:, :, ::1] Ez + ): """This function updates a temporary dispersive material array when disperisive materials (with multiple poles) are present. Args: @@ -190,7 +235,24 @@ cpdef void update_electric_dispersive_multipole_B(int nx, int ny, int nz, int nt Tz[pole, i, j, k] = Tz[pole, i, j, k] - updatecoeffsdispersive[material, 2 + (pole * 3)] * Ez[i, j, k] -cpdef void update_electric_dispersive_1pole_A(int nx, int ny, int nz, int nthreads, floattype_t[:, ::1] updatecoeffsE, complextype_t[:, ::1] updatecoeffsdispersive, np.uint32_t[:, :, :, ::1] ID, complextype_t[:, :, :, ::1] Tx, complextype_t[:, :, :, ::1] Ty, complextype_t[:, :, :, ::1] Tz, floattype_t[:, :, ::1] Ex, floattype_t[:, :, ::1] Ey, floattype_t[:, :, ::1] Ez, floattype_t[:, :, ::1] Hx, floattype_t[:, :, ::1] Hy, floattype_t[:, :, ::1] Hz): +cpdef void update_electric_dispersive_1pole_A( + int nx, + int ny, + int nz, + int nthreads, + floattype_t[:, ::1] updatecoeffsE, + complextype_t[:, ::1] updatecoeffsdispersive, + np.uint32_t[:, :, :, ::1] ID, + complextype_t[:, :, :, ::1] Tx, + complextype_t[:, :, :, ::1] Ty, + complextype_t[:, :, :, ::1] Tz, + floattype_t[:, :, ::1] Ex, + floattype_t[:, :, ::1] Ey, + floattype_t[:, :, ::1] Ez, + floattype_t[:, :, ::1] Hx, + floattype_t[:, :, ::1] Hy, + floattype_t[:, :, ::1] Hz + ): """This function updates the electric field components when dispersive materials (with 1 pole) are present. Args: @@ -234,7 +296,20 @@ cpdef void update_electric_dispersive_1pole_A(int nx, int ny, int nz, int nthrea Ez[i, j, k] = updatecoeffsE[material, 0] * Ez[i, j, k] + updatecoeffsE[material, 1] * (Hy[i, j, k] - Hy[i - 1, j, k]) - updatecoeffsE[material, 2] * (Hx[i, j, k] - Hx[i, j - 1, k]) - updatecoeffsE[material, 4] * phi -cpdef void update_electric_dispersive_1pole_B(int nx, int ny, int nz, int nthreads, complextype_t[:, ::1] updatecoeffsdispersive, np.uint32_t[:, :, :, ::1] ID, complextype_t[:, :, :, ::1] Tx, complextype_t[:, :, :, ::1] Ty, complextype_t[:, :, :, ::1] Tz, floattype_t[:, :, ::1] Ex, floattype_t[:, :, ::1] Ey, floattype_t[:, :, ::1] Ez): +cpdef void update_electric_dispersive_1pole_B( + int nx, + int ny, + int nz, + int nthreads, + complextype_t[:, ::1] updatecoeffsdispersive, + np.uint32_t[:, :, :, ::1] ID, + complextype_t[:, :, :, ::1] Tx, + complextype_t[:, :, :, ::1] Ty, + complextype_t[:, :, :, ::1] Tz, + floattype_t[:, :, ::1] Ex, + floattype_t[:, :, ::1] Ey, + floattype_t[:, :, ::1] Ez + ): """This function updates a temporary dispersive material array when disperisive materials (with 1 pole) are present. Args: @@ -274,7 +349,20 @@ cpdef void update_electric_dispersive_1pole_B(int nx, int ny, int nz, int nthrea ########################## # Magnetic field updates # ########################## -cpdef void update_magnetic(int nx, int ny, int nz, int nthreads, floattype_t[:, ::1] updatecoeffsH, np.uint32_t[:, :, :, ::1] ID, floattype_t[:, :, ::1] Ex, floattype_t[:, :, ::1] Ey, floattype_t[:, :, ::1] Ez, floattype_t[:, :, ::1] Hx, floattype_t[:, :, ::1] Hy, floattype_t[:, :, ::1] Hz): +cpdef void update_magnetic( + int nx, + int ny, + int nz, + int nthreads, + floattype_t[:, ::1] updatecoeffsH, + np.uint32_t[:, :, :, ::1] ID, + floattype_t[:, :, ::1] Ex, + floattype_t[:, :, ::1] Ey, + floattype_t[:, :, ::1] Ez, + floattype_t[:, :, ::1] Hx, + floattype_t[:, :, ::1] Hy, + floattype_t[:, :, ::1] Hz + ): """This function updates the magnetic field components. Args: diff --git a/gprMax/fractals.py b/gprMax/fractals.py index d4123bd8..13e61c0f 100644 --- a/gprMax/fractals.py +++ b/gprMax/fractals.py @@ -33,7 +33,8 @@ class FractalSurface(object): def __init__(self, xs, xf, ys, yf, zs, zf, dimension): """ Args: - xs, xf, ys, yf, zs, zf (float): Extent of the fractal surface (one pair of coordinates must be equal to correctly define a surface). + xs, xf, ys, yf, zs, zf (float): Extent of the fractal surface (one pair of + coordinates must be equal to correctly define a surface). dimension (float): Fractal dimension that controls the fractal distribution. """ @@ -100,7 +101,8 @@ class FractalSurface(object): fractalmin = np.amin(self.fractalsurface) fractalmax = np.amax(self.fractalsurface) fractalrange = fractalmax - fractalmin - self.fractalsurface = self.fractalsurface * ((self.fractalrange[1] - self.fractalrange[0]) / fractalrange) + self.fractalrange[0] - ((self.fractalrange[1] - self.fractalrange[0]) / fractalrange) * fractalmin + self.fractalsurface = self.fractalsurface * ((self.fractalrange[1] - self.fractalrange[0]) / fractalrange) \ + + self.fractalrange[0] - ((self.fractalrange[1] - self.fractalrange[0]) / fractalrange) * fractalmin class FractalVolume(object): @@ -189,7 +191,9 @@ class FractalVolume(object): self.fractalvolume[:, j, k] = np.digitize(self.fractalvolume[:, j, k], bins, right=True) def generate_volume_mask(self): - """Generate a 3D volume to use as a mask for adding rough surfaces, water and grass/roots. Zero signifies the mask is not set, one signifies the mask is set.""" + """ + Generate a 3D volume to use as a mask for adding rough surfaces, water and grass/roots. + Zero signifies the mask is not set, one signifies the mask is set.""" self.mask = np.zeros((self.nx, self.ny, self.nz), dtype=np.int8) maskxs = self.originalxs - self.xs diff --git a/gprMax/geometry_outputs.py b/gprMax/geometry_outputs.py index b1711195..e2357a87 100644 --- a/gprMax/geometry_outputs.py +++ b/gprMax/geometry_outputs.py @@ -41,7 +41,8 @@ class GeometryView(object): xs, xf, ys, yf, zs, zf (int): Extent of the volume in cells. dx, dy, dz (int): Spatial discretisation in cells. filename (str): Filename to save to. - fileext (str): File extension of VTK file - either '.vti' for a per cell geometry view, or '.vtp' for a per cell edge geometry view. + fileext (str): File extension of VTK file - either '.vti' for a per cell + geometry view, or '.vtp' for a per cell edge geometry view. """ self.xs = xs @@ -83,7 +84,8 @@ class GeometryView(object): self.datawritesize = np.dtype(np.float32).itemsize * self.vtk_numpoints * self.vtk_numpoint_components + np.dtype(np.uint32).itemsize * self.vtk_numlines * self.vtk_numline_components + np.dtype(np.uint32).itemsize * self.vtk_numlines + np.dtype(np.uint32).itemsize * self.vtk_numlines def set_filename(self, appendmodelnumber, G): - """Construct filename from user-supplied name and model run number. + """ + Construct filename from user-supplied name and model run number. Args: appendmodelnumber (str): Text to append to filename. @@ -94,7 +96,9 @@ class GeometryView(object): self.filename += self.fileext def write_vtk(self, G, pbar): - """Writes the geometry information to a VTK file. Either ImageData (.vti) for a per-cell geometry view, or PolygonalData (.vtp) for a per-cell-edge geometry view. + """ + Writes the geometry information to a VTK file. Either ImageData (.vti) for a + per-cell geometry view, or PolygonalData (.vtp) for a per-cell-edge geometry view. N.B. No Python 3 support for VTK at time of writing (03/2015) @@ -104,7 +108,8 @@ class GeometryView(object): """ if self.fileext == '.vti': - # Create arrays and add numeric IDs for PML, sources and receivers (0 is not set, 1 is PML, srcs and rxs numbered thereafter) + # Create arrays and add numeric IDs for PML, sources and receivers + # (0 is not set, 1 is PML, srcs and rxs numbered thereafter) self.srcs_pml = np.zeros((G.nx + 1, G.ny + 1, G.nz + 1), dtype=np.int8) self.rxs = np.zeros((G.nx + 1, G.ny + 1, G.nz + 1), dtype=np.int8) for pml in G.pmls: @@ -256,7 +261,9 @@ class GeometryView(object): self.write_gprmax_info(f, G, materialsonly=True) def write_gprmax_info(self, f, G, materialsonly=False): - """Writes gprMax specific information relating material, source, and receiver names to numeric identifiers. + """ + Writes gprMax specific information relating material, source, + and receiver names to numeric identifiers. Args: f (filehandle): VTK file. diff --git a/gprMax/geometry_primitives.pyx b/gprMax/geometry_primitives.pyx index 6f2cc73f..d434651d 100644 --- a/gprMax/geometry_primitives.pyx +++ b/gprMax/geometry_primitives.pyx @@ -34,7 +34,12 @@ from gprMax.yee_cell_setget_rigid cimport unset_rigid_H np.seterr(divide='raise') -cpdef bint are_clockwise(float v1x, float v1y, float v2x, float v2y): +cpdef bint are_clockwise( + float v1x, + float v1y, + float v2x, + float v2y + ): """Find if vector 2 is clockwise relative to vector 1. Args: @@ -47,7 +52,11 @@ cpdef bint are_clockwise(float v1x, float v1y, float v2x, float v2y): return -v1x*v2y + v1y*v2x > 0 -cpdef bint is_within_radius(float vx, float vy, float radius): +cpdef bint is_within_radius( + float vx, + float vy, + float radius + ): """Check if the point is within a given radius of the centre of the circle. Args: @@ -61,7 +70,15 @@ cpdef bint is_within_radius(float vx, float vy, float radius): return vx*vx + vy*vy <= radius*radius -cpdef bint is_inside_sector(float px, float py, float ctrx, float ctry, float sectorstartangle, float sectorangle, float radius): +cpdef bint is_inside_sector( + float px, + float py, + float ctrx, + float ctry, + float sectorstartangle, + float sectorangle, + float radius + ): """For a point to be inside a circular sector, it has to meet the following tests: It has to be positioned anti-clockwise from the start "arm" of the sector It has to be positioned clockwise from the end arm of the sector @@ -89,10 +106,16 @@ cpdef bint is_inside_sector(float px, float py, float ctrx, float ctry, float se relpoint1 = px - ctrx relpoint2 = py - ctry - return not are_clockwise(sectorstart1, sectorstart2, relpoint1, relpoint2) and are_clockwise(sectorend1, sectorend2, relpoint1, relpoint2) and is_within_radius(relpoint1, relpoint2, radius) + return (not are_clockwise(sectorstart1, sectorstart2, relpoint1, relpoint2) + and are_clockwise(sectorend1, sectorend2, relpoint1, relpoint2) + and is_within_radius(relpoint1, relpoint2, radius)) -cpdef bint point_in_polygon(float px, float py, list polycoords): +cpdef bint point_in_polygon( + float px, + float py, + list polycoords + ): """Calculates, using a ray casting algorithm, whether a point lies within a polygon. Args: @@ -141,7 +164,15 @@ cpdef bint point_in_polygon(float px, float py, list polycoords): return inside -cpdef void build_edge_x(int i, int j, int k, int numIDx, np.int8_t[:, :, :, ::1] rigidE, np.int8_t[:, :, :, ::1] rigidH, np.uint32_t[:, :, :, ::1] ID): +cpdef void build_edge_x( + int i, + int j, + int k, + int numIDx, + np.int8_t[:, :, :, ::1] rigidE, + np.int8_t[:, :, :, ::1] rigidH, + np.uint32_t[:, :, :, ::1] ID + ): """Set x-orientated edges in the rigid and ID arrays for a Yee voxel. Args: @@ -154,7 +185,15 @@ cpdef void build_edge_x(int i, int j, int k, int numIDx, np.int8_t[:, :, :, ::1] ID[0, i, j, k] = numIDx -cpdef void build_edge_y(int i, int j, int k, int numIDy, np.int8_t[:, :, :, ::1] rigidE, np.int8_t[:, :, :, ::1] rigidH, np.uint32_t[:, :, :, ::1] ID): +cpdef void build_edge_y( + int i, + int j, + int k, + int numIDy, + np.int8_t[:, :, :, ::1] rigidE, + np.int8_t[:, :, :, ::1] rigidH, + np.uint32_t[:, :, :, ::1] ID + ): """Set y-orientated edges in the rigid and ID arrays for a Yee voxel. Args: @@ -167,7 +206,15 @@ cpdef void build_edge_y(int i, int j, int k, int numIDy, np.int8_t[:, :, :, ::1] ID[1, i, j, k] = numIDy -cpdef void build_edge_z(int i, int j, int k, int numIDz, np.int8_t[:, :, :, ::1] rigidE, np.int8_t[:, :, :, ::1] rigidH, np.uint32_t[:, :, :, ::1] ID): +cpdef void build_edge_z( + int i, + int j, + int k, + int numIDz, + np.int8_t[:, :, :, ::1] rigidE, + np.int8_t[:, :, :, ::1] rigidH, + np.uint32_t[:, :, :, ::1] ID + ): """Set z-orientated edges in the rigid and ID arrays for a Yee voxel. Args: @@ -180,7 +227,16 @@ cpdef void build_edge_z(int i, int j, int k, int numIDz, np.int8_t[:, :, :, ::1] ID[2, i, j, k] = numIDz -cpdef void build_face_yz(int i, int j, int k, int numIDy, int numIDz, np.int8_t[:, :, :, ::1] rigidE, np.int8_t[:, :, :, ::1] rigidH, np.uint32_t[:, :, :, ::1] ID): +cpdef void build_face_yz( + int i, + int j, + int k, + int numIDy, + int numIDz, + np.int8_t[:, :, :, ::1] rigidE, + np.int8_t[:, :, :, ::1] rigidH, + np.uint32_t[:, :, :, ::1] ID + ): """Set the edges of the yz-plane face of a Yell cell in the rigid and ID arrays. Args: @@ -207,7 +263,16 @@ cpdef void build_face_yz(int i, int j, int k, int numIDy, int numIDz, np.int8_t[ ID[5, i, j + 1, k] = numIDz -cpdef void build_face_xz(int i, int j, int k, int numIDx, int numIDz, np.int8_t[:, :, :, ::1] rigidE, np.int8_t[:, :, :, ::1] rigidH, np.uint32_t[:, :, :, ::1] ID): +cpdef void build_face_xz( + int i, + int j, + int k, + int numIDx, + int numIDz, + np.int8_t[:, :, :, ::1] rigidE, + np.int8_t[:, :, :, ::1] rigidH, + np.uint32_t[:, :, :, ::1] ID + ): """Set the edges of the xz-plane face of a Yell cell in the rigid and ID arrays. Args: @@ -234,7 +299,16 @@ cpdef void build_face_xz(int i, int j, int k, int numIDx, int numIDz, np.int8_t[ ID[5, i + 1, j, k] = numIDz -cpdef void build_face_xy(int i, int j, int k, int numIDx, int numIDy, np.int8_t[:, :, :, ::1] rigidE, np.int8_t[:, :, :, ::1] rigidH, np.uint32_t[:, :, :, ::1] ID): +cpdef void build_face_xy( + int i, + int j, + int k, + int numIDx, + int numIDy, + np.int8_t[:, :, :, ::1] rigidE, + np.int8_t[:, :, :, ::1] rigidH, + np.uint32_t[:, :, :, ::1] ID + ): """Set the edges of the xy-plane face of a Yell cell in the rigid and ID arrays. Args: @@ -261,7 +335,20 @@ cpdef void build_face_xy(int i, int j, int k, int numIDx, int numIDy, np.int8_t[ ID[4, i + 1, j, k] = numIDy -cpdef void build_voxel(int i, int j, int k, int numID, int numIDx, int numIDy, int numIDz, bint averaging, np.uint32_t[:, :, ::1] solid, np.int8_t[:, :, :, ::1] rigidE, np.int8_t[:, :, :, ::1] rigidH, np.uint32_t[:, :, :, ::1] ID): +cpdef void build_voxel( + int i, + int j, + int k, + int numID, + int numIDx, + int numIDy, + int numIDz, + bint averaging, + np.uint32_t[:, :, ::1] solid, + np.int8_t[:, :, :, ::1] rigidE, + np.int8_t[:, :, :, ::1] rigidH, + np.uint32_t[:, :, :, ::1] ID + ): """Set values in the solid, rigid and ID arrays for a Yee voxel. Args: @@ -312,11 +399,38 @@ cpdef void build_voxel(int i, int j, int k, int numID, int numIDx, int numIDy, i ID[5, i, j + 1, k] = numIDz -cpdef void build_triangle(float x1, float y1, float z1, float x2, float y2, float z2, float x3, float y3, float z3, str normal, float thickness, float dx, float dy, float dz, int numID, int numIDx, int numIDy, int numIDz, bint averaging, np.uint32_t[:, :, ::1] solid, np.int8_t[:, :, :, ::1] rigidE, np.int8_t[:, :, :, ::1] rigidH, np.uint32_t[:, :, :, ::1] ID): - """Builds #triangle and #triangular_prism commands which sets values in the solid, rigid and ID arrays for a Yee voxel. +cpdef void build_triangle( + float x1, + float y1, + float z1, + float x2, + float y2, + float z2, + float x3, + float y3, + float z3, + str normal, + float thickness, + float dx, + float dy, + float dz, + int numID, + int numIDx, + int numIDy, + int numIDz, + bint averaging, + np.uint32_t[:, :, ::1] solid, + np.int8_t[:, :, :, ::1] rigidE, + np.int8_t[:, :, :, ::1] rigidH, + np.uint32_t[:, :, :, ::1] ID + ): + """ + Builds #triangle and #triangular_prism commands which sets values in the + solid, rigid and ID arrays for a Yee voxel. Args: - x1, y1, z1, x2, y2, z2, x3, y3, z3 (float): Coordinates of the vertices of the triangular prism. + x1, y1, z1, x2, y2, z2, x3, y3, z3 (float): Coordinates of the vertices + of the triangular prism. normal (char): Normal direction to the plane of the triangular prism. thickness (float): Thickness of the triangular prism. dx, dy, dz (float): Spatial discretisation. @@ -396,8 +510,35 @@ cpdef void build_triangle(float x1, float y1, float z1, float x2, float y2, floa build_voxel(i, j, k, numID, numIDx, numIDy, numIDz, averaging, solid, rigidE, rigidH, ID) -cpdef void build_cylindrical_sector(float ctr1, float ctr2, int level, float sectorstartangle, float sectorangle, float radius, str normal, float thickness, float dx, float dy, float dz, int numID, int numIDx, int numIDy, int numIDz, bint averaging, np.uint32_t[:, :, ::1] solid, np.int8_t[:, :, :, ::1] rigidE, np.int8_t[:, :, :, ::1] rigidH, np.uint32_t[:, :, :, ::1] ID): - """Builds #cylindrical_sector commands which sets values in the solid, rigid and ID arrays for a Yee voxel. It defines a sector of cylinder given by the direction of the axis of the coordinates of the cylinder face centre, depth coordinates, sector start point, sector angle, and sector radius. N.B Assumes sector start is always clockwise from sector end, i.e. sector defined in an anti-clockwise direction. +cpdef void build_cylindrical_sector( + float ctr1, + float ctr2, + int level, + float sectorstartangle, + float sectorangle, + float radius, + str normal, + float thickness, + float dx, + float dy, + float dz, + int numID, + int numIDx, + int numIDy, + int numIDz, + bint averaging, + np.uint32_t[:, :, ::1] solid, + np.int8_t[:, :, :, ::1] rigidE, + np.int8_t[:, :, :, ::1] rigidH, + np.uint32_t[:, :, :, ::1] ID + ): + """ + Builds #cylindrical_sector commands which sets values in the solid, rigid + and ID arrays for a Yee voxel. It defines a sector of cylinder given by the + direction of the axis of the coordinates of the cylinder face centre, depth + coordinates, sector start point, sector angle, and sector radius. N.B + Assumes sector start is always clockwise from sector end, i.e. sector + defined in an anti-clockwise direction. Args: ctr1, ctr2 (float): Coordinates of centre of circle. @@ -468,7 +609,23 @@ cpdef void build_cylindrical_sector(float ctr1, float ctr2, int level, float sec build_voxel(x, y, z, numID, numIDx, numIDy, numIDz, averaging, solid, rigidE, rigidH, ID) -cpdef void build_box(int xs, int xf, int ys, int yf, int zs, int zf, int numID, int numIDx, int numIDy, int numIDz, bint averaging, np.uint32_t[:, :, ::1] solid, np.int8_t[:, :, :, ::1] rigidE, np.int8_t[:, :, :, ::1] rigidH, np.uint32_t[:, :, :, ::1] ID): +cpdef void build_box( + int xs, + int xf, + int ys, + int yf, + int zs, + int zf, + int numID, + int numIDx, + int numIDy, + int numIDz, + bint averaging, + np.uint32_t[:, :, ::1] solid, + np.int8_t[:, :, :, ::1] rigidE, + np.int8_t[:, :, :, ::1] rigidH, + np.uint32_t[:, :, :, ::1] ID + ): """Builds #box commands which sets values in the solid, rigid and ID arrays. Args: @@ -526,7 +683,27 @@ cpdef void build_box(int xs, int xf, int ys, int yf, int zs, int zf, int numID, ID[5, i, j, k] = numIDz -cpdef void build_cylinder(float x1, float y1, float z1, float x2, float y2, float z2, float r, float dx, float dy, float dz, int numID, int numIDx, int numIDy, int numIDz, bint averaging, np.uint32_t[:, :, ::1] solid, np.int8_t[:, :, :, ::1] rigidE, np.int8_t[:, :, :, ::1] rigidH, np.uint32_t[:, :, :, ::1] ID): +cpdef void build_cylinder( + float x1, + float y1, + float z1, + float x2, + float y2, + float z2, + float r, + float dx, + float dy, + float dz, + int numID, + int numIDx, + int numIDy, + int numIDz, + bint averaging, + np.uint32_t[:, :, ::1] solid, + np.int8_t[:, :, :, ::1] rigidE, + np.int8_t[:, :, :, ::1] rigidH, + np.uint32_t[:, :, :, ::1] ID + ): """Builds #cylinder commands which sets values in the solid, rigid and ID arrays for a Yee voxel. Args: @@ -625,7 +802,24 @@ cpdef void build_cylinder(float x1, float y1, float z1, float x2, float y2, floa build_voxel(i, j, k, numID, numIDx, numIDy, numIDz, averaging, solid, rigidE, rigidH, ID) -cpdef void build_sphere(int xc, int yc, int zc, float r, float dx, float dy, float dz, int numID, int numIDx, int numIDy, int numIDz, bint averaging, np.uint32_t[:, :, ::1] solid, np.int8_t[:, :, :, ::1] rigidE, np.int8_t[:, :, :, ::1] rigidH, np.uint32_t[:, :, :, ::1] ID): +cpdef void build_sphere( + int xc, + int yc, + int zc, + float r, + float dx, + float dy, + float dz, + int numID, + int numIDx, + int numIDy, + int numIDz, + bint averaging, + np.uint32_t[:, :, ::1] solid, + np.int8_t[:, :, :, ::1] rigidE, + np.int8_t[:, :, :, ::1] rigidH, + np.uint32_t[:, :, :, ::1] ID + ): """Builds #sphere commands which sets values in the solid, rigid and ID arrays for a Yee voxel. Args: @@ -669,7 +863,18 @@ cpdef void build_sphere(int xc, int yc, int zc, float r, float dx, float dy, flo build_voxel(i, j, k, numID, numIDx, numIDy, numIDz, averaging, solid, rigidE, rigidH, ID) -cpdef void build_voxels_from_array(int xs, int ys, int zs, int numexistmaterials, bint averaging, np.int16_t[:, :, ::1] data, np.uint32_t[:, :, ::1] solid, np.int8_t[:, :, :, ::1] rigidE, np.int8_t[:, :, :, ::1] rigidH, np.uint32_t[:, :, :, ::1] ID): +cpdef void build_voxels_from_array( + int xs, + int ys, + int zs, + int numexistmaterials, + bint averaging, + np.int16_t[:, :, ::1] data, + np.uint32_t[:, :, ::1] solid, + np.int8_t[:, :, :, ::1] rigidE, + np.int8_t[:, :, :, ::1] rigidH, + np.uint32_t[:, :, :, ::1] ID + ): """Builds Yee voxels by reading integers from an array. Args: @@ -714,7 +919,20 @@ cpdef void build_voxels_from_array(int xs, int ys, int zs, int numexistmaterials build_voxel(i, j, k, numID, numID, numID, numID, averaging, solid, rigidE, rigidH, ID) -cpdef void build_voxels_from_array_mask(int xs, int ys, int zs, int waternumID, int grassnumID, bint averaging, np.int8_t[:, :, ::1] mask, np.int16_t[:, :, ::1] data, np.uint32_t[:, :, ::1] solid, np.int8_t[:, :, :, ::1] rigidE, np.int8_t[:, :, :, ::1] rigidH, np.uint32_t[:, :, :, ::1] ID): +cpdef void build_voxels_from_array_mask( + int xs, + int ys, + int zs, + int waternumID, + int grassnumID, + bint averaging, + np.int8_t[:, :, ::1] mask, + np.int16_t[:, :, ::1] data, + np.uint32_t[:, :, ::1] solid, + np.int8_t[:, :, :, ::1] rigidE, + np.int8_t[:, :, :, ::1] rigidH, + np.uint32_t[:, :, :, ::1] ID + ): """Builds Yee voxels by reading integers from an array. Args: diff --git a/gprMax/gprMax.py b/gprMax/gprMax.py index 464013bc..7256bfa2 100644 --- a/gprMax/gprMax.py +++ b/gprMax/gprMax.py @@ -67,7 +67,18 @@ def main(): run_main(args) -def api(inputfile, n=1, task=None, restart=None, mpi=False, benchmark=False, geometry_only=False, geometry_fixed=False, write_processed=False, opt_taguchi=False): +def api( + inputfile, + n=1, + task=None, + restart=None, + mpi=False, + benchmark=False, + geometry_only=False, + geometry_fixed=False, + write_processed=False, + opt_taguchi=False + ): """If installed as a module this is the entry point.""" # Print gprMax logo, version, and licencing/copyright information @@ -93,7 +104,8 @@ def api(inputfile, n=1, task=None, restart=None, mpi=False, benchmark=False, geo def run_main(args): - """Top-level function that controls what mode of simulation (standard/optimsation/benchmark etc...) is run. + """ + Top-level function that controls what mode of simulation (standard/optimsation/benchmark etc...) is run. Args: args (dict): Namespace with input arguments from command line or api. @@ -147,13 +159,17 @@ def run_main(args): def run_std_sim(args, inputfile, usernamespace, optparams=None): - """Run standard simulation - models are run one after another and each model is parallelised with OpenMP + """ + Run standard simulation - models are run one after another and each model + is parallelised with OpenMP Args: args (dict): Namespace with command line arguments inputfile (object): File object for the input file. - usernamespace (dict): Namespace that can be accessed by user in any Python code blocks in input file. - optparams (dict): Optional argument. For Taguchi optimisation it provides the parameters to optimise and their values. + usernamespace (dict): Namespace that can be accessed by user in any + Python code blocks in input file. + optparams (dict): Optional argument. For Taguchi optimisation it + provides the parameters to optimise and their values. """ # Set range for number of models to run @@ -171,7 +187,9 @@ def run_std_sim(args, inputfile, usernamespace, optparams=None): tsimstart = perf_counter() for currentmodelrun in range(modelstart, modelend): - if optparams: # If Taguchi optimistaion, add specific value for each parameter to optimise for each experiment to user accessible namespace + # If Taguchi optimistaion, add specific value for each parameter to + # optimise for each experiment to user accessible namespace + if optparams: tmp = {} tmp.update((key, value[currentmodelrun - 1]) for key, value in optparams.items()) modelusernamespace = usernamespace.copy() @@ -185,12 +203,15 @@ def run_std_sim(args, inputfile, usernamespace, optparams=None): def run_benchmark_sim(args, inputfile, usernamespace): - """Run standard simulation in benchmarking mode - models are run one after another and each model is parallelised with OpenMP + """ + Run standard simulation in benchmarking mode - models are run one + after another and each model is parallelised with OpenMP Args: args (dict): Namespace with command line arguments inputfile (object): File object for the input file. - usernamespace (dict): Namespace that can be accessed by user in any Python code blocks in input file. + usernamespace (dict): Namespace that can be accessed by user in any + Python code blocks in input file. """ # Get information about host machine @@ -244,13 +265,17 @@ def run_benchmark_sim(args, inputfile, usernamespace): def run_mpi_sim(args, inputfile, usernamespace, optparams=None): - """Run mixed mode MPI/OpenMP simulation - MPI task farm for models with each model parallelised with OpenMP + """ + Run mixed mode MPI/OpenMP simulation - MPI task farm for models with + each model parallelised with OpenMP Args: args (dict): Namespace with command line arguments inputfile (object): File object for the input file. - usernamespace (dict): Namespace that can be accessed by user in any Python code blocks in input file. - optparams (dict): Optional argument. For Taguchi optimisation it provides the parameters to optimise and their values. + usernamespace (dict): Namespace that can be accessed by user in any + Python code blocks in input file. + optparams (dict): Optional argument. For Taguchi optimisation it + provides the parameters to optimise and their values. """ from mpi4py import MPI @@ -318,7 +343,8 @@ def run_mpi_sim(args, inputfile, usernamespace, optparams=None): gpuinfo = '' print('MPI worker rank {} (PID {}) starting model {}/{}{} on {}'.format(rank, os.getpid(), currentmodelrun, numbermodelruns, gpuinfo, name)) - # If Taguchi optimistaion, add specific value for each parameter to optimise for each experiment to user accessible namespace + # If Taguchi optimistaion, add specific value for each parameter to + # optimise for each experiment to user accessible namespace if 'optparams' in work: tmp = {} tmp.update((key, value[currentmodelrun - 1]) for key, value in work['optparams'].items()) diff --git a/gprMax/grid.py b/gprMax/grid.py index df92b6d3..cc78b97a 100644 --- a/gprMax/grid.py +++ b/gprMax/grid.py @@ -75,7 +75,10 @@ class Grid(object): class FDTDGrid(Grid): - """Holds attributes associated with the entire grid. A convenient way for accessing regularly used parameters.""" + """ + Holds attributes associated with the entire grid. A convenient + way for accessing regularly used parameters. + """ def __init__(self): self.inputfilename = '' @@ -87,7 +90,8 @@ class FDTDGrid(Grid): # CPU - OpenMP threads self.nthreads = 0 - # Threshold (dB) down from maximum power (0dB) of main frequency used to calculate highest frequency for disperion analysis + # Threshold (dB) down from maximum power (0dB) of main frequency used + # to calculate highest frequency for disperion analysis self.highestfreqthres = 60 # Maximum allowable percentage physical phase-velocity phase error self.maxnumericaldisp = 2 @@ -105,7 +109,10 @@ class FDTDGrid(Grid): self.iterations = 0 self.timewindow = 0 - # Ordered dictionary required so that PMLs are always updated in the same order. The order itself does not matter, however, if must be the same from model to model otherwise the numerical precision from adding the PML corrections will be different. + # Ordered dictionary required so that PMLs are always updated in the + # same order. The order itself does not matter, however, if must be the + # same from model to model otherwise the numerical precision from adding + # the PML corrections will be different. self.pmlthickness = OrderedDict((key, 10) for key in PML.boundaryIDs) self.cfs = [] self.pmls = [] @@ -127,8 +134,12 @@ class FDTDGrid(Grid): self.snapshots = [] def initialise_geometry_arrays(self): - """Initialise an array for volumetric material IDs (solid); boolean arrays for specifying whether materials can have dielectric smoothing (rigid); - and an array for cell edge IDs (ID). Solid and ID arrays are initialised to free_space (one); rigid arrays to allow dielectric smoothing (zero). + """ + Initialise an array for volumetric material IDs (solid); + boolean arrays for specifying whether materials can have dielectric smoothing (rigid); + and an array for cell edge IDs (ID). + Solid and ID arrays are initialised to free_space (one); + rigid arrays to allow dielectric smoothing (zero). """ self.solid = np.ones((self.nx, self.ny, self.nz), dtype=np.uint32) self.rigidE = np.zeros((12, self.nx, self.ny, self.nz), dtype=np.int8) @@ -159,7 +170,9 @@ class FDTDGrid(Grid): def dispersion_analysis(G): - """Analysis of numerical dispersion (Taflove et al, 2005, p112) - worse case of maximum frequency and minimum wavelength + """ + Analysis of numerical dispersion (Taflove et al, 2005, p112) - + worse case of maximum frequency and minimum wavelength Args: G (class): Grid class instance - holds essential parameters describing the model. @@ -168,7 +181,8 @@ def dispersion_analysis(G): results (dict): Results from dispersion analysis """ - # Physical phase velocity error (percentage); grid sampling density; material with maximum permittivity; maximum significant frequency + # Physical phase velocity error (percentage); grid sampling density; + # material with maximum permittivity; maximum significant frequency results = {'deltavp': False, 'N': False, 'waveform': True, 'material': False, 'maxfreq': []} # Find maximum significant frequency diff --git a/gprMax/input_cmd_funcs.py b/gprMax/input_cmd_funcs.py index 071decae..28de45f4 100644 --- a/gprMax/input_cmd_funcs.py +++ b/gprMax/input_cmd_funcs.py @@ -19,7 +19,8 @@ import sys from collections import namedtuple -"""This module contains functional forms of some of the most commonly used gprMax commands. It can be useful to use these within Python scripting in an input file. +"""This module contains functional forms of some of the most commonly used gprMax +commands. It can be useful to use these within Python scripting in an input file. For convenience, x, y, z coordinates are lumped in a namedtuple Coordinate: >>> c = Coordinate(0.1, 0.2, 0.3) >>> c @@ -45,11 +46,14 @@ class Coordinate(Coordinate_tuple): def command(cmd, *parameters): - """Helper function. Prints the gprMax #: . None is ignored in the output. + """ + Helper function. Prints the gprMax #: . None is ignored + in the output. Args: cmd (str): the gprMax cmd string to be printed - *parameters: one or more strings as arguments, any None values are ignored + *parameters: one or more strings as arguments, any None values are + ignored Returns: s (str): the printed string @@ -112,7 +116,8 @@ def rotate90_edge(xs, ys, xf, yf, polarisation, rotate90origin): xsnew, ysnew = rotate90_point(xs, ys, rotate90origin) xfnew, yfnew = rotate90_point(xf, yf, rotate90origin) - # Swap coordinates for original y-directed edge, original x-directed edge does not require this. + # Swap coordinates for original y-directed edge, original x-directed + # edge does not require this. if polarisation == 'y': xs = xfnew xf = xsnew @@ -218,10 +223,13 @@ def geometry_view(xs, ys, zs, xf, yf, zf, dx, dy, dz, filename, type='n'): xs, ys, zs, xf, yf, zf (float): Start and finish coordinates. dx, dy, dz (float): Spatial discretisation of geometry view. filename (str): Filename where geometry file information will be stored. - type (str): Can be either n (normal) or f (fine) which specifies whether to output the geometry information on a per-cell basis (n) or a per-cell-edge basis (f). + type (str): Can be either n (normal) or f (fine) which specifies whether + to output the geometry information on a per-cell basis (n) or a + per-cell-edge basis (f). Returns: - s, f, d (tuple): 3 namedtuple Coordinate for the start, finish coordinates and spatial discretisation + s, f, d (tuple): 3 namedtuple Coordinate for the start, + finish coordinates and spatial discretisation """ s = Coordinate(xs, ys, zs) @@ -238,11 +246,14 @@ def snapshot(xs, ys, zs, xf, yf, zf, dx, dy, dz, time, filename): Args: xs, ys, zs, xf, yf, zf (float): Start and finish coordinates. dx, dy, dz (float): Spatial discretisation of geometry view. - time (float): Time in seconds (float) or the iteration number (integer) which denote the point in time at which the snapshot will be taken. + time (float): Time in seconds (float) or the iteration number + (integer) which denote the point in time at which the + snapshot will be taken. filename (str): Filename where geometry file information will be stored. Returns: - s, f, d (tuple): 3 namedtuple Coordinate for the start, finish coordinates and spatial discretisation + s, f, d (tuple): 3 namedtuple Coordinate for the start, + finish coordinates and spatial discretisation """ s = Coordinate(xs, ys, zs) @@ -401,16 +412,24 @@ def cylinder(x1, y1, z1, x2, y2, z2, radius, material, averaging='', rotate90ori return c1, c2 -def cylindrical_sector(axis, ctr1, ctr2, t1, t2, radius, startingangle, sweptangle, material, averaging=''): +def cylindrical_sector(axis, ctr1, ctr2, t1, t2, radius, + startingangle, sweptangle, material, averaging=''): """Prints the gprMax #cylindrical_sector command. Args: - axis (str): Axis of the cylinder from which the sector is defined and can be x, y, or z. + axis (str): Axis of the cylinder from which the sector is defined and + can be x, y, or z. ctr1, ctr2 (float): Coordinates of the centre of the cylindrical sector. - t1, t2 (float): Lower and higher coordinates of the axis of the cylinder from which the sector is defined (in effect they specify the thickness of the sector). + t1, t2 (float): Lower and higher coordinates of the axis of the cylinder + from which the sector is defined (in effect they specify the + thickness of the sector). radius (float): Radius. - startingangle (float): Starting angle (in degrees) for the cylindrical sector (with zero degrees defined on the positive first axis of the plane of the cylindrical sector). - sweptangle (float): Angle (in degrees) swept by the cylindrical sector (the finishing angle of the sector is always anti-clockwise from the starting angle). + startingangle (float): Starting angle (in degrees) for the cylindrical + sector (with zero degrees defined on the positive first axis of + the plane of the cylindrical sector). + sweptangle (float): Angle (in degrees) swept by the cylindrical sector + (the finishing angle of the sector is always anti-clockwise + from the starting angle). material (str): Material identifier(s). averaging (str): Turn averaging on or off. """ @@ -451,7 +470,8 @@ def waveform(shape, amplitude, frequency, identifier): return identifier -def hertzian_dipole(polarisation, f1, f2, f3, identifier, t0=None, t_remove=None, dxdy=None, rotate90origin=()): +def hertzian_dipole(polarisation, f1, f2, f3, identifier, + t0=None, t_remove=None, dxdy=None, rotate90origin=()): """Prints the #hertzian_dipole: polarisation, f1, f2, f3, identifier, [t0, t_remove] Args: @@ -487,7 +507,8 @@ def hertzian_dipole(polarisation, f1, f2, f3, identifier, t0=None, t_remove=None return c -def magnetic_dipole(polarisation, f1, f2, f3, identifier, t0=None, t_remove=None, dxdy=None, rotate90origin=()): +def magnetic_dipole(polarisation, f1, f2, f3, identifier, + t0=None, t_remove=None, dxdy=None, rotate90origin=()): """Prints the #magnetic_dipole: polarisation, f1, f2, f3, identifier, [t0, t_remove] Args: @@ -523,7 +544,8 @@ def magnetic_dipole(polarisation, f1, f2, f3, identifier, t0=None, t_remove=None return c -def voltage_source(polarisation, f1, f2, f3, resistance, identifier, t0=None, t_remove=None, dxdy=None, rotate90origin=()): +def voltage_source(polarisation, f1, f2, f3, resistance, identifier, + t0=None, t_remove=None, dxdy=None, rotate90origin=()): """Prints the #voltage_source: polarisation, f1, f2, f3, resistance, identifier, [t0, t_remove] Args: @@ -560,7 +582,8 @@ def voltage_source(polarisation, f1, f2, f3, resistance, identifier, t0=None, t_ return c -def transmission_line(polarisation, f1, f2, f3, resistance, identifier, t0=None, t_remove=None, dxdy=None, rotate90origin=()): +def transmission_line(polarisation, f1, f2, f3, resistance, identifier, + t0=None, t_remove=None, dxdy=None, rotate90origin=()): """Prints the #transmission_line: polarisation, f1, f2, f3, resistance, identifier, [t0, t_remove] Args: @@ -603,8 +626,10 @@ def rx(x, y, z, identifier=None, to_save=None, polarisation=None, dxdy=None, rot Args: x, y, z (float): are the coordinates (x,y,z) of the receiver in the model. identifier (str): is the optional identifier of the receiver - to_save (list): is a list of outputs with this receiver. It can be any selection from 'Ex', 'Ey', 'Ez', 'Hx', 'Hy', 'Hz', 'Ix', 'Iy', or 'Iz'. - polarisation (str): is the polarisation of the source and can be 'x', 'y', or 'z'. For rotation purposes only. + to_save (list): is a list of outputs with this receiver. It can be + any selection from 'Ex', 'Ey', 'Ez', 'Hx', 'Hy', 'Hz', 'Ix', 'Iy', or 'Iz'. + polarisation (str): is the polarisation of the source + and can be 'x', 'y', or 'z'. For rotation purposes only. dxdy (float): Tuple of x-y spatial resolutions. For rotation purposes only. rotate90origin (tuple): x, y origin for 90 degree CCW rotation in x-y plane. @@ -617,14 +642,16 @@ def rx(x, y, z, identifier=None, to_save=None, polarisation=None, dxdy=None, rot try: xf = x + dxdy[0] except Exception as e: - raise ValueError('With polarization = x, a dxdy[0] float values is required, got dxdy=%s' % dxdy) from e + raise ValueError('With polarization = x, a dxdy[0] float \ + values is required, got dxdy=%s' % dxdy) from e yf = y elif polarisation == 'y': xf = x try: yf = y + dxdy[1] except Exception as e: - raise ValueError('With polarization = y, a dxdy[1] float values is required, got dxdy=%s' % dxdy) from e + raise ValueError('With polarization = y, a dxdy[1] float \ + values is required, got dxdy=%s' % dxdy) from e x, y, xf, yf = rotate90_edge(x, y, xf, yf, polarisation, rotate90origin) @@ -642,7 +669,8 @@ def src_steps(dx=0, dy=0, dz=0): """Prints the #src_steps: dx, dy, dz command. Args: - dx, dy, dz (float): are the increments in (x, y, z) to move all simple sources or all receivers. + dx, dy, dz (float): are the increments in (x, y, z) to + move all simple sources or all receivers. Returns: coordinates (tuple): namedtuple Coordinate of the increments diff --git a/gprMax/input_cmds_file.py b/gprMax/input_cmds_file.py index 68ec3c6a..dfb8b6df 100644 --- a/gprMax/input_cmds_file.py +++ b/gprMax/input_cmds_file.py @@ -24,11 +24,17 @@ from gprMax.exceptions import CmdInputError def process_python_include_code(inputfile, usernamespace): - """Looks for and processes any Python code found in the input file. It will ignore any lines that are comments, i.e. begin with a double hash (##), and any blank lines. It will also ignore any lines that do not begin with a hash (#) after it has processed Python commands. It will also process any include file commands and insert the contents of the included file at that location. + """Looks for and processes any Python code found in the input file. + It will ignore any lines that are comments, i.e. begin with a + double hash (##), and any blank lines. It will also ignore any + lines that do not begin with a hash (#) after it has processed + Python commands. It will also process any include file commands + and insert the contents of the included file at that location. Args: inputfile (object): File object for input file. - usernamespace (dict): Namespace that can be accessed by user in any Python code blocks in input file. + usernamespace (dict): Namespace that can be accessed by user + in any Python code blocks in input file. Returns: processedlines (list): Input commands after Python processing. @@ -36,7 +42,7 @@ def process_python_include_code(inputfile, usernamespace): # Strip out any newline characters and comments that must begin with double hashes inputlines = [line.rstrip() for line in inputfile if(not line.startswith('##') and line.rstrip('\n'))] - + # Rewind input file in preparation for any subsequent reading function inputfile.seek(0) @@ -94,7 +100,7 @@ def process_python_include_code(inputfile, usernamespace): processedlines.append(inputlines[x]) x += 1 - + # Process any include file commands processedlines = process_include_files(processedlines, inputfile) @@ -102,14 +108,17 @@ def process_python_include_code(inputfile, usernamespace): def process_include_files(hashcmds, inputfile): - """Looks for and processes any include file commands and insert the contents of the included file at that location. + """ + Looks for and processes any include file commands and insert + the contents of the included file at that location. Args: hashcmds (list): Input commands. inputfile (object): File object for input file. Returns: - processedincludecmds (list): Input commands after processing any include file commands. + processedincludecmds (list): Input commands after processing + any include file commands. """ processedincludecmds = [] @@ -130,7 +139,7 @@ def process_include_files(hashcmds, inputfile): with open(includefile, 'r') as f: # Strip out any newline characters and comments that must begin with double hashes includelines = [includeline.rstrip() + '\n' for includeline in f if(not includeline.startswith('##') and includeline.rstrip('\n'))] - + # Add lines from include file processedincludecmds.extend(includelines) @@ -143,10 +152,13 @@ def process_include_files(hashcmds, inputfile): def write_processed_file(processedlines, appendmodelnumber, G): - """Writes an input file after any Python code and include commands in the original input file have been processed. + """ + Writes an input file after any Python code and include commands + in the original input file have been processed. Args: - processedlines (list): Input commands after after processing any Python code and include commands. + processedlines (list): Input commands after after processing any + Python code and include commands. appendmodelnumber (str): Text to append to filename. G (class): Grid class instance - holds essential parameters describing the model. """ @@ -161,7 +173,9 @@ def write_processed_file(processedlines, appendmodelnumber, G): def check_cmd_names(processedlines, checkessential=True): - """Checks the validity of commands, i.e. are they gprMax commands, and that all essential commands are present. + """ + Checks the validity of commands, i.e. are they gprMax commands, + and that all essential commands are present. Args: processedlines (list): Input commands after Python processing. @@ -183,12 +197,14 @@ def check_cmd_names(processedlines, checkessential=True): # 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']} - # Geometry object building commands that there can be multiple instances of in a model - these will be lists within the dictionary + # Geometry object building commands that there can be multiple instances + # of in a model - these will be lists within the dictionary geometrycmds = ['#geometry_objects_read', '#edge', '#plate', '#triangle', '#box', '#sphere', '#cylinder', '#cylindrical_sector', '#fractal_box', '#add_surface_roughness', '#add_surface_water', '#add_grass'] # List to store all geometry object commands in order from input file geometry = [] - # Check if command names are valid, if essential commands are present, and add command parameters to appropriate dictionary values or lists + # Check if command names are valid, if essential commands are present, and + # add command parameters to appropriate dictionary values or lists countessentialcmds = 0 lindex = 0 while(lindex < len(processedlines)): @@ -196,7 +212,9 @@ def check_cmd_names(processedlines, checkessential=True): cmdname = cmd[0] cmdparams = cmd[1] - # Check if there is space between command name and parameters, i.e. check first character of parameter string. Ignore case when there are no parameters for a command, e.g. for #taguchi: + # Check if there is space between command name and parameters, i.e. + # check first character of parameter string. Ignore case when there + # are no parameters for a command, e.g. for #taguchi: if ' ' not in cmdparams[0] and len(cmdparams.strip('\n')) != 0: raise CmdInputError('There must be a space between the command name and parameters in ' + processedlines[lindex]) diff --git a/gprMax/input_cmds_geometry.py b/gprMax/input_cmds_geometry.py index 4aef1091..354d05c6 100644 --- a/gprMax/input_cmds_geometry.py +++ b/gprMax/input_cmds_geometry.py @@ -49,13 +49,17 @@ from gprMax.utilities import get_terminal_width def process_geometrycmds(geometry, G): - """This function checks the validity of command parameters, creates instances of classes of parameters, and calls functions to directly set arrays solid, rigid and ID. + """ + This function checks the validity of command parameters, creates instances + of classes of parameters, and calls functions to directly set arrays + solid, rigid and ID. Args: geometry (list): Geometry commands in the model """ - # Disable progress bar if on Windows as it does not update properly when messages are printed + # Disable progress bar if on Windows as it does not update properly + # when messages are printed if sys.platform == 'win32': tqdmdisable = True else: @@ -112,11 +116,13 @@ def process_geometrycmds(geometry, G): data = f['/data'][:] - # Should be int16 to allow for -1 which indicates background, i.e. don't build anything, but AustinMan/Woman maybe uint16 + # Should be int16 to allow for -1 which indicates background, i.e. + # don't build anything, but AustinMan/Woman maybe uint16 if data.dtype != 'int16': data = data.astype('int16') - # Look to see if rigid and ID arrays are present (these should be present if the original geometry objects were written from gprMax) + # Look to see if rigid and ID arrays are present (these should be + # present if the original geometry objects were written from gprMax) try: rigidE = f['/rigidE'][:] rigidH = f['/rigidH'][:] diff --git a/gprMax/input_cmds_multiuse.py b/gprMax/input_cmds_multiuse.py index d4df19b9..9c5f2185 100644 --- a/gprMax/input_cmds_multiuse.py +++ b/gprMax/input_cmds_multiuse.py @@ -43,7 +43,9 @@ from gprMax.waveforms import Waveform def process_multicmds(multicmds, G): - """Checks the validity of command parameters and creates instances of classes of parameters. + """ + Checks the validity of command parameters and creates instances of + classes of parameters. Args: multicmds (dict): Commands that can have multiple instances in the model. diff --git a/gprMax/materials.py b/gprMax/materials.py index 2232419b..dd362b7f 100644 --- a/gprMax/materials.py +++ b/gprMax/materials.py @@ -85,10 +85,12 @@ class Material(object): """Calculates the electric update coefficients of the material. Args: - G (class): Grid class instance - holds essential parameters describing the model. + G (class): Grid class instance - holds essential parameters + describing the model. """ - # The implementation of the dispersive material modelling comes from the derivation in: http://dx.doi.org/10.1109/TAP.2014.2308549 + # The implementation of the dispersive material modelling comes from the + # derivation in: http://dx.doi.org/10.1109/TAP.2014.2308549 if self.maxpoles > 0: self.w = np.zeros(self.maxpoles, dtype=complextype) self.q = np.zeros(self.maxpoles, dtype=complextype) @@ -142,7 +144,9 @@ class Material(object): def process_materials(G): - """Process complete list of materials - calculate update coefficients, store in arrays, and build text list of materials/properties + """ + Process complete list of materials - calculate update coefficients, + store in arrays, and build text list of materials/properties Args: G (class): Grid class instance - holds essential parameters describing the model. @@ -207,7 +211,10 @@ def process_materials(G): class PeplinskiSoil(object): - """Soil objects that are characterised according to a mixing model by Peplinski (http://dx.doi.org/10.1109/36.387598).""" + """ + Soil objects that are characterised according to a mixing + model by Peplinski (http://dx.doi.org/10.1109/36.387598). + """ def __init__(self, ID, sandfraction, clayfraction, bulkdensity, sandpartdensity, watervolfraction): """ @@ -229,7 +236,9 @@ class PeplinskiSoil(object): self.startmaterialnum = 0 def calculate_debye_properties(self, nbins, G): - """Calculates the real and imaginery part of a Debye model for the soil as well as a conductivity. It uses a semi-empirical model (http://dx.doi.org/10.1109/36.387598). + """ + Calculates the real and imaginery part of a Debye model for the soil as + well as a conductivity. It uses a semi-empirical model (http://dx.doi.org/10.1109/36.387598). Args: nbins (int): Number of bins to use to create the different materials. diff --git a/gprMax/model_build_run.py b/gprMax/model_build_run.py index f8bb3069..b3e07413 100644 --- a/gprMax/model_build_run.py +++ b/gprMax/model_build_run.py @@ -69,7 +69,8 @@ def run_model(args, currentmodelrun, modelend, numbermodelruns, inputfile, usern modelend (int): Number of last model to run. numbermodelruns (int): Total number of model runs. inputfile (object): File object for the input file. - usernamespace (dict): Namespace that can be accessed by user in any Python code blocks in input file. + usernamespace (dict): Namespace that can be accessed by user + in any Python code blocks in input file. Returns: tsolve (int): Length of time (seconds) of main FDTD calculations @@ -95,7 +96,8 @@ def run_model(args, currentmodelrun, modelend, numbermodelruns, inputfile, usern inputfilestr = '\n--- Model {}/{}, input file: {}'.format(currentmodelrun, modelend, inputfile.name) print(Fore.GREEN + '{} {}\n'.format(inputfilestr, '-' * (get_terminal_width() - 1 - len(inputfilestr))) + Style.RESET_ALL) - # Add the current model run to namespace that can be accessed by user in any Python code blocks in input file + # Add the current model run to namespace that can be accessed by + # user in any Python code blocks in input file usernamespace['current_model_run'] = currentmodelrun # Read input file and process any Python and include file commands @@ -132,7 +134,8 @@ def run_model(args, currentmodelrun, modelend, numbermodelruns, inputfile, usern print() process_multicmds(multicmds, G) - # Initialise an array for volumetric material IDs (solid), boolean arrays for specifying materials not to be averaged (rigid), + # Initialise an array for volumetric material IDs (solid), boolean + # arrays for specifying materials not to be averaged (rigid), # an array for cell edge IDs (ID) G.initialise_geometry_arrays() @@ -162,7 +165,8 @@ def run_model(args, currentmodelrun, modelend, numbermodelruns, inputfile, usern build_pmls(G, pbar) pbar.close() - # Build the model, i.e. set the material properties (ID) for every edge of every Yee cell + # Build the model, i.e. set the material properties (ID) for every edge + # of every Yee cell print() pbar = tqdm(total=2, desc='Building main grid', ncols=get_terminal_width() - 1, file=sys.stdout, disable=G.tqdmdisable) build_electric_components(G.solid, G.rigidE, G.ID, G) @@ -171,18 +175,21 @@ def run_model(args, currentmodelrun, modelend, numbermodelruns, inputfile, usern pbar.update() pbar.close() - # Process any voltage sources (that have resistance) to create a new material at the source location + # Process any voltage sources (that have resistance) to create a new + # material at the source location for voltagesource in G.voltagesources: voltagesource.create_material(G) # Initialise arrays of update coefficients to pass to update functions G.initialise_std_update_coeff_arrays() - # Initialise arrays of update coefficients and temporary values if there are any dispersive materials + # Initialise arrays of update coefficients and temporary values if + # there are any dispersive materials if Material.maxpoles != 0: G.initialise_dispersive_arrays() - # Process complete list of materials - calculate update coefficients, store in arrays, and build text list of materials/properties + # Process complete list of materials - calculate update coefficients, + # store in arrays, and build text list of materials/properties materialsdata = process_materials(G) if G.messages: print('\nMaterials:') @@ -273,7 +280,8 @@ def run_model(args, currentmodelrun, modelend, numbermodelruns, inputfile, usern print('Memory (RAM) used: ~{}'.format(human_size(p.memory_info().rss))) print('Solving time [HH:MM:SS]: {}'.format(datetime.timedelta(seconds=tsolve))) - # If geometry information to be reused between model runs then FDTDGrid class instance must be global so that it persists + # If geometry information to be reused between model runs then FDTDGrid + # class instance must be global so that it persists if not args.geometry_fixed: del G @@ -281,7 +289,9 @@ def run_model(args, currentmodelrun, modelend, numbermodelruns, inputfile, usern def solve_cpu(currentmodelrun, modelend, G): - """Solving using FDTD method on CPU. Parallelised using Cython (OpenMP) for electric and magnetic field updates, and PML updates. + """ + Solving using FDTD method on CPU. Parallelised using Cython (OpenMP) for + electric and magnetic field updates, and PML updates. Args: currentmodelrun (int): Current model run number. @@ -318,9 +328,12 @@ def solve_cpu(currentmodelrun, modelend, G): source.update_magnetic(iteration, G.updatecoeffsH, G.ID, G.Hx, G.Hy, G.Hz, G) # Update electric field components - if Material.maxpoles == 0: # All materials are non-dispersive so do standard update + # All materials are non-dispersive so do standard update + if Material.maxpoles == 0: update_electric(G.nx, G.ny, G.nz, G.nthreads, G.updatecoeffsE, G.ID, G.Ex, G.Ey, G.Ez, G.Hx, G.Hy, G.Hz) - elif Material.maxpoles == 1: # If there are any dispersive materials do 1st part of dispersive update (it is split into two parts as it requires present and updated electric field values). + # If there are any dispersive materials do 1st part of dispersive update + # (it is split into two parts as it requires present and updated electric field values). + elif Material.maxpoles == 1: update_electric_dispersive_1pole_A(G.nx, G.ny, G.nz, G.nthreads, G.updatecoeffsE, G.updatecoeffsdispersive, G.ID, G.Tx, G.Ty, G.Tz, G.Ex, G.Ey, G.Ez, G.Hx, G.Hy, G.Hz) elif Material.maxpoles > 1: update_electric_dispersive_multipole_A(G.nx, G.ny, G.nz, G.nthreads, Material.maxpoles, G.updatecoeffsE, G.updatecoeffsdispersive, G.ID, G.Tx, G.Ty, G.Tz, G.Ex, G.Ey, G.Ez, G.Hx, G.Hy, G.Hz) @@ -333,7 +346,10 @@ def solve_cpu(currentmodelrun, modelend, G): for source in G.voltagesources + G.transmissionlines + G.hertziandipoles: source.update_electric(iteration, G.updatecoeffsE, G.ID, G.Ex, G.Ey, G.Ez, G) - # If there are any dispersive materials do 2nd part of dispersive update (it is split into two parts as it requires present and updated electric field values). Therefore it can only be completely updated after the electric field has been updated by the PML and source updates. + # If there are any dispersive materials do 2nd part of dispersive update + # (it is split into two parts as it requires present and updated electric + # field values). Therefore it can only be completely updated after the + # electric field has been updated by the PML and source updates. if Material.maxpoles == 1: update_electric_dispersive_1pole_B(G.nx, G.ny, G.nz, G.nthreads, G.updatecoeffsdispersive, G.ID, G.Tx, G.Ty, G.Tz, G.Ex, G.Ey, G.Ez) elif Material.maxpoles > 1: diff --git a/gprMax/optimisation_taguchi.py b/gprMax/optimisation_taguchi.py index 29ecf8e9..885d51f1 100644 --- a/gprMax/optimisation_taguchi.py +++ b/gprMax/optimisation_taguchi.py @@ -42,7 +42,8 @@ def run_opt_sim(args, inputfile, usernamespace): Args: args (dict): Namespace with command line arguments inputfile (object): File object for the input file. - usernamespace (dict): Namespace that can be accessed by user in any Python code blocks in input file. + usernamespace (dict): Namespace that can be accessed by user + in any Python code blocks in input file. """ tsimstart = perf_counter() @@ -52,10 +53,12 @@ def run_opt_sim(args, inputfile, usernamespace): inputfileparts = os.path.splitext(inputfile.name) - # Default maximum number of iterations of optimisation to perform (used if the stopping criterion is not achieved) + # Default maximum number of iterations of optimisation to perform (used + # if the stopping criterion is not achieved) maxiterations = 20 - # Process Taguchi code blocks in the input file; pass in ordered dictionary to hold parameters to optimise + # Process Taguchi code blocks in the input file; pass in ordered + # dictionary to hold parameters to optimise tmp = usernamespace.copy() tmp.update({'optparams': OrderedDict()}) taguchinamespace = taguchi_code_blocks(inputfile, tmp) @@ -110,9 +113,13 @@ def run_opt_sim(args, inputfile, usernamespace): optparams, levels, levelsdiff = calculate_ranges_experiments(optparams, optparamsinit, levels, levelsopt, levelsdiff, OA, N, k, s, iteration) # Run model for each experiment - if args.mpi: # Mixed mode MPI with OpenMP or CUDA - MPI task farm for models with each model parallelised with OpenMP (CPU) or CUDA (GPU) + # Mixed mode MPI with OpenMP or CUDA - MPI task farm for models with + # each model parallelised with OpenMP (CPU) or CUDA (GPU) + if args.mpi: run_mpi_sim(args, inputfile, usernamespace, optparams) - else: # Standard behaviour - models run serially with each model parallelised with OpenMP (CPU) or CUDA (GPU) + # Standard behaviour - models run serially with each model parallelised + # with OpenMP (CPU) or CUDA (GPU) + else: run_std_sim(args, inputfile, usernamespace, optparams) # Calculate fitness value for each experiment @@ -124,7 +131,8 @@ def run_opt_sim(args, inputfile, usernamespace): taguchistr = '\n--- Taguchi optimisation, iteration {}: {} initial experiments with fitness values {}.'.format(iteration + 1, N, fitnessvalues) print('{} {}\n'.format(taguchistr, '-' * (get_terminal_width() - 1 - len(taguchistr)))) - # Calculate optimal levels from fitness values by building a response table; update dictionary of parameters with optimal values + # Calculate optimal levels from fitness values by building a response + # table; update dictionary of parameters with optimal values optparams, levelsopt = calculate_optimal_levels(optparams, levels, levelsopt, fitnessvalues, OA, N, k) # Update dictionary with history of parameters with optimal values @@ -134,9 +142,13 @@ def run_opt_sim(args, inputfile, usernamespace): # Run a confirmation experiment with optimal values args.n = 1 usernamespace['number_model_runs'] = 1 - if args.mpi: # Mixed mode MPI with OpenMP or CUDA - MPI task farm for models with each model parallelised with OpenMP (CPU) or CUDA (GPU) + # Mixed mode MPI with OpenMP or CUDA - MPI task farm for models with + # each model parallelised with OpenMP (CPU) or CUDA (GPU) + if args.mpi: run_mpi_sim(args, inputfile, usernamespace, optparams) - else: # Standard behaviour - models run serially with each model parallelised with OpenMP (CPU) or CUDA (GPU) + # Standard behaviour - models run serially with each model parallelised + # with OpenMP (CPU) or CUDA (GPU) + else: run_std_sim(args, inputfile, usernamespace, optparams) # Calculate fitness value for confirmation experiment @@ -180,11 +192,15 @@ def run_opt_sim(args, inputfile, usernamespace): def taguchi_code_blocks(inputfile, taguchinamespace): - """Looks for and processes a Taguchi code block (containing Python code) in the input file. It will ignore any lines that are comments, i.e. begin with a double hash (##), and any blank lines. + """ + Looks for and processes a Taguchi code block (containing Python code) in + the input file. It will ignore any lines that are comments, i.e. begin + with a double hash (##), and any blank lines. Args: inputfile (object): File object for the input file. - taguchinamespace (dict): Namespace that can be accessed by user a Taguchi code block in input file. + taguchinamespace (dict): Namespace that can be accessed by user a + Taguchi code block in input file. Returns: processedlines (list): Input commands after Python processing. @@ -228,10 +244,13 @@ def taguchi_code_blocks(inputfile, taguchinamespace): def construct_OA(optparams): - """Load an orthogonal array (OA) from a numpy file. Configure and return OA and properties of OA. + """ + Load an orthogonal array (OA) from a numpy file. Configure and + return OA and properties of OA. Args: - optparams (dict): Dictionary containing name of parameters to optimise and their initial ranges + optparams (dict): Dictionary containing name of parameters to + optimise and their initial ranges Returns: OA (array): Orthogonal array @@ -292,7 +311,8 @@ def construct_OA(optparams): # Number of columns cols = int((N - 1) / (s - 1)) - # Algorithm to construct OA from: http://ieeexplore.ieee.org/xpl/articleDetails.jsp?reload=true&arnumber=6812898 + # Algorithm to construct OA from: + # http://ieeexplore.ieee.org/xpl/articleDetails.jsp?reload=true&arnumber=6812898 OA = np.zeros((N + 1, cols + 1), dtype=np.int8) # Construct basic columns @@ -308,7 +328,8 @@ def construct_OA(optparams): for kk in range(1, s): OA[:, col + (jj - 1) * (s - 1) + kk] = np.mod(OA[:, jj] * kk + OA[:, col], s) - # First row and first columns are unneccessary, only there to match algorithm, and cut down columns to number of parameters to optimise + # First row and first columns are unneccessary, only there to + # match algorithm, and cut down columns to number of parameters to optimise OA = OA[1:, 1:k + 1] return OA, N, cols, k, s, t @@ -423,7 +444,8 @@ def calculate_optimal_levels(optparams, levels, levelsopt, fitnessvalues, OA, N, # Calculate optimal level from table of responses optlevel = np.where(responses == np.amax(responses))[0] - # If 2 experiments produce the same fitness value pick first level (this shouldn't happen if the fitness function is designed correctly) + # If 2 experiments produce the same fitness value pick first level + # (this shouldn't happen if the fitness function is designed correctly) if len(optlevel) > 1: optlevel = optlevel[0] diff --git a/gprMax/pml.py b/gprMax/pml.py index 974df3db..69b5f5c5 100644 --- a/gprMax/pml.py +++ b/gprMax/pml.py @@ -169,7 +169,8 @@ class PML(object): self.ny = yf - ys self.nz = zf - zs - # Spatial discretisation and thickness (one extra cell of thickness required for interpolation of electric and magnetic scaling values) + # Spatial discretisation and thickness (one extra cell of thickness + # required for interpolation of electric and magnetic scaling values) if self.direction[0] == 'x': self.d = G.dx self.thickness = self.nx + 1 @@ -266,7 +267,9 @@ class PML(object): def build_pmls(G, pbar): - """This function builds instances of the PML and calculates the initial parameters and coefficients including setting profile + """ + This function builds instances of the PML and calculates the initial + parameters and coefficients including setting profile (based on underlying material er and mr from solid array). Args: diff --git a/gprMax/pml_1order_update.pyx b/gprMax/pml_1order_update.pyx index f81d06d0..a9f59d0f 100644 --- a/gprMax/pml_1order_update.pyx +++ b/gprMax/pml_1order_update.pyx @@ -27,7 +27,30 @@ from gprMax.constants cimport complextype_t ############################################ # 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): +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: @@ -71,7 +94,30 @@ cpdef void update_pml_1order_electric_xminus(int xs, int xf, int ys, int yf, int ########################################### # 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): +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: @@ -115,7 +161,30 @@ cpdef void update_pml_1order_electric_xplus(int xs, int xf, int ys, int yf, int ############################################ # 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): +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: @@ -159,7 +228,30 @@ cpdef void update_pml_1order_electric_yminus(int xs, int xf, int ys, int yf, int ########################################### # 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): +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: @@ -203,7 +295,30 @@ cpdef void update_pml_1order_electric_yplus(int xs, int xf, int ys, int yf, int ############################################ # 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): +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: @@ -247,7 +362,30 @@ cpdef void update_pml_1order_electric_zminus(int xs, int xf, int ys, int yf, int ########################################### # 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): +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: @@ -291,7 +429,30 @@ cpdef void update_pml_1order_electric_zplus(int xs, int xf, int ys, int yf, int ############################################ # 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): +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: @@ -335,7 +496,30 @@ cpdef void update_pml_1order_magnetic_xminus(int xs, int xf, int ys, int yf, int ########################################### # 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): +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: @@ -379,7 +563,30 @@ cpdef void update_pml_1order_magnetic_xplus(int xs, int xf, int ys, int yf, int ############################################ # 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): +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: @@ -423,7 +630,30 @@ cpdef void update_pml_1order_magnetic_yminus(int xs, int xf, int ys, int yf, int ########################################### # 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): +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: @@ -467,7 +697,30 @@ cpdef void update_pml_1order_magnetic_yplus(int xs, int xf, int ys, int yf, int ############################################ # 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): +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: @@ -511,7 +764,30 @@ cpdef void update_pml_1order_magnetic_zminus(int xs, int xf, int ys, int yf, int ########################################### # 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): +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: diff --git a/gprMax/pml_2order_update.pyx b/gprMax/pml_2order_update.pyx index d6da0a75..86ad4f3e 100644 --- a/gprMax/pml_2order_update.pyx +++ b/gprMax/pml_2order_update.pyx @@ -27,7 +27,30 @@ from gprMax.constants cimport complextype_t ############################################ # 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): +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: @@ -78,7 +101,30 @@ cpdef void update_pml_2order_electric_xminus(int xs, int xf, int ys, int yf, int ########################################### # 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): +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: @@ -129,7 +175,30 @@ cpdef void update_pml_2order_electric_xplus(int xs, int xf, int ys, int yf, int ############################################ # 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): +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: @@ -180,7 +249,30 @@ cpdef void update_pml_2order_electric_yminus(int xs, int xf, int ys, int yf, int ########################################### # 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): +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: @@ -231,7 +323,30 @@ cpdef void update_pml_2order_electric_yplus(int xs, int xf, int ys, int yf, int ############################################ # 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): +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: @@ -282,7 +397,30 @@ cpdef void update_pml_2order_electric_zminus(int xs, int xf, int ys, int yf, int ########################################### # 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): +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: @@ -333,7 +471,30 @@ cpdef void update_pml_2order_electric_zplus(int xs, int xf, int ys, int yf, int ############################################ # 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): +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: @@ -384,7 +545,30 @@ cpdef void update_pml_2order_magnetic_xminus(int xs, int xf, int ys, int yf, int ########################################### # 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): +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: @@ -435,7 +619,30 @@ cpdef void update_pml_2order_magnetic_xplus(int xs, int xf, int ys, int yf, int ############################################ # 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): +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: @@ -486,7 +693,30 @@ cpdef void update_pml_2order_magnetic_yminus(int xs, int xf, int ys, int yf, int ########################################### # 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): +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: @@ -537,7 +767,30 @@ cpdef void update_pml_2order_magnetic_yplus(int xs, int xf, int ys, int yf, int ############################################ # 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): +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: @@ -588,7 +841,30 @@ cpdef void update_pml_2order_magnetic_zminus(int xs, int xf, int ys, int yf, int ########################################### # 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): +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: diff --git a/gprMax/snapshots.py b/gprMax/snapshots.py index 73b976b0..be2180a8 100644 --- a/gprMax/snapshots.py +++ b/gprMax/snapshots.py @@ -136,7 +136,8 @@ class Snapshot(object): for j in range(self.ys, self.yf, self.dy): for i in range(self.xs, self.xf, self.dx): pbar.update(n=12) - # The magnetic field component value at a point comes from average of 2 magnetic field component values in that cell and the following cell + # The magnetic field component value at a point comes from average + # of 2 magnetic field component values in that cell and the following cell self.filehandle.write(pack(Snapshot.floatstring, (Hx[i, j, k] + Hx[i + 1, j, k]) / 2)) self.filehandle.write(pack(Snapshot.floatstring, (Hy[i, j, k] + Hy[i, j + 1, k]) / 2)) self.filehandle.write(pack(Snapshot.floatstring, (Hz[i, j, k] + Hz[i, j, k + 1]) / 2)) diff --git a/gprMax/sources.py b/gprMax/sources.py index b9709ed6..b0b461ef 100644 --- a/gprMax/sources.py +++ b/gprMax/sources.py @@ -69,7 +69,10 @@ class Source(object): class VoltageSource(Source): - """A voltage source can be a hard source if it's resistance is zero, i.e. the time variation of the specified electric field component is prescribed. If it's resistance is non-zero it behaves as a resistive voltage source.""" + """ + A voltage source can be a hard source if it's resistance is zero, i.e. the + time variation of the specified electric field component is prescribed. + If it's resistance is non-zero it behaves as a resistive voltage source.""" def __init__(self): super().__init__() @@ -111,7 +114,9 @@ class VoltageSource(Source): Ez[i, j, k] = -1 * self.waveformvaluesJ[iteration] / G.dz def create_material(self, G): - """Create a new material at the voltage source location that adds the voltage source conductivity to the underlying parameters. + """ + Create a new material at the voltage source location that adds the + voltage source conductivity to the underlying parameters. Args: G (class): Grid class instance - holds essential parameters describing the model. @@ -211,7 +216,10 @@ class MagneticDipole(Source): class TransmissionLine(Source): - """A transmission line source is a one-dimensional transmission line which is attached virtually to a grid cell.""" + """ + A transmission line source is a one-dimensional transmission + line which is attached virtually to a grid cell. + """ def __init__(self, G): """ @@ -226,10 +234,12 @@ class TransmissionLine(Source): self.abcv0 = 0 self.abcv1 = 0 - # Spatial step of transmission line (N.B if the magic time step is used it results in instabilities for certain impedances) + # Spatial step of transmission line (N.B if the magic time step is + # used it results in instabilities for certain impedances) self.dl = np.sqrt(3) * c * G.dt - # Number of cells in the transmission line (initially a long line to calculate incident voltage and current); consider putting ABCs/PML at end + # Number of cells in the transmission line (initially a long line to + # calculate incident voltage and current); consider putting ABCs/PML at end self.nl = round_value(0.667 * G.iterations) # Cell position of the one-way injector excitation in the transmission line @@ -246,7 +256,9 @@ class TransmissionLine(Source): self.Itotal = np.zeros(G.iterations, dtype=floattype) def calculate_incident_V_I(self, G): - """Calculates the incident voltage and current with a long length transmission line not connected to the main grid from: http://dx.doi.org/10.1002/mop.10415 + """ + Calculates the incident voltage and current with a long length + transmission line not connected to the main grid from: http://dx.doi.org/10.1002/mop.10415 Args: G (class): Grid class instance - holds essential parameters describing the model. diff --git a/gprMax/utilities.py b/gprMax/utilities.py index 28c85fae..1b41c5ca 100644 --- a/gprMax/utilities.py +++ b/gprMax/utilities.py @@ -83,7 +83,9 @@ def logo(version): @contextmanager def open_path_file(path_or_file): - """Accepts either a path as a string or a file object and returns a file object (http://stackoverflow.com/a/6783680). + """ + Accepts either a path as a string or a file object and returns a file + object (http://stackoverflow.com/a/6783680). Args: path_or_file: path as a string or a file object. @@ -157,7 +159,8 @@ def get_host_info(): """Get information about the machine, CPU, RAM, and OS. Returns: - hostinfo (dict): Manufacturer and model of machine; description of CPU type, speed, cores; RAM; name and version of operating system. + hostinfo (dict): Manufacturer and model of machine; description of CPU + type, speed, cores; RAM; name and version of operating system. """ # Default to 'unknown' if any of the detection fails @@ -275,7 +278,8 @@ def get_host_info(): hostinfo['hyperthreading'] = hyperthreading hostinfo['logicalcores'] = psutil.cpu_count() try: - hostinfo['physicalcores'] = psutil.cpu_count(logical=False) # Get number of physical CPU cores, i.e. avoid hyperthreading with OpenMP + # Get number of physical CPU cores, i.e. avoid hyperthreading with OpenMP + hostinfo['physicalcores'] = psutil.cpu_count(logical=False) except ValueError: hostinfo['physicalcores'] = hostinfo['logicalcores'] hostinfo['ram'] = psutil.virtual_memory().total