255 lines
8.8 KiB
Python
255 lines
8.8 KiB
Python
import numpy as np
|
|
|
|
class Simulation:
|
|
def __init__(self):
|
|
self.node_ids = {}
|
|
self.next_node_id = 0
|
|
|
|
self.branches = []
|
|
self.shunt_caps = []
|
|
|
|
self.branch_count = 0
|
|
self.cstate_count = 0
|
|
|
|
self.dirichlet = []
|
|
|
|
def add_node(self, label:str="") -> int:
|
|
if label in self.node_ids:
|
|
raise ValueError(f"Node {label} already exists")
|
|
idx = self.next_node_id
|
|
self.node_ids[label] = idx
|
|
self.next_node_id += 1
|
|
return idx
|
|
|
|
def _get_node_idx(self, label: str) -> int:
|
|
if label not in self.node_ids:
|
|
raise KeyError(f"Unknown node '{label}'")
|
|
return self.node_ids[label]
|
|
|
|
|
|
def add_branch(self, node1:str, node2:str, R:float=0.0, L:float=0.0, C:float=0.0, label:str="") -> int:
|
|
branch_idx = self.branch_count
|
|
self.branch_count += 1
|
|
|
|
if label == "":
|
|
label = f"b_{branch_idx}"
|
|
|
|
a, b = self._get_node_idx(node1), self._get_node_idx(node2)
|
|
|
|
c_idx = None
|
|
if C and C > 0.0:
|
|
c_idx = self.cstate_count
|
|
self.cstate_count += 1
|
|
self.branches.append((a, b, R, L, C, branch_idx, c_idx, label))
|
|
|
|
return branch_idx, c_idx
|
|
|
|
def add_shunt_C(self, node1: str, node2: str, C: float):
|
|
a, b = self._get_node_idx(node1), self._get_node_idx(node2)
|
|
self.shunt_caps.append((a, b, float(C)))
|
|
|
|
def build_matrices(self):
|
|
Nv = self.next_node_id
|
|
Ni = self.branch_count
|
|
Nc = self.cstate_count
|
|
N = Nv + Ni + Nc
|
|
|
|
E = np.zeros((N,N), dtype=float)
|
|
A = np.zeros((N,N), dtype=float)
|
|
|
|
for (a, b, R, L, C, branch_id, cap_id, _) in self.branches:
|
|
i = Nv + branch_id # Index of this branch's current
|
|
|
|
# KCL for the nodes that have this current. Current leaves a and enters b
|
|
A[a, i] += 1.0
|
|
A[b, i] -= 1.0
|
|
|
|
# KVL for branch: v(a)-v(b) - R*i - L di/dt - v_c = 0
|
|
A[i, a] += 1.0
|
|
A[i, b] -= 1.0
|
|
if R: A[i, i] += -R
|
|
if L: E[i, i] += L
|
|
if C and cap_id is not None:
|
|
vc = Nv + Ni + cap_id
|
|
A[i, vc] += 1.0 # subtract v_c in KVL
|
|
E[vc, vc] += C # C * d v_c/dt = i
|
|
A[vc, i] -= 1.0
|
|
|
|
# Shunt capacitors into E(v,v) block
|
|
for (a, b, C) in self.shunt_caps:
|
|
if a != b:
|
|
E[a, a] -= C; E[a, b] += C
|
|
E[b, a] += C; E[b, b] -= C
|
|
|
|
return E, A, N
|
|
|
|
|
|
def step_BE(self, E, A, x_n, h, dirichlet_now):
|
|
# dirichlet_now: dict label->value at t_{n+1}
|
|
Nv = self.next_node_id
|
|
Ni = self.branch_count
|
|
Nc = self.cstate_count
|
|
N = Nv + Ni + Nc
|
|
|
|
# fixed voltage indices and values
|
|
fixed_v_idx = np.array([self._get_node_idx(lbl) for lbl in dirichlet_now.keys()], dtype=int)
|
|
v_d = np.array([float(dirichlet_now[lbl]) for lbl in dirichlet_now.keys()], dtype=float)
|
|
|
|
# free sets
|
|
all_v_idx = np.arange(Nv, dtype=int)
|
|
free_v_idx = np.array(sorted(set(all_v_idx) - set(fixed_v_idx)), dtype=int)
|
|
free_i_idx = np.arange(Nv, Nv+Ni+Nc, dtype=int)
|
|
free_idx = np.concatenate([free_v_idx, free_i_idx])
|
|
|
|
# Partitioned matrices
|
|
Ef_f = E[np.ix_(free_idx, free_idx)]
|
|
Af_f = A[np.ix_(free_idx, free_idx)]
|
|
Ef_d = E[np.ix_(free_idx, fixed_v_idx)]
|
|
Af_d = A[np.ix_(free_idx, fixed_v_idx)]
|
|
|
|
# Left matrix and RHS for BE
|
|
LHS = Ef_f - h * Af_f
|
|
# print(np.linalg.cond(LHS))
|
|
rhs = (E[np.ix_(free_idx, np.arange(N))] @ x_n) - (Ef_d - h * Af_d) @ v_d
|
|
|
|
x_free_next = np.linalg.solve(LHS, rhs)
|
|
x_next = np.zeros_like(x_n)
|
|
x_next[free_idx] = x_free_next
|
|
x_next[fixed_v_idx]= v_d
|
|
return x_next
|
|
|
|
def build_trace_with_split(sim, N=20, Rseg=0.05, Lseg=8e-9, Csh0=1e-12, Cshi=1e-14,
|
|
Rsrc=10.0, Rs_gnd=120.0, Rload=1e5,
|
|
Rsplit=0.3, Lsplit=50e-9, Cgap=100e-12,Rsplit_parallel=0.1):
|
|
# nodes: Vs, V0..V{N}, GNDL, GNDR
|
|
sim.add_node("Vs")
|
|
for k in range(N+1):
|
|
sim.add_node(f"V{k}")
|
|
sim.add_node("GNDL"); sim.add_node("GNDR")
|
|
|
|
# source
|
|
sim.add_branch("Vs","V0", R=Rsrc)
|
|
sim.add_branch("Vs","GNDL", R=Rs_gnd)
|
|
|
|
# series ladder
|
|
for k in range(N):
|
|
Lk = Lseg if (k==0 or k==N-1) else (Lseg*(N/(N-1))) # tweak if you want end-sections different
|
|
sim.add_branch(f"V{k}", f"V{k+1}", R=Rseg, L=Lk)
|
|
|
|
# shunt caps to local grounds
|
|
sim.add_shunt_C("V0","GNDL", Csh0)
|
|
for k in range(1, N):
|
|
g = "GNDL" if k <= N//2 else "GNDR"
|
|
sim.add_shunt_C(f"V{k}", g, Cshi)
|
|
sim.add_shunt_C(f"V{N}", "GNDR", Csh0)
|
|
|
|
# load
|
|
sim.add_branch(f"V{N}", "GNDR", R=Rload)
|
|
|
|
# ground tie: lossy via in parallel with gap capacitance
|
|
sim.add_branch("GNDL","GNDR", R=Rsplit, L=Lsplit) # via
|
|
sim.add_shunt_C("GNDL","GNDR", Cgap) # gap
|
|
|
|
sim.add_branch("GNDL","GNDR", R=Rsplit_parallel) # Parallel across the gap
|
|
|
|
if __name__ == "__main__":
|
|
import matplotlib.pyplot as plt
|
|
|
|
N = 50
|
|
sim = Simulation()
|
|
# build_trace_with_split(sim, N=N, Rseg=0.02, Lseg=8e-9, Csh0=1e-12, Cshi=10e-15,
|
|
# Rsrc=50.0, Rs_gnd=120.0, Rload=120,
|
|
# Rsplit=0.5, Lsplit=50e-9, Cgap=50e-12, Rsplit_parallel=0.1)
|
|
build_trace_with_split(sim, N=N, Rseg=0.02, Lseg=8e-9, Csh0=1e-12, Cshi=1e-12,
|
|
Rsrc=50.0, Rs_gnd=120.0, Rload=120,
|
|
Rsplit=0.0, Lsplit=0.0, Cgap=0.0, Rsplit_parallel=0.1)
|
|
|
|
E,A,NN = sim.build_matrices()
|
|
dt = 1e-12 # 1 ps
|
|
t_final = 50e-9 # 50 ns
|
|
steps = int(t_final / dt)
|
|
x = np.zeros(NN)
|
|
drive = {"Vs":3.0, "GNDL":0.0}
|
|
|
|
nodes = ["Vs"]+[f"V{k}" for k in range(N+1)]+["GNDL","GNDR"]
|
|
idx = [sim._get_node_idx(s) for s in nodes]
|
|
tvec = np.arange(steps+1)*dt
|
|
traces = np.zeros((len(nodes), steps+1))
|
|
traces[:,0] = x[idx]
|
|
|
|
for n in range(steps):
|
|
x = sim.step_BE(E, A, x, dt, drive)
|
|
traces[:,n+1] = x[idx]
|
|
|
|
# local-ground spatial profile
|
|
def spatial_profile(n):
|
|
v = np.zeros(N+1)
|
|
for k in range(N+1):
|
|
if k <= N//2:
|
|
v[k] = traces[nodes.index(f"V{k}"), n] - traces[nodes.index("GNDL"), n]
|
|
else:
|
|
v[k] = traces[nodes.index(f"V{k}"), n] - traces[nodes.index("GNDR"), n]
|
|
return v
|
|
|
|
# time plots for 4 evenly spaced taps (include V0 and VN vs local ground) ---
|
|
# pick indices: 0, ~N/3, ~2N/3, N
|
|
tap_idx = sorted(set([0, N//3, (2*N)//3, N]))
|
|
tap_labels = [f"V{k}" for k in tap_idx]
|
|
|
|
def v_local(k, n):
|
|
# local ground reference by side
|
|
if k <= N//2:
|
|
return traces[nodes.index(f"V{k}"), n] - traces[nodes.index("GNDL"), n]
|
|
else:
|
|
return traces[nodes.index(f"V{k}"), n] - traces[nodes.index("GNDR"), n]
|
|
|
|
Vt = {k: np.array([v_local(k, n) for n in range(traces.shape[1])]) for k in tap_idx}
|
|
|
|
plt.figure(figsize=(8,4))
|
|
for k in tap_idx:
|
|
plt.plot(tvec*1e9, Vt[k], label=f"{'V0' if k==0 else ('V_N' if k==N else f'V{k}')}")
|
|
plt.xlabel("Time [ns]")
|
|
plt.ylabel("Voltage vs local ground [V]")
|
|
plt.title("Selected nodes over time")
|
|
plt.xlim(0, 10)
|
|
plt.grid(True)
|
|
plt.legend()
|
|
plt.tight_layout()
|
|
plt.savefig("reflection_wave.png")
|
|
plt.show()
|
|
|
|
# GIF
|
|
from PIL import Image
|
|
import io, numpy as np, matplotlib.pyplot as plt
|
|
|
|
xpos = np.arange(N+1)
|
|
Nt = traces.shape[1]
|
|
target_fps = 5
|
|
max_frames = 300
|
|
stride = max(1, Nt // max_frames)
|
|
frame_idx = np.arange(0, Nt, stride)
|
|
vals = np.vstack([spatial_profile(i) for i in frame_idx])
|
|
ymin, ymax = 1.1*vals.min(), 1.1*vals.max()
|
|
duration_ms = int(1000/target_fps)
|
|
|
|
frames = []
|
|
for i in frame_idx:
|
|
v = spatial_profile(i)
|
|
fig, ax = plt.subplots(figsize=(6,3.5))
|
|
ax.plot(xpos, v, marker=None)
|
|
ax.set_xlim(0, N)
|
|
ax.set_ylim(ymin, ymax)
|
|
ax.set_xlabel("Position index (0 … N)")
|
|
ax.set_ylabel("V vs local ground [V]")
|
|
ax.set_title(f"t = {tvec[i]*1e9:.2f} ns")
|
|
ax.grid(True)
|
|
buf = io.BytesIO()
|
|
plt.tight_layout()
|
|
fig.savefig(buf, format='png', dpi=120)
|
|
plt.close(fig)
|
|
buf.seek(0)
|
|
frames.append(Image.open(buf).convert("P"))
|
|
|
|
out_path = "reflection_wave_N24.gif"
|
|
frames[0].save(out_path, save_all=True, append_images=frames[1:], duration=duration_ms, loop=0, optimize=False)
|
|
print("Saved", out_path) |