Using lumped elements in Tidy3D simulations¶

# Now we create the simulation and include the definition of the lumped element.
sim = td.Simulation(
    center=(0, strip_height/2, 0),
    size=(simx, simy, simz),
    boundary_spec=boundary_spec,
    grid_spec=td.GridSpec.auto(min_steps_per_wvl=80, wavelength=3000),
    lumped_elements=[lumped_element],
    structures=[substrate, strip, strip2],
    # Add monitor around lumped element to enable computation of voltage, current, and impedance.
    monitors=[prop_mon, lumped_element.to_monitor(freqs=freqs)],
    # The source will be added by the TerminalComponentModeler
    sources=[],
    run_time= 50*(tline_length/td.C_0),
    plot_length_units="mm", # This option will make plots default to units of millimeters.
)
# The next step is to setup a LumpedPort for the TerminalComponentModeler, which is needed for extracting scattering parameters
# as well as providing a source for the simulation.
port1 = smatrix.LumpedPort(
    voltage_axis=1,
    size=(0,strip_height,strip_width),
    center=(-tline_length/2, strip_height/2, 0),
    num_grid_cells=5,
    name="port_1",
)
# Finally, we combine the base simulation with the LumpedPort in the TerminalComponentModeler, which will enable use
# to compute scattering parameters as well as retrieve any other desired data from the simulation.
modeler = smatrix.TerminalComponentModeler(
    simulation = sim,
    freqs = freqs,
    ports = [port1],
    verbose = True,
    remove_dc_component=False, # Include DC component for more accuracy at low frequencies
)
# Before running the solver, we plot the created structure along its propagation direction. We also plot the cross-section 
# of the structure at the terminated end of the transmission line, which includes the location of the lumped element.
sim_temp = list(modeler.sim_dict.values())[0]
f, (ax1, ax2) = plt.subplots(2, 1, tight_layout=True, figsize=(15, 12))
sim_temp.plot(z = 0*mill, ax=ax1)
sim_temp.plot(x = tline_length/2, ax=ax2)
sim_temp.plot_grid(x = tline_length/2, ax=ax2, vlim=[-2*strip_width, 2*strip_width], hlim=[-strip_height, 2*strip_height])
# Add arrows that indicate locations of the Port and the load
ax1.annotate("Port 1", xy=(-tline_length/2 - 10*mill, strip_height + 10*mill), xytext=(-tline_length/1.2,strip_width*10),
            arrowprops=dict(facecolor='black', shrink=0.05))
ax1.annotate("Load", xy=(tline_length/2 + 10*mill, strip_height + 10*mill), xytext=(tline_length/1.5,strip_width*10),
            arrowprops=dict(facecolor='black', shrink=0.05))
plt.show()

s_matrix = modeler.run()
batch_data = modeler.batch_data

sim_data = batch_data["smatrix_port_1"]

f, (ax1, ax2) = plt.subplots(2, 1, tight_layout=True, figsize=(6, 6))
sim_data.plot_field(
    field_monitor_name="prop", field_name="Ey", val="abs", f=freq0, ax=ax1
)
sim_data.plot_field(
    field_monitor_name="prop", field_name="Hz", val="abs", f=freq0, ax=ax2
)
plt.show()

# The eigenmode associated with the tranmission line is computed at the central frequency 25 GHz
mode_spec = td.ModeSpec(num_modes=1, target_neff=2.2)
mode_plane = td.Box(center = sim.center, size = (0, simy, simz))

mode_solver = mode.ModeSolver(simulation=sim_temp, plane=mode_plane, mode_spec=mode_spec, freqs=[freq0])
mode_data = web.run(mode_solver)

# Convert the effective refractive index into an effective permittivity
er_eff = (np.real(mode_data.n_eff.values[0]))**2

# Tidy3D uses the physics convention for time-harmonic fields,
# so we take the conjugate in order to switch to the engineering convention.
load_impedance = np.conj(lumped_element.impedance(freqs))
Gamma = compute_Gamma_from_transmission_line_theory(load_impedance, characteristic_impedance, er_eff)

# Get S11 data
freq = s_matrix.f / 1e9
S11_RL_term = s_matrix.isel(port_out=0, port_in=0).values.flatten()

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(11, 3))
fig.suptitle("Transmission line terminated by RL series network")
ax1.plot(
    s_matrix.f / 1e9,
    20 * np.log10(np.abs(S11_RL_term)),
    "--b",
    label="Tidy3D"
)
ax1.plot(
    freqs / 1e9,
    20 * np.log10(np.abs(Gamma)),
    "-k",
    label="Transmission line theory"
)
ax1.set_xlabel("Frequency (GHz)")
ax1.set_ylabel(r"$|S_{11}|$ (dB)")
ax1.set_xlim([0, 50])
ax1.grid()
ax1.legend()

ax2.plot(
    s_matrix.f / 1e9,
    -np.angle(S11_RL_term), #Follow engineering convention exp(j omega t)
    "--b",
    label="Tidy3D"
)
ax2.plot(
    freqs / 1e9,
    np.angle(Gamma),
    "-k",
    label="Transmission line theory"
)
ax2.set_xlabel("Frequency (GHz)")
ax2.set_ylabel(r"$\angle~S_{11}$ (rad)")
ax2.set_xlim([0, 50])
ax2.grid()
ax2.legend()
plt.show()

# Get the monitor data associated the lumped element
vi_mon_data = sim_data[lumped_element.monitor_name]
# Setup voltage and current path integrals automatically around the lumped element. 
voltage_integral, current_integral = mw.path_integrals_from_lumped_element(lumped_element, sim_data.simulation.grid)
# Compute the voltage and current using the path integral tools and the frequency-domain field data.
V = voltage_integral.compute_voltage(vi_mon_data)
I = current_integral.compute_current(vi_mon_data)
# Apply Ohm's law to retrieve impedance
impedance = V/I
resistance = np.real(impedance).values
reactance = -np.imag(impedance).values # Change of sign, since Tidy3D uses the physics convention for time-harmonic fields (exp(-j omega t))
# Here we use the convenience methods from the lumped element to compute the exact impedance of the modeled network.
ideal_impedance = lumped_element.impedance(freqs=freqs)
ideal_resistance = np.real(ideal_impedance)
ideal_reactance = -np.imag(ideal_impedance) # Change of sign, since Tidy3D uses the physics convention for time-harmonic fields (exp(-j omega t))

# Plot the results for visual comparison
f, (ax1, ax2) = plt.subplots(2, 1, tight_layout=True, figsize=(8, 6))
f.suptitle("RL series network")
ax1.plot(freqs / 1e9, resistance, label="Computed")
ax1.plot(freqs / 1e9, ideal_resistance,"--", label="Desired")
ax1.set_ylabel(r"Resistance ($\Omega$)")
ax1.set_xlabel("Frequency (GHz)")
ax1.set_xlim(0, 50)
ax1.set_ylim(0, 100)
ax1.legend(loc="lower right")
ax2.plot(freqs / 1e9, reactance, label="Computed")
ax2.plot(freqs / 1e9, ideal_reactance,"--", label="Desired")
ax2.set_ylabel(r"Reactance ($\Omega$)")
ax2.set_xlabel("Frequency (GHz)")
ax2.set_xlim(0, 50)
ax2.legend(loc="lower right")
plt.show()

# Example network from Sec 15.9.6 [1] 
L1 = 1e-9
L2 = 1.5e-9
C1 = 0.2e-12
C2 = 0.2e-12
R1 = 10
R2 = 250
R3 = 50

# If you have installed lcapy, the transfer function will be automatically generated using a simple description of the network.
# Otherwise, the precomputed transfer function coefficients will be used.
try:
    from lcapy import R, C, L, s
    print("Generating transfer function using 'lcapy' ...")
    # Create the 4th order network using 'lcapy'
    # The network is described by adding each of the elements in series and parallel
    network = L(L1) + (C(C1) | (R(R1) + L(L2) + (R(R2) | C(C2)))) + R(R3)
    # Can also get `lcapy` to draw the network
    network.draw(layout="horizontal", scale=0.4)
    lca_a = network.Y(s).N.coeffs()
    lca_b = network.Y(s).D.coeffs()
    # Reverse order of coefficients
    a = tuple(np.flip(lca_a.fval))
    b = tuple(np.flip(lca_b.fval))

except ModuleNotFoundError:
    print("Skipping steps involving 'lcapy' package ...")
    # Manually derived admittance transfer function for network
    a = [
        1,
        C1 * R1 + C1 * R2 + C2 * R2,
        C1 * C2 * R1 * R2 + C1 * L2,
        C1 * C2 * L2 * R2,
        0,
    ]
    b = [
        R1 + R2 + R3,
        C1 * R1 * R3 + C1 * R2 * R3 + C2 * R1 * R2 + C2 * R2 * R3 + L1 + L2,
        C1 * C2 * R1 * R2 * R3 + C1 * L1 * R1 + C1 * L1 * R2
        + C1 * L2 * R3 + C2 * L1 * R2 + C2 * L2 * R2,
        C1 * C2 * L1 * R1 * R2 + C1 * C2 * L2 * R2 * R3 + C1 * L1 * L2,
        C1 * C2 * L1 * L2 * R2,
    ]

# With the admittance transfer function in hand, we can construct the AdmittanceNetwork
fourth_order_network = td.AdmittanceNetwork(a=a, b=b)
# We create an updated copy of the lumped element from the previous simulation.
lumped_element = lumped_element.updated_copy(network=fourth_order_network)
# Instead of repeating all of the previous steps to build the Simulation and TerminalComponentModeler,
# we create an updated copy of the modeler which includes the newly created lumped element.
modeler = modeler.updated_copy(lumped_elements=[lumped_element], path="simulation")

s_matrix = modeler.run(modeler.path_dir)
batch_data = modeler.batch_data

# Tidy3D uses the physics convention for time-harmonic fields,
# so we take the conjugate in order to switch to the engineering convention.
load_impedance = np.conj(lumped_element.impedance(freqs))
Gamma = compute_Gamma_from_transmission_line_theory(load_impedance, characteristic_impedance, er_eff)

# Get the S11 data
freq = s_matrix.f / 1e9
S11_ntwk = s_matrix.isel(port_out=0, port_in=0).values.flatten()

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(11, 3))
fig.suptitle("Transmission line terminated by the seven-element 4th-order network")
ax1.plot(
    s_matrix.f / 1e9,
    20 * np.log10(np.abs(S11_ntwk)),
    "--b",
    label="Tidy3D"
)
ax1.plot(
    freqs / 1e9,
    20 * np.log10(np.abs(Gamma)),
    "-k",
    label="Transmission line theory"
)
ax1.set_xlabel("Frequency (GHz)")
ax1.set_ylabel(r"$|S_{11}|$ (dB)")
ax1.set_xlim([0, 50])
ax1.grid()
ax1.legend()
ax2.plot(
    s_matrix.f / 1e9,
    -np.angle(S11_ntwk), #Follow engineering convention exp(j omega t)
    "--b",
    label="Tidy3D"
)
ax2.plot(
    freqs / 1e9,
    np.angle(Gamma),
    "-k",
    label="Transmission line theory"
)
ax2.set_xlabel("Frequency (GHz)")
ax2.set_ylabel(r"$\angle~S_{11}$ (rad)")
ax2.set_xlim([0, 50])
ax2.grid()
ax2.legend()
plt.show()

# Retrieve the simulation data
sim_data = batch_data["smatrix_port_1"]
# Get the monitor data associated the lumped element
vi_mon_data = sim_data[lumped_element.monitor_name]
# Setup voltage and current path integrals automatically around the lumped element. 
voltage_integral, current_integral = mw.path_integrals_from_lumped_element(lumped_element, sim_data.simulation.grid)
# Compute the voltage and current using the path integral tools and the frequency-domain field data.
V = voltage_integral.compute_voltage(vi_mon_data)
I = current_integral.compute_current(vi_mon_data)
# Apply Ohm's law to retrieve impedance
impedance = V/I
resistance = np.real(impedance).values
reactance = -np.imag(impedance).values # Change of sign, since Tidy3D uses the physics convention for time-harmonic fields (exp(-j omega t))
# Here we use the convenience methods from the lumped element to compute the intended impedance of the modeled network.
ideal_impedance = lumped_element.impedance(freqs=freqs)
ideal_resistance = np.real(ideal_impedance)
ideal_reactance = -np.imag(ideal_impedance) # Change of sign, since Tidy3D uses the physics convention for time-harmonic fields (exp(-j omega t))

# Plot the results for visual comparison
f, (ax1, ax2) = plt.subplots(2, 1, tight_layout=True, figsize=(8, 6))
f.suptitle("Seven-element 4th-order network")
ax1.plot(freqs / 1e9, resistance, label="Computed")
ax1.plot(freqs / 1e9, ideal_resistance,"--", label="Desired")
ax1.set_ylabel(r"Resistance ($\Omega$)")
ax1.set_xlabel("Frequency (GHz)")
ax1.set_xlim(0, 50)
ax1.set_ylim(0, 250)
ax1.legend(loc="upper right")
ax2.plot(freqs / 1e9, reactance, label="Computed")
ax2.plot(freqs / 1e9, ideal_reactance,"--", label="Desired")
ax2.set_ylabel(r"Reactance ($\Omega$)")
ax2.set_xlabel("Frequency (GHz)")
ax2.set_xlim(0, 50)
ax2.legend(loc="lower right")
plt.show()

# We create an updated copy of the lumped element from the previous simulation,
# where Distribute Type is turned 'off' and the width is set to 0.
lumped_element = lumped_element.updated_copy(
    dist_type="off", size=(0, strip_height, 0), network=RLC
)
# Use a convenience method for estimating the parasitic inductance and capacitance 
# associated with the wire connections, which are used in the lumped element.
Lp, Cp = lumped_element.estimate_parasitic_elements(
    grid=modeler.sim_dict["smatrix_port_1"].grid
)
print(
    f"The parasitic inductance of the wire connections is {Lp*1e9:.2f} nH, "
    f"while the parasitic capacitance is {Cp*1e15:.2f} fF."
)

# We create an updated copy of the lumped element that compensates for the parasitic inductance by subtracting it.
lumped_element_comp = lumped_element.updated_copy(path="network", inductance=Lval - Lp)
modeler = modeler.updated_copy(lumped_elements=[lumped_element], path="simulation")
s_matrix_not_comp = modeler.run(modeler.path_dir)
# Run again with the compensated version
modeler = modeler.updated_copy(lumped_elements=[lumped_element_comp], path="simulation")
s_matrix_comp = modeler.run(modeler.path_dir)

# Tidy3D uses the physics convention for time-harmonic fields,
# so we take the conjugate in order to switch to the engineering convention.
# Get S11 data
freq = s_matrix.f / 1e9
S11_not_comp = s_matrix_not_comp.isel(port_out=0, port_in=0).values.flatten()
S11_comp = s_matrix_comp.isel(port_out=0, port_in=0).values.flatten()

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(11, 3))
fig.suptitle("Transmission line terminated by RL series network")
ax1.plot(
    s_matrix.f / 1e9,
    20 * np.log10(np.abs(S11_RL_term)),
    "-b",
    label="Distributed"
)
ax1.plot(
    s_matrix.f / 1e9,
    20 * np.log10(np.abs(S11_not_comp)),
    "--r",
    label="Single-cell"
)
ax1.plot(
    s_matrix.f / 1e9,
    20 * np.log10(np.abs(S11_comp)),
    "--g",
    label="Compensated-single-cell"
)
ax1.set_xlabel("Frequency (GHz)")
ax1.set_ylabel(r"$|S_{11}|$ (dB)")
ax1.set_xlim([0, 50])
ax1.grid()
ax1.legend()

ax2.plot(
    s_matrix.f / 1e9,
    -np.angle(S11_RL_term), #Follow engineering convention exp(j omega t)
    "-b",
    label="Distributed"
)
ax2.plot(
    s_matrix.f / 1e9,
    -np.angle(S11_not_comp), #Follow engineering convention exp(j omega t)
    "--r",
    label="Single-cell"
)
ax2.plot(
    s_matrix.f / 1e9,
    -np.angle(S11_comp), #Follow engineering convention exp(j omega t)
    "--g",
    label="Compensated-single-cell"
)
ax2.set_xlabel("Frequency (GHz)")
ax2.set_ylabel(r"$\angle~S_{11}$ (rad)")
ax2.set_xlim([0, 50])
ax2.grid()
ax2.legend()
plt.show()