EDA Skillpack

"""Coupled copper-sheet resistor network. Grid and contact approximations are reported."""
from pathlib import Path
import sys,os,json,pickle,time,math
D=Path(__file__).resolve().parent;sys.path.insert(0,str(D/'vendor'))
import numpy as np,shapely
from shapely.geometry import Point,box
from scipy import sparse
from scipy.sparse.csgraph import connected_components
from scipy.sparse.linalg import cg
import pyamg
rho=1.724e-8;layers=['F.Cu','In1.Cu','In2.Cu','B.Cu'];th=np.array([.035,.0152,.0152,.035])*1e-3;zz=np.array([.0175,.253,1.3332,1.5683])*1e-3
geo=pickle.load(open(D/'copper.pkl','rb'));h=float(sys.argv[1]) if len(sys.argv)>1 else .15;plating=float(sys.argv[2]) if len(sys.argv)>2 else .025
x=np.arange(104+h/2,168,h);y=np.arange(58+h/2,132,h);xx,yy=np.meshgrid(x,y);xy=np.column_stack([xx.ravel(),yy.ravel()]);shape=xx.shape;nxy=len(xy)
padmap={p['name']:p for p in geo['pads']};results=[]
def network(net):
 masks=[];starts=[];ends=[];cond=[];grids=[];count=0
 for li,l in enumerate(layers):
  g=geo['copper'].get((net,l),Point(0,0).buffer(0));shapely.prepare(g);mask=shapely.contains_xy(g,xx,yy);grid=np.full(shape,-1,dtype=np.int32);grid[mask]=np.arange(count,count+mask.sum());count+=int(mask.sum());grids.append(grid);masks.append(mask)
  for axis in [0,1]:
   a,b=(grid[:-1,:],grid[1:,:]) if axis==0 else (grid[:,:-1],grid[:,1:]);hit=(a>=0)&(b>=0);iy,ix=np.nonzero(hit);j=iy*shape[1]+ix;k=j+(shape[1] if axis==0 else 1)
   # Exact segment containment prevents hopping over a narrow slot or drill void.
   valid=shapely.covers(g,shapely.linestrings(np.stack([xy[j],xy[k]],axis=1)))
   starts.extend(a[hit][valid]);ends.extend(b[hit][valid]);cond.extend(np.full(valid.sum(),th[li]/rho))
 viaedges=0
 for v in geo['bars']:
  if v['net']!=net:continue
  dist=np.hypot(xx-v['x'],yy-v['y']);ann=(dist>=v['drill']/2)&(dist<v['diameter']/2)
  area=np.pi*((v['drill']/2+plating)**2-(v['drill']/2)**2)*1e-6
  for l in range(3):
   hit=ann & (grids[l]>=0)&(grids[l+1]>=0);num=hit.sum()
   if num:
    starts.extend(grids[l][hit]);ends.extend(grids[l+1][hit]);cond.extend(np.full(num,area/(rho*(zz[l+1]-zz[l]))/num));viaedges+=num
 u=np.array(starts,dtype=np.int32);v=np.array(ends,dtype=np.int32);c=np.array(cond);A=sparse.coo_matrix((np.r_[c,c,-c,-c],(np.r_[u,v,u,v],np.r_[u,v,v,u])),shape=(count,count)).tocsr();_,labels=connected_components(A,directed=False)
 print('mesh',net,h,count,len(c),flush=True)
 return A,grids,labels,u,v,c

def solve(net,name,src,dst,N):
 A,grids,labels,u,v,c=N
 def ids(refs):
  ids=[]
  for ref in refs:
   p=padmap[ref];mask=shapely.contains_xy(p['g'],xx,yy)
   # Source/sink on the actual front pad surface. Plated barrel handled separately.
   a=grids[0][mask];ids.extend(a[a>=0])
  return np.unique(ids).astype(np.int32)
 a,b=ids(src),ids(dst)
 if not len(a) or not len(b):raise RuntimeError(('terminal unresolved',name,len(a),len(b)))
 lab=labels[a[0]];a=a[labels[a]==lab];b=b[labels[b]==lab]
 if not len(b):raise RuntimeError(('grid disconnected',name,h))
 active=np.flatnonzero(labels==lab);free=np.setdiff1d(active,np.r_[a,b]);fixed=np.r_[a,b];U=np.zeros(A.shape[0]);U[a]=1
 M=A[free][:,free];rhs=-A[free][:,a].sum(axis=1).A.ravel();ml=pyamg.smoothed_aggregation_solver(M);sol,info=cg(M,rhs,M=ml.aspreconditioner(),rtol=2e-9,maxiter=2000);U[free]=sol
 if info:raise RuntimeError(('solver',name,info))
 cur=A@U;I=cur[a].sum();drop=U[u]-U[v];power=c*drop**2;R=1/I;energy=power.sum()/I;balance=abs(cur[a].sum()+cur[b].sum())/I
 out=dict(case=name,net=net,source=src,sink=dst,resistanceOhm=R,lossAt58_5A_W=58.5**2*R,powerBalanceRelative=abs(energy-1),currentBalanceRelative=balance,gridMm=h,platingMm=plating,nodes=len(active),freeResidualMaxA=float(np.max(np.abs(cur[free]))))
 print(json.dumps(out),flush=True);results.append(out)
 # Unit-current heating at each layer cell, including half of each via/edge loss at either end.
 Q=np.bincount(np.r_[u,v],weights=np.r_[power,power]/(2*I*I),minlength=A.shape[0]);maps=[]
 for grid in grids:
  q=np.zeros(shape);valid=grid>=0;q[valid]=Q[grid[valid]];maps.append(q)
 np.savez_compressed(D/f'dc-{name}-{h}-{plating}.npz',heatWPerA2=np.array(maps),x=x,y=y)

if __name__=='__main__':
 specs={'+VBAT':[(f'vbat-{q}',['MC1.1'],[f'{q}.5']) for q in ['Q1','Q3','Q5']], 'GND':[(f'return-{q}',[f'{q}.{i}' for i in [1,2,3]],['R32.2']) for q in ['Q2','Q4','Q6']], 'GND_OUT':[('battery-return',['R32.1'],['MC3.1'])]}
 for phase,hi,lo,contact in [('A','Q1','Q2','MC5.1'),('B','Q3','Q4','MC9.1'),('C','Q5','Q6','MC11.1')]:specs['/DRV_SH'+phase]=[(f'phase-{phase}-high',[f'{hi}.{i}' for i in [1,2,3]],[contact]),(f'phase-{phase}-low',[f'{lo}.5'],[contact])]
 for net,cases in specs.items():
  if os.environ.get('DC_NET') and net!=os.environ['DC_NET']:continue
  N=network(net)
  for name,src,dst in cases:solve(net,name,src,dst,N)
 (D/f'dc-results-{h}-{plating}.json').write_text(json.dumps(results,indent=2))