From 3698a0a6f0a67dcdb2089536c8692eb2c45bb705 Mon Sep 17 00:00:00 2001 From: Craig Warren Date: Wed, 30 Aug 2023 09:25:45 +0100 Subject: [PATCH 01/35] Removed 0.5mm & 1mm resolutions for 400MHz model as not valid anymore. --- toolboxes/GPRAntennaModels/GSSI.py | 617 +++++++---------------------- 1 file changed, 152 insertions(+), 465 deletions(-) diff --git a/toolboxes/GPRAntennaModels/GSSI.py b/toolboxes/GPRAntennaModels/GSSI.py index 641cc439..6273d61a 100644 --- a/toolboxes/GPRAntennaModels/GSSI.py +++ b/toolboxes/GPRAntennaModels/GSSI.py @@ -420,12 +420,12 @@ def antenna_like_GSSI_1500(x, y, z, resolution=0.001, **kwargs): return scene_objects -def antenna_like_GSSI_400(x, y, z, resolution=0.001, **kwargs): +def antenna_like_GSSI_400(x, y, z, resolution=0.002, **kwargs): """Inserts a description of an antenna similar to the GSSI 400MHz antenna. This model represents an update to the previous model of the GSSI 400MHz antenna and was created and optimised by Stadler et al. (2022) in: https://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=9686638. - Can be used with 0.5mm, 1mm (default) or 2mm spatial resolution. + Can be used with 2mm spatial resolution. The external dimensions of the antenna are 300x300x178mm. One output point is defined between the arms of the receiver bowtie. The bowties are aligned with the y axis so the output is the y component @@ -451,7 +451,6 @@ def antenna_like_GSSI_400(x, y, z, resolution=0.001, **kwargs): casesize = (0.3, 0.3, 0.178) # original casethickness = 0.002 shieldthickness = 0.002 - foamsurroundthickness = 0.003 pcbthickness = 0.002 bowtiebase = 0.06 bowtieheight = 0.06 # original 0.056 @@ -488,25 +487,12 @@ def antenna_like_GSSI_400(x, y, z, resolution=0.001, **kwargs): # Coordinates of source excitation point in antenna tx = x + 0.01 + 0.005 + 0.056, y + casethickness + 0.005 + 0.143, z + skidthickness - if resolution == 0.0005: - dx = 0.0005 - dy = 0.0005 - dz = 0.0005 - tx = x + 0.01 + 0.005 + 0.056, y + casethickness + 0.005 + 0.1435, z + skidthickness - elif resolution == 0.001: - dx = 0.001 - dy = 0.001 - dz = 0.001 - elif resolution == 0.002: - dx = 0.002 - dy = 0.002 - dz = 0.002 - foamsurroundthickness = 0.002 - metalboxheight = 0.088 - tx = x + 0.01 + 0.004 + 0.056, y + casethickness + 0.005 + 0.143 - 0.002, z + skidthickness - 0.002 - else: - logger.exception("This antenna module can only be used with a spatial discretisation of 0.5mm, 1mm, 2mm") - raise ValueError + dx = 0.002 + dy = 0.002 + dz = 0.002 + foamsurroundthickness = 0.002 + metalboxheight = 0.088 + tx = x + 0.01 + 0.004 + 0.056, y + casethickness + 0.005 + 0.143 - 0.002, z + skidthickness - 0.002 # Material definitions absorber = gprMax.Material(er=absorberEr, se=absorbersig, mr=1, sm=0, id="absorber") @@ -565,68 +551,25 @@ def antenna_like_GSSI_400(x, y, z, resolution=0.001, **kwargs): scene_objects.extend((b1, b2, b3, b4, b5)) # PCB - if resolution == 0.0005: - b6 = gprMax.Box( - p1=(x + 0.01 + 0.005 + 0.018, y + casethickness + 0.005 + 0.0235, z + skidthickness), - p2=( - x + 0.01 + 0.005 + 0.034 + bowtiebase, - y + casethickness + 0.005 + 0.204 + patchheight, - z + skidthickness + pcbthickness, - ), - material_id="pcb", - ) - b7 = gprMax.Box( - p1=(x + 0.01 + 0.005 + 0.178, y + casethickness + 0.005 + 0.0235, z + skidthickness), - p2=( - x + 0.01 + 0.005 + 0.194 + bowtiebase, - y + casethickness + 0.005 + 0.204 + patchheight, - z + skidthickness + pcbthickness, - ), - material_id="pcb", - ) - scene_objects.extend((b6, b7)) - - elif resolution == 0.001: - b6 = gprMax.Box( - p1=(x + 0.01 + 0.005 + 0.018, y + casethickness + 0.005 + 0.023, z + skidthickness), - p2=( - x + 0.01 + 0.005 + 0.034 + bowtiebase, - y + casethickness + 0.005 + 0.204 + patchheight, - z + skidthickness + pcbthickness, - ), - material_id="pcb", - ) - b7 = gprMax.Box( - p1=(x + 0.01 + 0.005 + 0.178, y + casethickness + 0.005 + 0.023, z + skidthickness), - p2=( - x + 0.01 + 0.005 + 0.194 + bowtiebase, - y + casethickness + 0.005 + 0.204 + patchheight, - z + skidthickness + pcbthickness, - ), - material_id="pcb", - ) - scene_objects.extend((b6, b7)) - - elif resolution == 0.002: - b6 = gprMax.Box( - p1=(x + 0.01 + 0.005 + 0.017, y + casethickness + 0.005 + 0.021, z + skidthickness - 0.002), - p2=( - x + 0.01 + 0.005 + 0.033 + bowtiebase, - y + casethickness + 0.006 + 0.202 + patchheight, - z + skidthickness + pcbthickness - 0.002, - ), - material_id="pcb", - ) - b7 = gprMax.Box( - p1=(x + 0.01 + 0.005 + 0.179, y + casethickness + 0.005 + 0.021, z + skidthickness - 0.002), - p2=( - x + 0.01 + 0.005 + 0.195 + bowtiebase, - y + casethickness + 0.006 + 0.202 + patchheight, - z + skidthickness + pcbthickness - 0.002, - ), - material_id="pcb", - ) - scene_objects.extend((b6, b7)) + b6 = gprMax.Box( + p1=(x + 0.01 + 0.005 + 0.017, y + casethickness + 0.005 + 0.021, z + skidthickness - 0.002), + p2=( + x + 0.01 + 0.005 + 0.033 + bowtiebase, + y + casethickness + 0.006 + 0.202 + patchheight, + z + skidthickness + pcbthickness - 0.002, + ), + material_id="pcb", + ) + b7 = gprMax.Box( + p1=(x + 0.01 + 0.005 + 0.179, y + casethickness + 0.005 + 0.021, z + skidthickness - 0.002), + p2=( + x + 0.01 + 0.005 + 0.195 + bowtiebase, + y + casethickness + 0.006 + 0.202 + patchheight, + z + skidthickness + pcbthickness - 0.002, + ), + material_id="pcb", + ) + scene_objects.extend((b6, b7)) elif smooth_dec == "no": # Plastic case @@ -682,366 +625,124 @@ def antenna_like_GSSI_400(x, y, z, resolution=0.001, **kwargs): scene_objects.extend((b8, b9, b10, b11, b12)) # PCB - if resolution == 0.0005: - b13 = gprMax.Box( - p1=(x + 0.01 + 0.005 + 0.018, y + casethickness + 0.005 + 0.0235, z + skidthickness), - p2=( - x + 0.01 + 0.005 + 0.034 + bowtiebase, - y + casethickness + 0.005 + 0.204 + patchheight, - z + skidthickness + pcbthickness, - ), - material_id="pcb", - averaging="n", - ) - b14 = gprMax.Box( - p1=(x + 0.01 + 0.005 + 0.178, y + casethickness + 0.005 + 0.0235, z + skidthickness), - p2=( - x + 0.01 + 0.005 + 0.194 + bowtiebase, - y + casethickness + 0.005 + 0.204 + patchheight, - z + skidthickness + pcbthickness, - ), - material_id="pcb", - averaging="n", - ) - scene_objects.extend((b13, b14)) - - elif resolution == 0.001: - b13 = gprMax.Box( - p1=(x + 0.01 + 0.005 + 0.018, y + casethickness + 0.005 + 0.023, z + skidthickness), - p2=( - x + 0.01 + 0.005 + 0.034 + bowtiebase, - y + casethickness + 0.005 + 0.204 + patchheight, - z + skidthickness + pcbthickness, - ), - material_id="pcb", - averaging="n", - ) - b14 = gprMax.Box( - p1=(x + 0.01 + 0.005 + 0.178, y + casethickness + 0.005 + 0.023, z + skidthickness), - p2=( - x + 0.01 + 0.005 + 0.194 + bowtiebase, - y + casethickness + 0.005 + 0.204 + patchheight, - z + skidthickness + pcbthickness, - ), - material_id="pcb", - averaging="n", - ) - scene_objects.extend((b13, b14)) - - elif resolution == 0.002: - b13 = gprMax.Box( - p1=(x + 0.01 + 0.005 + 0.017, y + casethickness + 0.005 + 0.021, z + skidthickness - 0.002), - p2=( - x + 0.01 + 0.005 + 0.033 + bowtiebase, - y + casethickness + 0.006 + 0.202 + patchheight, - z + skidthickness + pcbthickness - 0.002, - ), - material_id="pcb", - averaging="n", - ) - b14 = gprMax.Box( - p1=(x + 0.01 + 0.005 + 0.179, y + casethickness + 0.005 + 0.021, z + skidthickness), - p2=( - x + 0.01 + 0.005 + 0.195 + bowtiebase, - y + casethickness + 0.006 + 0.202 + patchheight, - z + skidthickness + pcbthickness, - ), - material_id="pcb", - averaging="n", - ) - scene_objects.extend((b13, b14)) + b13 = gprMax.Box( + p1=(x + 0.01 + 0.005 + 0.017, y + casethickness + 0.005 + 0.021, z + skidthickness - 0.002), + p2=( + x + 0.01 + 0.005 + 0.033 + bowtiebase, + y + casethickness + 0.006 + 0.202 + patchheight, + z + skidthickness + pcbthickness - 0.002, + ), + material_id="pcb", + averaging="n", + ) + b14 = gprMax.Box( + p1=(x + 0.01 + 0.005 + 0.179, y + casethickness + 0.005 + 0.021, z + skidthickness), + p2=( + x + 0.01 + 0.005 + 0.195 + bowtiebase, + y + casethickness + 0.006 + 0.202 + patchheight, + z + skidthickness + pcbthickness, + ), + material_id="pcb", + averaging="n", + ) + scene_objects.extend((b13, b14)) # PCB components # My own bowties with triangle commands - if resolution == 0.0005: - # "left" side - # extension plates - p1 = gprMax.Plate( - p1=(x + 0.01 + 0.005 + 0.026, y + casethickness + 0.005 + 0.0235, z + skidthickness), - p2=( - x + 0.01 + 0.005 + 0.026 + bowtiebase, - y + casethickness + 0.005 + 0.0235 + patchheight, - z + skidthickness, - ), - material_id="pec", - ) - p2 = gprMax.Plate( - p1=(x + 0.01 + 0.005 + 0.026, y + casethickness + 0.005 + 0.204, z + skidthickness), - p2=( - x + 0.01 + 0.005 + 0.026 + bowtiebase, - y + casethickness + 0.005 + 0.204 + patchheight, - z + skidthickness, - ), - material_id="pec", - ) - # triangles - t1 = gprMax.Triangle( - p1=(x + 0.01 + 0.005 + 0.026, y + casethickness + 0.005 + 0.0835, z + skidthickness), - p2=(x + 0.01 + 0.005 + 0.026 + bowtiebase, y + casethickness + 0.005 + 0.0835, z + skidthickness), - p3=( - x + 0.01 + 0.005 + 0.026 + (bowtiebase / 2), - y + casethickness + 0.005 + 0.0835 + bowtieheight, - z + skidthickness, - ), - thickness=0, - material_id="pec", - ) - t2 = gprMax.Triangle( - p1=(x + 0.01 + 0.005 + 0.026, y + casethickness + 0.005 + 0.204, z + skidthickness), - p2=(x + 0.01 + 0.005 + 0.026 + bowtiebase, y + casethickness + 0.005 + 0.204, z + skidthickness), - p3=( - x + 0.01 + 0.005 + 0.026 + (bowtiebase / 2), - y + casethickness + 0.005 + 0.204 - bowtieheight, - z + skidthickness, - ), - thickness=0, - material_id="pec", - ) + # "left" side + # extension plates + p1 = gprMax.Plate( + p1=(x + 0.01 + 0.005 + 0.025, y + casethickness + 0.005 + 0.021, z + skidthickness - 0.002), + p2=( + x + 0.01 + 0.005 + 0.025 + bowtiebase, + y + casethickness + 0.005 + 0.021 + patchheight, + z + skidthickness - 0.002, + ), + material_id="pec", + ) + p2 = gprMax.Plate( + p1=(x + 0.01 + 0.005 + 0.025, y + casethickness + 0.005 + 0.203, z + skidthickness - 0.002), + p2=( + x + 0.01 + 0.005 + 0.025 + bowtiebase, + y + casethickness + 0.005 + 0.203 + patchheight, + z + skidthickness - 0.002, + ), + material_id="pec", + ) + # triangles + t1 = gprMax.Triangle( + p1=(x + 0.01 + 0.005 + 0.025, y + casethickness + 0.005 + 0.081, z + skidthickness - 0.002), + p2=(x + 0.01 + 0.005 + 0.025 + bowtiebase, y + casethickness + 0.005 + 0.081, z + skidthickness - 0.002), + p3=( + x + 0.01 + 0.005 + 0.025 + (bowtiebase / 2), + y + casethickness + 0.005 + 0.081 + bowtieheight, + z + skidthickness - 0.002, + ), + thickness=0, + material_id="pec", + ) + t2 = gprMax.Triangle( + p1=(x + 0.01 + 0.005 + 0.025, y + casethickness + 0.005 + 0.203, z + skidthickness - 0.002), + p2=(x + 0.01 + 0.005 + 0.025 + bowtiebase, y + casethickness + 0.005 + 0.203, z + skidthickness - 0.002), + p3=( + x + 0.01 + 0.005 + 0.025 + (bowtiebase / 2), + y + casethickness + 0.005 + 0.203 - bowtieheight, + z + skidthickness - 0.002, + ), + thickness=0, + material_id="pec", + ) + # "right" side + p3 = gprMax.Plate( + p1=(x + 0.01 + 0.005 + 0.187, y + casethickness + 0.005 + 0.021, z + skidthickness - 0.002), + p2=( + x + 0.01 + 0.005 + 0.187 + bowtiebase, + y + casethickness + 0.005 + 0.021 + patchheight, + z + skidthickness - 0.002, + ), + material_id="pec", + ) + p4 = gprMax.Plate( + p1=(x + 0.01 + 0.005 + 0.187, y + casethickness + 0.005 + 0.203, z + skidthickness - 0.002), + p2=( + x + 0.01 + 0.005 + 0.187 + bowtiebase, + y + casethickness + 0.005 + 0.203 + patchheight, + z + skidthickness - 0.002, + ), + material_id="pec", + ) + # triangles + t3 = gprMax.Triangle( + p1=(x + 0.01 + 0.005 + 0.187, y + casethickness + 0.005 + 0.081, z + skidthickness - 0.002), + p2=(x + 0.01 + 0.005 + 0.187 + bowtiebase, y + casethickness + 0.005 + 0.081, z + skidthickness - 0.002), + p3=( + x + 0.01 + 0.005 + 0.187 + (bowtiebase / 2), + y + casethickness + 0.005 + 0.081 + bowtieheight, + z + skidthickness - 0.002, + ), + thickness=0, + material_id="pec", + ) + t4 = gprMax.Triangle( + p1=(x + 0.01 + 0.005 + 0.187, y + casethickness + 0.005 + 0.203, z + skidthickness - 0.002), + p2=(x + 0.01 + 0.005 + 0.187 + bowtiebase, y + casethickness + 0.005 + 0.203, z + skidthickness - 0.002), + p3=( + x + 0.01 + 0.005 + 0.187 + (bowtiebase / 2), + y + casethickness + 0.005 + 0.203 - bowtieheight, + z + skidthickness - 0.002, + ), + thickness=0, + material_id="pec", + ) - # "right" side - p3 = gprMax.Plate( - p1=(x + 0.01 + 0.005 + 0.186, y + casethickness + 0.005 + 0.0235, z + skidthickness), - p2=( - x + 0.01 + 0.005 + 0.186 + bowtiebase, - y + casethickness + 0.005 + 0.0235 + patchheight, - z + skidthickness, - ), - material_id="pec", - ) - p4 = gprMax.Plate( - p1=(x + 0.01 + 0.005 + 0.186, y + casethickness + 0.005 + 0.204, z + skidthickness), - p2=( - x + 0.01 + 0.005 + 0.186 + bowtiebase, - y + casethickness + 0.005 + 0.204 + patchheight, - z + skidthickness, - ), - material_id="pec", - ) - # triangles - t3 = gprMax.Triangle( - p1=(x + 0.01 + 0.005 + 0.186, y + casethickness + 0.005 + 0.0835, z + skidthickness), - p2=(x + 0.01 + 0.005 + 0.186 + bowtiebase, y + casethickness + 0.005 + 0.0835, z + skidthickness), - p3=( - x + 0.01 + 0.005 + 0.186 + (bowtiebase / 2), - y + casethickness + 0.005 + 0.0835 + bowtieheight, - z + skidthickness, - ), - thickness=0, - material_ID="pec", - ) - t4 = gprMax.Triangle( - p1=(x + 0.01 + 0.005 + 0.186, y + casethickness + 0.005 + 0.204, z + skidthickness), - p2=(x + 0.01 + 0.005 + 0.186 + bowtiebase, y + casethickness + 0.005 + 0.204, z + skidthickness), - p3=( - x + 0.01 + 0.005 + 0.186 + (bowtiebase / 2), - y + casethickness + 0.005 + 0.204 - bowtieheight, - z + skidthickness, - ), - thickness=0, - material_id="pec", - ) - - # Edges that represent wire between bowtie halves in 1mm model - e1 = gprMax.Edge(p1=(tx[0] + 0.16, tx[1] - dy, tx[2]), p2=(tx[0] + 0.16, tx[1], tx[2]), material_id="pec") - e2 = gprMax.Edge( - p1=(tx[0] + 0.16, tx[1] + dy, tx[2]), p2=(tx[0] + 0.16, tx[1] + 2 * dy, tx[2]), material_id="pec" - ) - e3 = gprMax.Edge(p1=(tx[0], tx[1] - dy, tx[2]), p2=(tx[0], tx[1], tx[2]), material_id="pec") - e4 = gprMax.Edge(p1=(tx[0], tx[1] + dy, tx[2]), p2=(tx[0], tx[1] + 2 * dy, tx[2]), material_id="pec") - scene_objects.extend((p1, p2, t1, t2, p3, p4, t3, t4, e1, e2, e3, e4)) - - elif resolution == 0.001: - # "left" side - # extension plates - p1 = gprMax.Plate( - p1=(x + 0.01 + 0.005 + 0.026, y + casethickness + 0.005 + 0.023, z + skidthickness), - p2=( - x + 0.01 + 0.005 + 0.026 + bowtiebase, - y + casethickness + 0.005 + 0.023 + patchheight, - z + skidthickness, - ), - material_id="pec", - ) - p2 = gprMax.Plate( - p1=(x + 0.01 + 0.005 + 0.026, y + casethickness + 0.005 + 0.204, z + skidthickness), - p2=( - x + 0.01 + 0.005 + 0.026 + bowtiebase, - y + casethickness + 0.005 + 0.204 + patchheight, - z + skidthickness, - ), - material_id="pec", - ) - # triangles - t1 = gprMax.Triangle( - p1=(x + 0.01 + 0.005 + 0.026, y + casethickness + 0.005 + 0.083, z + skidthickness), - p2=(x + 0.01 + 0.005 + 0.026 + bowtiebase, y + casethickness + 0.005 + 0.083, z + skidthickness), - p3=( - x + 0.01 + 0.005 + 0.026 + (bowtiebase / 2), - y + casethickness + 0.005 + 0.083 + bowtieheight, - z + skidthickness, - ), - thickness=0, - material_id="pec", - ) - t2 = gprMax.Triangle( - p1=(x + 0.01 + 0.005 + 0.026, y + casethickness + 0.005 + 0.204, z + skidthickness), - p2=(x + 0.01 + 0.005 + 0.026 + bowtiebase, y + casethickness + 0.005 + 0.204, z + skidthickness), - p3=( - x + 0.01 + 0.005 + 0.026 + (bowtiebase / 2), - y + casethickness + 0.005 + 0.204 - bowtieheight, - z + skidthickness, - ), - thickness=0, - material_id="pec", - ) - - # "right" side - p3 = gprMax.Plate( - p1=(x + 0.01 + 0.005 + 0.186, y + casethickness + 0.005 + 0.023, z + skidthickness), - p2=( - x + 0.01 + 0.005 + 0.186 + bowtiebase, - y + casethickness + 0.005 + 0.023 + patchheight, - z + skidthickness, - ), - material_id="pec", - ) - p4 = gprMax.Plate( - p1=(x + 0.01 + 0.005 + 0.186, y + casethickness + 0.005 + 0.204, z + skidthickness), - p2=( - x + 0.01 + 0.005 + 0.186 + bowtiebase, - y + casethickness + 0.005 + 0.204 + patchheight, - z + skidthickness, - ), - material_id="pec", - ) - # triangles - t3 = gprMax.Triangle( - p1=(x + 0.01 + 0.005 + 0.186, y + casethickness + 0.005 + 0.083, z + skidthickness), - p2=(x + 0.01 + 0.005 + 0.186 + bowtiebase, y + casethickness + 0.005 + 0.083, z + skidthickness), - p3=( - x + 0.01 + 0.005 + 0.186 + (bowtiebase / 2), - y + casethickness + 0.005 + 0.083 + bowtieheight, - z + skidthickness, - ), - thickness=0, - material_id="pec", - ) - t4 = gprMax.Triangle( - p1=(x + 0.01 + 0.005 + 0.186, y + casethickness + 0.005 + 0.204, z + skidthickness), - p2=(x + 0.01 + 0.005 + 0.186 + bowtiebase, y + casethickness + 0.005 + 0.204, z + skidthickness), - p3=( - x + 0.01 + 0.005 + 0.186 + (bowtiebase / 2), - y + casethickness + 0.005 + 0.204 - bowtieheight, - z + skidthickness, - ), - thickness=0, - material_id="pec", - ) - - # Edges that represent wire between bowtie halves in 1mm model - e1 = gprMax.Edge(p1=(tx[0] + 0.16, tx[1] - dy, tx[2]), p2=(tx[0] + 0.16, tx[1], tx[2]), material_id="pec") - e2 = gprMax.Edge( - p1=(tx[0] + 0.16, tx[1] + dy, tx[2]), p2=(tx[0] + 0.16, tx[1] + 2 * dy, tx[2]), material_id="pec" - ) - e3 = gprMax.Edge(p1=(tx[0], tx[1] - dy, tx[2]), p2=(tx[0], tx[1], tx[2]), material_id="pec") - e4 = gprMax.Edge(p1=(tx[0], tx[1] + dy, tx[2]), p2=(tx[0], tx[1] + 2 * dy, tx[2]), material_id="pec") - scene_objects.extend((p1, p2, t1, t2, p3, p4, t3, t4, e1, e2, e3, e4)) - - elif resolution == 0.002: - # "left" side - # extension plates - p1 = gprMax.Plate( - p1=(x + 0.01 + 0.005 + 0.025, y + casethickness + 0.005 + 0.021, z + skidthickness - 0.002), - p2=( - x + 0.01 + 0.005 + 0.025 + bowtiebase, - y + casethickness + 0.005 + 0.021 + patchheight, - z + skidthickness - 0.002, - ), - material_id="pec", - ) - p2 = gprMax.Plate( - p1=(x + 0.01 + 0.005 + 0.025, y + casethickness + 0.005 + 0.203, z + skidthickness - 0.002), - p2=( - x + 0.01 + 0.005 + 0.025 + bowtiebase, - y + casethickness + 0.005 + 0.203 + patchheight, - z + skidthickness - 0.002, - ), - material_id="pec", - ) - # triangles - t1 = gprMax.Triangle( - p1=(x + 0.01 + 0.005 + 0.025, y + casethickness + 0.005 + 0.081, z + skidthickness - 0.002), - p2=(x + 0.01 + 0.005 + 0.025 + bowtiebase, y + casethickness + 0.005 + 0.081, z + skidthickness - 0.002), - p3=( - x + 0.01 + 0.005 + 0.025 + (bowtiebase / 2), - y + casethickness + 0.005 + 0.081 + bowtieheight, - z + skidthickness - 0.002, - ), - thickness=0, - material_id="pec", - ) - t2 = gprMax.Triangle( - p1=(x + 0.01 + 0.005 + 0.025, y + casethickness + 0.005 + 0.203, z + skidthickness - 0.002), - p2=(x + 0.01 + 0.005 + 0.025 + bowtiebase, y + casethickness + 0.005 + 0.203, z + skidthickness - 0.002), - p3=( - x + 0.01 + 0.005 + 0.025 + (bowtiebase / 2), - y + casethickness + 0.005 + 0.203 - bowtieheight, - z + skidthickness - 0.002, - ), - thickness=0, - material_id="pec", - ) - # "right" side - p3 = gprMax.Plate( - p1=(x + 0.01 + 0.005 + 0.187, y + casethickness + 0.005 + 0.021, z + skidthickness - 0.002), - p2=( - x + 0.01 + 0.005 + 0.187 + bowtiebase, - y + casethickness + 0.005 + 0.021 + patchheight, - z + skidthickness - 0.002, - ), - material_id="pec", - ) - p4 = gprMax.Plate( - p1=(x + 0.01 + 0.005 + 0.187, y + casethickness + 0.005 + 0.203, z + skidthickness - 0.002), - p2=( - x + 0.01 + 0.005 + 0.187 + bowtiebase, - y + casethickness + 0.005 + 0.203 + patchheight, - z + skidthickness - 0.002, - ), - material_id="pec", - ) - # triangles - t3 = gprMax.Triangle( - p1=(x + 0.01 + 0.005 + 0.187, y + casethickness + 0.005 + 0.081, z + skidthickness - 0.002), - p2=(x + 0.01 + 0.005 + 0.187 + bowtiebase, y + casethickness + 0.005 + 0.081, z + skidthickness - 0.002), - p3=( - x + 0.01 + 0.005 + 0.187 + (bowtiebase / 2), - y + casethickness + 0.005 + 0.081 + bowtieheight, - z + skidthickness - 0.002, - ), - thickness=0, - material_id="pec", - ) - t4 = gprMax.Triangle( - p1=(x + 0.01 + 0.005 + 0.187, y + casethickness + 0.005 + 0.203, z + skidthickness - 0.002), - p2=(x + 0.01 + 0.005 + 0.187 + bowtiebase, y + casethickness + 0.005 + 0.203, z + skidthickness - 0.002), - p3=( - x + 0.01 + 0.005 + 0.187 + (bowtiebase / 2), - y + casethickness + 0.005 + 0.203 - bowtieheight, - z + skidthickness - 0.002, - ), - thickness=0, - material_id="pec", - ) - - # Edges that represent wire between bowtie halves in 2mm model - e1 = gprMax.Edge(p1=(tx[0] + 0.162, tx[1] - dy, tx[2]), p2=(tx[0] + 0.162, tx[1], tx[2]), material_id="pec") - e2 = gprMax.Edge( - p1=(tx[0] + 0.162, tx[1] + dy, tx[2]), p2=(tx[0] + 0.162, tx[1] + 2 * dy, tx[2]), material_id="pec" - ) - e3 = gprMax.Edge(p1=(tx[0], tx[1] - dy, tx[2]), p2=(tx[0], tx[1], tx[2]), material_id="pec") - e4 = gprMax.Edge(p1=(tx[0], tx[1] + dy, tx[2]), p2=(tx[0], tx[1] + 2 * dy, tx[2]), material_id="pec") - scene_objects.extend((p1, p2, t1, t2, p3, p4, t3, t4, e1, e2, e3, e4)) + # Edges that represent wire between bowtie halves in 2mm model + e1 = gprMax.Edge(p1=(tx[0] + 0.162, tx[1] - dy, tx[2]), p2=(tx[0] + 0.162, tx[1], tx[2]), material_id="pec") + e2 = gprMax.Edge( + p1=(tx[0] + 0.162, tx[1] + dy, tx[2]), p2=(tx[0] + 0.162, tx[1] + 2 * dy, tx[2]), material_id="pec" + ) + e3 = gprMax.Edge(p1=(tx[0], tx[1] - dy, tx[2]), p2=(tx[0], tx[1], tx[2]), material_id="pec") + e4 = gprMax.Edge(p1=(tx[0], tx[1] + dy, tx[2]), p2=(tx[0], tx[1] + 2 * dy, tx[2]), material_id="pec") + scene_objects.extend((p1, p2, t1, t2, p3, p4, t3, t4, e1, e2, e3, e4)) # Metallic plate extension b15 = gprMax.Box( @@ -1098,31 +799,17 @@ def antenna_like_GSSI_400(x, y, z, resolution=0.001, **kwargs): w2 = gprMax.Waveform(wave_type="gaussian", amp=0, freq=excitationfreq, id="my_zero_wave") scene_objects.append(w2) - if resolution == 0.001 or resolution == 0.0005: - if src_type == "transmission_line": - tl2 = gprMax.TransmissionLine( - polarisation="y", - p1=(tx[0] + 0.16, tx[1], tx[2]), - resistance=receiverresistance, - waveform_id="my_zero_wave", - ) - scene_objects.append(tl2) - elif src_type == "voltage_source": - r1 = gprMax.Rx(p1=(tx[0] + 0.16, tx[1], tx[2]), id="rxbowtie", outputs="Ey") - scene_objects.append(r1) - - elif resolution == 0.002: - if src_type == "transmission_line": - tl2 = gprMax.TransmissionLine( - polarisation="y", - p1=(tx[0] + 0.162, tx[1], tx[2]), - resistance=receiverresistance, - waveform_id="my_zero_wave", - ) - scene_objects.append(tl2) - elif src_type == "voltage_source": - r1 = gprMax.Rx(p1=(tx[0] + 0.162, tx[1], tx[2]), id="rxbowtie", outputs="Ey") - scene_objects.append(r1) + if src_type == "transmission_line": + tl2 = gprMax.TransmissionLine( + polarisation="y", + p1=(tx[0] + 0.162, tx[1], tx[2]), + resistance=receiverresistance, + waveform_id="my_zero_wave", + ) + scene_objects.append(tl2) + elif src_type == "voltage_source": + r1 = gprMax.Rx(p1=(tx[0] + 0.162, tx[1], tx[2]), id="rxbowtie", outputs="Ey") + scene_objects.append(r1) # Geometry views gv1 = gprMax.GeometryView( From ab922a97ebbdd60654888f8bb037a6601b7bd14a Mon Sep 17 00:00:00 2001 From: Craig Warren Date: Wed, 30 Aug 2023 09:28:45 +0100 Subject: [PATCH 02/35] Updated email & web links --- toolboxes/GPRAntennaModels/README.rst | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/toolboxes/GPRAntennaModels/README.rst b/toolboxes/GPRAntennaModels/README.rst index a5fcd241..e15b8c0d 100644 --- a/toolboxes/GPRAntennaModels/README.rst +++ b/toolboxes/GPRAntennaModels/README.rst @@ -14,7 +14,7 @@ Manufacturer/Model Dimensions Resolution(s) Author/Contact ======================== ============= ============= ========================================================================================================================================================================================================================= ================ GSSI 1.5GHz (Model 5100) 170x108x45mm 1, 2mm Craig Warren (craig.warren@northumbria.ac.uk), Northumbria University, UK 1,2 MALA 1.2GHz 184x109x46mm 1, 2mm Craig Warren (craig.warren@northumbria.ac.uk), Northumbria University, UK 1 -GSSI 400MHz 300x300x170mm 0.5, 1, 2mm Sam Stadler (Sam.Stadler@liag-hannover.de), `Leibniz Institute for Applied Geophysics `_, Germany 3 +GSSI 400MHz 300x300x170mm 2mm Sam Stadler (sam.stadler@leibniz-liag.de), `Leibniz Institute for Applied Geophysics `_, Germany 3 ======================== ============= ============= ========================================================================================================================================================================================================================= ================ **License**: `Creative Commons Attribution-ShareAlike 4.0 International License `_ @@ -36,7 +36,7 @@ Descriptions of how the models were created can be found in the aforementioned a How to use the package ====================== -The antenna models can be accessed from within a block of Python code in an input file. The models are inserted at location x,y,z. The coordinates are relative to the geometric centre of the antenna in the x-y plane and the bottom of the antenna skid in the z direction. The models must be used with cubic spatial resolutions of either 0.5mm (GSSI 400MHz antenna only), 1mm (default), or 2mm by setting the keyword argument, e.g. ``resolution=0.002``. +The antenna models can be accessed from within a block of Python code in an input file. The models are inserted at location x,y,z. The coordinates are relative to the geometric centre of the antenna in the x-y plane and the bottom of the antenna skid in the z direction. The models must be used with cubic spatial resolutions of either 1mm or 2mm by setting the keyword argument, e.g. ``resolution=0.002``. .. note:: From c91ab8331533b17ac4690f7da9fc619470c6e3fd Mon Sep 17 00:00:00 2001 From: Craig Warren Date: Wed, 30 Aug 2023 10:08:53 +0100 Subject: [PATCH 03/35] Improved compiler opts for macOS --- setup.py | 45 ++++++++++++++++++--------------------------- 1 file changed, 18 insertions(+), 27 deletions(-) diff --git a/setup.py b/setup.py index df082040..7c201976 100644 --- a/setup.py +++ b/setup.py @@ -18,6 +18,7 @@ import glob import os +import platform import re import shutil import subprocess @@ -187,35 +188,25 @@ else: ) cpuID = " ".join(cpuID.split()) if "Apple" in cpuID: - gccpath = glob.glob("/opt/homebrew/bin/gcc-[4-9]*") - gccpath += glob.glob("/opt/homebrew/bin/gcc-[10-11]*") - if gccpath: - # Use newest gcc found - os.environ["CC"] = gccpath[-1].split(os.sep)[-1] - rpath = "/opt/homebrew/opt/gcc/lib/gcc/" + gccpath[-1].split(os.sep)[-1][-1] + "/" - else: - raise ( - "Cannot find gcc in /opt/homebrew/bin. gprMax requires gcc " - + "to be installed - easily done through the Homebrew package " - + "manager (http://brew.sh). Note: gcc with OpenMP support " - + "is required." - ) + gccbasepath = "/opt/homebrew/bin/" else: - gccpath = glob.glob("/usr/local/bin/gcc-[4-9]*") - gccpath += glob.glob("/usr/local/bin/gcc-[10-11]*") - if gccpath: - # Use newest gcc found - os.environ["CC"] = gccpath[-1].split(os.sep)[-1] - else: - raise ( - "Cannot find gcc in /usr/local/bin. gprMax requires gcc " - + "to be installed - easily done through the Homebrew package " - + "manager (http://brew.sh). Note: gcc with OpenMP support " - + "is required." - ) + gccbasepath = "/usr/local/bin/" + gccpath = glob.glob(gccbasepath + "gcc-[0-9][0-9]") + if gccpath: + # Use newest gcc found + os.environ["CC"] = gccpath[-1].split(os.sep)[-1] + if "Apple" in cpuID: + rpath = "/opt/homebrew/opt/gcc/lib/gcc/" + gccpath[-1].split(os.sep)[-1][-1] + "/" + else: + raise ( + f"Cannot find gcc in {gccbasepath}. gprMax requires gcc " + + "to be installed - easily done through the Homebrew package " + + "manager (http://brew.sh). Note: gcc with OpenMP support " + + "is required." + ) - # Minimum supported macOS deployment target - MIN_MACOS_VERSION = "10.13" + # Set minimum supported macOS deployment target to installed macOS version + MIN_MACOS_VERSION = platform.mac_ver()[0] try: os.environ["MACOSX_DEPLOYMENT_TARGET"] del os.environ["MACOSX_DEPLOYMENT_TARGET"] From c3e357e7b736efff3af93c6a9604dcca5b11be2c Mon Sep 17 00:00:00 2001 From: Craig Warren Date: Wed, 30 Aug 2023 10:50:48 +0100 Subject: [PATCH 04/35] Fixed typo bug with custom waveforms --- gprMax/cmds_singleuse.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/gprMax/cmds_singleuse.py b/gprMax/cmds_singleuse.py index 2570bb86..b2560548 100644 --- a/gprMax/cmds_singleuse.py +++ b/gprMax/cmds_singleuse.py @@ -433,7 +433,7 @@ class ExcitationFile(UserObjectSingle): waveformtime = np.arange(0, G.timewindow + G.dt, G.dt) timestr = "simulation time array" - for waveformID in waveformIDs: + for i, waveformID in enumerate(waveformIDs): if any(x.ID == waveformID for x in G.waveforms): logger.exception(f"Waveform with ID {waveformID} already exists") raise ValueError @@ -442,7 +442,7 @@ class ExcitationFile(UserObjectSingle): w.type = "user" # Select correct column of waveform values depending on array shape - singlewaveformvalues = waveformvalues[:] if len(waveformvalues.shape) == 1 else waveformvalues[:, waveform] + singlewaveformvalues = waveformvalues[:] if len(waveformvalues.shape) == 1 else waveformvalues[:, i] # Truncate waveform array if it is longer than time array if len(singlewaveformvalues) > len(waveformtime): From 8891e89de5d660f55625ed7e715787bf6b199a27 Mon Sep 17 00:00:00 2001 From: Craig Warren Date: Wed, 30 Aug 2023 10:51:26 +0100 Subject: [PATCH 05/35] Corrected some bugs with 400MHz model --- examples/antenna_like_GSSI_400_fs.py | 15 +++++++++------ toolboxes/GPRAntennaModels/GSSI.py | 20 ++++++++------------ 2 files changed, 17 insertions(+), 18 deletions(-) diff --git a/examples/antenna_like_GSSI_400_fs.py b/examples/antenna_like_GSSI_400_fs.py index 1477c12f..56c52b6b 100644 --- a/examples/antenna_like_GSSI_400_fs.py +++ b/examples/antenna_like_GSSI_400_fs.py @@ -30,7 +30,7 @@ scene = gprMax.Scene() title = gprMax.Title(name=fn.with_suffix("").name) domain = gprMax.Domain(p1=(x, y, z)) dxdydz = gprMax.Discretisation(p1=(dl, dl, dl)) -time_window = gprMax.TimeWindow(time=6e-9) +time_window = gprMax.TimeWindow(time=15e-9) scene.add(title) scene.add(domain) @@ -44,17 +44,20 @@ for obj in gssi_objects: scene.add(obj) gv1 = gprMax.GeometryView( - p1=(0, 0, 0), p2=(x, y, z), dl=(dl, dl, dl), filename="antenna_like_GSSI_400", output_type="n" + p1=(0, 0, 0), + p2=(x, y, z), + dl=(dl, dl, dl), + filename="antenna_like_GSSI_400", output_type="n" ) gv2 = gprMax.GeometryView( - p1=(ant_pos[0] - 0.150 / 2, ant_pos[1] - 0.150 / 2, ant_pos[2] - 0.050), - p2=(ant_pos[0] + 0.150 / 2, ant_pos[1] + 0.150 / 2, ant_pos[2] + 0.010), + p1=(ant_pos[0] - 0.150, ant_pos[1] - 0.150, ant_pos[2]), + p2=(ant_pos[0] + 0.150, ant_pos[1] + 0.150, ant_pos[2] + 0.050), dl=(dl, dl, dl), filename="antenna_like_GSSI_400_pcb", output_type="f", ) -scene.add(gv1) -scene.add(gv2) +#scene.add(gv1) +#scene.add(gv2) # Run model gprMax.run(scenes=[scene], geometry_only=False, outputfile=fn, gpu=None) diff --git a/toolboxes/GPRAntennaModels/GSSI.py b/toolboxes/GPRAntennaModels/GSSI.py index 6273d61a..58bf1c9e 100644 --- a/toolboxes/GPRAntennaModels/GSSI.py +++ b/toolboxes/GPRAntennaModels/GSSI.py @@ -770,20 +770,18 @@ def antenna_like_GSSI_400(x, y, z, resolution=0.002, **kwargs): scene_objects.extend((b15, b16)) # Source - # Excitation - Gaussian pulse - w1 = gprMax.Waveform(wave_type="gaussian", amp=1, freq=excitationfreq, id="my_gaussian") - scene_objects.append(w1) - if src_type == "voltage_source": + w1 = gprMax.Waveform(wave_type="gaussian", amp=1, freq=excitationfreq, id="my_gaussian") vs1 = gprMax.VoltageSource( polarisation="y", p1=(tx[0], tx[1], tx[2]), resistance=sourceresistance, waveform_id="my_gaussian" ) - scene_objects.append(vs1) + scene_objects.extend((w1, vs1)) elif src_type == "transmission_line": + w1 = gprMax.Waveform(wave_type="gaussian", amp=1, freq=excitationfreq, id="my_gaussian") tl1 = gprMax.TransmissionLine( polarisation="y", p1=(tx[0], tx[1], tx[2]), resistance=sourceresistance, waveform_id="my_gaussian" ) - scene_objects.append(tl1) + scene_objects.extend((w1, tl1)) else: # Optimised custom pulse exc1 = gprMax.ExcitationFile( @@ -795,19 +793,17 @@ def antenna_like_GSSI_400(x, y, z, resolution=0.002, **kwargs): scene_objects.extend((exc1, vs1)) # Receiver - # Zero waveform to use with transmission line at receiver output - w2 = gprMax.Waveform(wave_type="gaussian", amp=0, freq=excitationfreq, id="my_zero_wave") - scene_objects.append(w2) - if src_type == "transmission_line": + # Zero waveform to use with transmission line at receiver output + w2 = gprMax.Waveform(wave_type="gaussian", amp=0, freq=excitationfreq, id="my_zero_wave") tl2 = gprMax.TransmissionLine( polarisation="y", p1=(tx[0] + 0.162, tx[1], tx[2]), resistance=receiverresistance, waveform_id="my_zero_wave", ) - scene_objects.append(tl2) - elif src_type == "voltage_source": + scene_objects.extend((w2, tl2)) + else: r1 = gprMax.Rx(p1=(tx[0] + 0.162, tx[1], tx[2]), id="rxbowtie", outputs="Ey") scene_objects.append(r1) From 9763280b61efa9aa6998417f64fe203646cd902d Mon Sep 17 00:00:00 2001 From: Craig Warren Date: Thu, 31 Aug 2023 12:10:06 +0100 Subject: [PATCH 06/35] Fixed typo in voltage src and mag dipole. --- gprMax/cuda_opencl/knl_source_updates.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/gprMax/cuda_opencl/knl_source_updates.py b/gprMax/cuda_opencl/knl_source_updates.py index 5e1e6372..6e2a5c41 100644 --- a/gprMax/cuda_opencl/knl_source_updates.py +++ b/gprMax/cuda_opencl/knl_source_updates.py @@ -129,7 +129,7 @@ update_magnetic_dipole = { __global const int* restrict rcinfo1, __global const $REAL* restrict rcinfo2, __global const $REAL* restrict rcwaveforms, - __global const unsigned int* estrict ID, + __global const unsigned int* restrict ID, __global $REAL *Hx, __global $REAL *Hy, __global $REAL *Hz @@ -211,7 +211,7 @@ update_voltage_source = { __global const int* restrict rcinfo1, __global const $REAL* restrict rcinfo2, __global const $REAL* restrict rcwaveforms, - __global const unsigned int* estrict ID, + __global const unsigned int* restrict ID, __global $REAL *Ex, __global $REAL *Ey, __global $REAL *Ez From 6bab9e32375aebaad80813bf2bb6b7182afd15cf Mon Sep 17 00:00:00 2001 From: Craig Warren Date: Thu, 31 Aug 2023 12:14:03 +0100 Subject: [PATCH 07/35] Fixed further typos with voltage src and mag dipole opencl --- gprMax/cuda_opencl/knl_source_updates.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/gprMax/cuda_opencl/knl_source_updates.py b/gprMax/cuda_opencl/knl_source_updates.py index 6e2a5c41..85203c32 100644 --- a/gprMax/cuda_opencl/knl_source_updates.py +++ b/gprMax/cuda_opencl/knl_source_updates.py @@ -126,9 +126,9 @@ update_magnetic_dipole = { $REAL dx, $REAL dy, $REAL dz, - __global const int* restrict rcinfo1, - __global const $REAL* restrict rcinfo2, - __global const $REAL* restrict rcwaveforms, + __global const int* restrict srcinfo1, + __global const $REAL* restrict srcinfo2, + __global const $REAL* restrict srcwaveforms, __global const unsigned int* restrict ID, __global $REAL *Hx, __global $REAL *Hy, @@ -208,9 +208,9 @@ update_voltage_source = { $REAL dx, $REAL dy, $REAL dz, - __global const int* restrict rcinfo1, - __global const $REAL* restrict rcinfo2, - __global const $REAL* restrict rcwaveforms, + __global const int* restrict srcinfo1, + __global const $REAL* restrict srcinfo2, + __global const $REAL* restrict srcwaveforms, __global const unsigned int* restrict ID, __global $REAL *Ex, __global $REAL *Ey, From 2d08accca39a1094abf5fcbb553a9c5324385ba9 Mon Sep 17 00:00:00 2001 From: Craig Warren Date: Thu, 31 Aug 2023 13:41:54 +0100 Subject: [PATCH 08/35] Removed unused max outputs variable --- gprMax/receivers.py | 5 ----- 1 file changed, 5 deletions(-) diff --git a/gprMax/receivers.py b/gprMax/receivers.py index f24c8263..64193636 100644 --- a/gprMax/receivers.py +++ b/gprMax/receivers.py @@ -26,9 +26,7 @@ class Rx: allowableoutputs = ["Ex", "Ey", "Ez", "Hx", "Hy", "Hz", "Ix", "Iy", "Iz"] defaultoutputs = allowableoutputs[:-3] - allowableoutputs_dev = allowableoutputs[:-3] - maxnumoutputs_dev = 0 def __init__(self): self.ID = None @@ -61,9 +59,6 @@ def htod_rx_arrays(G, queue=None): rxcoords[i, 0] = rx.xcoord rxcoords[i, 1] = rx.ycoord rxcoords[i, 2] = rx.zcoord - # Store maximum number of output components - if len(rx.outputs) > Rx.maxnumoutputs_dev: - Rx.maxnumoutputs_dev = len(rx.outputs) # Array to store field components for receivers on compute device - # rows are field components; columns are iterations; pages are receivers From 596dc6e3301d9d1b1d41f3e799f987789b05ed97 Mon Sep 17 00:00:00 2001 From: Craig Warren Date: Thu, 31 Aug 2023 13:42:20 +0100 Subject: [PATCH 09/35] Removed refs to PML as these are tested separately --- testing/test_models.py | 7 +------ 1 file changed, 1 insertion(+), 6 deletions(-) diff --git a/testing/test_models.py b/testing/test_models.py index 196617b9..18c0e13d 100644 --- a/testing/test_models.py +++ b/testing/test_models.py @@ -40,14 +40,12 @@ if sys.platform == "linux": python -m testing.test_models """ -# Specify directoty with set of models to test +# Specify directory with set of models to test modelset = "models_basic" # modelset += 'models_advanced' -# modelset += 'models_pmls' basepath = Path(__file__).parents[0] / modelset - # List of available basic test models testmodels = [ "hertzian_dipole_fs_analytical", @@ -64,9 +62,6 @@ testmodels = [ # List of available advanced test models # testmodels = ['antenna_GSSI_1500_fs', 'antenna_MALA_1200_fs'] -# List of available PML models -# testmodels = ['pml_x0', 'pml_y0', 'pml_z0', 'pml_xmax', 'pml_ymax', 'pml_zmax', 'pml_3D_pec_plate'] - # Select a specific model if desired # testmodels = [testmodels[0]] testresults = dict.fromkeys(testmodels) From 631d9a8414b2942d9036838ede2d1b8b1c2dc84b Mon Sep 17 00:00:00 2001 From: Craig Warren Date: Thu, 31 Aug 2023 14:47:47 +0100 Subject: [PATCH 10/35] Added more options for checking PMLs --- testing/diff_output_files.py | 93 ++++++++++++ testing/models_pmls/plot_pml_comparison.py | 135 ------------------ .../pml_3D_pec_plate/plot_pml_comparison.py | 80 +++++++++++ .../pml_basic/plot_pml_comparison.py | 75 ++++++++++ .../models_pmls/{ => pml_basic}/pml_basic.py | 0 5 files changed, 248 insertions(+), 135 deletions(-) create mode 100644 testing/diff_output_files.py delete mode 100644 testing/models_pmls/plot_pml_comparison.py create mode 100644 testing/models_pmls/pml_3D_pec_plate/plot_pml_comparison.py create mode 100644 testing/models_pmls/pml_basic/plot_pml_comparison.py rename testing/models_pmls/{ => pml_basic}/pml_basic.py (100%) diff --git a/testing/diff_output_files.py b/testing/diff_output_files.py new file mode 100644 index 00000000..b2a3da6b --- /dev/null +++ b/testing/diff_output_files.py @@ -0,0 +1,93 @@ +# Copyright (C) 2015-2023: The University of Edinburgh, United Kingdom +# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley +# +# This file is part of gprMax. +# +# gprMax is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# gprMax is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with gprMax. If not, see . + +import logging +from pathlib import Path + +import h5py +import numpy as np + +logger = logging.getLogger(__name__) + +def diff_output_files(filename1, filename2): + """Calculates differences between two output files. + + Args: + filename1: string of filename (including path) of output file 1. + filename2: string of filename (including path) of output file 2. + + Returns: + time: numpy array containing time. + datadiffs: numpy array containing power (dB) of differences. + """ + + file1 = h5py.File(Path(filename1), "r") + file2 = h5py.File(Path(filename2), "r") + # Path to receivers in files + path = "rxs/rx1/" + + # Get available field output component names + outputs1 = list(file1[path].keys()) + outputs2 = list(file2[path].keys()) + if outputs1 != outputs2: + logger.exception("Field output components are not the same in each file") + raise ValueError + + # Check that type of float used to store fields matches + floattype1 = file1[path + outputs1[0]].dtype + floattype2 = file2[path + outputs2[0]].dtype + if floattype1 != floattype2: + logger.warning( + f"Type of floating point number in test model ({file1[path + outputs1[0]].dtype}) " + f"does not match type in reference solution ({file2[path + outputs2[0]].dtype})\n" + ) + + # Arrays for storing time + time1 = np.zeros((file1.attrs["Iterations"]), dtype=floattype1) + time1 = np.linspace(0, (file1.attrs["Iterations"] - 1), num=file1.attrs["Iterations"]) + time2 = np.zeros((file2.attrs["Iterations"]), dtype=floattype2) + time2 = np.linspace(0, (file2.attrs["Iterations"] - 1), num=file2.attrs["Iterations"]) + + # Arrays for storing field data + data1 = np.zeros((file1.attrs["Iterations"], len(outputs1)), dtype=floattype1) + data2 = np.zeros((file2.attrs["Iterations"], len(outputs2)), dtype=floattype2) + for ID, name in enumerate(outputs1): + data1[:, ID] = file1[path + str(name)][:] + data2[:, ID] = file2[path + str(name)][:] + if np.any(np.isnan(data1[:, ID])) or np.any(np.isnan(data2[:, ID])): + logger.exception("Data contains NaNs") + raise ValueError + + file1.close() + file2.close() + + # Diffs + datadiffs = np.zeros(data1.shape, dtype=np.float64) + for i in range(len(outputs2)): + maxi = np.amax(np.abs(data1[:, i])) + datadiffs[:, i] = np.divide( + np.abs(data2[:, i] - data1[:, i]), maxi, out=np.zeros_like(data1[:, i]), where=maxi != 0 + ) # Replace any division by zero with zero + + # Calculate power (ignore warning from taking a log of any zero values) + with np.errstate(divide="ignore"): + datadiffs[:, i] = 20 * np.log10(datadiffs[:, i]) + # Replace any NaNs or Infs from zero division + datadiffs[:, i][np.invert(np.isfinite(datadiffs[:, i]))] = 0 + + return time1, datadiffs diff --git a/testing/models_pmls/plot_pml_comparison.py b/testing/models_pmls/plot_pml_comparison.py deleted file mode 100644 index 95c5a79f..00000000 --- a/testing/models_pmls/plot_pml_comparison.py +++ /dev/null @@ -1,135 +0,0 @@ -# Copyright (C) 2015-2023: The University of Edinburgh, United Kingdom -# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley -# -# This file is part of gprMax. -# -# gprMax is free software: you can redistribute it and/or modify -# it under the terms of the GNU General Public License as published by -# the Free Software Foundation, either version 3 of the License, or -# (at your option) any later version. -# -# gprMax is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU General Public License for more details. -# -# You should have received a copy of the GNU General Public License -# along with gprMax. If not, see . - -import itertools -import logging -from operator import add -from pathlib import Path - -import h5py -import matplotlib.pyplot as plt -import numpy as np - -logger = logging.getLogger(__name__) - -# Create/setup plot figure -# colors = ['#E60D30', '#5CB7C6', '#A21797', '#A3B347'] # Plot colours from http://tools.medialab.sciences-po.fr/iwanthue/index.php -# colorIDs = ["#62a85b", "#9967c7", "#b3943f", "#6095cd", "#cb5c42", "#c95889"] -colorIDs = ["#79c72e", "#5774ff", "#ff7c2c", "#4b4e80", "#d7004e", "#007545", "#ff83ec"] -# colorIDs = ["#ba0044", "#b2d334", "#470055", "#185300", "#ff96b1", "#3e2700", "#0162a9", "#fdb786"] -colors = itertools.cycle(colorIDs) -# for i in range(2): -# next(colors) -lines = itertools.cycle(("--", ":", "-.", "-")) -markers = ["o", "d", "^", "s", "*"] - -parts = Path(__file__).parts -path = "rxs/rx1/" -basename = "pml_3D_pec_plate" -PMLIDs = ["CFS-PML", "HORIPML-1", "HORIPML-2", "MRIPML-1", "MRIPML-2"] -maxerrors = [] -testmodels = ["pml_3D_pec_plate_" + s for s in PMLIDs] - -fig, ax = plt.subplots( - subplot_kw=dict(xlabel="Iterations", ylabel="Error [dB]"), figsize=(20, 10), facecolor="w", edgecolor="w" -) - -for x, model in enumerate(testmodels): - # Open output file and read iterations - fileref = h5py.File(Path(*parts[:-1], basename, basename + "_ref.h5"), "r") - filetest = h5py.File(Path(*parts[:-1], basename, basename + str(x + 1) + ".h5"), "r") - - # Get available field output component names - outputsref = list(fileref[path].keys()) - outputstest = list(filetest[path].keys()) - if outputsref != outputstest: - logger.exception("Field output components do not match reference solution") - raise ValueError - - # Check that type of float used to store fields matches - if filetest[path + outputstest[0]].dtype != fileref[path + outputsref[0]].dtype: - logger.warning( - f"Type of floating point number in test model ({filetest[path + outputstest[0]].dtype}) " - f"does not match type in reference solution ({fileref[path + outputsref[0]].dtype})\n" - ) - floattyperef = fileref[path + outputsref[0]].dtype - floattypetest = filetest[path + outputstest[0]].dtype - # logger.info(f'Data type: {floattypetest}') - - # Arrays for storing time - # timeref = np.zeros((fileref.attrs['Iterations']), dtype=floattyperef) - # timeref = np.linspace(0, (fileref.attrs['Iterations'] - 1) * fileref.attrs['dt'], num=fileref.attrs['Iterations']) / 1e-9 - # timetest = np.zeros((filetest.attrs['Iterations']), dtype=floattypetest) - # timetest = np.linspace(0, (filetest.attrs['Iterations'] - 1) * filetest.attrs['dt'], num=filetest.attrs['Iterations']) / 1e-9 - timeref = np.zeros((fileref.attrs["Iterations"]), dtype=floattyperef) - timeref = np.linspace(0, (fileref.attrs["Iterations"] - 1), num=fileref.attrs["Iterations"]) - timetest = np.zeros((filetest.attrs["Iterations"]), dtype=floattypetest) - timetest = np.linspace(0, (filetest.attrs["Iterations"] - 1), num=filetest.attrs["Iterations"]) - - # Arrays for storing field data - dataref = np.zeros((fileref.attrs["Iterations"], len(outputsref)), dtype=floattyperef) - datatest = np.zeros((filetest.attrs["Iterations"], len(outputstest)), dtype=floattypetest) - for ID, name in enumerate(outputsref): - dataref[:, ID] = fileref[path + str(name)][:] - datatest[:, ID] = filetest[path + str(name)][:] - if np.any(np.isnan(datatest[:, ID])): - logger.exception("Test data contains NaNs") - raise ValueError - - fileref.close() - filetest.close() - - # Diffs - datadiffs = np.zeros(datatest.shape, dtype=np.float64) - for i in range(len(outputstest)): - maxi = np.amax(np.abs(dataref[:, i])) - datadiffs[:, i] = np.divide( - np.abs(datatest[:, i] - dataref[:, i]), maxi, out=np.zeros_like(dataref[:, i]), where=maxi != 0 - ) # Replace any division by zero with zero - - # Calculate power (ignore warning from taking a log of any zero values) - with np.errstate(divide="ignore"): - datadiffs[:, i] = 20 * np.log10(datadiffs[:, i]) - # Replace any NaNs or Infs from zero division - datadiffs[:, i][np.invert(np.isfinite(datadiffs[:, i]))] = 0 - - # Print maximum error value - start = 210 - maxerrors.append(f": {np.amax(datadiffs[start::, 1]):.1f} [dB]") - logger.info(f"{model}: Max. error {maxerrors[x]}") - - # Plot diffs (select column to choose field component, 0-Ex, 1-Ey etc..) - ax.plot(timeref[start::], datadiffs[start::, 1], color=next(colors), lw=2, ls=next(lines), label=model) - ax.set_xticks(np.arange(0, 2200, step=100)) - ax.set_xlim([0, 2100]) - ax.set_yticks(np.arange(-160, 0, step=20)) - ax.set_ylim([-160, -20]) - ax.set_axisbelow(True) - ax.grid(color=(0.75, 0.75, 0.75), linestyle="dashed") - -mylegend = list(map(add, PMLIDs, maxerrors)) -legend = ax.legend(mylegend, loc=1, fontsize=14) -frame = legend.get_frame() -frame.set_edgecolor("white") -frame.set_alpha(0) - -plt.show() - -# Save a PDF/PNG of the figure -# fig.savefig(basepath + '.pdf', dpi=None, format='pdf', bbox_inches='tight', pad_inches=0.1) -# fig.savefig(savename + '.png', dpi=150, format='png', bbox_inches='tight', pad_inches=0.1) diff --git a/testing/models_pmls/pml_3D_pec_plate/plot_pml_comparison.py b/testing/models_pmls/pml_3D_pec_plate/plot_pml_comparison.py new file mode 100644 index 00000000..52968ef1 --- /dev/null +++ b/testing/models_pmls/pml_3D_pec_plate/plot_pml_comparison.py @@ -0,0 +1,80 @@ +# Copyright (C) 2015-2023: The University of Edinburgh, United Kingdom +# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley +# +# This file is part of gprMax. +# +# gprMax is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# gprMax is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with gprMax. If not, see . + +import itertools +import logging +from operator import add +from pathlib import Path + +import matplotlib.pyplot as plt +import numpy as np + +from testing.diff_output_files import diff_output_files + +logger = logging.getLogger(__name__) + +# Create/setup plot figure +# colors = ['#E60D30', '#5CB7C6', '#A21797', '#A3B347'] # Plot colours from http://tools.medialab.sciences-po.fr/iwanthue/index.php +# colorIDs = ["#62a85b", "#9967c7", "#b3943f", "#6095cd", "#cb5c42", "#c95889"] +colorIDs = ["#79c72e", "#5774ff", "#ff7c2c", "#4b4e80", "#d7004e", "#007545", "#ff83ec"] +# colorIDs = ["#ba0044", "#b2d334", "#470055", "#185300", "#ff96b1", "#3e2700", "#0162a9", "#fdb786"] +colors = itertools.cycle(colorIDs) +# for i in range(2): +# next(colors) +lines = itertools.cycle(("--", ":", "-.", "-")) +markers = ["o", "d", "^", "s", "*"] + +fn = Path(__file__) +basename = "pml_3D_pec_plate" +PMLIDs = ["CFS-PML", "HORIPML-1", "HORIPML-2", "MRIPML-1", "MRIPML-2"] +maxerrors = [] +testmodels = [basename + "_" + s for s in PMLIDs] + +fig, ax = plt.subplots( + subplot_kw=dict(xlabel="Iterations", ylabel="Error [dB]"), figsize=(20, 10), facecolor="w", edgecolor="w" +) + +for x, model in enumerate(testmodels): + time, datadiffs = diff_output_files(fn.parent.joinpath(basename + "_ref.h5"), + fn.parent.joinpath(basename + str(x + 1) + ".h5")) + + # Print maximum error value + start = 210 + maxerrors.append(f": {np.amax(datadiffs[start::, 1]):.1f} [dB]") + logger.info(f"{model}: Max. error {maxerrors[x]}") + + # Plot diffs (select column to choose field component, 0-Ex, 1-Ey etc..) + ax.plot(time[start::], datadiffs[start::, 1], color=next(colors), lw=2, ls=next(lines), label=model) + ax.set_xticks(np.arange(0, 2200, step=100)) + ax.set_xlim([0, 2100]) + ax.set_yticks(np.arange(-160, 0, step=20)) + ax.set_ylim([-160, -20]) + ax.set_axisbelow(True) + ax.grid(color=(0.75, 0.75, 0.75), linestyle="dashed") + +mylegend = list(map(add, PMLIDs, maxerrors)) +legend = ax.legend(mylegend, loc=1, fontsize=14) +frame = legend.get_frame() +frame.set_edgecolor("white") +frame.set_alpha(0) + +plt.show() + +# Save a PDF/PNG of the figure +# fig.savefig(basepath + '.pdf', dpi=None, format='pdf', bbox_inches='tight', pad_inches=0.1) +# fig.savefig(savename + '.png', dpi=150, format='png', bbox_inches='tight', pad_inches=0.1) diff --git a/testing/models_pmls/pml_basic/plot_pml_comparison.py b/testing/models_pmls/pml_basic/plot_pml_comparison.py new file mode 100644 index 00000000..7429dca2 --- /dev/null +++ b/testing/models_pmls/pml_basic/plot_pml_comparison.py @@ -0,0 +1,75 @@ +# Copyright (C) 2015-2023: The University of Edinburgh, United Kingdom +# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley +# +# This file is part of gprMax. +# +# gprMax is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# gprMax is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with gprMax. If not, see . + +import itertools +import logging +from operator import add +from pathlib import Path + +import matplotlib.pyplot as plt +import numpy as np + +from testing.diff_output_files import diff_output_files + +logger = logging.getLogger(__name__) + +# Create/setup plot figure +# Plot colours from http://tools.medialab.sciences-po.fr/iwanthue/index.php +colorIDs = ["#79c72e", "#5774ff", "#ff7c2c", "#4b4e80", "#d7004e", "#007545", "#ff83ec"] +colors = itertools.cycle(colorIDs) +lines = itertools.cycle(("--", ":", "-.", "-")) +markers = ["o", "d", "^", "s", "*"] + +fn = Path(__file__) +basename = "pml_basic" +PMLIDs = ["off", "x0", "y0", "z0", "xmax", "ymax", "zmax"] +maxerrors = [] + +for x in range(len(PMLIDs)): + file1 = fn.parent.joinpath(basename + str(x + 1) + "_CPU.h5") + file2 = fn.parent.joinpath(basename + str(x + 1) + "_GPU.h5") + time, datadiffs = diff_output_files(file1, file2) + + # Print maximum error value + start = 0 + maxerrors.append(f": {np.amax(datadiffs[start::, 1]):.1f} [dB]") + logger.info(f"{file1.name} - {file2.name}: Max. error {maxerrors[x]}") + + fig, ax = plt.subplots( + subplot_kw=dict(xlabel="Iterations", ylabel="Error [dB]"), figsize=(20, 10), facecolor="w", edgecolor="w") + + # Plot diffs (select column to choose field component, 0-Ex, 1-Ey etc..) + ax.plot(time[start::], datadiffs[start::, 1], color=next(colors), lw=2, ls=next(lines), label=f"{file1.name} - {file2.name}") + ax.set_xticks(np.arange(0, 2200, step=100)) + ax.set_xlim([0, 2100]) + ax.set_yticks(np.arange(-160, 0, step=20)) + ax.set_ylim([-160, -20]) + ax.set_axisbelow(True) + ax.grid(color=(0.75, 0.75, 0.75), linestyle="dashed") + +mylegend = list(map(add, PMLIDs, maxerrors)) +legend = ax.legend(mylegend, loc=1, fontsize=14) +frame = legend.get_frame() +frame.set_edgecolor("white") +frame.set_alpha(0) + +plt.show() + +# Save a PDF/PNG of the figure +# fig.savefig(basepath + '.pdf', dpi=None, format='pdf', bbox_inches='tight', pad_inches=0.1) +# fig.savefig(savename + '.png', dpi=150, format='png', bbox_inches='tight', pad_inches=0.1) diff --git a/testing/models_pmls/pml_basic.py b/testing/models_pmls/pml_basic/pml_basic.py similarity index 100% rename from testing/models_pmls/pml_basic.py rename to testing/models_pmls/pml_basic/pml_basic.py From b5fe96ae216c161bbb84d976ce4621303b85cbfa Mon Sep 17 00:00:00 2001 From: Craig Warren Date: Thu, 31 Aug 2023 20:59:24 +0100 Subject: [PATCH 11/35] Update to plot PML diffs --- .../pml_basic/plot_pml_comparison.py | 55 +++++++++++-------- 1 file changed, 31 insertions(+), 24 deletions(-) diff --git a/testing/models_pmls/pml_basic/plot_pml_comparison.py b/testing/models_pmls/pml_basic/plot_pml_comparison.py index 7429dca2..86786b07 100644 --- a/testing/models_pmls/pml_basic/plot_pml_comparison.py +++ b/testing/models_pmls/pml_basic/plot_pml_comparison.py @@ -21,6 +21,7 @@ import logging from operator import add from pathlib import Path +import matplotlib.gridspec as gridspec import matplotlib.pyplot as plt import numpy as np @@ -30,7 +31,7 @@ logger = logging.getLogger(__name__) # Create/setup plot figure # Plot colours from http://tools.medialab.sciences-po.fr/iwanthue/index.php -colorIDs = ["#79c72e", "#5774ff", "#ff7c2c", "#4b4e80", "#d7004e", "#007545", "#ff83ec"] +colorIDs = ["#79c72e", "#5774ff", "#ff7c2c", "#4b4e80", "#d7004e", "#007545"] colors = itertools.cycle(colorIDs) lines = itertools.cycle(("--", ":", "-.", "-")) markers = ["o", "d", "^", "s", "*"] @@ -40,36 +41,42 @@ basename = "pml_basic" PMLIDs = ["off", "x0", "y0", "z0", "xmax", "ymax", "zmax"] maxerrors = [] -for x in range(len(PMLIDs)): +for x, PMLID in enumerate(PMLIDs): file1 = fn.parent.joinpath(basename + str(x + 1) + "_CPU.h5") file2 = fn.parent.joinpath(basename + str(x + 1) + "_GPU.h5") time, datadiffs = diff_output_files(file1, file2) # Print maximum error value - start = 0 - maxerrors.append(f": {np.amax(datadiffs[start::, 1]):.1f} [dB]") - logger.info(f"{file1.name} - {file2.name}: Max. error {maxerrors[x]}") + start = 150 + maxerrors.append(f": {np.amax(np.amax(datadiffs[start::, :])):.1f} [dB]") + print(f"{PMLID}: Max. error {maxerrors[x]}") - fig, ax = plt.subplots( - subplot_kw=dict(xlabel="Iterations", ylabel="Error [dB]"), figsize=(20, 10), facecolor="w", edgecolor="w") + # Plot diffs + gs = gridspec.GridSpec(3, 2, hspace=0.3, wspace=0.3) + fig, ax = plt.subplots(figsize=(20, 10), facecolor="w", edgecolor="w") + ax.remove() + fig.suptitle(f"{PMLID}") + + outputs = ["Ex", "Ey", "Ez", "Hx", "Hy", "Hz"] + for i, output in enumerate(outputs): + if i < 3: + ax = plt.subplot(gs[i, 0]) + else: + ax = plt.subplot(gs[i - 3, 1]) + ax.plot(time[start::], datadiffs[start::, i], color=next(colors), lw=2) + ax.set_xticks(np.arange(0, 1800, step=200)) + ax.set_xlim([0, 1600]) + ax.set_yticks(np.arange(-400, 80, step=40)) + ax.set_ylim([-400, 40]) + ax.set_axisbelow(True) + ax.grid(color=(0.75, 0.75, 0.75), linestyle="dashed") + ax.set_xlabel("Time [iterations]") + ax.set_ylabel(f"{output} error [dB]") - # Plot diffs (select column to choose field component, 0-Ex, 1-Ey etc..) - ax.plot(time[start::], datadiffs[start::, 1], color=next(colors), lw=2, ls=next(lines), label=f"{file1.name} - {file2.name}") - ax.set_xticks(np.arange(0, 2200, step=100)) - ax.set_xlim([0, 2100]) - ax.set_yticks(np.arange(-160, 0, step=20)) - ax.set_ylim([-160, -20]) - ax.set_axisbelow(True) - ax.grid(color=(0.75, 0.75, 0.75), linestyle="dashed") - -mylegend = list(map(add, PMLIDs, maxerrors)) -legend = ax.legend(mylegend, loc=1, fontsize=14) -frame = legend.get_frame() -frame.set_edgecolor("white") -frame.set_alpha(0) + # Save a PDF/PNG of the figure + fig.savefig(basename + "_diffs_" + PMLID + ".pdf", dpi=None, format='pdf', bbox_inches='tight', pad_inches=0.1) + # fig.savefig(basename + "_diffs_" + PMLID + ".png", dpi=150, format='png', bbox_inches='tight', pad_inches=0.1) plt.show() -# Save a PDF/PNG of the figure -# fig.savefig(basepath + '.pdf', dpi=None, format='pdf', bbox_inches='tight', pad_inches=0.1) -# fig.savefig(savename + '.png', dpi=150, format='png', bbox_inches='tight', pad_inches=0.1) + From 8dbb673d679781f0b1e22d07b3ae0d754d0bcb8c Mon Sep 17 00:00:00 2001 From: Craig Warren Date: Fri, 1 Sep 2023 14:49:01 +0100 Subject: [PATCH 12/35] Ensure blank scenes list created even from CLI --- gprMax/gprMax.py | 60 +++++++++++++++++++++++++----------------------- 1 file changed, 31 insertions(+), 29 deletions(-) diff --git a/gprMax/gprMax.py b/gprMax/gprMax.py index ba61add0..aecf22b0 100644 --- a/gprMax/gprMax.py +++ b/gprMax/gprMax.py @@ -25,7 +25,7 @@ from .utilities.logging import logging_config # Arguments (used for API) and their default values (used for API and CLI) args_defaults = { - "scenes": None, + "scenes": [], "inputfile": None, "outputfile": None, "n": 1, @@ -44,37 +44,39 @@ args_defaults = { # Argument help messages (used for CLI argparse) help_msg = { - "scenes": "(list, req): Scenes to run the model. " - "Multiple scene objects can given in order to run multiple " - "simulation runs. Each scene must contain the essential " - "simulation objects", - "inputfile": "(str, opt): Input file path. Can also run simulation " "by providing an input file.", + "scenes": "(list, req): Scenes to run the model. Multiple scene objects " + "can given in order to run multiple simulation runs. Each scene " + "must contain the essential simulation objects", + "inputfile": "(str, opt): Input file path. Can also run simulation by " + "providing an input file.", "outputfile": "(str, req): File path to the output data file.", "n": "(int, req): Number of required simulation runs.", - "i": "(int, opt): Model number to start/restart simulation " - "from. It would typically be used to restart a series of " - "models from a specific model number, with the n argument, " - "e.g. to restart from A-scan 45 when creating a B-scan " - "with 60 traces.", - "mpi": "(bool, opt): Flag to use Message Passing Interface (MPI) " - "task farm. This option is most usefully combined with n to " - "allow individual models to be farmed out using a MPI task " - "farm, e.g. to create a B-scan with 60 traces and use MPI to " - "farm out each trace. For further details see the parallel " - "performance section of the User Guide.", - "gpu": "(list/bool, opt): Flag to use NVIDIA GPU or list of NVIDIA " "GPU device ID(s) for specific GPU card(s).", - "opencl": "(list/bool, opt): Flag to use OpenCL or list of OpenCL " "device ID(s) for specific compute device(s).", + "i": "(int, opt): Model number to start/restart simulation from. It would " + "typically be used to restart a series of models from a specific " + "model number, with the n argument, e.g. to restart from A-scan 45 " + "when creating a B-scan with 60 traces.", + "mpi": "(bool, opt): Flag to use Message Passing Interface (MPI) task farm. " + "This option is most usefully combined with n to allow individual " + "models to be farmed out using a MPI task farm, e.g. to create a " + "B-scan with 60 traces and use MPI to farm out each trace. For " + "further details see the performance section of the User Guide.", + "gpu": "(list/bool, opt): Flag to use NVIDIA GPU or list of NVIDIA GPU " + "device ID(s) for specific GPU card(s).", + "opencl": "(list/bool, opt): Flag to use OpenCL or list of OpenCL device " + "ID(s) for specific compute device(s).", "subgrid": "(bool, opt): Flag to use sub-gridding.", - "autotranslate": "(bool, opt): For sub-gridding - auto translate " - "objects with main grid coordinates to their " - "equivalent local grid coordinate within the " - "subgrid. If this option is off users must specify " - "sub-grid object point within the global subgrid space.", - "geometry_only": "(bool, opt): Build a model and produce any " "geometry views but do not run the simulation.", - "geometry_fixed": "(bool, opt): Run a series of models where the " "geometry does not change between models.", - "write_processed": "(bool, opt): Writes another input file after " - "any Python code (#python blocks) and in the " - "original input file has been processed.", + "autotranslate": "(bool, opt): For sub-gridding - auto translate objects " + "with main grid coordinates to their equivalent local " + "grid coordinate within the subgrid. If this option is " + "off users must specify sub-grid object point within the " + "global subgrid space.", + "geometry_only": "(bool, opt): Build a model and produce any geometry " + "views but do not run the simulation.", + "geometry_fixed": "(bool, opt): Run a series of models where the geometry " + "does not change between models.", + "write_processed": "(bool, opt): Writes another input file after any " + "Python code (#python blocks) and in the original input " + "file has been processed.", "log_level": "(int, opt): Level of logging to use.", "log_file": "(bool, opt): Write logging information to file.", } From 23bb2633f8aa5123e35d51e76b7fd99dac3cccb7 Mon Sep 17 00:00:00 2001 From: Craig Warren Date: Fri, 1 Sep 2023 16:22:59 +0100 Subject: [PATCH 13/35] Formatting cleanup --- gprMax/gprMax.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/gprMax/gprMax.py b/gprMax/gprMax.py index aecf22b0..d3b31521 100644 --- a/gprMax/gprMax.py +++ b/gprMax/gprMax.py @@ -211,9 +211,7 @@ def run_main(args): """ results = {} - logging_config(level=args.log_level, log_file=args.log_file) - config.sim_config = config.SimulationConfig(args) # MPI running with (OpenMP/CUDA/OpenCL) @@ -224,4 +222,5 @@ def run_main(args): context = Context() results = context.run() + return results From 4dbc5c4d6a3d32565083274d82ccfc9ba521827a Mon Sep 17 00:00:00 2001 From: Craig Warren Date: Fri, 1 Sep 2023 16:23:14 +0100 Subject: [PATCH 14/35] Corrected path for plotting --- docs/source/examples_simple_2D.rst | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/docs/source/examples_simple_2D.rst b/docs/source/examples_simple_2D.rst index 0ed67e0d..c42ed962 100644 --- a/docs/source/examples_simple_2D.rst +++ b/docs/source/examples_simple_2D.rst @@ -112,7 +112,7 @@ You should have produced an output file ``cylinder_Ascan_2D.h5``. You can view t .. code-block:: none - python -m tools.plot_Ascan examples/cylinder_Ascan_2D.h5 + python -m toolboxes.Plotting.plot_Ascan examples/cylinder_Ascan_2D.h5 :numref:`cylinder_Ascan_results` shows the time history of the electric and magnetic field components and currents at the receiver location. The :math:`E_z` field component can be converted to a voltage that represents the A-scan (trace). The initial part of the signal (~0.5-1.5 ns) represents the direct wave from transmitter to receiver. Then comes the reflected wavelet (~1.8-2.6 ns), which has opposite polarity, from the metal cylinder. @@ -153,7 +153,7 @@ You should have produced 60 output files, one for each A-scan, with names ``cyli .. code-block:: none - python -m tools.outputfiles_merge examples/cylinder_Bscan_2D + python -m toolboxes.Utilities.outputfiles_merge examples/cylinder_Bscan_2D You should see a combined output file ``cylinder_Bscan_2D_merged.h5``. You can add the optional argument ``--remove-files`` if you want to automatically delete the original single A-scan output files. @@ -161,7 +161,7 @@ You can now view an image of the B-scan using the command: .. code-block:: none - python -m tools.plot_Bscan examples/cylinder_Bscan_2D_merged.h5 Ez + python -m toolboxes.Plotting.plot_Bscan examples/cylinder_Bscan_2D_merged.h5 Ez :numref:`cylinder_Bscan_results` shows the B-scan (of the :math:`E_z` field component). Again, the initial part of the signal (~0.5-1.5 ns) represents the direct wave from transmitter to receiver. Then comes the reflected wave (~2-3 ns) from the metal cylinder which creates the hyperbolic shape. From e3890bc2c58c22ba5887df1bd5a49e19003a0b96 Mon Sep 17 00:00:00 2001 From: Craig Warren Date: Fri, 1 Sep 2023 16:23:34 +0100 Subject: [PATCH 15/35] Constants import from scipy --- testing/analytical_solutions.py | 19 ++++++++++--------- 1 file changed, 10 insertions(+), 9 deletions(-) diff --git a/testing/analytical_solutions.py b/testing/analytical_solutions.py index 76656dd0..877826d5 100644 --- a/testing/analytical_solutions.py +++ b/testing/analytical_solutions.py @@ -17,8 +17,9 @@ # along with gprMax. If not, see . import numpy as np +from scipy.constants import c +from scipy.constants import epsilon_0 as e0 -import gprMax.config as config from gprMax.waveforms import Waveform @@ -78,35 +79,35 @@ def hertzian_dipole_fs(iterations, dt, dxdydz, rx): Ex_y = y Ex_z = z - 0.5 * dz Er_x = np.sqrt((Ex_x**2 + Ex_y**2 + Ex_z**2)) - tau_Ex = Er_x / config.sim_config.em_consts["c"] + tau_Ex = Er_x / c # Coordinates of Rx for Ey FDTD component Ey_x = x Ey_y = y + 0.5 * dy Ey_z = z - 0.5 * dz Er_y = np.sqrt((Ey_x**2 + Ey_y**2 + Ey_z**2)) - tau_Ey = Er_y / config.sim_config.em_consts["c"] + tau_Ey = Er_y / c # Coordinates of Rx for Ez FDTD component Ez_x = x Ez_y = y Ez_z = z Er_z = np.sqrt((Ez_x**2 + Ez_y**2 + Ez_z**2)) - tau_Ez = Er_z / config.sim_config.em_consts["c"] + tau_Ez = Er_z / c # Coordinates of Rx for Hx FDTD component Hx_x = x Hx_y = y + 0.5 * dy Hx_z = z Hr_x = np.sqrt((Hx_x**2 + Hx_y**2 + Hx_z**2)) - tau_Hx = Hr_x / config.sim_config.em_consts["c"] + tau_Hx = Hr_x / c # Coordinates of Rx for Hy FDTD component Hy_x = x + 0.5 * dx Hy_y = y Hy_z = z Hr_y = np.sqrt((Hy_x**2 + Hy_y**2 + Hy_z**2)) - tau_Hy = Hr_y / config.sim_config.em_consts["c"] + tau_Hy = Hr_y / c # Initialise fields fields = np.zeros((iterations, 6)) @@ -134,7 +135,7 @@ def hertzian_dipole_fs(iterations, dt, dxdydz, rx): fdot_Hy = wdot.calculate_value((timestep * dt) - tau_Hy, dt) * dl # Ex - fields[timestep, 0] = ((Ex_x * Ex_z) / (4 * np.pi * config.sim_config.em_consts["e0"] * Er_x**5)) * ( + fields[timestep, 0] = ((Ex_x * Ex_z) / (4 * np.pi * e0 * Er_x**5)) * ( 3 * (fint_Ex + (tau_Ex * f_Ex)) + (tau_Ex**2 * fdot_Ex) ) @@ -145,12 +146,12 @@ def hertzian_dipole_fs(iterations, dt, dxdydz, rx): tmp = 0 fields[timestep, 1] = ( tmp - * ((Ey_x * Ey_z) / (4 * np.pi * config.sim_config.em_consts["e0"] * Er_y**5)) + * ((Ey_x * Ey_z) / (4 * np.pi * e0 * Er_y**5)) * (3 * (fint_Ey + (tau_Ey * f_Ey)) + (tau_Ey**2 * fdot_Ey)) ) # Ez - fields[timestep, 2] = (1 / (4 * np.pi * config.sim_config.em_consts["e0"] * Er_z**5)) * ( + fields[timestep, 2] = (1 / (4 * np.pi * e0 * Er_z**5)) * ( (2 * Ez_z**2 - (Ez_x**2 + Ez_y**2)) * (fint_Ez + (tau_Ez * f_Ez)) - (Ez_x**2 + Ez_y**2) * tau_Ez**2 * fdot_Ez ) From 7230a8bf1f520338df225d2c677625909769fb9f Mon Sep 17 00:00:00 2001 From: Craig Warren Date: Fri, 1 Sep 2023 20:43:20 +0100 Subject: [PATCH 16/35] Cleaned up function arg formatting for CUDA --- .../knl_pml_updates_electric_HORIPML.py | 36 +++++++++++++++---- .../knl_pml_updates_electric_MRIPML.py | 36 +++++++++++++++---- .../knl_pml_updates_magnetic_HORIPML.py | 32 +++++++++++++---- .../knl_pml_updates_magnetic_MRIPML.py | 32 +++++++++++++---- 4 files changed, 110 insertions(+), 26 deletions(-) diff --git a/gprMax/cuda_opencl/knl_pml_updates_electric_HORIPML.py b/gprMax/cuda_opencl/knl_pml_updates_electric_HORIPML.py index 8974a3a6..8d6a357e 100644 --- a/gprMax/cuda_opencl/knl_pml_updates_electric_HORIPML.py +++ b/gprMax/cuda_opencl/knl_pml_updates_electric_HORIPML.py @@ -35,11 +35,19 @@ x_args = { int NZ_PHI2, int NY_R, const unsigned int* __restrict__ ID, - const $REAL* __restrict__ Ex, $REAL *Ey, + const $REAL* __restrict__ Ex, + $REAL *Ey, $REAL *Ez, - const $REAL* __restrict__ Hx, const $REAL* __restrict__ Hy, const $REAL* __restrict__ Hz, $REAL *PHI1, + const $REAL* __restrict__ Hx, + const $REAL* __restrict__ Hy, + const $REAL* __restrict__ Hz, + $REAL *PHI1, $REAL *PHI2, - const $REAL* __restrict__ RA, const $REAL* __restrict__ RB, const $REAL* __restrict__ RE, const $REAL* __restrict__ RF, $REAL d) + const $REAL* __restrict__ RA, + const $REAL* __restrict__ RB, + const $REAL* __restrict__ RE, + const $REAL* __restrict__ RF, + $REAL d) """ ), "opencl": Template( @@ -95,9 +103,16 @@ y_args = { $REAL *Ex, const $REAL* __restrict__ Ey, $REAL *Ez, - const $REAL* __restrict__ Hx, const $REAL* __restrict__ Hy, const $REAL* __restrict__ Hz, $REAL *PHI1, + const $REAL* __restrict__ Hx, + const $REAL* __restrict__ Hy, + const $REAL* __restrict__ Hz, + $REAL *PHI1, $REAL *PHI2, - const $REAL* __restrict__ RA, const $REAL* __restrict__ RB, const $REAL* __restrict__ RE, const $REAL* __restrict__ RF, $REAL d) + const $REAL* __restrict__ RA, + const $REAL* __restrict__ RB, + const $REAL* __restrict__ RE, + const $REAL* __restrict__ RF, + $REAL d) """ ), "opencl": Template( @@ -153,9 +168,16 @@ z_args = { $REAL *Ex, $REAL *Ey, const $REAL* __restrict__ Ez, - const $REAL* __restrict__ Hx, const $REAL* __restrict__ Hy, const $REAL* __restrict__ Hz, $REAL *PHI1, + const $REAL* __restrict__ Hx, + const $REAL* __restrict__ Hy, + const $REAL* __restrict__ Hz, + $REAL *PHI1, $REAL *PHI2, - const $REAL* __restrict__ RA, const $REAL* __restrict__ RB, const $REAL* __restrict__ RE, const $REAL* __restrict__ RF, $REAL d) + const $REAL* __restrict__ RA, + const $REAL* __restrict__ RB, + const $REAL* __restrict__ RE, + const $REAL* __restrict__ RF, + $REAL d) """ ), "opencl": Template( diff --git a/gprMax/cuda_opencl/knl_pml_updates_electric_MRIPML.py b/gprMax/cuda_opencl/knl_pml_updates_electric_MRIPML.py index 1effe788..33adee2d 100644 --- a/gprMax/cuda_opencl/knl_pml_updates_electric_MRIPML.py +++ b/gprMax/cuda_opencl/knl_pml_updates_electric_MRIPML.py @@ -35,11 +35,19 @@ x_args = { int NZ_PHI2, int NY_R, const unsigned int* __restrict__ ID, - const $REAL* __restrict__ Ex, $REAL *Ey, + const $REAL* __restrict__ Ex, + $REAL *Ey, $REAL *Ez, - const $REAL* __restrict__ Hx, const $REAL* __restrict__ Hy, const $REAL* __restrict__ Hz, $REAL *PHI1, + const $REAL* __restrict__ Hx, + const $REAL* __restrict__ Hy, + const $REAL* __restrict__ Hz, + $REAL *PHI1, $REAL *PHI2, - const $REAL* __restrict__ RA, const $REAL* __restrict__ RB, const $REAL* __restrict__ RE, const $REAL* __restrict__ RF, $REAL d) + const $REAL* __restrict__ RA, + const $REAL* __restrict__ RB, + const $REAL* __restrict__ RE, + const $REAL* __restrict__ RF, + $REAL d) """ ), "opencl": Template( @@ -95,9 +103,16 @@ y_args = { $REAL *Ex, const $REAL* __restrict__ Ey, $REAL *Ez, - const $REAL* __restrict__ Hx, const $REAL* __restrict__ Hy, const $REAL* __restrict__ Hz, $REAL *PHI1, + const $REAL* __restrict__ Hx, + const $REAL* __restrict__ Hy, + const $REAL* __restrict__ Hz, + $REAL *PHI1, $REAL *PHI2, - const $REAL* __restrict__ RA, const $REAL* __restrict__ RB, const $REAL* __restrict__ RE, const $REAL* __restrict__ RF, $REAL d) + const $REAL* __restrict__ RA, + const $REAL* __restrict__ RB, + const $REAL* __restrict__ RE, + const $REAL* __restrict__ RF, + $REAL d) """ ), "opencl": Template( @@ -153,9 +168,16 @@ z_args = { $REAL *Ex, $REAL *Ey, const $REAL* __restrict__ Ez, - const $REAL* __restrict__ Hx, const $REAL* __restrict__ Hy, const $REAL* __restrict__ Hz, $REAL *PHI1, + const $REAL* __restrict__ Hx, + const $REAL* __restrict__ Hy, + const $REAL* __restrict__ Hz, + $REAL *PHI1, $REAL *PHI2, - const $REAL* __restrict__ RA, const $REAL* __restrict__ RB, const $REAL* __restrict__ RE, const $REAL* __restrict__ RF, $REAL d) + const $REAL* __restrict__ RA, + const $REAL* __restrict__ RB, + const $REAL* __restrict__ RE, + const $REAL* __restrict__ RF, + $REAL d) """ ), "opencl": Template( diff --git a/gprMax/cuda_opencl/knl_pml_updates_magnetic_HORIPML.py b/gprMax/cuda_opencl/knl_pml_updates_magnetic_HORIPML.py index 3f8bc823..6fccb4f6 100644 --- a/gprMax/cuda_opencl/knl_pml_updates_magnetic_HORIPML.py +++ b/gprMax/cuda_opencl/knl_pml_updates_magnetic_HORIPML.py @@ -35,11 +35,19 @@ x_args = { int NZ_PHI2, int NY_R, const unsigned int* __restrict__ ID, - const $REAL* __restrict__ Ex, const $REAL* __restrict__ Ey, const $REAL* __restrict__ Ez, const $REAL* __restrict__ Hx, $REAL *Hy, + const $REAL* __restrict__ Ex, + const $REAL* __restrict__ Ey, + const $REAL* __restrict__ Ez, + const $REAL* __restrict__ Hx, + $REAL *Hy, $REAL *Hz, $REAL *PHI1, $REAL *PHI2, - const $REAL* __restrict__ RA, const $REAL* __restrict__ RB, const $REAL* __restrict__ RE, const $REAL* __restrict__ RF, $REAL d) + const $REAL* __restrict__ RA, + const $REAL* __restrict__ RB, + const $REAL* __restrict__ RE, + const $REAL* __restrict__ RF, + $REAL d) """ ), "opencl": Template( @@ -92,13 +100,19 @@ y_args = { int NZ_PHI2, int NY_R, const unsigned int* __restrict__ ID, - const $REAL* __restrict__ Ex, const $REAL* __restrict__ Ey, const $REAL* __restrict__ Ez, + const $REAL* __restrict__ Ex, + const $REAL* __restrict__ Ey, + const $REAL* __restrict__ Ez, $REAL *Hx, const $REAL* __restrict__ Hy, $REAL *Hz, $REAL *PHI1, $REAL *PHI2, - const $REAL* __restrict__ RA, const $REAL* __restrict__ RB, const $REAL* __restrict__ RE, const $REAL* __restrict__ RF, $REAL d) + const $REAL* __restrict__ RA, + const $REAL* __restrict__ RB, + const $REAL* __restrict__ RE, + const $REAL* __restrict__ RF, + $REAL d) """ ), "opencl": Template( @@ -151,13 +165,19 @@ z_args = { int NZ_PHI2, int NY_R, const unsigned int* __restrict__ ID, - const $REAL* __restrict__ Ex, const $REAL* __restrict__ Ey, const $REAL* __restrict__ Ez, + const $REAL* __restrict__ Ex, + const $REAL* __restrict__ Ey, + const $REAL* __restrict__ Ez, $REAL *Hx, $REAL *Hy, const $REAL* __restrict__ Hz, $REAL *PHI1, $REAL *PHI2, - const $REAL* __restrict__ RA, const $REAL* __restrict__ RB, const $REAL* __restrict__ RE, const $REAL* __restrict__ RF, $REAL d) + const $REAL* __restrict__ RA, + const $REAL* __restrict__ RB, + const $REAL* __restrict__ RE, + const $REAL* __restrict__ RF, + $REAL d) """ ), "opencl": Template( diff --git a/gprMax/cuda_opencl/knl_pml_updates_magnetic_MRIPML.py b/gprMax/cuda_opencl/knl_pml_updates_magnetic_MRIPML.py index 05ca8b9e..89df47ae 100644 --- a/gprMax/cuda_opencl/knl_pml_updates_magnetic_MRIPML.py +++ b/gprMax/cuda_opencl/knl_pml_updates_magnetic_MRIPML.py @@ -35,11 +35,19 @@ x_args = { int NZ_PHI2, int NY_R, const unsigned int* __restrict__ ID, - const $REAL* __restrict__ Ex, const $REAL* __restrict__ Ey, const $REAL* __restrict__ Ez, const $REAL* __restrict__ Hx, $REAL *Hy, + const $REAL* __restrict__ Ex, + const $REAL* __restrict__ Ey, + const $REAL* __restrict__ Ez, + const $REAL* __restrict__ Hx, + $REAL *Hy, $REAL *Hz, $REAL *PHI1, $REAL *PHI2, - const $REAL* __restrict__ RA, const $REAL* __restrict__ RB, const $REAL* __restrict__ RE, const $REAL* __restrict__ RF, $REAL d) + const $REAL* __restrict__ RA, + const $REAL* __restrict__ RB, + const $REAL* __restrict__ RE, + const $REAL* __restrict__ RF, + $REAL d) """ ), "opencl": Template( @@ -92,13 +100,19 @@ y_args = { int NZ_PHI2, int NY_R, const unsigned int* __restrict__ ID, - const $REAL* __restrict__ Ex, const $REAL* __restrict__ Ey, const $REAL* __restrict__ Ez, + const $REAL* __restrict__ Ex, + const $REAL* __restrict__ Ey, + const $REAL* __restrict__ Ez, $REAL *Hx, const $REAL* __restrict__ Hy, $REAL *Hz, $REAL *PHI1, $REAL *PHI2, - const $REAL* __restrict__ RA, const $REAL* __restrict__ RB, const $REAL* __restrict__ RE, const $REAL* __restrict__ RF, $REAL d) + const $REAL* __restrict__ RA, + const $REAL* __restrict__ RB, + const $REAL* __restrict__ RE, + const $REAL* __restrict__ RF, + $REAL d) """ ), "opencl": Template( @@ -151,13 +165,19 @@ z_args = { int NZ_PHI2, int NY_R, const unsigned int* __restrict__ ID, - const $REAL* __restrict__ Ex, const $REAL* __restrict__ Ey, const $REAL* __restrict__ Ez, + const $REAL* __restrict__ Ex, + const $REAL* __restrict__ Ey, + const $REAL* __restrict__ Ez, $REAL *Hx, $REAL *Hy, const $REAL* __restrict__ Hz, $REAL *PHI1, $REAL *PHI2, - const $REAL* __restrict__ RA, const $REAL* __restrict__ RB, const $REAL* __restrict__ RE, const $REAL* __restrict__ RF, $REAL d) + const $REAL* __restrict__ RA, + const $REAL* __restrict__ RB, + const $REAL* __restrict__ RE, + const $REAL* __restrict__ RF, + $REAL d) """ ), "opencl": Template( From b3e8b983e0fb8a0b60eed20233fdd7299b0dfc54 Mon Sep 17 00:00:00 2001 From: Craig Warren Date: Fri, 1 Sep 2023 21:23:34 +0100 Subject: [PATCH 17/35] Cleaned up formatting and fixed an issue with double being set incorrectly --- gprMax/config.py | 21 +++++++-------------- 1 file changed, 7 insertions(+), 14 deletions(-) diff --git a/gprMax/config.py b/gprMax/config.py index fbbb7628..8d589a56 100644 --- a/gprMax/config.py +++ b/gprMax/config.py @@ -209,33 +209,29 @@ class SimulationConfig: self.args = args - if args.mpi and args.geometry_fixed: + if self.args.mpi and self.args.geometry_fixed: logger.exception("The geometry fixed option cannot be used with MPI.") raise ValueError - if args.gpu and args.opencl: + if self.args.gpu and self.args.opencl: logger.exception("You cannot use both CUDA and OpenCl simultaneously.") raise ValueError # General settings for the simulation - # inputfilepath: path to inputfile location. - # outputfilepath: path to outputfile location. - # progressbars: whether to show progress bars on stdoout or not. # solver: cpu, cuda, opencl. - # subgrid: whether the simulation uses sub-grids. # precision: data type for electromagnetic field output (single/double). # progressbars: progress bars on stdoout or not - switch off # progressbars when logging level is greater than # info (20) - self.general = {"solver": "cpu", "subgrid": False, "precision": "single", "progressbars": args.log_level <= 20} + self.general = {"solver": "cpu", "precision": "single", "progressbars": args.log_level <= 20} self.em_consts = { "c": c, # Speed of light in free space (m/s) "e0": e0, # Permittivity of free space (F/m) "m0": m0, # Permeability of free space (H/m) - "z0": np.sqrt(m0 / e0), - } # Impedance of free space (Ohms) + "z0": np.sqrt(m0 / e0), # Impedance of free space (Ohms) + } # Store information about host machine self.hostinfo = get_host_info() @@ -266,9 +262,8 @@ class SimulationConfig: # Add pyopencl available device(s) self.devices["devs"] = detect_opencl() - # Subgrid parameter may not exist if user enters via CLI - try: - self.general["subgrid"] = self.args.subgrid + self.general["subgrid"] = self.args.subgrid + if self.args.subgrid: # Double precision should be used with subgrid for best accuracy self.general["precision"] = "double" if (self.general["subgrid"] and self.general["solver"] == "cuda") or ( @@ -278,8 +273,6 @@ class SimulationConfig: "You cannot currently use CUDA or OpenCL-based " "solvers with models that contain sub-grids." ) raise ValueError - except AttributeError: - self.general["subgrid"] = False # Scenes parameter may not exist if user enters via CLI try: From 3e17dc649961f51b2fd9177cc4ddc66b0b4a73bd Mon Sep 17 00:00:00 2001 From: Craig Warren Date: Mon, 4 Sep 2023 10:49:20 +0100 Subject: [PATCH 18/35] Add new geo primitives --- gprMax/__init__.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/gprMax/__init__.py b/gprMax/__init__.py index 7c48639f..4e572e47 100644 --- a/gprMax/__init__.py +++ b/gprMax/__init__.py @@ -14,9 +14,11 @@ from .cmds_geometry.add_grass import AddGrass from .cmds_geometry.add_surface_roughness import AddSurfaceRoughness from .cmds_geometry.add_surface_water import AddSurfaceWater from .cmds_geometry.box import Box +from .cmds_geometry.cone import Cone from .cmds_geometry.cylinder import Cylinder from .cmds_geometry.cylindrical_sector import CylindricalSector from .cmds_geometry.edge import Edge +from .cmds_geometry.ellipsoid import Ellipsoid from .cmds_geometry.fractal_box import FractalBox from .cmds_geometry.geometry_objects_read import GeometryObjectsRead from .cmds_geometry.plate import Plate @@ -38,6 +40,7 @@ from .cmds_multiuse import ( RxArray, Snapshot, SoilPeplinski, + Subgrid, TransmissionLine, VoltageSource, Waveform, @@ -47,6 +50,7 @@ from .cmds_singleuse import ( Domain, ExcitationFile, OMPThreads, + OutputDir, PMLProps, RxSteps, SrcSteps, From f52e1d63cffa697b03bc5fbc13c5b610d53c077b Mon Sep 17 00:00:00 2001 From: Craig Warren Date: Mon, 4 Sep 2023 10:49:58 +0100 Subject: [PATCH 19/35] Reinstate PML CFS and formulation via hash cmds --- docs/source/input_hash_cmds.rst | 46 ++++++++++++++++++++++++++++++++- gprMax/cmds_multiuse.py | 7 +---- gprMax/hash_cmds_file.py | 2 ++ gprMax/hash_cmds_multiuse.py | 27 ++++++++++++++++++- gprMax/hash_cmds_singleuse.py | 19 ++++++++++++-- 5 files changed, 91 insertions(+), 10 deletions(-) diff --git a/docs/source/input_hash_cmds.rst b/docs/source/input_hash_cmds.rst index c4710f32..e6ebf26a 100644 --- a/docs/source/input_hash_cmds.rst +++ b/docs/source/input_hash_cmds.rst @@ -965,7 +965,7 @@ For example to save a snapshot of the electromagnetic fields in the model at a s PML commands ============ -The default behaviour for the absorbing boundary conditions (ABC) is first order Complex Frequency Shifted (CFS) Perfectly Matched Layers (PML), with thicknesses of 10 cells on each of the six sides of the model domain. The thickness of the PML can be altered by using the following command (further customisation of the PML is possible using our Python API): +The default behaviour for the absorbing boundary conditions (ABC) is first order Complex Frequency Shifted (CFS) Perfectly Matched Layers (PML), with thicknesses of 10 cells on each of the six sides of the model domain. #pml_cells: ------------ @@ -989,3 +989,47 @@ For example to use a PML with 20 cells (thicker than the default 10 cells) on on .. code-block:: none #pml_cells: 10 10 20 10 10 20 + +#pml_formulation: +----------------- + +Allows you to alter the formulation used for the PML. The current options are to use the Higher Order RIPML (HORIPML) - https://doi.org/10.1109/TAP.2011.2180344, or Multipole RIPML (MRIPML) - https://doi.org/10.1109/TAP.2018.2823864. The syntax of the command is: + +.. code-block:: none + + #pml_formulation: str + +* ``str`` can be either 'HORIPML' or 'MRIPML' + +For example to use the Multipole RIPML: + +.. code-block:: none + + #pml_formulation: MRIPML + +#pml_cfs: +--------- + +Allows you (advanced) control of the parameters that are used to build each order of the PML. Up to a second order PML can currently be specified, i.e. by using two ``#pml_cfs`` commands. The syntax of the command is: + +.. code-block:: none + + #pml_cfs: str1 str2 f1 f2 str3 str4 f3 f4 str5 str6 f5 f6 + +* ``str1`` is the type of scaling to use for the CFS :math:`\alpha` parameter. It can be ``constant``, ``linear``, ``quadratic``, ``cubic``, ``quartic``, ``quintic`` and ``sextic``. +* ``str2`` is the direction of the scaling to use for the CFS :math:`\alpha` parameter. It can be ``forward`` or ``reverse``. +* ``f1 f2`` are the minimum and maximum values for the CFS :math:`\alpha` parameter. +* ``str3`` is the type of scaling to use for the CFS :math:`\kappa` parameter. It can be ``constant``, ``linear``, ``quadratic``, ``cubic``, ``quartic``, ``quintic`` and ``sextic``. +* ``str4`` is the direction of the scaling to use for the CFS :math:`\kappa` parameter. It can be ``forward`` or ``reverse``. +* ``f3 f4`` are the minimum and maximum values for the CFS :math:`\kappa` parameter. The minimum value for the CFS :math:`\kappa` parameter is one. +* ``str5`` is the type of scaling to use for the CFS :math:`\sigma` parameter. It can be ``constant``, ``linear``, ``quadratic``, ``cubic``, ``quartic``, ``quintic`` and ``sextic``. +* ``str6`` is the direction of the scaling to use for the CFS :math:`\sigma` parameter. It can be ``forward`` or ``reverse``. +* ``f5 f6`` are the minimum and maximum values for the CFS :math:`\sigma` parameter. + +The CFS values (which are internally specified) used for the default standard first order PML are: ``#pml_cfs: constant forward 0 0 constant forward 1 1 quartic forward 0 None``. Specifying 'None' for the maximum value of :math:`\sigma` forces gprMax to calculate it internally based on the relative permittivity and permeability of the underlying materials in the model. + +The parameters will be applied to all slabs of the PML that are switched on. + +.. tip:: + + ``forward`` direction implies minimum parameter value at the inner boundary of the PML and maximum parameter value at the edge of computational domain, ``reverse`` is the opposite. diff --git a/gprMax/cmds_multiuse.py b/gprMax/cmds_multiuse.py index 7d3a4868..50cee93c 100644 --- a/gprMax/cmds_multiuse.py +++ b/gprMax/cmds_multiuse.py @@ -1737,11 +1737,6 @@ class PMLCFS(UserObjectMulti): ): logger.exception(f"{self.params_str()} minimum and maximum scaling values must be greater than zero.") raise ValueError - # TODO: Fix handling of kappa for 2nd order PMLs - # if float(kappamin) < 1: - # logger.exception(f'{self.params_str()} minimum scaling value for ' - # 'kappa must be greater than or equal to one.') - # raise ValueError cfsalpha = CFSParameter() cfsalpha.ID = "alpha" @@ -1800,4 +1795,4 @@ class Subgrid(UserObjectMulti): self.children_geometry.append(node) else: logger.exception("This object is unknown to gprMax.") - raise ValueError + raise ValueError diff --git a/gprMax/hash_cmds_file.py b/gprMax/hash_cmds_file.py index b5ba0bcb..90539a92 100644 --- a/gprMax/hash_cmds_file.py +++ b/gprMax/hash_cmds_file.py @@ -217,6 +217,7 @@ def check_cmd_names(processedlines, checkessential=True): "#title", "#omp_threads", "#time_step_stability_factor", + "#pml_formulation", "#pml_cells", "#excitation_file", "#src_steps", @@ -248,6 +249,7 @@ def check_cmd_names(processedlines, checkessential=True): "#rx", "#rx_array", "#snapshot", + "#pml_cfs", "#include_file", ] } diff --git a/gprMax/hash_cmds_multiuse.py b/gprMax/hash_cmds_multiuse.py index 18e71c06..f3e1f5ff 100644 --- a/gprMax/hash_cmds_multiuse.py +++ b/gprMax/hash_cmds_multiuse.py @@ -29,6 +29,7 @@ from .cmds_multiuse import ( Material, MaterialList, MaterialRange, + PMLCFS, Rx, RxArray, Snapshot, @@ -375,7 +376,7 @@ def process_multicmds(multicmds): tmp = cmdinstance.split() if len(tmp) < 2: - logger.exception("'" + cmdname + ": " + " ".join(tmp) + "'" + " requires at least 2 parameters") + logger.exception("'" + cmdname + ": " + " ".join(tmp) + "'" + " requires at least two parameters") raise ValueError tokens = len(tmp) @@ -386,4 +387,28 @@ def process_multicmds(multicmds): material_list = MaterialList(list_of_materials=lmats, id=tmp[tokens - 1]) scene_objects.append(material_list) + cmdname = '#pml_cfs' + if multicmds[cmdname] is not None: + for cmdinstance in multicmds[cmdname]: + tmp = cmdinstance.split() + + if len(tmp) != 12: + logger.exception("'" + cmdname + ": " + " ".join(tmp) + "'" + " requires exactly twelve parameters") + raise ValueError + + pml_cfs = PMLCFS(alphascalingprofile=tmp[0], + alphascalingdirection=tmp[1], + alphamin=tmp[2], + alphamax=tmp[3], + kappascalingprofile=tmp[4], + kappascalingdirection=tmp[5], + kappamin=tmp[6], + kappamax=tmp[7], + sigmascalingprofile=tmp[8], + sigmascalingdirection=tmp[9], + sigmamin=tmp[10], + sigmamax=tmp[11]) + + scene_objects.append(pml_cfs) + return scene_objects diff --git a/gprMax/hash_cmds_singleuse.py b/gprMax/hash_cmds_singleuse.py index 023e78a1..6dd3765b 100644 --- a/gprMax/hash_cmds_singleuse.py +++ b/gprMax/hash_cmds_singleuse.py @@ -120,17 +120,32 @@ def process_singlecmds(singlecmds): scene_objects.append(tw) + cmd = "#pml_formulation" + if singlecmds[cmd] is not None: + tmp = singlecmds[cmd].split() + if len(tmp) != 1: + logger.exception(f"{cmd} requires one parameter") + raise ValueError + else: + pml_formulation = tmp[0] + cmd = "#pml_cells" if singlecmds[cmd] is not None: tmp = singlecmds[cmd].split() if len(tmp) not in [1, 6]: logger.exception(f"{cmd} requires either one or six parameter(s)") raise ValueError + if pml_formulation: + formulation = pml_formulation + else: + formulation = None if len(tmp) == 1: - pml_cells = PMLProps(thickness=int(tmp[0])) + pml_cells = PMLProps(formulation=formulation, thickness=int(tmp[0])) else: pml_cells = PMLProps( - x0=int(tmp[0]), y0=int(tmp[1]), z0=int(tmp[2]), xmax=int(tmp[3]), ymax=int(tmp[4]), zmax=int(tmp[5]) + formulation=formulation, + x0=int(tmp[0]), y0=int(tmp[1]), z0=int(tmp[2]), + xmax=int(tmp[3]), ymax=int(tmp[4]), zmax=int(tmp[5]) ) scene_objects.append(pml_cells) From ce2ac508a41a30fc41f6b37473929f61bd04fa8e Mon Sep 17 00:00:00 2001 From: Craig Warren Date: Mon, 4 Sep 2023 10:54:32 +0100 Subject: [PATCH 20/35] Clarified some comments --- gprMax/grid.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/gprMax/grid.py b/gprMax/grid.py index bb0ab22e..c6610190 100644 --- a/gprMax/grid.py +++ b/gprMax/grid.py @@ -56,10 +56,11 @@ class FDTDGrid: self.pmls["formulation"] = "HORIPML" self.pmls["cfs"] = [] self.pmls["slabs"] = [] - # 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 *updating* the PMLs always follows the + # same order (the order for *building* PMLs does not matter). 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.pmls["thickness"] = OrderedDict((key, 10) for key in PML.boundaryIDs) self.materials = [] From 2e891861a2b514d0f6b19e4c6e700fbdef76cf20 Mon Sep 17 00:00:00 2001 From: Craig Warren Date: Mon, 4 Sep 2023 11:03:02 +0100 Subject: [PATCH 21/35] Fix subgrid arg for CLI --- gprMax/config.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/gprMax/config.py b/gprMax/config.py index 8d589a56..f9b1dcf8 100644 --- a/gprMax/config.py +++ b/gprMax/config.py @@ -262,8 +262,9 @@ class SimulationConfig: # Add pyopencl available device(s) self.devices["devs"] = detect_opencl() - self.general["subgrid"] = self.args.subgrid - if self.args.subgrid: + # Subgrids + try: + self.general["subgrid"] = self.args.subgrid # Double precision should be used with subgrid for best accuracy self.general["precision"] = "double" if (self.general["subgrid"] and self.general["solver"] == "cuda") or ( @@ -273,6 +274,8 @@ class SimulationConfig: "You cannot currently use CUDA or OpenCL-based " "solvers with models that contain sub-grids." ) raise ValueError + except AttributeError: + self.general["subgrid"] = False # Scenes parameter may not exist if user enters via CLI try: From f9a940182b3ddb5447e245bd5b461e35802ed81c Mon Sep 17 00:00:00 2001 From: Craig Warren Date: Mon, 4 Sep 2023 11:17:53 +0100 Subject: [PATCH 22/35] Updates to handle sigmamax None via CLI --- gprMax/cmds_multiuse.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/gprMax/cmds_multiuse.py b/gprMax/cmds_multiuse.py index 50cee93c..fc62fe2a 100644 --- a/gprMax/cmds_multiuse.py +++ b/gprMax/cmds_multiuse.py @@ -1755,8 +1755,11 @@ class PMLCFS(UserObjectMulti): cfssigma.scalingprofile = sigmascalingprofile cfssigma.scalingdirection = sigmascalingdirection cfssigma.min = float(sigmamin) + if sigmamax == "None": + sigmamax = None if sigmamax is not None: - cfssigma.max = float(sigmamax) + sigmamax = float(sigmamax) + cfssigma.max = sigmamax cfs = CFS() cfs.alpha = cfsalpha cfs.kappa = cfskappa @@ -1770,7 +1773,7 @@ class PMLCFS(UserObjectMulti): f"{cfskappa.scalingdirection}, min: {cfskappa.min:g}, max: " f"{cfskappa.max:g}), sigma (scaling: {cfssigma.scalingprofile}, " f"scaling direction: {cfssigma.scalingdirection}, min: " - f"{cfssigma.min:g}, max: {cfssigma.max:g}) created." + f"{cfssigma.min:g}, max: {cfssigma.max}) created." ) grid.pmls["cfs"].append(cfs) From fb14166e1dfbdcf085ee8ace60ce9c459ec83933 Mon Sep 17 00:00:00 2001 From: Craig Warren Date: Mon, 4 Sep 2023 11:26:06 +0100 Subject: [PATCH 23/35] Fixed issue with PML formulation and CLI --- gprMax/hash_cmds_singleuse.py | 12 +++++------- 1 file changed, 5 insertions(+), 7 deletions(-) diff --git a/gprMax/hash_cmds_singleuse.py b/gprMax/hash_cmds_singleuse.py index 6dd3765b..76456174 100644 --- a/gprMax/hash_cmds_singleuse.py +++ b/gprMax/hash_cmds_singleuse.py @@ -121,6 +121,7 @@ def process_singlecmds(singlecmds): scene_objects.append(tw) cmd = "#pml_formulation" + pml_formulation = None if singlecmds[cmd] is not None: tmp = singlecmds[cmd].split() if len(tmp) != 1: @@ -135,19 +136,16 @@ def process_singlecmds(singlecmds): if len(tmp) not in [1, 6]: logger.exception(f"{cmd} requires either one or six parameter(s)") raise ValueError - if pml_formulation: - formulation = pml_formulation - else: - formulation = None + if len(tmp) == 1: - pml_cells = PMLProps(formulation=formulation, thickness=int(tmp[0])) + pml_cells = PMLProps(formulation=pml_formulation, thickness=int(tmp[0])) else: pml_cells = PMLProps( - formulation=formulation, + formulation=pml_formulation, x0=int(tmp[0]), y0=int(tmp[1]), z0=int(tmp[2]), xmax=int(tmp[3]), ymax=int(tmp[4]), zmax=int(tmp[5]) ) - + scene_objects.append(pml_cells) cmd = "#src_steps" From e33163bbf9d18818cc9bd5b9c30b3ee389341534 Mon Sep 17 00:00:00 2001 From: Craig Warren Date: Mon, 4 Sep 2023 11:38:15 +0100 Subject: [PATCH 24/35] Further tweak to handle pml formulation and CLI --- gprMax/hash_cmds_singleuse.py | 25 ++++++++++++++++--------- 1 file changed, 16 insertions(+), 9 deletions(-) diff --git a/gprMax/hash_cmds_singleuse.py b/gprMax/hash_cmds_singleuse.py index 76456174..61145e13 100644 --- a/gprMax/hash_cmds_singleuse.py +++ b/gprMax/hash_cmds_singleuse.py @@ -121,7 +121,6 @@ def process_singlecmds(singlecmds): scene_objects.append(tw) cmd = "#pml_formulation" - pml_formulation = None if singlecmds[cmd] is not None: tmp = singlecmds[cmd].split() if len(tmp) != 1: @@ -137,16 +136,24 @@ def process_singlecmds(singlecmds): logger.exception(f"{cmd} requires either one or six parameter(s)") raise ValueError - if len(tmp) == 1: - pml_cells = PMLProps(formulation=pml_formulation, thickness=int(tmp[0])) + if pml_formulation in locals(): + if len(tmp) == 1: + pml_props = PMLProps(formulation=pml_formulation, thickness=int(tmp[0])) + else: + pml_props = PMLProps( + formulation=pml_formulation, + x0=int(tmp[0]), y0=int(tmp[1]), z0=int(tmp[2]), + xmax=int(tmp[3]), ymax=int(tmp[4]), zmax=int(tmp[5]) + ) else: - pml_cells = PMLProps( - formulation=pml_formulation, - x0=int(tmp[0]), y0=int(tmp[1]), z0=int(tmp[2]), - xmax=int(tmp[3]), ymax=int(tmp[4]), zmax=int(tmp[5]) - ) + if len(tmp) == 1: + pml_props = PMLProps(thickness=int(tmp[0])) + else: + pml_props = PMLProps( + x0=int(tmp[0]), y0=int(tmp[1]), z0=int(tmp[2]), xmax=int(tmp[3]), ymax=int(tmp[4]), zmax=int(tmp[5]) + ) - scene_objects.append(pml_cells) + scene_objects.append(pml_props) cmd = "#src_steps" if singlecmds[cmd] is not None: From af374440d9a915123ccb704690882c3063d9c454 Mon Sep 17 00:00:00 2001 From: Craig Warren Date: Mon, 4 Sep 2023 11:40:32 +0100 Subject: [PATCH 25/35] Further tweaks for pml formulation and CLI --- gprMax/hash_cmds_singleuse.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/gprMax/hash_cmds_singleuse.py b/gprMax/hash_cmds_singleuse.py index 61145e13..21fea96d 100644 --- a/gprMax/hash_cmds_singleuse.py +++ b/gprMax/hash_cmds_singleuse.py @@ -136,7 +136,7 @@ def process_singlecmds(singlecmds): logger.exception(f"{cmd} requires either one or six parameter(s)") raise ValueError - if pml_formulation in locals(): + if 'pml_formulation' in locals(): if len(tmp) == 1: pml_props = PMLProps(formulation=pml_formulation, thickness=int(tmp[0])) else: From da7be0133d65e2331fbec60a28a1e1ba047f0ef0 Mon Sep 17 00:00:00 2001 From: Craig Warren Date: Mon, 4 Sep 2023 14:50:36 +0100 Subject: [PATCH 26/35] Removed unused funcs --- gprMax/pml.py | 19 ------------------- 1 file changed, 19 deletions(-) diff --git a/gprMax/pml.py b/gprMax/pml.py index fba68ddd..6ce24e5b 100644 --- a/gprMax/pml.py +++ b/gprMax/pml.py @@ -433,19 +433,6 @@ class CUDAPML(PML): 1, ) - def get_update_funcs(self, kernelselectric, kernelsmagnetic): - """Gets update functions from PML kernels. - - Args: - kernelselectric: pycuda SourceModule containing PML kernels for - electric updates. - kernelsmagnetic: pycuda SourceModule containing PML kernels for - magnetic updates. - """ - - self.update_electric_dev = kernelselectric.get_function("order" + str(len(self.CFS)) + "_" + self.direction) - self.update_magnetic_dev = kernelsmagnetic.get_function("order" + str(len(self.CFS)) + "_" + self.direction) - def update_electric(self): """Updates electric field components with the PML correction on the GPU.""" self.update_electric_dev( @@ -549,12 +536,6 @@ class OpenCLPML(PML): self.HPhi1_dev = clarray.to_device(self.queue, self.HPhi1) self.HPhi2_dev = clarray.to_device(self.queue, self.HPhi2) - def set_blocks_per_grid(): - pass - - def get_update_funcs(): - pass - def update_electric(self): """Updates electric field components with the PML correction on the compute device. From 833e86b2853348309d46d0525319120e504127e9 Mon Sep 17 00:00:00 2001 From: Craig Warren Date: Mon, 4 Sep 2023 14:50:49 +0100 Subject: [PATCH 27/35] Formatting --- gprMax/grid.py | 1 + 1 file changed, 1 insertion(+) diff --git a/gprMax/grid.py b/gprMax/grid.py index c6610190..79141a45 100644 --- a/gprMax/grid.py +++ b/gprMax/grid.py @@ -354,6 +354,7 @@ class CUDAGrid(FDTDGrid): self.Hx_dev = gpuarray.to_gpu(self.Hx) self.Hy_dev = gpuarray.to_gpu(self.Hy) self.Hz_dev = gpuarray.to_gpu(self.Hz) + elif config.sim_config.general["solver"] == "opencl": import pyopencl.array as clarray From 119d2d35e0700d21daf6e8bee90472ac9677f3ed Mon Sep 17 00:00:00 2001 From: Craig Warren Date: Mon, 4 Sep 2023 16:15:15 +0100 Subject: [PATCH 28/35] Removed inherited funcs --- gprMax/updates.py | 20 -------------------- 1 file changed, 20 deletions(-) diff --git a/gprMax/updates.py b/gprMax/updates.py index 9adece69..296360ca 100644 --- a/gprMax/updates.py +++ b/gprMax/updates.py @@ -595,11 +595,6 @@ class CUDAUpdates: grid=self.grid.bpg, ) - def update_magnetic_pml(self): - """Updates magnetic field components with the PML correction.""" - for pml in self.grid.pmls["slabs"]: - pml.update_magnetic() - def update_magnetic_sources(self): """Updates magnetic field components from sources.""" if self.grid.magneticdipoles: @@ -662,11 +657,6 @@ class CUDAUpdates: grid=self.grid.bpg, ) - def update_electric_pml(self): - """Updates electric field components with the PML correction.""" - for pml in self.grid.pmls["slabs"]: - pml.update_electric() - def update_electric_sources(self): """Updates electric field components from sources - update any Hertzian dipole sources last. @@ -1169,11 +1159,6 @@ class OpenCLUpdates: ) event.wait() - def update_magnetic_pml(self): - """Updates magnetic field components with the PML correction.""" - for pml in self.grid.pmls["slabs"]: - pml.update_magnetic() - def update_magnetic_sources(self): """Updates magnetic field components from sources.""" if self.grid.magneticdipoles: @@ -1233,11 +1218,6 @@ class OpenCLUpdates: event.wait() - def update_electric_pml(self): - """Updates electric field components with the PML correction.""" - for pml in self.grid.pmls["slabs"]: - pml.update_electric() - def update_electric_sources(self): """Updates electric field components from sources - update any Hertzian dipole sources last. From 1d4f91b9121e49919288cb38b9a159087084bffb Mon Sep 17 00:00:00 2001 From: Craig Warren Date: Mon, 4 Sep 2023 16:21:21 +0100 Subject: [PATCH 29/35] Revert previous change --- gprMax/updates.py | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/gprMax/updates.py b/gprMax/updates.py index 296360ca..9adece69 100644 --- a/gprMax/updates.py +++ b/gprMax/updates.py @@ -595,6 +595,11 @@ class CUDAUpdates: grid=self.grid.bpg, ) + def update_magnetic_pml(self): + """Updates magnetic field components with the PML correction.""" + for pml in self.grid.pmls["slabs"]: + pml.update_magnetic() + def update_magnetic_sources(self): """Updates magnetic field components from sources.""" if self.grid.magneticdipoles: @@ -657,6 +662,11 @@ class CUDAUpdates: grid=self.grid.bpg, ) + def update_electric_pml(self): + """Updates electric field components with the PML correction.""" + for pml in self.grid.pmls["slabs"]: + pml.update_electric() + def update_electric_sources(self): """Updates electric field components from sources - update any Hertzian dipole sources last. @@ -1159,6 +1169,11 @@ class OpenCLUpdates: ) event.wait() + def update_magnetic_pml(self): + """Updates magnetic field components with the PML correction.""" + for pml in self.grid.pmls["slabs"]: + pml.update_magnetic() + def update_magnetic_sources(self): """Updates magnetic field components from sources.""" if self.grid.magneticdipoles: @@ -1218,6 +1233,11 @@ class OpenCLUpdates: event.wait() + def update_electric_pml(self): + """Updates electric field components with the PML correction.""" + for pml in self.grid.pmls["slabs"]: + pml.update_electric() + def update_electric_sources(self): """Updates electric field components from sources - update any Hertzian dipole sources last. From fa87b6c57adb0ed281ddf4f35c972315225e1fc4 Mon Sep 17 00:00:00 2001 From: Craig Warren Date: Mon, 4 Sep 2023 20:23:55 +0100 Subject: [PATCH 30/35] More work on setting precision with subgrids --- gprMax/config.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/gprMax/config.py b/gprMax/config.py index f9b1dcf8..387888fb 100644 --- a/gprMax/config.py +++ b/gprMax/config.py @@ -263,7 +263,7 @@ class SimulationConfig: self.devices["devs"] = detect_opencl() # Subgrids - try: + if hasattr(self.args, "subgrid") and self.args.subgrid: self.general["subgrid"] = self.args.subgrid # Double precision should be used with subgrid for best accuracy self.general["precision"] = "double" @@ -274,7 +274,7 @@ class SimulationConfig: "You cannot currently use CUDA or OpenCL-based " "solvers with models that contain sub-grids." ) raise ValueError - except AttributeError: + else: self.general["subgrid"] = False # Scenes parameter may not exist if user enters via CLI From 31d39e67f546b6057a5ae4e83ac0f02fc77b03bc Mon Sep 17 00:00:00 2001 From: Craig Warren Date: Mon, 4 Sep 2023 21:15:39 +0100 Subject: [PATCH 31/35] Improved syntax --- testing/test_models.py | 19 ++++++++++--------- 1 file changed, 10 insertions(+), 9 deletions(-) diff --git a/testing/test_models.py b/testing/test_models.py index 18c0e13d..c4a85348 100644 --- a/testing/test_models.py +++ b/testing/test_models.py @@ -75,7 +75,7 @@ for i, model in enumerate(testmodels): # Run model file = basepath / model / model - gprMax.run(inputfile=file.with_suffix(".in"), gpu=None) + gprMax.run(inputfile=file.with_suffix(".in"), gpu=None, opencl=None) # Special case for analytical comparison if model == "hertzian_dipole_fs_analytical": @@ -111,6 +111,7 @@ for i, model in enumerate(testmodels): dataref = hertzian_dipole_fs( filetest.attrs["Iterations"], filetest.attrs["dt"], filetest.attrs["dx_dy_dz"], rxposrelative ) + filetest.close() else: # Get output for model and reference files @@ -129,15 +130,15 @@ for i, model in enumerate(testmodels): raise ValueError # Check that type of float used to store fields matches - if filetest[path + outputstest[0]].dtype != fileref[path + outputsref[0]].dtype: - logger.warning( - f"Type of floating point number in test model " - f"({filetest[path + outputstest[0]].dtype}) does not " - f"match type in reference solution ({fileref[path + outputsref[0]].dtype})\n" - ) float_or_doubleref = fileref[path + outputsref[0]].dtype float_or_doubletest = filetest[path + outputstest[0]].dtype - + if float_or_doubleref != float_or_doubletest: + logger.warning( + f"Type of floating point number in test model " + f"({float_or_doubletest}) does not " + f"match type in reference solution ({float_or_doubleref})\n" + ) + # Arrays for storing time timeref = np.zeros((fileref.attrs["Iterations"]), dtype=float_or_doubleref) timeref = ( @@ -161,7 +162,7 @@ for i, model in enumerate(testmodels): raise ValueError fileref.close() - filetest.close() + filetest.close() # Diffs datadiffs = np.zeros(datatest.shape, dtype=np.float64) From 8d477704059057ac42a21c16e4a2fb7d319440db Mon Sep 17 00:00:00 2001 From: Craig Warren Date: Mon, 4 Sep 2023 21:16:33 +0100 Subject: [PATCH 32/35] Fixed PML bug with CUDA solver --- gprMax/updates.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/gprMax/updates.py b/gprMax/updates.py index 9adece69..5691c3b3 100644 --- a/gprMax/updates.py +++ b/gprMax/updates.py @@ -407,7 +407,7 @@ class CUDAUpdates: knlH = self.source_module(bld, options=config.sim_config.devices["nvcc_opts"]) pml.update_magnetic_dev = knlH.get_function(knl_name) - self._copy_mat_coeffs(knlE, knlH) + self._copy_mat_coeffs(knlE, knlH) def _set_rx_knl(self): """Receivers - initialises arrays on GPU, prepares kernel and gets kernel From 77428a300eeb398a52e9ed32a0bcaaccb100aa74 Mon Sep 17 00:00:00 2001 From: Craig Warren Date: Mon, 4 Sep 2023 21:16:52 +0100 Subject: [PATCH 33/35] Tidied some comments --- gprMax/config.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/gprMax/config.py b/gprMax/config.py index 387888fb..950072aa 100644 --- a/gprMax/config.py +++ b/gprMax/config.py @@ -242,7 +242,7 @@ class SimulationConfig: # Both single and double precision are possible on GPUs, but single # provides best performance. self.general["precision"] = "single" - self.devices = {"devs": [], "nvcc_opts": None} # pycuda device objects # nvcc compiler options + self.devices = {"devs": [], "nvcc_opts": None} # pycuda device objects; nvcc compiler options # Suppress nvcc warnings on Microsoft Windows if sys.platform == "win32": self.devices["nvcc_opts"] = ["-w"] @@ -254,7 +254,7 @@ class SimulationConfig: if self.args.opencl is not None: self.general["solver"] = "opencl" self.general["precision"] = "single" - self.devices = {"devs": [], "compiler_opts": None} # pyopencl available device(s) + self.devices = {"devs": [], "compiler_opts": None} # pyopencl device device(s); compiler options # Suppress unused variable warnings on gcc # if sys.platform != 'win32': self.devices['compiler_opts'] = ['-w'] From 1dd3b8d247f2abfa5dd67c29f35cd5066e2e56cc Mon Sep 17 00:00:00 2001 From: Craig Warren Date: Tue, 5 Sep 2023 08:41:28 +0100 Subject: [PATCH 34/35] Added comment on coeff array copying to GPU for PMLs --- gprMax/updates.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/gprMax/updates.py b/gprMax/updates.py index 5691c3b3..c6d103ed 100644 --- a/gprMax/updates.py +++ b/gprMax/updates.py @@ -407,6 +407,8 @@ class CUDAUpdates: knlH = self.source_module(bld, options=config.sim_config.devices["nvcc_opts"]) pml.update_magnetic_dev = knlH.get_function(knl_name) + # Copy material coefficient arrays to constant memory of GPU - must + # be done for each kernel self._copy_mat_coeffs(knlE, knlH) def _set_rx_knl(self): From 290a397880b06cfdaab0ab3fa58f49137b93a6c6 Mon Sep 17 00:00:00 2001 From: Craig Warren Date: Tue, 5 Sep 2023 09:12:50 +0100 Subject: [PATCH 35/35] Fixing Win path backslashes --- README.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.rst b/README.rst index 07b3c31b..8b39624d 100644 --- a/README.rst +++ b/README.rst @@ -117,7 +117,7 @@ Microsoft Windows ^^^^^^^^^^^^^^^^^ * Download and install Microsoft `Build Tools for Visual Studio 2022 `_ (direct link). You can also find it on the `Microsoft Visual Studio downloads page `_ by scrolling down to the 'All Downloads' section, clicking the disclosure triangle by 'Tools for Visual Studio 2022', then clicking the download button next to 'Build Tools for Visual Studio 2022'. When installing, choose the 'Desktop development with C++' Workload and select only 'MSVC v143' and 'Windows 10 SDK' or 'Windows 11 SDK options. -* Set the Path and Environment Variables - this can be done by following the `instructions from Microsoft `_, or manually by adding a form of :code:`C:\Program Files (x86)\Microsoft Visual Studio\2019\BuildTools\VC\Tools\MSVC\14.23.28105\bin\Hostx64\x64` (this may vary according to your exact machine and installation) to your system Path environment variable. +* Set the Path and Environment Variables - this can be done by following the `instructions from Microsoft `_, or manually by adding a form of :code:``C:\Program Files (x86)\Microsoft Visual Studio\2019\BuildTools\VC\Tools\MSVC\14.23.28105\bin\Hostx64\x64`` (this may vary according to your exact machine and installation) to your system Path environment variable. Alternatively, if you are using Windows 10/11 you can install the `Windows Subsystem for Linux `_ and then follow the Linux install instructions for gprMax. Note however that currently, WSL does not aim to support GUI desktops or applications, e.g. Gnome, KDE, etc....