Add subgrid and antenna reframe tests

reframe_tests/base_tests.py

@@ -212,5 +212,9 @@ class GprMaxRegressionTest(GprMaxBaseTest):
             )
 
 
+class GprMaxAPIRegressionTest(GprMaxRegressionTest):
+    executable = "time -p python"
+
+
 class GprMaxMpiTest(GprMaxBaseTest):
     pass

reframe_tests/reframe_tests.py

@@ -1,5 +1,5 @@
 import reframe as rfm
-from base_tests import GprMaxRegressionTest
+from base_tests import GprMaxAPIRegressionTest, GprMaxRegressionTest
 from reframe.core.builtins import parameter, run_after
 
 """ReFrame tests for basic functionality
@@ -82,3 +82,71 @@ class BasicModelsTest(GprMaxRegressionTest):
         self.executable_opts = [self.input_file, "-o", self.output_file]
         self.postrun_cmds = [f"python -m toolboxes.Plotting.plot_Ascan -save {self.output_file}"]
         self.keep_files = [self.input_file, self.output_file, f"{self.model}.pdf"]
+
+
+@rfm.simple_test
+class AntennaModelsTest(GprMaxRegressionTest):
+    tags = {"test", "serial", "regression", "antenna"}
+
+    # List of available antenna test models
+    model = parameter(
+        [
+            "antenna_wire_dipole_fs",
+        ]
+    )
+    num_cpus_per_task = 16
+
+    @run_after("init")
+    def set_filenames(self):
+        self.input_file = f"{self.model}.in"
+        self.output_file = f"{self.model}.h5"
+        self.executable_opts = [self.input_file, "-o", self.output_file]
+        self.postrun_cmds = [f"python -m toolboxes.Plotting.plot_Ascan -save {self.output_file}"]
+        self.postrun_cmds = [
+            f"python -m toolboxes.Plotting.plot_antenna_params -save {self.output_file}"
+        ]
+
+        antenna_t1_params = f"{self.model}_t1_params.pdf"
+        antenna_ant_params = f"{self.model}_ant_params.pdf"
+        plot_ascan_output = f"{self.model}.pdf"
+        geometry_view = f"{self.model}.vtu"
+        self.keep_files = [
+            self.input_file,
+            self.output_file,
+            antenna_t1_params,
+            antenna_ant_params,
+            plot_ascan_output,
+            geometry_view,
+        ]
+
+
+@rfm.simple_test
+class SubgridTest(GprMaxAPIRegressionTest):
+    tags = {"test", "api", "serial", "regression", "subgrid"}
+
+    # List of available subgrid test models
+    model = parameter(
+        [
+            "cylinder_fs",
+            # "gssi_400_over_fractal_subsurface",  # Takes ~1hr 30m on ARCHER2
+        ]
+    )
+    num_cpus_per_task = 16
+
+    @run_after("init")
+    def set_filenames(self):
+        self.input_file = f"{self.model}.py"
+        self.output_file = f"{self.model}.h5"
+        self.executable_opts = [self.input_file, "-o", self.output_file]
+        self.postrun_cmds = [f"python -m toolboxes.Plotting.plot_Ascan -save {self.output_file}"]
+
+        geometry_view = f"{self.model}.vti"
+        subgrid_geometry_view = f"{self.model}_sg.vti"
+        plot_ascan_output = f"{self.model}.pdf"
+        self.keep_files = [
+            self.input_file,
+            self.output_file,
+            geometry_view,
+            subgrid_geometry_view,
+            plot_ascan_output,
+        ]

reframe_tests/src/antenna_wire_dipole_fs.in

@@ -0,0 +1,15 @@
+#title: Wire antenna - half-wavelength dipole in free-space
+#domain: 0.050 0.050 0.200
+#dx_dy_dz: 0.001 0.001 0.001
+#time_window: 60e-9
+
+#waveform: gaussian 1 1e9 mypulse
+#transmission_line: z 0.025 0.025 0.100 73 mypulse
+
+## 150mm length
+#edge: 0.025 0.025 0.025 0.025 0.025 0.175 pec
+
+## 1mm gap at centre of dipole
+#edge: 0.025 0.025 0.100 0.025 0.025 0.101 free_space
+
+#geometry_view: 0.020 0.020 0.020 0.030 0.030 0.180 0.001 0.001 0.001 antenna_wire_dipole_fs f

reframe_tests/src/cylinder_fs.py

@@ -0,0 +1,115 @@
+"""Cylinder in freespace
+
+This example model demonstrates how to use subgrids at a basic level.
+
+The geometry is 3D (required for any use of subgrids) and is of a water-filled
+cylindrical object in freespace. The subgrid encloses the cylinderical object
+using a fine spatial discretisation (1mm), and a courser spatial discretisation
+(5mm) is used in the rest of the model (main grid). A simple Hertzian dipole
+source is used with a waveform shaped as the first derivative of a gaussian.
+"""
+
+from pathlib import Path
+
+import gprMax
+from gprMax.materials import calculate_water_properties
+
+# File path - used later to specify name of output files
+fn = Path(__file__)
+parts = fn.parts
+
+# Subgrid spatial discretisation in x, y, z directions
+dl_sg = 1e-3
+
+# Subgrid ratio - must always be an odd integer multiple
+ratio = 5
+dl = dl_sg * ratio
+
+# Domain extent
+x = 0.500
+y = 0.500
+z = 0.500
+
+# Time window
+tw = 6e-9
+
+scene = gprMax.Scene()
+
+title = gprMax.Title(name=fn.name)
+dxdydz = gprMax.Discretisation(p1=(dl, dl, dl))
+domain = gprMax.Domain(p1=(x, y, z))
+time_window = gprMax.TimeWindow(time=tw)
+
+wf = gprMax.Waveform(wave_type="gaussiandot", amp=1, freq=1.5e9, id="mypulse")
+hd = gprMax.HertzianDipole(polarisation="z", p1=(0.205, 0.400, 0.250), waveform_id="mypulse")
+rx = gprMax.Rx(p1=(0.245, 0.400, 0.250))
+
+scene.add(title)
+scene.add(dxdydz)
+scene.add(domain)
+scene.add(time_window)
+scene.add(wf)
+scene.add(hd)
+scene.add(rx)
+
+# Cylinder parameters
+c1 = (0.225, 0.250, 0.100)
+c2 = (0.225, 0.250, 0.400)
+r = 0.010
+sg1 = (c1[0] - r, c1[1] - r, c1[2])
+sg2 = (c2[0] + r, c2[1] + r, c2[2])
+
+# Create subgrid
+subgrid = gprMax.SubGridHSG(p1=sg1, p2=sg2, ratio=ratio, id="sg")
+scene.add(subgrid)
+
+# Create water material
+eri, er, tau, sig = calculate_water_properties()
+water = gprMax.Material(er=eri, se=sig, mr=1, sm=0, id="water")
+subgrid.add(water)
+water = gprMax.AddDebyeDispersion(poles=1, er_delta=[er - eri], tau=[tau], material_ids=["water"])
+subgrid.add(water)
+
+# Add cylinder to subgrid
+cylinder = gprMax.Cylinder(p1=c1, p2=c2, r=r, material_id="water")
+subgrid.add(cylinder)
+
+# Create some geometry views for both subgrid and main grid
+gvsg = gprMax.GeometryView(
+    p1=sg1,
+    p2=sg2,
+    dl=(dl_sg, dl_sg, dl_sg),
+    filename=fn.with_suffix("").parts[-1] + "_sg",
+    output_type="n",
+)
+subgrid.add(gvsg)
+
+gv1 = gprMax.GeometryView(
+    p1=(0, 0, 0),
+    p2=(x, y, z),
+    dl=(dl, dl, dl),
+    filename=fn.with_suffix("").parts[-1],
+    output_type="n",
+)
+scene.add(gv1)
+
+# Create some snapshots of entire domain
+for i in range(5):
+    s = gprMax.Snapshot(
+        p1=(0, 0, 0),
+        p2=(x, y, z),
+        dl=(dl, dl, dl),
+        time=(i + 0.5) * 1e-9,
+        filename=fn.with_suffix("").parts[-1] + "_" + str(i + 1),
+    )
+    scene.add(s)
+
+gprMax.run(
+    scenes=[scene],
+    n=1,
+    geometry_only=False,
+    outputfile=fn,
+    subgrid=True,
+    autotranslate=True,
+    log_level=25,
+)

reframe_tests/src/gssi_400_over_fractal_subsurface.py

@@ -0,0 +1,166 @@
+"""GPR antenna model (like a GSSI 400MHz antenna) over layered media with a
+    rough subsurface interface.
+
+This example model demonstrates how to use subgrids at a more advanced level -
+    combining use of an imported antenna model and rough subsurface interface.
+
+The geometry is 3D (required for any use of subgrids) and is of a 2 layered
+subsurface. The top layer in a sandy soil and the bottom layer a soil with
+higher permittivity (both have some simple conductive loss). There is a rough
+interface between the soil layers. A GPR antenna model (like a GSSI 400MHz
+antenna) is imported and placed on the surface of the layered media. The antenna
+is meshed using a subgrid with a fine spatial discretisation (1mm), and a
+courser spatial discretisation (9mm) is used in the rest of the model (main
+grid).
+"""
+
+from pathlib import Path
+
+import numpy as np
+
+import gprMax
+from toolboxes.GPRAntennaModels.GSSI import antenna_like_GSSI_400
+
+# File path - used later to specify name of output files
+fn = Path(__file__)
+parts = fn.parts
+
+# Subgrid spatial discretisation in x, y, z directions
+dl_sg = 1e-3
+
+# Subgrid ratio - must always be an odd integer multiple
+ratio = 9
+dl = dl_sg * ratio
+
+# Domain extent
+x = 3
+y = 1
+z = 2
+
+# Time window
+# Estimated two way travel time over 1 metre in material with highest
+# permittivity, slowest velocity.
+tw = 2 / 3e8 * (np.sqrt(3.2) + np.sqrt(9))
+
+scene = gprMax.Scene()
+
+title = gprMax.Title(name=fn.name)
+dxdydz = gprMax.Discretisation(p1=(dl, dl, dl))
+domain = gprMax.Domain(p1=(x, y, z))
+time_window = gprMax.TimeWindow(time=tw)
+
+scene.add(title)
+scene.add(dxdydz)
+scene.add(domain)
+scene.add(time_window)
+
+# Dimensions of antenna case
+antenna_case = (0.3, 0.3, 0.178)
+
+# Position of antenna
+antenna_p = (x / 2, y / 2, 170 * dl)
+
+# Extra distance surrounding antenna for subgrid
+bounding_box = 2 * dl
+
+# Subgrid extent
+sg_x0 = antenna_p[0] - antenna_case[0] / 2 - bounding_box
+sg_y0 = antenna_p[1] - antenna_case[1] / 2 - bounding_box
+sg_z0 = antenna_p[2] - bounding_box
+sg_x1 = antenna_p[0] + antenna_case[0] / 2 + bounding_box
+sg_y1 = antenna_p[1] + antenna_case[1] / 2 + bounding_box
+sg_z1 = antenna_p[2] + antenna_case[2] + bounding_box
+
+# Create subgrid
+sg = gprMax.SubGridHSG(p1=[sg_x0, sg_y0, sg_z0], p2=[sg_x1, sg_y1, sg_z1], ratio=ratio, id="sg")
+scene.add(sg)
+
+# Create and add a box of homogeneous material to main grid - sandy_soil
+sandy_soil = gprMax.Material(er=3.2, se=0.397e-3, mr=1, sm=0, id="sandy_soil")
+scene.add(sandy_soil)
+b1 = gprMax.Box(p1=(0, 0, 0), p2=(x, y, antenna_p[2]), material_id="sandy_soil")
+scene.add(b1)
+
+# Position box of sandy_soil in the subgrid.
+#   It has to be positioned manually because it traverses the main grid/subgrid
+#   interface. Grid traversal is when objects extend beyond the outer surface.
+#   Setting autotranslate to false allows you to place objects beyond the outer
+#   surface.
+
+# PML separation from the outer surface
+ps = ratio // 2 + 2
+# Number of PML cells in the subgrid
+pc = 6
+# Inner surface/outer surface separation
+isos = 3 * ratio
+
+# Calculate maximum z-coordinate (height) for box of sandy_soil in subgrid
+h = antenna_p[2] - sg_z0 + (ps + pc + isos) * dl_sg
+
+# Create and add a box of homogeneous material to subgrid - sandy_soil
+sg.add(sandy_soil)
+b2 = gprMax.Box(p1=(0, 0, 0), p2=(411 * dl_sg, 411 * dl_sg, h), material_id="sandy_soil")
+# Set autotranslate for the box object to false
+b2.autotranslate = False
+sg.add(b2)
+
+# Import antenna model and add components to subgrid
+gssi_objects = antenna_like_GSSI_400(*antenna_p, resolution=dl_sg)
+for obj in gssi_objects:
+    sg.add(obj)
+
+# Create and add a homogeneous material with a rough surface
+soil = gprMax.Material(er=9, se=0.397e-3, mr=1, sm=0, id="soil")
+scene.add(soil)
+
+fb = gprMax.FractalBox(
+    p1=(0, 0, 0),
+    p2=(3, 1, 1),
+    frac_dim=1.5,
+    weighting=(1, 1, 1),
+    n_materials=1,
+    mixing_model_id="soil",
+    id="fbox",
+    seed=1,
+)
+scene.add(fb)
+
+rough_surf = gprMax.AddSurfaceRoughness(
+    p1=(0, 0, 1),
+    p2=(3, 1, 1),
+    frac_dim=1.5,
+    weighting=(1, 1),
+    limits=(0.4, 1.2),
+    fractal_box_id="fbox",
+    seed=1,
+)
+scene.add(rough_surf)
+
+# Create some snapshots and geometry views
+for i in range(1, 51):
+    snap = gprMax.Snapshot(
+        p1=(0, y / 2, 0),
+        p2=(x, y / 2 + dl, z),
+        dl=(dl, dl, dl),
+        filename=Path(*parts[:-1], f"{parts[-1]}_{str(i)}").name,
+        time=i * tw / 50,
+    )
+    scene.add(snap)
+
+gvsg = gprMax.GeometryView(
+    p1=(sg_x0, sg_y0, sg_z0),
+    p2=(sg_x1, sg_y1, sg_z1),
+    dl=(dl_sg, dl_sg, dl_sg),
+    filename=fn.with_suffix("").parts[-1] + "_sg",
+    output_type="n",
+)
+sg.add(gvsg)
+
+gv1 = gprMax.GeometryView(
+    p1=(0, 0, 0), p2=domain.props.p1, dl=dl, filename=fn.with_suffix("").parts[-1], output_type="n"
+)
+scene.add(gv1)
+
+gprMax.run(
+    scenes=[scene], n=1, geometry_only=False, outputfile=fn, subgrid=True, autotranslate=True
+)

toolboxes/Plotting/plot_antenna_params.py

@@ -28,7 +28,9 @@ import numpy as np
 logger = logging.getLogger(__name__)
 
 
-def calculate_antenna_params(filename, tltxnumber=1, tlrxnumber=None, rxnumber=None, rxcomponent=None):
+def calculate_antenna_params(
+    filename, tltxnumber=1, tlrxnumber=None, rxnumber=None, rxcomponent=None
+):
     """Calculates antenna parameters - incident, reflected and total volatges
         and currents; s11, (s21) and input impedance.
 
@@ -224,7 +226,10 @@ def mpl_plot(
 
     # Print some useful values from s11, and input impedance
     s11minfreq = np.where(s11[pltrange] == np.amin(s11[pltrange]))[0][0]
-    logger.info(f"s11 minimum: {np.amin(s11[pltrange]):g} dB at " + f"{freqs[s11minfreq + pltrangemin]:g} Hz")
+    logger.info(
+        f"s11 minimum: {np.amin(s11[pltrange]):g} dB at "
+        + f"{freqs[s11minfreq + pltrangemin]:g} Hz"
+    )
     logger.info(f"At {freqs[s11minfreq + pltrangemin]:g} Hz...")
     logger.info(
         f"Input impedance: {np.abs(zin[s11minfreq + pltrangemin]):.1f}"
@@ -236,7 +241,10 @@ def mpl_plot(
     # Figure 1
     # Plot incident voltage
     fig1, ax = plt.subplots(
-        num="Transmitter transmission line parameters", figsize=(20, 12), facecolor="w", edgecolor="w"
+        num="Transmitter transmission line parameters",
+        figsize=(20, 12),
+        facecolor="w",
+        edgecolor="w",
     )
     gs1 = gridspec.GridSpec(4, 2, hspace=0.7)
     ax = plt.subplot(gs1[0, 0])
@@ -368,7 +376,9 @@ def mpl_plot(
 
     # Figure 2
     # Plot frequency spectra of s11
-    fig2, ax = plt.subplots(num="Antenna parameters", figsize=(20, 12), facecolor="w", edgecolor="w")
+    fig2, ax = plt.subplots(
+        num="Antenna parameters", figsize=(20, 12), facecolor="w", edgecolor="w"
+    )
     gs2 = gridspec.GridSpec(2, 2, hspace=0.3)
     ax = plt.subplot(gs2[0, 0])
     markerline, stemlines, baseline = ax.stem(freqs[pltrange], s11[pltrange], "-.")
@@ -458,13 +468,26 @@ def mpl_plot(
     # ax.grid(which='both', axis='both', linestyle='-.')
 
     if save:
-        savename1 = filename.stem + "_tl_params"
-        savename1 = filename.parent / savename1
-        savename2 = filename.stem + "_ant_params"
-        savename2 = filename.parent / savename2
+        filepath = Path(filename)
+        savename1 = filepath.stem + "_tl_params"
+        savename1 = filepath.parent / savename1
+        savename2 = filepath.stem + "_ant_params"
+        savename2 = filepath.parent / savename2
         # Save a PDF of the figure
-        fig1.savefig(savename1.with_suffix(".pdf"), dpi=None, format="pdf", bbox_inches="tight", pad_inches=0.1)
-        fig2.savefig(savename2.with_suffix(".pdf"), dpi=None, format="pdf", bbox_inches="tight", pad_inches=0.1)
+        fig1.savefig(
+            savename1.with_suffix(".pdf"),
+            dpi=None,
+            format="pdf",
+            bbox_inches="tight",
+            pad_inches=0.1,
+        )
+        fig2.savefig(
+            savename2.with_suffix(".pdf"),
+            dpi=None,
+            format="pdf",
+            bbox_inches="tight",
+            pad_inches=0.1,
+        )
         # Save a PNG of the figure
         # fig1.savefig(savename1.with_suffix('.png'), dpi=150, format='png',
         #             bbox_inches='tight', pad_inches=0.1)
@@ -485,7 +508,9 @@ if __name__ == "__main__":
         usage="cd gprMax; python -m toolboxes.Plotting.plot_antenna_params outputfile",
     )
     parser.add_argument("outputfile", help="name of output file including path")
-    parser.add_argument("--tltx-num", default=1, type=int, help="transmitter antenna - transmission line number")
+    parser.add_argument(
+        "--tltx-num", default=1, type=int, help="transmitter antenna - transmission line number"
+    )
     parser.add_argument("--tlrx-num", type=int, help="receiver antenna - transmission line number")
     parser.add_argument("--rx-num", type=int, help="receiver antenna - output number")
     parser.add_argument(
@@ -495,7 +520,10 @@ if __name__ == "__main__":
         choices=["Ex", "Ey", "Ez"],
     )
     parser.add_argument(
-        "-save", action="store_true", default=False, help="save plot directly to file, i.e. do not display"
+        "-save",
+        action="store_true",
+        default=False,
+        help="save plot directly to file, i.e. do not display",
     )
     args = parser.parse_args()