From 4a610faebb9d888bcfe1376cd343ede897fe6210 Mon Sep 17 00:00:00 2001 From: Craig Warren Date: Fri, 9 Jun 2023 14:56:34 +0100 Subject: [PATCH] Cleaned up with better comments --- .../gssi_400_over_fractal_subsurface.py | 154 +++++++++--------- 1 file changed, 76 insertions(+), 78 deletions(-) diff --git a/examples/subgrids/gssi_400_over_fractal_subsurface.py b/examples/subgrids/gssi_400_over_fractal_subsurface.py index 9a4546a8..358f271b 100644 --- a/examples/subgrids/gssi_400_over_fractal_subsurface.py +++ b/examples/subgrids/gssi_400_over_fractal_subsurface.py @@ -5,102 +5,117 @@ from pathlib import Path +import numpy as np + import gprMax from toolboxes.GPRAntennaModels.GSSI import antenna_like_GSSI_400 -import numpy as np - -# file path step +# File path - used later to specify name of output files fn = Path(__file__) parts = fn.parts -# subgrid Discretisation is 1 mm in x, y, z directions. This allows us -# to model the geometry of the antenna -dl_s = 1e-3 -# subgridding ratio. This must always be an odd integer multiple. In this case -# the main grid discrestisation is 9e-3 m. +# 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_s * 9 +dl = dl_sg * ratio -# estimated return time for signal to propagate 1 metre and back -tw = 2 / 3e8 * (np.sqrt(3.2) + np.sqrt(9)) - -# domain extent +# 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_gpr = gprMax.Title(name=fn.name) +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(domain) -scene.add(title_gpr) +scene.add(title) scene.add(dxdydz) +scene.add(domain) scene.add(time_window) -# half space material -halfspace_m = gprMax.Material(er=3.2, se=0.397e-3, mr=1, sm=0, id='soil') -scene.add(halfspace_m) - +# Dimensions of antenna case antenna_case = (0.3, 0.3, 0.178) -bounding_box = 2 * dl - -# pml + boundary_cells + is_os + subgrid+boundary + half antenna -x0 = (10 + 15 + 5 + 2) * dl + antenna_case[0] / 2 - -#antenna_p = (x / 2, y / 2, z - 30 * dl - bounding_box - antenna_case[2]) - -# Position of antenna base +# 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 -sg = gprMax.SubGridHSG(p1=[sg_x0, sg_y0, sg_z0], p2=[sg_x1, sg_y1, sg_z1], ratio=ratio, id='subgrid1') +# 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) -# half space box in main grid -halfspace = gprMax.Box(p1=(0, 0, 0), p2=(x, y, antenna_p[2]), material_id='soil') -scene.add(halfspace) +# 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 half space box in the subgrid. The halfspace has to be positioned -# manually because it traverses the grid. Grid traversal requires that objects extend -# beyond the OS. Turning off autotranslate allows you to place objects beyond the OS. +# 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 seperation from the OS +# PML separation from the outer surface ps = ratio // 2 + 2 -# number of pml cells in the subgrid +# Number of PML cells in the subgrid pc = 6 -# is os seperation +# Inner surface/outer surface separation isos = 3 * ratio -h = antenna_p[2] - sg_z0 + (ps + pc + isos) * dl_s +# Calculate maximum z-coordinate (height) for box of sandy_soil in subgrid +h = antenna_p[2] - sg_z0 + (ps + pc + isos) * dl_sg -# half space box -halfspace = gprMax.Box(p1=(0, 0, 0), p2=(411 * dl_s, 411 * dl_s, h), material_id='soil') -# position the box using local coordinates3e8 / 400e6 -halfspace.autotranslate = False -sg.add(halfspace) +# 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 the antenna model and add components to subgrid -gssi_objects = antenna_like_GSSI_400(*antenna_p, resolution=dl_s) +# 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) -# half space box -halfspace = gprMax.Box(p1=(0, 0, 0), p2=(x, y, antenna_p[2]), material_id='soil') -scene.add(halfspace) +# 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), @@ -109,34 +124,17 @@ for i in range(1, 51): time=i * tw / 50,) scene.add(snap) -# create a geometry view of the main grid and the sub grid stitched together -gv1 = gprMax.GeometryView(p1=(sg_x0, sg_y0, sg_z0), p2=(sg_x1, sg_y1, sg_z1), - dl=(dl_s, dl_s, dl_s), - filename=fn.with_suffix('').parts[-1] + '_sg', - output_type='n') -#sg.add(gv1) +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) -gv3 = gprMax.GeometryView(p1=(sg_x0, sg_y0, 1.512), p2=(sg_x1, sg_y1, 1.513), - dl=(dl_s, dl_s, dl_s), - filename=fn.with_suffix('').parts[-1] + '_sg', - output_type='f') -sg.add(gv3) - -gv2 = gprMax.GeometryView(p1=(0, 0, 0), - p2=domain.props.p1, - dl=dl, - filename=fn.with_suffix('').parts[-1], - output_type='n') -scene.add(gv2) - -# half space material -layer_m = gprMax.Material(er=9, se=0.397e-3, mr=1, sm=0, id='soil_2') -scene.add(layer_m) - -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_2', id='fbox', seed=1) -scene.add(fb) - -sr = 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(sr) +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=True, outputfile=fn, subgrid=True, autotranslate=True)