diff --git a/Simulation.py b/Simulation.py new file mode 100644 index 0000000..b4186d9 --- /dev/null +++ b/Simulation.py @@ -0,0 +1,119 @@ +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 \ No newline at end of file diff --git a/__pycache__/Simulation.cpython-312.pyc b/__pycache__/Simulation.cpython-312.pyc new file mode 100644 index 0000000..c221c85 Binary files /dev/null and b/__pycache__/Simulation.cpython-312.pyc differ diff --git a/__pycache__/helpers.cpython-312.pyc b/__pycache__/helpers.cpython-312.pyc new file mode 100644 index 0000000..da037f1 Binary files /dev/null and b/__pycache__/helpers.cpython-312.pyc differ diff --git a/helpers.py b/helpers.py new file mode 100644 index 0000000..7ec709e --- /dev/null +++ b/helpers.py @@ -0,0 +1,45 @@ +from Simulation import Simulation + +def build_trace_with_split(sim: 'Simulation', + 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))) + 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) + + sim.add_branch("GNDL","GNDR", R=Rsplit, L=Lsplit) + sim.add_shunt_C("GNDL","GNDR", Cgap) + + sim.add_branch("GNDL","GNDR", R=Rsplit_parallel) \ No newline at end of file diff --git a/main.py b/main.py index 4239943..8dac255 100644 --- a/main.py +++ b/main.py @@ -1,171 +1,42 @@ 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 +from Simulation import Simulation +from helpers import build_trace_with_split 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) + # 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 + dt = 1e-12 # 1 ps t_final = 50e-9 # 50 ns steps = int(t_final / dt) x = np.zeros(NN) @@ -181,7 +52,7 @@ if __name__ == "__main__": x = sim.step_BE(E, A, x, dt, drive) traces[:,n+1] = x[idx] - # local-ground spatial profile + # Node voltages referenced to their local ground def spatial_profile(n): v = np.zeros(N+1) for k in range(N+1): @@ -191,13 +62,11 @@ if __name__ == "__main__": 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 + # Time plots for 4 evenly spaced taps (include V0 and VN vs local ground) --- 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: