Possible bug in v4: erroneous asymmetry in TWTL sim compared with v3

[cstarkjp]

Running a sim of EM wave propagation along a two-wire transmission line (TWTL), I found a possible bug in gprMax v.4.0.0b0 (or possibly some misuse by me of the new API). The v4 sim exhibits non-physical asymmetry, while the same sim in v.3.1.7 is perfectly symmetric as expected.

The main sim domain is ~150mm long by about 40x40mm and has 1mm cells. The ambient medium has a relative permittivity of 81 (simplified from 1-Debye model of water). A pair of parallel 1mm-thick PEC lines 11mm apart models the TWTL; these lines are bent at 90º at each end to bring them together at 1mm-wide transmission line inputs. One input provides a continuous sine ±1V source at 400MHz (50Ω); the other dummy input is set at 0V. (An alternative setup using voltage sources and receivers behaves in a similar manner). Therefore the sim generates a ~TEM wave that initiates at the source, propagates along the parallel "wires", and should set up a steadily varying EM field perfectly symmetric about the centre line between them.

The sim runs for 500ns with ~0.5ns time steps. As you can see from the pic here (at time 250ns), the Ez field component reveals unexpected asymmetry in the v4 sim, whereas the v3 sim behaves as we would expect (visualization in ParaView 5.12). There's also a strong E-field at places within the PEC in v4 but not in v3.

It's possible that I have set up the v4 sim incorrectly, but it's not at all obvious to me how. I'm attaching here the v3 .in file and the v4 .py file (which exploits the new pure-Python API). If I have accidentally introduced an asymmetry, I don't know how.

Video

twtl_stackedvideo_lowres.mov

HDF5 output

v3.1.7

v4.0.0b0

Input files

v3.1.7 input file:

#title: TWTL

#output_dir: Data_v3

#python:

import numpy as np
from gprMax.input_cmd_funcs import (
    domain,
    dx_dy_dz,
    time_window,
    waveform,
    transmission_line,
    box,
    geometry_view,
    snapshot,
    material,
)

###########################################################################


n_frames = 1000
frequency_Hz = 400e6
rod_length_mm = 150
conductivity = 0.0
waveform_type = "contsine"

cell_size_mm = 1
PML_thickness = cell_size_mm*10
tank_width_mm = 40 + cell_size_mm
tank_length_mm = rod_length_mm
do_snapshots = True
do_save_geometry = True


###########################################################################



domain_x_mm = PML_thickness*2 + tank_length_mm
domain_y_mm = PML_thickness*2 + tank_width_mm
domain_z_mm = PML_thickness*2 + tank_width_mm
x_rod_start_mm = PML_thickness

rod_axes_spacing_mm = 6+6+cell_size_mm  # = 13
T_per_frame = 50e-12
T_simulation = T_per_frame*n_frames
dt_snapshot = (T_simulation / n_frames)
time_window(T_simulation)



###########################################################################



rod_width_mm = cell_size_mm if cell_size_mm>=0.2 else 1
source_impedance = 50

cell_size = np.round(cell_size_mm*0.001, 6)
cell = dx_dy_dz(cell_size, cell_size, cell_size,)
set_cell_size = lambda dl, dimension: np.round(dl * dimension/cell_size_mm, 6)
dimensions = [
    set_cell_size(cell_,dimension_) 
    for cell_, dimension_ in zip(cell, [domain_x_mm, domain_y_mm, domain_z_mm])
]
print(f"Dimensions:  {dimensions}", flush=True,)
domain_ = domain(*dimensions)

# Assume the whole domain has a PML buffer layer 10 cells thick
PML_width = cell_size*10



###########################################################################


material(
    permittivity=81,
    conductivity=conductivity,
    permeability=1,
    magconductivity=0,
    name="unspecified_high_dielectric",
)


###########################################################################



tank_interior_indent = PML_width
box(
    tank_interior_indent, tank_interior_indent, tank_interior_indent, 
    domain_.x-tank_interior_indent, domain_.y-tank_interior_indent, domain_.z-tank_interior_indent, 
    material="unspecified_high_dielectric",
)



###########################################################################



## Probe built from two bent rods using PEC
rod_width = rod_width_mm/1000
y_mid = (dimensions[1] + cell[1])/2
z_mid = (dimensions[2] + cell[2])/2
x_rod_start = x_rod_start_mm/1000
x_rod_end = (x_rod_start_mm + rod_length_mm)/1000

rod_axes_spacing = rod_axes_spacing_mm/1000
rod_z_separation = (rod_axes_spacing-(2*rod_width)/2)
z_rod1 = y_mid - rod_z_separation/2 - rod_width
z_rod2 = y_mid + rod_z_separation/2 - cell.z

## Rod 1 
box(
    x_rod_start, y_mid, z_rod1,  
    x_rod_end, y_mid+rod_width, z_rod1+rod_width, 
    "pec",
)
## Rod 2
box(
    x_rod_start, y_mid, z_rod2, 
    x_rod_end, y_mid+rod_width, z_rod2+rod_width, 
    "pec",
)
## Linking (bent) rod
box(
    x_rod_start, y_mid,  z_rod1, 
    x_rod_start+rod_width, y_mid+rod_width, z_rod2+rod_width, 
    "pec",
)
## Gap between bent rods
box(
    x_rod_start, y_mid, z_mid-cell.z,  
    x_rod_start+cell.x, y_mid+rod_width, z_mid,  
    "unspecified_high_dielectric",
)
## Linking (bent) rod
box(
    x_rod_end-rod_width, y_mid, z_rod1, 
    x_rod_end, y_mid+rod_width, z_rod2+rod_width, 
    "pec",
)
## Gap between bent rods
box(
    x_rod_end-rod_width, y_mid, z_mid-cell.z, 
    x_rod_end, y_mid+rod_width, z_mid, 
    "unspecified_high_dielectric",
)


###########################################################################



tx_waveform = waveform(
    waveform_type, 
    amplitude=1, 
    frequency=frequency_Hz, 
    identifier="contsine_in"
)
print(f"Source waveform = {type(tx_waveform)}: {tx_waveform}")

rx_waveform = waveform(
    waveform_type, 
    amplitude=1e-40, 
    frequency=1e-40, 
    identifier="dummy_out"
)



## Transmission line source across gap between probe rods in z direction with impedance=50ohms
transmission_line(
    polarisation="z",
    f1=x_rod_start, f2=y_mid, f3=z_mid-cell.z,
    resistance=source_impedance,
    identifier="contsine_in",
)
transmission_line(
    polarisation="z",
    f1=x_rod_end-cell.x, f2=y_mid, f3=z_mid-cell.z,
    resistance=source_impedance,
    identifier="dummy_out",
)



###########################################################################



# Geometry - for use in ParaView
if do_save_geometry:
    geometry_view(
        0, 0, 0, 
        dimensions[0], dimensions[1], dimensions[2], 
        cell.x, cell.y, cell.z, 
        "twtl" \
            +f"_f{frequency_Hz/1000000:3.0f}MHz".replace(".","p") \
            +f"_rodsep{rod_axes_spacing*1000:0.4g}mm".replace(".","p") \
            +f"_res{cell_size*1000:0.2g}mm".replace(".","p") \
            +f"_{"unspecified_high_dielectric"}".replace(".","p"), 
        "n",
    )



# Snapshot just a thin vertical slice of the domain through the probe
if do_snapshots:
    for i in range(1, n_frames+1):
        _ = snapshot(
            PML_width, y_mid-cell.y, PML_width, 
            dimensions[0]-PML_width, y_mid+cell.y, dimensions[2]-PML_width,
            cell.x, cell.y, cell.z, 
            i * dt_snapshot,
            f"xzslice_{i}",
        )

#end_python:

v4.0.0b0 file:

"""
Two-Wire Transmission Line probe
"""

import numpy as np
from pathlib import Path
import gprMax

script_name = Path(__file__)

scene = gprMax.Scene()

###########################################################################


n_frames = 1000
frequency_Hz = 400e6
rod_length_mm = 150
conductivity = 0.0
waveform_type = "contsine"

cell_size_mm = 1
PML_thickness = cell_size_mm*10
unspecified_high_dielectric_tank_width_mm = 40 + cell_size_mm
unspecified_high_dielectric_tank_length_mm = rod_length_mm
do_snapshots = True
do_save_geometry = True


###########################################################################


title_ = f"{script_name.with_suffix('').name}" \
            +f"_f{frequency_Hz/1000000:3.0f}MHz".replace(".","p") \
            +f"_res{cell_size_mm:0.2g}mm".replace(".","p") \
            +f"_{"unspecified_high_dielectric"}".replace(".","p"), 
title = gprMax.Title(name=title_)
output_dir = gprMax.OutputDir(dir="Data_v4",)
scene.add(title)
scene.add(output_dir)



###########################################################################



domain_x_mm = PML_thickness*2 + unspecified_high_dielectric_tank_length_mm
domain_y_mm = PML_thickness*2 + unspecified_high_dielectric_tank_width_mm
domain_z_mm = PML_thickness*2 + unspecified_high_dielectric_tank_width_mm
x_rod_start_mm = PML_thickness

rod_axes_spacing_mm = 6+6+cell_size_mm  # = 13
T_per_frame = 50e-12
T_simulation = T_per_frame*n_frames
dt_snapshot = (T_simulation / n_frames)
time_window = gprMax.TimeWindow(time=T_simulation)
scene.add(time_window)



###########################################################################



rod_width_mm = cell_size_mm if cell_size_mm>=0.2 else 1
source_impedance = 50

cell_size = np.round(cell_size_mm*0.001, 6)
cell = (cell_size, cell_size, cell_size,)
set_cell_size = lambda dl, dimension: np.round(dl * dimension/cell_size_mm, 6)
dimensions = [
    set_cell_size(cell_,dimension_) 
    for cell_, dimension_ in zip(cell, [domain_x_mm, domain_y_mm, domain_z_mm])
]
gprMax.config.logger.warning(f"Dimensions:  {dimensions}")
domain = gprMax.Domain(p1=dimensions)
dxdydz = gprMax.Discretisation(p1=cell)
scene.add(domain)
scene.add(dxdydz)

# Assume the whole domain has a PML buffer layer 10 cells thick
PML_width = cell_size*10



###########################################################################



unspecified_high_dielectric = gprMax.Material(
    er=81,
    se=0,
    mr=1,
    sm=0,
    id="unspecified_high_dielectric",
)
scene.add(unspecified_high_dielectric)



###########################################################################



tank_interior_indent = PML_width
unspecified_high_dielectric_box = gprMax.Box(
    p1=(tank_interior_indent, tank_interior_indent, tank_interior_indent,), 
    p2=(dimensions[0]-tank_interior_indent-cell_size, dimensions[1]-tank_interior_indent-cell_size, dimensions[2]-tank_interior_indent,), 
    material_id="unspecified_high_dielectric",
)    
scene.add(unspecified_high_dielectric_box)



###########################################################################



## Probe built from two bent rods using PEC
rod_width = rod_width_mm/1000
y_mid = (dimensions[1] + cell[1])/2
z_mid = (dimensions[2] + cell[2])/2
x_rod_start = x_rod_start_mm/1000
x_rod_end = (x_rod_start_mm + rod_length_mm)/1000

rod_axes_spacing = rod_axes_spacing_mm/1000
rod_z_separation = (rod_axes_spacing-(2*rod_width)/2)
z_rod1 = y_mid - rod_z_separation/2 - rod_width
z_rod2 = y_mid + rod_z_separation/2 - cell[2]

## Rod 1 
rod1_box = gprMax.Box(
    p1=(x_rod_start, y_mid, z_rod1,),  
    p2=(x_rod_end, y_mid+rod_width, z_rod1+rod_width,), 
    material_id="pec",
)
scene.add(rod1_box)

## Rod 2
rod2_box = gprMax.Box(
    p1=(x_rod_start, y_mid, z_rod2,), 
    p2=(x_rod_end, y_mid+rod_width, z_rod2+rod_width,), 
    material_id="pec",
)
scene.add(rod2_box)

## Linking (bent) rod
linkrod1_box = gprMax.Box(
    p1=(x_rod_start, y_mid,  z_rod1,), 
    p2=(x_rod_start+rod_width, y_mid+rod_width, z_rod2+rod_width,), 
    material_id="pec",
)
## Gap between bent rods
gap1_box = gprMax.Box(
    p1=(x_rod_start, y_mid, z_mid-cell[2],),  
    p2=(x_rod_start+cell[0], y_mid+rod_width, z_mid,),  
    material_id="unspecified_high_dielectric",
)
scene.add(linkrod1_box)
scene.add(gap1_box)

## Linking (bent) rod
linkrod2_box = gprMax.Box(
    p1=(x_rod_end-rod_width, y_mid, z_rod1,), 
    p2=(x_rod_end, y_mid+rod_width, z_rod2+rod_width,), 
    material_id="pec",
)
## Gap between bent rods
gap2_box = gprMax.Box(
    p1=(x_rod_end-rod_width, y_mid, z_mid-cell[2],), 
    p2=(x_rod_end, y_mid+rod_width, z_mid,), 
    material_id="unspecified_high_dielectric",
)
scene.add(linkrod2_box)
scene.add(gap2_box)




###########################################################################



waveform_in = gprMax.Waveform(
    wave_type=waveform_type, 
    amp=1, 
    freq=frequency_Hz, 
    id="contsine_in"
)
gprMax.config.logger.warning(f"Source waveform = {type(waveform_in)}: {waveform_in}")

waveform_out = gprMax.Waveform(
    wave_type=waveform_type, 
    amp=1e-40, 
    freq=1e-40, 
    id="dummy_out"
)
scene.add(waveform_in)
scene.add(waveform_out)


## Transmission line source across gap between probe rod in z direction with impedance=50ohms
transmission_line_in = gprMax.TransmissionLine(
    polarisation="z",
    p1=(x_rod_start, y_mid, z_mid-cell[2],),
    resistance=source_impedance,
    waveform_id="contsine_in",
)
transmission_line_out = gprMax.TransmissionLine(
    polarisation="z",
    p1=(x_rod_end-cell[0], y_mid, z_mid-cell[2],),
    resistance=source_impedance,
    waveform_id="dummy_out",
)
scene.add(transmission_line_in)
scene.add(transmission_line_out)

# ## Voltage source across gap between probe rod in z direction with impedance=50ohms
# voltage_source = gprMax.VoltageSource(
#     polarisation="z",
#     p1=(x_rod_start, y_mid, z_mid-cell[2],),
#     resistance=source_impedance,
#     waveform_id="contsine_in",
# )
# scene.add(voltage_source)

# ## Receivers to record electric field across input & output cells => voltages
# port1_rx = gprMax.Rx(
#     p1=(x_rod_start, y_mid, z_mid-cell[2],),
#     id="port1_rx",
#     outputs=("Ex", "Ey", "Ez", "Ix", "Iy", "Iz",),
# )
# port2_rx = gprMax.Rx(
#     p1=(x_rod_end-cell[0], y_mid, z_mid-cell[2],),
#     id="port2_rx",
#     outputs=("Ex", "Ey", "Ez", "Ix", "Iy", "Iz",),
# )
# scene.add(port1_rx)
# scene.add(port2_rx)



###########################################################################



# Geometry - for use in ParaView
if do_save_geometry:
    geometry_view = gprMax.GeometryView(
        p1=(0, 0, 0,),
        p2=(dimensions[0], dimensions[1], dimensions[2],),
        dl=(cell[0], cell[1], cell[2],),
        output_type="n",
        filename="twtl" \
            +f"_f{frequency_Hz/1000000:3.0f}MHz".replace(".","p") \
            +f"_rodsep{rod_axes_spacing*1000:0.4g}mm".replace(".","p") \
            +f"_res{cell_size*1000:0.2g}mm".replace(".","p") \
            +f"_{"unspecified_high_dielectric"}".replace(".","p"), 
    )
    scene.add(geometry_view)



# Snapshot just a thin vertical slice of the domain through the probe
if do_snapshots:
    for i in range(1, n_frames+1):
        snapshot_ = gprMax.Snapshot(
            p1=(PML_width, y_mid-cell[1], PML_width,),
            p2=(dimensions[0]-PML_width, y_mid+cell[1], dimensions[2]-PML_width,),
            dl=(cell[0], cell[1], cell[2],),
            filename=f"xzslice_{i}",
            time=dt_snapshot*i,
            outputs=("Ex", "Ey", "Ez",)
        )
        scene.add(snapshot_)



###########################################################################


# Run the simulation
gprMax.run(scenes=[scene], n=1, outputfile=script_name,)

[cstarkjp]

If confirmed, this should be marked as a v4 bug

[MikaPhysics]

Is v4 reliable for simpler scenarios?