From b54dc1b1afaae9f5cb476050e485d6809a34348c Mon Sep 17 00:00:00 2001 From: judsonupchurch Date: Sun, 18 Jan 2026 15:07:40 -0600 Subject: [PATCH] Comment changes --- Simulation.py | 119 +++++++++++++++ __pycache__/Simulation.cpython-312.pyc | Bin 0 -> 6278 bytes __pycache__/helpers.cpython-312.pyc | Bin 0 -> 1806 bytes helpers.py | 45 ++++++ main.py | 193 ++++--------------------- 5 files changed, 195 insertions(+), 162 deletions(-) create mode 100644 Simulation.py create mode 100644 __pycache__/Simulation.cpython-312.pyc create mode 100644 __pycache__/helpers.cpython-312.pyc create mode 100644 helpers.py 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 0000000000000000000000000000000000000000..c221c85e4d09f8891af7e8abc8d5306e1047e3a2 GIT binary patch literal 6278 zcmb^#OKcm*b(hPfmfs~QQ6JlqDA}@1e`4F>o^x@f*E!2XLyD)!`}jbjmp7Ufy!K(17jf*X9G$LXNMiEoCBb3f%Ul+ z%W!ylHWG-3W6>a4_TuQR0GUYVra( z!Ax2k!!u10#IZ?4i`T4lN*j1P@ArK8t@SjpF=`grj<0B;k8<{h$e$zIsiDjV0Cmd#>!$1cFKn?-DCv}LDk{Iz2nC<%p*1jhJ? zViNfHtPl;BgbiazocN%!&{!cP5b-3ms>>W9uyOsdkT*eDKcqe@9^+#msVR}SaY9PwVfz9+#8PeyrW&F3hkY;@rHfTp5pV& zt^#A*U|Q}mEt#P_({`U}DX?zYBA}&fpT;PRJmxgWr#KB*p}I(93!p%b)L_Cf5}k_1 zu0_4bIK16IlHERoDm~^G`~sF8AWC8rkv4 z{e64h+V#M_BXcEpI(Phb$Lh{!*WAZeZO001{f*>eGBc8Acds(Lv1k+;@z$$WM-=kB zLY`B|{vbM1tspmx_4h?k0FYQs?M>EUJ=17(HAfcptb~)mT~6%-WuZArs=G;iXh@o1 z)spU)CW5VK+W&}ggEpQtCmD`NS|nQ7E76>B>J+Ts0aUdt#_~v;(0DD`FRV~Y%k{HP>QIb)jr8>zXnL>=l8%ZEc2w4|BC^P_pLMVnY|35WSA1xpu3p)|8BESGp zOo5phKFay5s*dbJ1FZ`_?eCX+WpdNmAk*JiZ7El3_%}_3hUSfit~{bE*qj@-)_b&T7X@!X|XTVGjc&TO(RnPaQ0ca`y~rmjD3 zq{LWpe+vz)xJ6eA-BPJgv0CG06{C|zP=}Igl$IbSOA)h-Ly8bDnUPd-F0va;8@SUfneDyB=lVXGSZ(jQ&m2}|t(O+!9?SYaAfHD%2^7L=$PYqT zZ9}CXSQ=0UCQ%3`jb({WA`_`9Aai5L5J$Fz%jcw)=OuM#@e?VLYqDu8mg`f*<=+#; zo4`e~F3EuVu_UcYThg9%0M9neEi0M8@Y#@I9)2h+a0cuOsh|M|W)rxK%F!kWR^?Gv zs(TDAD=`(LdRnvKguPk)>XO+?RP1gMSa4O9b1W4ztRlH2O7^$OrMc2M*zuRCt2-<9 zN6Cc7J1sdS-Kb=(Sl0xjkNobX&y6U=psb;)fW!wt;-F<%;#QC?Gj$yFnOw4Cr>m0| z?5LLqa>+VrRq;c^IaJbCehX7Yt6O>LRp6EeSsOzX4*H*x3eCkAX83WH`>mIJ3{H?# z509L&Vu;Qtbb=RRqHqxU*x3unfP5}&@`_<_4*p^I2NiNqAx|r0U)k`h+`nQjd49zN z*)}+Rh4d>Vkk1806dN?xU#eAq;TW3@M>zjJ|Z4!Lj3MA#j2Amx&|ch986#ld&a6m0fN#M8|A zhIOVoQhlkbsZ-0{nLy@z`p|NBC2}ZzaJgH7d2i%B-33omrZasZGne;tr3?k93#iYh zx|W$tf0|up5TY|RmfE{)&J3j;%Vrf3&}O{HY<0kwZS5s?w0+e6e)~g$G9O$#tA1ty z_UL!@mkySPVNAr3A{Es|a1NQbj(xCDeotzC0Rhi=E3d{&-x#Tv!@ zJ<00UK-ZJ1A~_^>yqr`yB~mwkN7ZELqEMoRaU}7iOEUaK3ing&VXfGs{wemTq-c;_ zkbbxmkP0FxCp}NzQ}-3MWNp??O3xW2LSVU4J?zz>?FH!u;P9@E@shPqqWeNsnM9X{ zsv>!$Ms0^=la42E)kNr))7WIQ)Fd@aPMpZtpk>9a=3%fz5?kXnT?;)Nq0|!s97>7H zu~0x)SA`A;LSR8LPVoyOPKMzHAsWZQ5iUXSpft`aCNU<&c}}53KCW1U@E8V9NKyF1 z+A;mX(S0#NCq{{UiDQk+uE{U$&C6{o?vt!SzN?wmxt;Z%iys+_}7Q^7Z`5*Vj%Cecg2acJKo_Go9C{UGB?$e}3t!VE*SrdB>1UecRZ)5>AI#V(C~m zm^-@G*aO0I)vvUsTQlOCYqxAFwC}`wtI$3GXHH2c_w} zh%O^vI+_O#Fn?+A3>;<#eWcI{IJijefl3E3){dS(qk(DO?^i5-|8$I-ji9#8@4q@5 zh^Qk1dhZ6?iVMQaXbluk)N3iuQ1_$mAq3FfQMe4?QvzofSP=zqKLiFS28izfxK2DY zSm<5fH=F4C%m4i*7x;?Uus)tMm-LySGdg(@$O~j#(f_JXj>6r!W wbyTLu5a%dj;}7|fIibquMS%kxy5NfFy=IFdN%HT6``^YmNj5$t5X4mf1^?nS4FCWD literal 0 HcmV?d00001 diff --git a/__pycache__/helpers.cpython-312.pyc b/__pycache__/helpers.cpython-312.pyc new file mode 100644 index 0000000000000000000000000000000000000000..da037f13eb6a216415d91eac7c5d8d34c5c7f651 GIT binary patch literal 1806 zcmaJ?-A@}w5Z^uTnZw2=B#{$EY6T@Bamq&!khW2PAiJchEF&CBE274|crWLJ^9Of! zBhVp%6tz$y^%|+l0jheufru3HShf8N`dHJ@ZL3D=OX(AWwyN^f-LoybDDGI=nVp@T z`Tb^h&)@lc9zf;R+bQJ~55RA1QWjUqIBcfILx2H>DoB9|071R;!wB=3yT?z!8CTG~ z$HIeV2Is@`%&Fi2`^-|UL1#vy6U)U zhUF@h+ZBhK4+~Ql0MG@)=kTWo#pW7kxY8cOWx%p`)(gNatvBbdQMc3&Om(r3tUCs0 z_#f48p&}dFqpMDj*eM+m?8eAI6?PTWm^0uOJxwMQ_t_$Dipw)gy;ElxT8z>V^K*?9 z+lz9Suo;a8U#W1vttJ-Ao1w>iHKLll%dtXr*K0)6%WTvep6Z^|z|11O6}118uc5l` z3+$klDa1ISkI)p z!je9f%VfJ^kx)lS(Gqzj*D;h9)frvUJH#~Z7|G_8lrrmlEV`1Cnw2!YGdn{9xATgM z<2g-CNb$Q$E*W<~X|1#-du943w1x9vUjV%DabzK4o-Cli&%Wbkw>kP*&+^0pfc(pg za98kt+OX7MPFqcu%Y1i7cxm577oMRyXH%enPFt;Cwy(6?9~RDZ?z>sG=t9)&wUCYG zqXl$smxdNSZo7?nf`7V!|Ag7S{N4fJ_%~>cz$vTg0k%ZD$?jf1v3h4|Y$LqgvOT&z z{p7}z^PAB{j=}IR`6-BQtJS<~x7b6rxE`|m?D~zd?bxD!M|gdor`h&B48h34h#9s5 zHs#b-yN^`44Tsdpm~) 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: