6. \(S\)-matrix extraction of a transmission line with the driven Maxwell solver

6.1. Requirements

6.1.1. Software components

• QTCAD

• Gmsh

6.1.2. Geometry file

• qtcad/examples/tutorials/meshes/tx_line.py

6.1.3. Python script

• qtcad/examples/tutorials/maxwell_driven_tx_line.py

6.1.4. References

• Getting Started

• Maxwell eigenmode extraction of a coplanar waveguide resonator with adaptive meshing

• Driven Maxwell Solver

6.2. Briefing

This tutorial is a simple example demonstrating extraction of the \(S\) parameters using the driven Maxwell Solver with adaptive mesh refinement. In this example, we extract the scattering matrix of a parallel plate transmission line with three ports. The results are compared against the analytical solution.

6.3. Geometry of the problem

The geometry of the transmission line is shown in Fig. 6.3.3. The transmission line has three ports: one on each end and one in the middle. Those ports will be modelled via the current sources. The following table shows the physical groups set up in the geometry file:

Gmsh physical group	Dimension	Purpose
"domain"	3D	The silicon dielectric medium separating the parallel plates
"y=0"	2D	The bottom plate of the transmission line, modelled as a PEC [1]
"y=height"	2D	The top plate of the transmission line, modelled as a PEC
"x=-lenx/2"	2D	Port 1 at the left end of the transmission line, modelled as a current source
"x=0"	2D	Port 2 in the middle of the transmission line, modelled as a current source
"x=lenx/2"	2D	Port 3 at the right end of the transmission line, modelled as a current source

Fig. 6.3.3 Geometry of the transmission line.

6.4. Setting up the device and running the driven Maxwell solver

The simulation process is analogous to the one presented in the Maxwell eigenmode extraction of a coplanar waveguide resonator with adaptive meshing tutorial for the Maxwell eigenmode solver. The following sections show the full simulation workflow for the driven Solver, highlighting notable differences with respect to the eigenmode extraction process.

6.4.1. Header, input parameters, and file paths

Similarly to Maxwell eigenmode extraction of a coplanar waveguide resonator with adaptive meshing, we start by importing the relevant modules and setting up the physical dimensions and file paths.

"""
Extraction of the S-matrix for a parallel plate transmission line with three ports.
This tutorial demonstrates the use of the driven Maxwell solver.
"""

import os
from time import time
from pathlib import Path
import numpy as np

from qtcad.device.maxwell_driven import Solver
from qtcad.device.maxwell_driven import SolverParams
from qtcad.device.device import Device
from qtcad.device.mesh3d import Mesh
from qtcad.device import materials as mt
from qtcad.device import constants as ct

# Scale in the Gmsh files
scale = 1e-3

# Dimensions of the transmission line in meters (defined in the mesh generation script).
lenx = 4.0 * scale
height = 0.5 * scale

# Directories and file paths
script_dir = Path(__file__).parent.resolve()
input_dir = script_dir / "meshes"
result_dir = script_dir / "output" / Path(__file__).stem
fpath_mesh = input_dir / "tx_line.msh"
fpath_xao = input_dir / "tx_line.xao"

# Check if the mesh and raw geometry files exist.
if not os.path.isfile(fpath_mesh) or not os.path.isfile(fpath_xao):
    raise Exception(
        f"Please run {input_dir}/tx_line.py to generate the mesh and raw geometry files."
    )

6.4.2. Loading the mesh and defining the device

Next, we load the mesh into the Device object, assign the material properties, and define the boundary conditions. In this case, we define PEC boundaries for the two conductors of the transmission line and current source boundaries for the three ports.

# Setup the device.
mesh = Mesh(scale, fpath_mesh)
dvc = Device(mesh)

# Extract the width in meters from the mesh.
width = np.amax(mesh.glob_nodes[:, 2])

# Assign media to regions.
dvc.new_region("domain", mt.Si)

# Assign Perfect Electric Conductor (PEC) boundaries.
dvc.new_pec_bnd("y=0")
dvc.new_pec_bnd("y=height")

# Define current sources at the three ports.
I0 = 1.0

# Port 1 at x = -lenx/2
dvc.new_current_source("x=-lenx/2", I0, "y", height, width=width)
# Port 2 at x = 0
dvc.new_current_source("x=0", I0, "y", height, width=width)
# Port 3 at x = lenx/2
dvc.new_current_source("x=lenx/2", I0, "y", height, width=width)

print(f"\nMesh file: {dvc.mesh.filename}")
print(f"Node count in the initial mesh: {dvc.mesh.node_number}")

6.4.3. Creating and running the driven Maxwell solver

Next, we set up the solver parameters that will be used by the driven Maxwell Solver.

# Setup the solver.
freq = 1e9

params = SolverParams()
params.freq = freq
# Split excitations to evaluate multi-port excitations.
params.split_excitations = True
params.name = Path(__file__).stem
params.output_dir = result_dir

# Enable adaptive meshing using S-parameter convergence.
params.converg_criterion = "S"
params.tol = 0.001

In the code above:

• The parameter freq specifies the driving frequency of the problem.

• We set the parameter split_excitations to True, which allows the solver to evaluate multi-port excitations to compute the \(S\)-matrix.

• The parameter convergence_criterion specifies whether the convergence of the adaptive meshing algorithm is determined using the error in energy (value "energy") or the error in the scattering parameters (value "S").

• The parameter tol specifies the adaptive meshing tolerance.

In this example, we choose "S" as the convergence criterion, since we are interested in extracting the \(S\) parameters and controlling the error associated with their calculation. The adaptive meshing algorithm is then considered to have converged when the \(S\) parameters agree with several previous iterations within the specified tolerance tol. Next, we use the solver parameters params to initialize and run the solver.

# Run the solver.
slv = Solver(dvc, params, geo_file=str(fpath_xao))
t0 = time()
slv.solve()
dt = time() - t0

print(f"Solution completed in {dt:.2f} s")
print(f"Node count in the final mesh: {dvc.mesh.node_number}")

6.5. Analytical comparison and \(S\)-matrix extraction

Once the solution is complete, we calculate the theoretical \(S\)-matrix for comparison. The methodology for finding these theoretical values can be found in [Poz12].

# Analytical Comparison
# See Pozar, David M. "Microwave engineering." Fourth Edition, John Wiley & Sons, Inc (2012).
L_pul = ct.mu0 * height / width
C_pul = mt.Si.eps * width / height
# Characteristic impedance (in this case, equals to the reference impedance)
Z0 = np.sqrt(L_pul / C_pul)

omega = 2 * np.pi * freq
beta = omega * np.sqrt(ct.mu0 * mt.Si.eps)

def Gamma(ZL):
    """Reflection coefficient of a load `ZL` on a transmission line with the
    characteric impedance `Z0`.
    """
    return (ZL - Z0) / (ZL + Z0)

def ztrans(ZL, len):
    """Input impedance of a transmission line of length `len` terminated with a
    load `ZL`.
    """
    return Z0 * (ZL + 1j * Z0 * np.tan(beta * len)) / (Z0 + 1j * ZL * np.tan(beta * len))

S11_analyt = Gamma(ztrans(Z0 / 2, lenx / 2))
S21_analyt = np.exp(-1j * beta * lenx / 2) * (1 + Gamma(Z0 / 2))
S31_analyt = (1 + Gamma(Z0 / 2)) * np.exp(-1j * beta * lenx)

S33_analyt = Gamma(ztrans(Z0 / 2, lenx / 2))
S23_analyt = np.exp(-1j * beta * lenx / 2) * (1 + Gamma(Z0 / 2))
S13_analyt = (1 + Gamma(Z0 / 2)) * np.exp(-1j * beta * lenx)

S22_analyt = Gamma(Z0 / 2)

S_analyt = np.array(
    [
        [S11_analyt, S21_analyt, S31_analyt],
        [S21_analyt, S22_analyt, S23_analyt],
        [S31_analyt, S23_analyt, S33_analyt],
    ]
)

Finally, we extract the computed \(S\)-matrix from the device using the method get_scattering_matrix, and compute the relative error compared to the analytical solution.

# Retrieve and output results.
ports = ["x=-lenx/2", "x=0", "x=lenx/2"]

# Retrieve the numerical S-matrix.
S_computed = dvc.get_scattering_matrix(ports, z0=Z0)

print("\nComputed S-matrix:")
print(np.array2string(S_computed, precision=4, separator=", "))

print("\nAnalytical S-matrix:")
print(np.array2string(S_analyt, precision=4, separator=", "))

# Compute and print the relative error.
err = np.linalg.norm(S_computed - S_analyt) / np.linalg.norm(S_analyt)
print(f"\nThe relative error compared to the analytical solution is {err:.4e}.")

The results are shown below:

Computed S-matrix:
[[-0.3196+9.4675e-02j,  0.6598-9.5661e-02j,  0.6392-1.8934e-01j],
 [ 0.6598-9.5661e-02j, -0.3333+2.0596e-06j,  0.6598-9.5661e-02j],
 [ 0.6392-1.8934e-01j,  0.6598-9.5661e-02j, -0.3196+9.4675e-02j]]

Analytical S-matrix:
[[-0.3196+0.0947j,  0.6598-0.0957j,  0.6392-0.1893j],
 [ 0.6598-0.0957j, -0.3333+0.j    ,  0.6598-0.0957j],
 [ 0.6392-0.1893j,  0.6598-0.0957j, -0.3196+0.0947j]]

The relative error compared to the analytical solution is 3.8204e-06.

6.6. Full code

__copyright__ = "Copyright 2022-2026, Nanoacademic Technologies Inc."

"""
Extraction of the S-matrix for a parallel plate transmission line with three ports.
This tutorial demonstrates the use of the driven Maxwell solver.
"""
import os
from time import time
from pathlib import Path
import numpy as np
from qtcad.device.maxwell_driven import Solver
from qtcad.device.maxwell_driven import SolverParams
from qtcad.device.device import Device
from qtcad.device.mesh3d import Mesh
from qtcad.device import materials as mt
from qtcad.device import constants as ct
# Scale in the Gmsh files
scale = 1e-3
# Dimensions of the transmission line in meters (defined in the mesh generation script).
lenx = 4.0 * scale
height = 0.5 * scale
# Directories and file paths
script_dir = Path(__file__).parent.resolve()
input_dir = script_dir / "meshes"
result_dir = script_dir / "output" / Path(__file__).stem
fpath_mesh = input_dir / "tx_line.msh"
fpath_xao = input_dir / "tx_line.xao"
# Check if the mesh and raw geometry files exist.
if not os.path.isfile(fpath_mesh) or not os.path.isfile(fpath_xao):
    raise Exception(
        f"Please run {input_dir}/tx_line.py to generate the mesh and raw geometry files."
    )

# Setup the device.
mesh = Mesh(scale, fpath_mesh)
dvc = Device(mesh)

# Extract the width in meters from the mesh.
width = np.amax(mesh.glob_nodes[:, 2])

# Assign media to regions.
dvc.new_region("domain", mt.Si)

# Assign Perfect Electric Conductor (PEC) boundaries.
dvc.new_pec_bnd("y=0")
dvc.new_pec_bnd("y=height")

# Define current sources at the three ports.
I0 = 1.0

# Port 1 at x = -lenx/2
dvc.new_current_source("x=-lenx/2", I0, "y", height, width=width)
# Port 2 at x = 0
dvc.new_current_source("x=0", I0, "y", height, width=width)
# Port 3 at x = lenx/2
dvc.new_current_source("x=lenx/2", I0, "y", height, width=width)

print(f"\nMesh file: {dvc.mesh.filename}")
print(f"Node count in the initial mesh: {dvc.mesh.node_number}")


# Setup the solver.
freq = 1e9

params = SolverParams()
params.freq = freq
# Split excitations to evaluate multi-port excitations.
params.split_excitations = True
params.name = Path(__file__).stem
params.output_dir = result_dir

# Enable adaptive meshing using S-parameter convergence.
params.converg_criterion = "S"
params.tol = 0.001

# Run the solver.
slv = Solver(dvc, params, geo_file=str(fpath_xao))
t0 = time()
slv.solve()
dt = time() - t0

print(f"Solution completed in {dt:.2f} s")
print(f"Node count in the final mesh: {dvc.mesh.node_number}")

# Analytical Comparison
# See Pozar, David M. "Microwave engineering." Fourth Edition, John Wiley & Sons, Inc (2012).
L_pul = ct.mu0 * height / width
C_pul = mt.Si.eps * width / height
# Characteristic impedance (in this case, equals to the reference impedance)
Z0 = np.sqrt(L_pul / C_pul)

omega = 2 * np.pi * freq
beta = omega * np.sqrt(ct.mu0 * mt.Si.eps)

def Gamma(ZL):
    """Reflection coefficient of a load `ZL` on a transmission line with the
    characteric impedance `Z0`.
    """
    return (ZL - Z0) / (ZL + Z0)

def ztrans(ZL, len):
    """Input impedance of a transmission line of length `len` terminated with a
    load `ZL`.
    """
    return Z0 * (ZL + 1j * Z0 * np.tan(beta * len)) / (Z0 + 1j * ZL * np.tan(beta * len))

S11_analyt = Gamma(ztrans(Z0 / 2, lenx / 2))
S21_analyt = np.exp(-1j * beta * lenx / 2) * (1 + Gamma(Z0 / 2))
S31_analyt = (1 + Gamma(Z0 / 2)) * np.exp(-1j * beta * lenx)

S33_analyt = Gamma(ztrans(Z0 / 2, lenx / 2))
S23_analyt = np.exp(-1j * beta * lenx / 2) * (1 + Gamma(Z0 / 2))
S13_analyt = (1 + Gamma(Z0 / 2)) * np.exp(-1j * beta * lenx)

S22_analyt = Gamma(Z0 / 2)

S_analyt = np.array(
    [
        [S11_analyt, S21_analyt, S31_analyt],
        [S21_analyt, S22_analyt, S23_analyt],
        [S31_analyt, S23_analyt, S33_analyt],
    ]
)


# Retrieve and output results.
ports = ["x=-lenx/2", "x=0", "x=lenx/2"]

# Retrieve the numerical S-matrix.
S_computed = dvc.get_scattering_matrix(ports, z0=Z0)

print("\nComputed S-matrix:")
print(np.array2string(S_computed, precision=4, separator=", "))

print("\nAnalytical S-matrix:")
print(np.array2string(S_analyt, precision=4, separator=", "))

# Compute and print the relative error.
err = np.linalg.norm(S_computed - S_analyt) / np.linalg.norm(S_analyt)
print(f"\nThe relative error compared to the analytical solution is {err:.4e}.")

[1]

Perfect electric conductor.