kopia lustrzana https://gitlab.com/gerbolyze/gerbonara
Move coil stuff to separate repo
jaseg
rodzic
a35125b123
commit
9624e46147
File diff suppressed because one or more lines are too long
318
coil_gen.py
318
coil_gen.py
|
|
@ -1,318 +0,0 @@
|
|||
#!/usr/bin/env python3
|
||||
|
||||
import subprocess
|
||||
import sys
|
||||
import os
|
||||
from math import *
|
||||
from pathlib import Path
|
||||
from itertools import cycle
|
||||
|
||||
from gerbonara.cad.kicad import pcb as kicad_pcb
|
||||
from gerbonara.cad.kicad import footprints as kicad_fp
|
||||
from gerbonara.cad.kicad import graphical_primitives as kicad_gr
|
||||
import click
|
||||
|
||||
|
||||
__version__ = '1.0.0'
|
||||
|
||||
|
||||
def point_line_distance(p, l1, l2):
|
||||
x0, y0 = p
|
||||
x1, y1 = l1
|
||||
x2, y2 = l2
|
||||
# https://en.wikipedia.org/wiki/Distance_from_a_point_to_a_line
|
||||
return abs((x2-x1)*(y1-y0) - (x1-x0)*(y2-y1)) / sqrt((x2-x1)**2 + (y2-y1)**2)
|
||||
|
||||
def line_line_intersection(l1, l2):
|
||||
p1, p2 = l1
|
||||
p3, p4 = l2
|
||||
x1, y1 = p1
|
||||
x2, y2 = p2
|
||||
x3, y3 = p3
|
||||
x4, y4 = p4
|
||||
|
||||
# https://en.wikipedia.org/wiki/Line%E2%80%93line_intersection
|
||||
px = ((x1*y2-y1*x2)*(x3-x4)-(x1-x2)*(x3*y4-y3*x4))/((x1-x2)*(y3-y4)-(y1-y2)*(x3-x4))
|
||||
py = ((x1*y2-y1*x2)*(y3-y4)-(y1-y2)*(x3*y4-y3*x4))/((x1-x2)*(y3-y4)-(y1-y2)*(x3-x4))
|
||||
return px, py
|
||||
|
||||
def angle_between_vectors(va, vb):
|
||||
angle = atan2(vb[1], vb[0]) - atan2(va[1], va[0])
|
||||
if angle < 0:
|
||||
angle += 2*pi
|
||||
return angle
|
||||
|
||||
|
||||
@click.command()
|
||||
@click.argument('infile', required=False, type=click.Path(exists=True, dir_okay=False, path_type=Path))
|
||||
@click.argument('outfile', required=False, type=click.Path(writable=True, dir_okay=False, path_type=Path))
|
||||
@click.option('--footprint-name', help="Name for the generated footprint. Default: Output file name sans extension.")
|
||||
@click.option('--target-layer', default='F.Cu', help="Target KiCad layer for the generated footprint. Default: F.Cu.")
|
||||
@click.option('--polygon', type=int, default=0, help="Use n'th polygon instead of first one. 0-based index.")
|
||||
@click.option('--start-angle', type=float, default=0, help='Angle for the start at the outermost layer of the spiral in degree')
|
||||
@click.option('--windings', type=int, default=5, help='Number of windings to generate')
|
||||
@click.option('--trace-width', type=float, default=0.15)
|
||||
@click.option('--clearance', type=float, default=0.15)
|
||||
@click.option('--clipboard/--no-clipboard', help='Use clipboard integration (requires wl-clipboard)')
|
||||
@click.option('--counter-clockwise/--clockwise', help='Direction of generated spiral. Default: clockwise when wound from the inside.')
|
||||
def generate(infile, outfile, polygon, start_angle, windings, trace_width, clearance, footprint_name, target_layer, clipboard, counter_clockwise):
|
||||
if 'WAYLAND_DISPLAY' in os.environ:
|
||||
copy, paste, cliputil = ['wl-copy'], ['wl-paste'], 'xclip'
|
||||
else:
|
||||
copy, paste, cliputil = ['xclip', '-i', '-sel', 'clipboard'], ['xclip', '-o', '-sel' 'clipboard'], 'wl-clipboard'
|
||||
|
||||
if clipboard:
|
||||
try:
|
||||
proc = subprocess.run(paste, capture_output=True, text=True, check=True)
|
||||
except FileNotFoundError:
|
||||
print(f'Error: --clipboard requires the {copy[0]} and {paste[0]} utilities from {cliputil} to be installed.', file=sys.stderr)
|
||||
board = kicad_pcb.Board.load(proc.stdout)
|
||||
elif not infile:
|
||||
board = kicad_pcb.Board.load(sys.stdin.read())
|
||||
else:
|
||||
board = kicad_pcb.Board.open(infile)
|
||||
|
||||
objs = [obj for obj in board.objects() if isinstance(obj, kicad_gr.Polygon)]
|
||||
print(f'Found {len(objs)} polygon(s).', file=sys.stderr)
|
||||
poly = objs[polygon]
|
||||
xy = [(pt.x, pt.y) for pt in poly.pts.xy]
|
||||
|
||||
if counter_clockwise:
|
||||
xy = [(-x, y) for x, y in xy]
|
||||
|
||||
segments = list(zip(xy, xy[1:] + xy[:1]))
|
||||
|
||||
# normalize orientation, make xy counter-clockwise
|
||||
if sum((x2 - x1) * (y2 + y1) for (x1, y1), (x2, y2) in segments) < 0:
|
||||
print(f'Reversing polygon direction.', file=sys.stderr)
|
||||
xy = xy[::-1]
|
||||
segments = list(zip(xy, xy[1:] + xy[:1]))
|
||||
|
||||
vbx, vby = min(x for x, y in xy), min(y for x, y, in xy)
|
||||
vbw, vbh = max(x for x, y in xy), max(y for x, y, in xy)
|
||||
vbw, vbh = vbw-vbx, vbh-vby
|
||||
|
||||
vbx -= 5
|
||||
vby -= 5
|
||||
vbw += 10
|
||||
vbh += 10
|
||||
|
||||
cx, cy = 0, 0
|
||||
ls = 0
|
||||
for (x1, y1), (x2, y2) in segments:
|
||||
l = dist((x1, y1), (x2, y2))
|
||||
cx += x1*l/2 + x2*l/2
|
||||
cy += y1*l/2 + y2*l/2
|
||||
ls += l
|
||||
cx /= ls
|
||||
cy /= ls
|
||||
|
||||
segment_angles = [(atan2(y1-cy, x1-cx) - atan2(y2-cy, x2-cx) + 2*pi) % (2*pi) for (x1, y1), (x2, y2) in segments]
|
||||
angle_strs = [f'{degrees(a):.2f}' for a in segment_angles]
|
||||
|
||||
segment_heights = [point_line_distance((cx, cy), (x1, y1), (x2, y2)) for (x1, y1), (x2, y2) in segments]
|
||||
segment_foo = list(zip(segment_heights, segments))
|
||||
|
||||
midpoints = []
|
||||
for h, ((x1, y1), (x2, y2)) in segment_foo:
|
||||
xb = (x1 + x2) / 2
|
||||
yb = (y1 + y2) / 2
|
||||
midpoints.append((xb, yb))
|
||||
|
||||
normals = []
|
||||
for h, ((x1, y1), (x2, y2)) in segment_foo:
|
||||
d12 = dist((x1, y1), (x2, y2))
|
||||
dx = x2 - x1
|
||||
dy = y2 - y1
|
||||
normals.append((-dy/d12, dx/d12))
|
||||
|
||||
smallest_radius = min(segment_heights)
|
||||
#trace_radius = smallest_radius - stop_radius
|
||||
trace_radius = smallest_radius
|
||||
|
||||
segment_foo = list(zip(segment_heights, segments, segment_angles, midpoints, normals))
|
||||
|
||||
dbg_lines1, dbg_lines2 = [], []
|
||||
spiral_points = []
|
||||
dr_tot = trace_width/2
|
||||
for n in range(windings):
|
||||
for (ha, (pa1, pa2), aa, ma, na), (hb, (pb1, pb2), ab, mb, nb) in zip(segment_foo[-1:] + segment_foo[:-1], segment_foo):
|
||||
pitch = clearance + trace_width
|
||||
dr_tot_a = dr_tot
|
||||
dr_tot_b = n * pitch + trace_width/2
|
||||
|
||||
xma, yma = ma
|
||||
xna, yna = na
|
||||
xmb, ymb = mb
|
||||
xnb, ynb = nb
|
||||
|
||||
xa1, ya1 = pa1
|
||||
xa2, ya2 = pa2
|
||||
xb1, yb1 = pb1
|
||||
xb2, yb2 = pb2
|
||||
|
||||
dma = dist(pa2, ma)
|
||||
dmb = dist(pb1, mb)
|
||||
|
||||
x_cons_a, y_cons_a = p_cons_a = line_line_intersection((pa2, (cx, cy)), (ma, (xma-xna, yma-yna)))
|
||||
d_cons_a = dist(p_cons_a, ma)
|
||||
qa = dma * dr_tot_a / d_cons_a
|
||||
dra = hypot(qa, dr_tot_a)
|
||||
|
||||
nrax = (xa2 - cx) / dist((cx, cy), pa2)
|
||||
nray = (ya2 - cy) / dist((cx, cy), pa2)
|
||||
|
||||
xea = xa2 - nrax*dra
|
||||
yea = ya2 - nray*dra
|
||||
|
||||
x_cons_b, y_cons_b = p_cons_b = line_line_intersection((pb1, (cx, cy)), (mb, (xmb-xnb, ymb-ynb)))
|
||||
d_cons_b = dist(p_cons_b, mb)
|
||||
qb = dmb * dr_tot_b / d_cons_b
|
||||
drb = hypot(qb, dr_tot_b)
|
||||
|
||||
nrbx = (xb1 - cx) / dist((cx, cy), pb1)
|
||||
nrby = (yb1 - cy) / dist((cx, cy), pb1)
|
||||
|
||||
xeb = xb1 - nrbx*drb
|
||||
yeb = yb1 - nrby*drb
|
||||
|
||||
xsa = xma - xna*dr_tot_a
|
||||
ysa = yma - yna*dr_tot_a
|
||||
|
||||
xsb = xmb - xnb*dr_tot_b
|
||||
ysb = ymb - ynb*dr_tot_b
|
||||
|
||||
l1 = (xsa, ysa), (xea, yea)
|
||||
l2 = (xsb, ysb), (xeb, yeb)
|
||||
|
||||
dbg_lines1.append(l1)
|
||||
dbg_lines2.append(l2)
|
||||
|
||||
pic = line_line_intersection(l1, l2)
|
||||
spiral_points.append(pic)
|
||||
|
||||
dr_tot = dr_tot_b
|
||||
|
||||
#spiral_points = []
|
||||
#r_now = 0
|
||||
#for winding in range(num_windings):
|
||||
# for angle, ((x1, y1), (x2, y2)) in zip(segment_angles, segments):
|
||||
# angle_frac = angle/(2*pi)
|
||||
# d_r = angle_frac * (clearance + trace_width)
|
||||
# r_pt = dist((cx, cy), (x1, y1)) * (num_windings - winding) / num_windings
|
||||
#
|
||||
# x1, y1 = x1-cx, y1-cy
|
||||
# x2, y2 = x2-cx, y2-cy
|
||||
# l1, l2 = hypot(x1, y1), hypot(x2, y2)
|
||||
# x1, y1 = x1/l1, y1/l1
|
||||
# x2, y2 = x2/l2, y2/l2
|
||||
#
|
||||
# r_now += d_r
|
||||
# spiral_points.append((cx + x1*r_pt, cy + y1*r_pt))
|
||||
|
||||
out = [spiral_points[0]]
|
||||
ndrop = 0
|
||||
for i, (pa, pb, pc) in enumerate(zip(spiral_points, spiral_points[1:], spiral_points[2:])):
|
||||
xa, ya = pa
|
||||
xb, yb = pb
|
||||
xc, yc = pc
|
||||
if ndrop:
|
||||
ndrop -= 1
|
||||
continue
|
||||
|
||||
angle = angle_between_vectors((xa-xb, ya-yb), (xc-xb, yc-yb))
|
||||
if angle > pi:
|
||||
ndrop += 1
|
||||
for pd, pe in zip(spiral_points[i+2:], spiral_points[i+3:]):
|
||||
xd, yd = pd
|
||||
xe, ye = pe
|
||||
angle = angle_between_vectors((xa-xb, ya-yb), (xe-xd, ye-yd))
|
||||
if angle > pi:
|
||||
ndrop += 1
|
||||
else:
|
||||
out.append(line_line_intersection((pa, pb), (pd, pe)))
|
||||
break
|
||||
|
||||
else:
|
||||
out.append(pb)
|
||||
spiral_points = out
|
||||
|
||||
path_d = ' '.join([f'M {xy[0][0]} {xy[0][1]}', *[f'L {x} {y}' for x, y in xy[1:]], 'Z'])
|
||||
path_d2 = ' '.join(f'M {cx} {cy} L {x} {y}' for x, y in xy)
|
||||
path_d3 = ' '.join([f'M {spiral_points[0][0]} {spiral_points[0][1]}', *[f'L {x} {y}' for x, y in spiral_points[1:]]])
|
||||
|
||||
with open('/tmp/test.svg', 'w') as f:
|
||||
f.write('<?xml version="1.0" standalone="no"?>\n')
|
||||
f.write('<!DOCTYPE svg PUBLIC "-//W3C//DTD SVG 1.1//EN" "http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd">\n')
|
||||
f.write(f'<svg version="1.1" width="200mm" height="200mm" viewBox="{vbx} {vby} {vbw} {vbh}" xmlns="http://www.w3.org/2000/svg" xmlns:xlink="http://www.w3.org/1999/xlink">>\n')
|
||||
f.write(f'<path fill="none" stroke="#303030" stroke-width="0.05" d="{path_d}"/>\n')
|
||||
f.write(f'<path fill="none" stroke="#a0a0a0" stroke-width="0.05" d="{path_d2}"/>\n')
|
||||
f.write(f'<path fill="none" stroke="#ff00ff" opacity="0.5" stroke-width="{trace_width}" d="{path_d3}"/>\n')
|
||||
|
||||
for (x1, y1), (x2, y2) in dbg_lines1:
|
||||
f.write(f'<path fill="none" stroke="#ff0000" opacity="0.2" stroke-width="0.05" d="M {x1} {y1} L {x2} {y2}"/>')
|
||||
|
||||
for (x1, y1), (x2, y2) in dbg_lines2:
|
||||
f.write(f'<path fill="none" stroke="#0000ff" opacity="0.2" stroke-width="0.05" d="M {x1} {y1} L {x2} {y2}"/>')
|
||||
|
||||
for x, y in midpoints:
|
||||
f.write(f'<path fill="none" stroke="#a0a0ff" stroke-width="0.05" d="M {cx} {cy} L {x} {y}"/>')
|
||||
f.write(f'<circle r="0.1" fill="blue" stroke="none" cx="{x}" cy="{y}"/>\n')
|
||||
|
||||
f.write(f'<circle r="0.1" fill="red" stroke="none" cx="{cx}" cy="{cy}"/>\n')
|
||||
f.write('</svg>\n')
|
||||
|
||||
if counter_clockwise:
|
||||
spiral_points = [(-x, y) for x, y in spiral_points]
|
||||
|
||||
fp_lines = [
|
||||
kicad_fp.Line(
|
||||
start=kicad_fp.XYCoord(x=x1, y=y1),
|
||||
end=kicad_fp.XYCoord(x=x2, y=y2),
|
||||
layer=target_layer,
|
||||
stroke=kicad_fp.Stroke(width=trace_width))
|
||||
for (x1, y1), (x2, y2) in zip(spiral_points, spiral_points[1:])]
|
||||
|
||||
make_pad = lambda num, x, y: kicad_fp.Pad(
|
||||
number=str(num),
|
||||
type=kicad_fp.Atom.smd,
|
||||
shape=kicad_fp.Atom.circle,
|
||||
at=kicad_fp.AtPos(x=x, y=y),
|
||||
size=kicad_fp.XYCoord(x=trace_width, y=trace_width),
|
||||
layers=[target_layer],
|
||||
clearance=clearance,
|
||||
zone_connect=0,
|
||||
)
|
||||
|
||||
if footprint_name:
|
||||
name = footprint_name
|
||||
elif outfile:
|
||||
name = outfile.stem,
|
||||
else:
|
||||
name = 'generated_coil'
|
||||
|
||||
fp = kicad_fp.Footprint(
|
||||
name=name,
|
||||
generator=kicad_fp.Atom('GerbonaraCoilGenV1'),
|
||||
layer='F.Cu',
|
||||
descr=f"{windings} winding coil footprint generated by gerbonara'c Coil generator, version {__version__}",
|
||||
clearance=clearance,
|
||||
zone_connect=0,
|
||||
lines=fp_lines,
|
||||
pads=[make_pad(1, *spiral_points[0]), make_pad(2, *spiral_points[-1])],
|
||||
)
|
||||
|
||||
if clipboard:
|
||||
try:
|
||||
print(f'Running {copy[0]}.')
|
||||
proc = subprocess.Popen(copy, stdin=subprocess.PIPE, text=True)
|
||||
proc.communicate(fp.serialize())
|
||||
except FileNotFoundError:
|
||||
print(f'Error: --clipboard requires the {copy[0]} and {paste[0]} utilities from {cliputil} to be installed.', file=sys.stderr)
|
||||
elif not outfile:
|
||||
print(fp.serialize())
|
||||
else:
|
||||
fp.write(outfile)
|
||||
|
||||
if __name__ == '__main__':
|
||||
generate()
|
||||
|
|
@ -1,38 +0,0 @@
|
|||
air:
|
||||
Density: 1.1885 # 20°C
|
||||
Electric Conductivity: 0.0
|
||||
Heat Capacity: 1006.4 # 20°C
|
||||
Heat Conductivity: 0.025873 # 20°C
|
||||
Relative Permeability: 1
|
||||
Relative Permittivity: 1
|
||||
ro4003c:
|
||||
Density: 1790 # 23°C
|
||||
Relative Permeability: 1
|
||||
Relative Permittivity: 3.55
|
||||
fr4:
|
||||
Density: 1850 # 23°C
|
||||
Relative Permeability: 1
|
||||
Relative Permittivity: 4.4
|
||||
Heat Conductivity: 0.81 # in-plane
|
||||
ideal:
|
||||
Relative Permittivity: 1
|
||||
copper:
|
||||
Density: 8960.0 # 0°C
|
||||
Electric Conductivity: 59600000 # 20°C
|
||||
Emissivity: 0.012 # 327°C
|
||||
Heat Capacity: 415.0 # 200°C
|
||||
Heat Conductivity: 401.0 # 0°C
|
||||
Relative Permeability: 1
|
||||
Relative Permittivity: 1
|
||||
steel_1.4541:
|
||||
Density: 7900.0 # 20°C
|
||||
Electric Conductivity: 1370
|
||||
Emissivity: 0.111 # 200°C
|
||||
Heat Capacity: 470.0 # 20°C
|
||||
Heat Conductivity: 15.0 # 20°C
|
||||
Relative Permeability: 1
|
||||
Relative Permittivity: 1
|
||||
water:
|
||||
Density: 1000.0
|
||||
Heat Capacity: 4182.0
|
||||
Heat Conductivity: 0.6
|
||||
|
|
@ -1,14 +0,0 @@
|
|||
3D_steady:
|
||||
Mesh Levels: 1
|
||||
Max Output Level: 7
|
||||
Coordinate System: Cartesian
|
||||
Coordinate Mapping(3): 1 2 3
|
||||
Simulation Type: Steady state
|
||||
Steady State Max Iterations: 1
|
||||
Output Intervals: 1
|
||||
Timestepping Method: BDF
|
||||
Simulation Timing: True
|
||||
BDF Order: 1
|
||||
Solver Input File: case.sif
|
||||
Post File: case.vtu
|
||||
Output File: case.result
|
||||
|
|
@ -1,55 +0,0 @@
|
|||
Static_Current_Conduction:
|
||||
Equation: Static Current Conduction
|
||||
Variable: Potential
|
||||
Variable DOFs: 1
|
||||
Procedure: '"StatCurrentSolve" "StatCurrentSolver"'
|
||||
Calculate Volume Current: True
|
||||
Calculate Joule Heating: False
|
||||
Optimize Bandwidth: True
|
||||
Nonlinear System Max Iterations: 1
|
||||
Linear System Solver: Iterative
|
||||
Linear System Iterative Method: CG
|
||||
Linear System Max Iterations: 10000
|
||||
Linear System Convergence Tolerance: 1.0e-10
|
||||
Linear System Preconditioning: ILU3
|
||||
Linear System ILUT Tolerance: 1.0e-3
|
||||
Linear System Abort Not Converged: False
|
||||
Linear System Residual Output: 20
|
||||
Linear System Precondition Recompute: 1
|
||||
|
||||
Magneto_Dynamics:
|
||||
Equation: MGDynamics
|
||||
Procedure: '"MagnetoDynamics" "WhitneyAVSolver"'
|
||||
Newton-Raphson Iteration: True
|
||||
Nonlinear System Max Iterations: 1
|
||||
Fix Input Current Density: True
|
||||
Linear System Solver: Iterative
|
||||
Linear System Preconditioning: none
|
||||
Linear System Convergence Tolerance: 1e-8
|
||||
Linear System Abort Not Converged: False
|
||||
Linear System Residual Output: 20
|
||||
Linear System Max Iterations: 5000
|
||||
Linear System Iterative Method: TFQMR
|
||||
BiCGStabL Polynomial Degree: 4
|
||||
"Jfix: Linear System Max Iterations": 5000
|
||||
Idrs Parameter: 6
|
||||
Solver Timing: True
|
||||
|
||||
Magneto_Dynamics_Calculations:
|
||||
Equation: MGDynamicsCalc
|
||||
Procedure: '"MagnetoDynamics" "MagnetoDynamicsCalcFields"'
|
||||
Linear System Symmetric: True
|
||||
Calculate Magnetic Field Strength: True
|
||||
Calculate JxB: False
|
||||
Calculate Current Density: True
|
||||
Calculate Electric Field: True
|
||||
Calculate Nodal Fields: False
|
||||
Calculate Elemental Fields: True
|
||||
Linear System Solver: Iterative
|
||||
Linear System Preconditioning: ILU0
|
||||
Linear System Abort Not Converged: False
|
||||
Linear System Residual Output: 0
|
||||
Linear System Max Iterations: 5000
|
||||
Linear System Iterative Method: CG
|
||||
Linear System Convergence Tolerance: 1.0e-8
|
||||
|
||||
|
|
@ -1,595 +0,0 @@
|
|||
#!/usr/bin/env python3
|
||||
|
||||
import math
|
||||
from pathlib import Path
|
||||
import multiprocessing
|
||||
import re
|
||||
import tempfile
|
||||
import fnmatch
|
||||
import shutil
|
||||
import numpy as np
|
||||
|
||||
from pyelmer import elmer
|
||||
import subprocess_tee
|
||||
import click
|
||||
from scipy import constants
|
||||
|
||||
def enumerate_mesh_bodies(msh_file):
|
||||
with open(msh_file, 'r') as f:
|
||||
for line in f:
|
||||
if line.startswith('$PhysicalNames'):
|
||||
break
|
||||
else:
|
||||
raise ValueError('No physcial bodies found in mesh file.')
|
||||
|
||||
_num_names = next(f)
|
||||
|
||||
for line in f:
|
||||
if line.startswith('$EndPhysicalNames'):
|
||||
break
|
||||
|
||||
dim, _, line = line.strip().partition(' ')
|
||||
tag, _, name = line.partition(' ')
|
||||
yield name.strip().strip('"'), (int(dim), int(tag))
|
||||
|
||||
# https://en.wikipedia.org/wiki/Metric_prefix
|
||||
SI_PREFIX = 'QRYZEPTGMk mµnpfazyrq'
|
||||
|
||||
def format_si(value, unit='', fractional_digits=1):
|
||||
mag = int(math.log10(abs(value))//3)
|
||||
value /= 1000**mag
|
||||
prefix = SI_PREFIX[SI_PREFIX.find(' ') - mag].strip()
|
||||
value = f'{{:.{fractional_digits}f}}'.format(value)
|
||||
return f'{value} {prefix}{unit}'
|
||||
|
||||
|
||||
INPUT_EXT_MAP = {
|
||||
'.grd': 1,
|
||||
'.mesh*': 2,
|
||||
'.ep': 3,
|
||||
'.ansys': 4,
|
||||
'.inp': 5,
|
||||
'.fil': 6,
|
||||
'.FDNEUT': 7,
|
||||
'.unv': 8,
|
||||
'.mphtxt': 9,
|
||||
'.dat': 10,
|
||||
'.node': 11,
|
||||
'.ele': 11,
|
||||
'.mesh': 12,
|
||||
'.msh': 14,
|
||||
'.ep.i': 15,
|
||||
'.2dm': 16}
|
||||
|
||||
OUTPUT_EXT_MAP = {
|
||||
'.grd': 1,
|
||||
'.mesh*': 2,
|
||||
'.ep': 3,
|
||||
'.msh': 4,
|
||||
'.vtu': 5}
|
||||
|
||||
def elmer_grid(infile, outfile=None, intype=None, outtype=None, cwd=None, stdout_log=None, stderr_log=None, **kwargs):
|
||||
infile = Path(infile)
|
||||
if outfile is not None:
|
||||
outfile = Path(outfile)
|
||||
|
||||
if intype is None:
|
||||
intype = str(INPUT_EXT_MAP[infile.suffix])
|
||||
|
||||
if outtype is None:
|
||||
if outfile is not None and outfile.suffix:
|
||||
outtype = str(OUTPUT_EXT_MAP[outfile.suffix])
|
||||
else:
|
||||
outtype = '2'
|
||||
|
||||
if outfile is not None:
|
||||
kwargs['out'] = str(outfile)
|
||||
|
||||
args = ['ElmerGrid', intype, outtype, str(infile)]
|
||||
for key, value in kwargs.items():
|
||||
args.append(f'-{key}')
|
||||
if isinstance(value, (tuple, list)):
|
||||
args.extend(str(v) for v in value)
|
||||
else:
|
||||
args.append(str(value))
|
||||
|
||||
result = subprocess_tee.run(args, cwd=cwd, check=True)
|
||||
if stdout_log:
|
||||
Path(stdout_log).write_text(result.stdout or '')
|
||||
if stderr_log:
|
||||
Path(stderr_log).write_text(result.stderr or '')
|
||||
|
||||
def elmer_solver(cwd, stdout_log=None, stderr_log=None):
|
||||
result = subprocess_tee.run(['ElmerSolver'], cwd=cwd, check=True)
|
||||
if stdout_log:
|
||||
Path(stdout_log).write_text(result.stdout or '')
|
||||
if stderr_log:
|
||||
Path(stderr_log).write_text(result.stderr or '')
|
||||
return result
|
||||
|
||||
|
||||
@click.group()
|
||||
def cli():
|
||||
pass
|
||||
|
||||
|
||||
@cli.command()
|
||||
@click.option('-d', '--sim-dir', type=click.Path(dir_okay=True, file_okay=False, path_type=Path))
|
||||
@click.option('-o', '--output', type=click.Path(dir_okay=False, writable=True, path_type=Path), help='Capacitance matrix output file')
|
||||
@click.argument('mesh_file', type=click.Path(dir_okay=False, path_type=Path))
|
||||
def capacitance_matrix(mesh_file, sim_dir, output):
|
||||
physical = dict(enumerate_mesh_bodies(mesh_file))
|
||||
if sim_dir is not None:
|
||||
sim_dir = Path(sim_dir)
|
||||
sim_dir.mkdir(exist_ok=True)
|
||||
|
||||
sim = elmer.load_simulation('3D_steady', 'coil_parasitics_sim.yml')
|
||||
mesh_dir = '.'
|
||||
mesh_fn = 'mesh'
|
||||
sim.header['Mesh DB'] = f'"{mesh_dir}" "{mesh_fn}"'
|
||||
sim.constants.update({
|
||||
'Permittivity of Vacuum': str(constants.epsilon_0),
|
||||
'Gravity(4)': f'0 -1 0 {constants.g}',
|
||||
'Boltzmann Constant': str(constants.Boltzmann),
|
||||
'Unit Charge': str(constants.elementary_charge)})
|
||||
|
||||
air = elmer.load_material('air', sim, 'coil_parasitics_materials.yml')
|
||||
fr4 = elmer.load_material('fr4', sim, 'coil_parasitics_materials.yml')
|
||||
|
||||
solver_electrostatic = elmer.load_solver('Electrostatics_Capacitance', sim, 'coil_parasitics_solvers.yml')
|
||||
solver_electrostatic.data['Potential Difference'] = '1.0'
|
||||
eqn = elmer.Equation(sim, 'main', [solver_electrostatic])
|
||||
|
||||
bdy_sub = elmer.Body(sim, 'substrate', [physical['substrate'][1]])
|
||||
bdy_sub.material = fr4
|
||||
bdy_sub.equation = eqn
|
||||
|
||||
bdy_ab = elmer.Body(sim, 'airbox', [physical['airbox'][1]])
|
||||
bdy_ab.material = air
|
||||
bdy_ab.equation = eqn
|
||||
|
||||
max_num = -1
|
||||
|
||||
# boundaries
|
||||
for name, identity in physical.items():
|
||||
if (m := re.fullmatch(r'trace([0-9]+)', name)):
|
||||
num = int(m.group(1))
|
||||
max_num = max(num, max_num)
|
||||
|
||||
bndry_m2 = elmer.Boundary(sim, name, [identity[1]])
|
||||
bndry_m2.data['Capacitance Body'] = str(num)
|
||||
|
||||
if (tr := physical.get('trace')):
|
||||
bndry_m2 = elmer.Boundary(sim, 'trace', [tr[1]])
|
||||
bndry_m2.data['Capacitance Body'] = f'{max_num+1}'
|
||||
|
||||
boundary_airbox = elmer.Boundary(sim, 'FarField', [physical['airbox_surface'][1]])
|
||||
boundary_airbox.data['Electric Infinity BC'] = 'True'
|
||||
|
||||
with tempfile.TemporaryDirectory() as tmpdir:
|
||||
tmpdir = sim_dir if sim_dir else Path(tmpdir)
|
||||
|
||||
sim.write_startinfo(tmpdir)
|
||||
sim.write_sif(tmpdir)
|
||||
# Convert mesh from gmsh to elemer formats. Also scale it from 1 unit = 1 mm to 1 unit = 1 m (SI units)
|
||||
elmer_grid(mesh_file.absolute(), 'mesh', cwd=tmpdir, scale=[1e-3, 1e-3, 1e-3],
|
||||
stdout_log=(tmpdir / 'ElmerGrid_stdout.log'),
|
||||
stderr_log=(tmpdir / 'ElmerGrid_stderr.log'))
|
||||
elmer_solver(tmpdir,
|
||||
stdout_log=(tmpdir / 'ElmerSolver_stdout.log'),
|
||||
stderr_log=(tmpdir / 'ElmerSolver_stderr.log'))
|
||||
|
||||
capacitance_matrix = np.loadtxt(tmpdir / 'capacitance.txt')
|
||||
np.savetxt(output, capacitance_matrix)
|
||||
|
||||
|
||||
@cli.command()
|
||||
@click.option('-d', '--sim-dir', type=click.Path(dir_okay=True, file_okay=False, path_type=Path))
|
||||
@click.option('--solver-method')
|
||||
@click.argument('mesh_file', type=click.Path(dir_okay=False, path_type=Path))
|
||||
def inductance(mesh_file, sim_dir, solver_method):
|
||||
physical = dict(enumerate_mesh_bodies(mesh_file))
|
||||
|
||||
if sim_dir is not None:
|
||||
sim_dir = Path(sim_dir)
|
||||
sim_dir.mkdir(exist_ok=True)
|
||||
|
||||
sim = elmer.load_simulation('3D_steady', 'coil_mag_sim.yml')
|
||||
mesh_dir = '.'
|
||||
mesh_fn = 'mesh'
|
||||
sim.header['Mesh DB'] = f'"{mesh_dir}" "{mesh_fn}"'
|
||||
sim.constants.update({
|
||||
'Permittivity of Vacuum': str(constants.epsilon_0),
|
||||
'Gravity(4)': f'0 -1 0 {constants.g}',
|
||||
'Boltzmann Constant': str(constants.Boltzmann),
|
||||
'Unit Charge': str(constants.elementary_charge)})
|
||||
|
||||
air = elmer.load_material('air', sim, 'coil_mag_materials.yml')
|
||||
fr4 = elmer.load_material('fr4', sim, 'coil_mag_materials.yml')
|
||||
copper = elmer.load_material('copper', sim, 'coil_mag_materials.yml')
|
||||
|
||||
solver_current = elmer.load_solver('Static_Current_Conduction', sim, 'coil_mag_solvers.yml')
|
||||
solver_magdyn = elmer.load_solver('Magneto_Dynamics', sim, 'coil_mag_solvers.yml')
|
||||
if solver_method:
|
||||
solver_magdyn.data['Linear System Iterative Method'] = solver_method
|
||||
solver_magdyn_calc = elmer.load_solver('Magneto_Dynamics_Calculations', sim, 'coil_mag_solvers.yml')
|
||||
|
||||
copper_eqn = elmer.Equation(sim, 'copperEqn', [solver_current, solver_magdyn, solver_magdyn_calc])
|
||||
air_eqn = elmer.Equation(sim, 'airEqn', [solver_magdyn, solver_magdyn_calc])
|
||||
|
||||
bdy_trace = elmer.Body(sim, 'trace', [physical['trace'][1]])
|
||||
bdy_trace.material = copper
|
||||
bdy_trace.equation = copper_eqn
|
||||
|
||||
bdy_sub = elmer.Body(sim, 'substrate', [physical['substrate'][1]])
|
||||
bdy_sub.material = fr4
|
||||
bdy_sub.equation = air_eqn
|
||||
|
||||
bdy_ab = elmer.Body(sim, 'airbox', [physical['airbox'][1]])
|
||||
bdy_ab.material = air
|
||||
bdy_ab.equation = air_eqn
|
||||
|
||||
bdy_if_top = elmer.Body(sim, 'interface_top', [physical['interface_top'][1]])
|
||||
bdy_if_top.material = copper
|
||||
bdy_if_top.equation = copper_eqn
|
||||
|
||||
bdy_if_bottom = elmer.Body(sim, 'interface_bottom', [physical['interface_bottom'][1]])
|
||||
bdy_if_bottom.material = copper
|
||||
bdy_if_bottom.equation = copper_eqn
|
||||
|
||||
potential_force = elmer.BodyForce(sim, 'electric_potential', {'Electric Potential': 'Equals "Potential"'})
|
||||
bdy_trace.body_force = potential_force
|
||||
|
||||
# boundaries
|
||||
boundary_airbox = elmer.Boundary(sim, 'FarField', [physical['airbox_surface'][1]])
|
||||
boundary_airbox.data['Electric Infinity BC'] = 'True'
|
||||
|
||||
boundary_vplus = elmer.Boundary(sim, 'Vplus', [physical['interface_top'][1]])
|
||||
boundary_vplus.data['Potential'] = 1.0
|
||||
boundary_vplus.data['Save Scalars'] = True
|
||||
|
||||
boundary_vminus = elmer.Boundary(sim, 'Vminus', [physical['interface_bottom'][1]])
|
||||
boundary_vminus.data['Potential'] = 0.0
|
||||
|
||||
with tempfile.TemporaryDirectory() as tmpdir:
|
||||
tmpdir = sim_dir if sim_dir else Path(tmpdir)
|
||||
|
||||
sim.write_startinfo(tmpdir)
|
||||
sim.write_sif(tmpdir)
|
||||
# Convert mesh from gmsh to elemer formats. Also scale it from 1 unit = 1 mm to 1 unit = 1 m (SI units)
|
||||
elmer_grid(mesh_file.absolute(), 'mesh', cwd=tmpdir, scale=[1e-3, 1e-3, 1e-3],
|
||||
stdout_log=(tmpdir / 'ElmerGrid_stdout.log'),
|
||||
stderr_log=(tmpdir / 'ElmerGrid_stderr.log'))
|
||||
solver_stdout, solver_stderr = (tmpdir / 'ElmerSolver_stdout.log'), (tmpdir / 'ElmerSolver_stderr.log')
|
||||
res = elmer_solver(tmpdir,
|
||||
stdout_log=solver_stdout,
|
||||
stderr_log=solver_stderr)
|
||||
|
||||
P, R, U_mag = None, None, None
|
||||
solver_error = False
|
||||
for l in res.stdout.splitlines():
|
||||
if (m := re.fullmatch(r'StatCurrentSolve:\s*Total Heating Power\s*:\s*([0-9.+-Ee]+)\s*', l)):
|
||||
P = float(m.group(1))
|
||||
elif (m := re.fullmatch(r'StatCurrentSolve:\s*Effective Resistance\s*:\s*([0-9.+-Ee]+)\s*', l)):
|
||||
R = float(m.group(1))
|
||||
elif (m := re.fullmatch(r'MagnetoDynamicsCalcFields:\s*ElectroMagnetic Field Energy\s*:\s*([0-9.+-Ee]+)\s*', l)):
|
||||
U_mag = float(m.group(1))
|
||||
elif re.fullmatch(r'IterSolve: Linear iteration did not converge to tolerance', l):
|
||||
solver_error = True
|
||||
|
||||
if solver_error:
|
||||
raise click.ClickException(f'Error: One of the solvers did not converge. See log files for details:\n{solver_stdout.absolute()}\n{solver_stderr.absolute()}')
|
||||
elif P is None or R is None or U_mag is None:
|
||||
raise click.ClickException(f'Error during solver execution. Electrical parameters could not be calculated. See log files for details:\n{solver_stdout.absolute()}\n{solver_stderr.absolute()}')
|
||||
|
||||
V = math.sqrt(P*R)
|
||||
I = math.sqrt(P/R)
|
||||
L = 2*U_mag / (I**2)
|
||||
|
||||
assert math.isclose(V, 1.0, abs_tol=1e-3)
|
||||
|
||||
print(f'Total magnetic field energy: {format_si(U_mag, "J")}')
|
||||
print(f'Reference coil current: {format_si(I, "Ω")}')
|
||||
print(f'Coil resistance calculated by solver: {format_si(R, "Ω")}')
|
||||
print(f'Inductance calucated from field: {format_si(L, "H")}')
|
||||
|
||||
|
||||
@cli.command()
|
||||
@click.option('-r', '--reference-field', type=float, required=True)
|
||||
@click.option('-d', '--sim-dir', type=click.Path(dir_okay=True, file_okay=False, path_type=Path))
|
||||
@click.argument('mesh_file', type=click.Path(dir_okay=False, path_type=Path))
|
||||
def mutual_inductance(mesh_file, sim_dir, reference_field):
|
||||
physical = dict(enumerate_mesh_bodies(mesh_file))
|
||||
|
||||
if sim_dir is not None:
|
||||
sim_dir = Path(sim_dir)
|
||||
sim_dir.mkdir(exist_ok=True)
|
||||
|
||||
sim = elmer.load_simulation('3D_steady', 'coil_mag_sim.yml')
|
||||
mesh_dir = '.'
|
||||
mesh_fn = 'mesh'
|
||||
sim.header['Mesh DB'] = f'"{mesh_dir}" "{mesh_fn}"'
|
||||
sim.constants.update({
|
||||
'Permittivity of Vacuum': str(constants.epsilon_0),
|
||||
'Gravity(4)': f'0 -1 0 {constants.g}',
|
||||
'Boltzmann Constant': str(constants.Boltzmann),
|
||||
'Unit Charge': str(constants.elementary_charge)})
|
||||
|
||||
air = elmer.load_material('air', sim, 'coil_mag_materials.yml')
|
||||
fr4 = elmer.load_material('fr4', sim, 'coil_mag_materials.yml')
|
||||
copper = elmer.load_material('copper', sim, 'coil_mag_materials.yml')
|
||||
|
||||
solver_current = elmer.load_solver('Static_Current_Conduction', sim, 'coil_mag_solvers.yml')
|
||||
solver_magdyn = elmer.load_solver('Magneto_Dynamics', sim, 'coil_mag_solvers.yml')
|
||||
solver_magdyn_calc = elmer.load_solver('Magneto_Dynamics_Calculations', sim, 'coil_mag_solvers.yml')
|
||||
|
||||
copper_eqn = elmer.Equation(sim, 'copperEqn', [solver_current, solver_magdyn, solver_magdyn_calc])
|
||||
air_eqn = elmer.Equation(sim, 'airEqn', [solver_magdyn, solver_magdyn_calc])
|
||||
|
||||
bdy_trace1 = elmer.Body(sim, 'trace1', [physical['trace1'][1]])
|
||||
bdy_trace1.material = copper
|
||||
bdy_trace1.equation = copper_eqn
|
||||
|
||||
bdy_trace2 = elmer.Body(sim, 'trace2', [physical['trace2'][1]])
|
||||
bdy_trace2.material = copper
|
||||
bdy_trace2.equation = copper_eqn
|
||||
|
||||
bdy_sub1 = elmer.Body(sim, 'substrate1', [physical['substrate1'][1]])
|
||||
bdy_sub1.material = fr4
|
||||
bdy_sub1.equation = air_eqn
|
||||
|
||||
bdy_sub2 = elmer.Body(sim, 'substrate2', [physical['substrate2'][1]])
|
||||
bdy_sub2.material = fr4
|
||||
bdy_sub2.equation = air_eqn
|
||||
|
||||
|
||||
bdy_ab = elmer.Body(sim, 'airbox', [physical['airbox'][1]])
|
||||
bdy_ab.material = air
|
||||
bdy_ab.equation = air_eqn
|
||||
|
||||
bdy_if_top1 = elmer.Body(sim, 'interface_top1', [physical['interface_top1'][1]])
|
||||
bdy_if_top1.material = copper
|
||||
bdy_if_top1.equation = copper_eqn
|
||||
|
||||
bdy_if_bottom1 = elmer.Body(sim, 'interface_bottom1', [physical['interface_bottom1'][1]])
|
||||
bdy_if_bottom1.material = copper
|
||||
bdy_if_bottom1.equation = copper_eqn
|
||||
|
||||
bdy_if_top2 = elmer.Body(sim, 'interface_top2', [physical['interface_top2'][1]])
|
||||
bdy_if_top2.material = copper
|
||||
bdy_if_top2.equation = copper_eqn
|
||||
|
||||
bdy_if_bottom2 = elmer.Body(sim, 'interface_bottom2', [physical['interface_bottom2'][1]])
|
||||
bdy_if_bottom2.material = copper
|
||||
bdy_if_bottom2.equation = copper_eqn
|
||||
|
||||
potential_force = elmer.BodyForce(sim, 'electric_potential', {'Electric Potential': 'Equals "Potential"'})
|
||||
bdy_trace1.body_force = potential_force
|
||||
bdy_trace2.body_force = potential_force
|
||||
|
||||
# boundaries
|
||||
boundary_airbox = elmer.Boundary(sim, 'FarField', [physical['airbox_surface'][1]])
|
||||
boundary_airbox.data['Electric Infinity BC'] = 'True'
|
||||
|
||||
boundary_vplus1 = elmer.Boundary(sim, 'Vplus1', [physical['interface_top1'][1]])
|
||||
boundary_vplus1.data['Potential'] = 1.0
|
||||
boundary_vplus1.data['Save Scalars'] = True
|
||||
|
||||
boundary_vminus1 = elmer.Boundary(sim, 'Vminus1', [physical['interface_bottom1'][1]])
|
||||
boundary_vminus1.data['Potential'] = 0.0
|
||||
|
||||
boundary_vplus2 = elmer.Boundary(sim, 'Vplus2', [physical['interface_top2'][1]])
|
||||
boundary_vplus2.data['Potential'] = 1.0
|
||||
boundary_vplus2.data['Save Scalars'] = True
|
||||
|
||||
boundary_vminus2 = elmer.Boundary(sim, 'Vminus2', [physical['interface_bottom2'][1]])
|
||||
boundary_vminus2.data['Potential'] = 0.0
|
||||
|
||||
with tempfile.TemporaryDirectory() as tmpdir:
|
||||
tmpdir = sim_dir if sim_dir else Path(tmpdir)
|
||||
|
||||
sim.write_startinfo(tmpdir)
|
||||
sim.write_sif(tmpdir)
|
||||
# Convert mesh from gmsh to elemer formats. Also scale it from 1 unit = 1 mm to 1 unit = 1 m (SI units)
|
||||
elmer_grid(mesh_file.absolute(), 'mesh', cwd=tmpdir, scale=[1e-3, 1e-3, 1e-3],
|
||||
stdout_log=(tmpdir / 'ElmerGrid_stdout.log'),
|
||||
stderr_log=(tmpdir / 'ElmerGrid_stderr.log'))
|
||||
solver_stdout, solver_stderr = (tmpdir / 'ElmerSolver_stdout.log'), (tmpdir / 'ElmerSolver_stderr.log')
|
||||
res = elmer_solver(tmpdir,
|
||||
stdout_log=solver_stdout,
|
||||
stderr_log=solver_stderr)
|
||||
|
||||
P, R, U_mag = None, None, None
|
||||
solver_error = False
|
||||
for l in res.stdout.splitlines():
|
||||
if (m := re.fullmatch(r'StatCurrentSolve:\s*Total Heating Power\s*:\s*([0-9.+-Ee]+)\s*', l)):
|
||||
P = float(m.group(1))
|
||||
elif (m := re.fullmatch(r'StatCurrentSolve:\s*Effective Resistance\s*:\s*([0-9.+-Ee]+)\s*', l)):
|
||||
R = float(m.group(1))
|
||||
elif (m := re.fullmatch(r'MagnetoDynamicsCalcFields:\s*ElectroMagnetic Field Energy\s*:\s*([0-9.+-Ee]+)\s*', l)):
|
||||
U_mag = float(m.group(1))
|
||||
elif re.fullmatch(r'IterSolve: Linear iteration did not converge to tolerance', l):
|
||||
solver_error = True
|
||||
|
||||
if solver_error:
|
||||
raise click.ClickException(f'Error: One of the solvers did not converge. See log files for details:\n{solver_stdout.absolute()}\n{solver_stderr.absolute()}')
|
||||
elif P is None or R is None or U_mag is None:
|
||||
raise click.ClickException(f'Error during solver execution. Electrical parameters could not be calculated. See log files for details:\n{solver_stdout.absolute()}\n{solver_stderr.absolute()}')
|
||||
|
||||
V = math.sqrt(P*R)
|
||||
I = math.sqrt(P/R)
|
||||
Lm = (U_mag - 2*reference_field) / ((I/2)**2)
|
||||
|
||||
assert math.isclose(V, 1.0, abs_tol=1e-3)
|
||||
|
||||
print(f'Mutual inductance calucated from field: {format_si(Lm, "H")}')
|
||||
|
||||
|
||||
@cli.command()
|
||||
@click.option('-d', '--sim-dir', type=click.Path(dir_okay=True, file_okay=False, path_type=Path))
|
||||
@click.argument('mesh_file', type=click.Path(dir_okay=False, path_type=Path))
|
||||
def self_capacitance(mesh_file, sim_dir):
|
||||
physical = dict(enumerate_mesh_bodies(mesh_file))
|
||||
|
||||
if sim_dir is not None:
|
||||
sim_dir = Path(sim_dir)
|
||||
sim_dir.mkdir(exist_ok=True)
|
||||
|
||||
sim = elmer.load_simulation('3D_steady', 'self_capacitance_sim.yml')
|
||||
mesh_dir = '.'
|
||||
mesh_fn = 'mesh'
|
||||
sim.header['Mesh DB'] = f'"{mesh_dir}" "{mesh_fn}"'
|
||||
sim.constants.update({
|
||||
'Permittivity of Vacuum': str(constants.epsilon_0),
|
||||
'Gravity(4)': f'0 -1 0 {constants.g}',
|
||||
'Boltzmann Constant': str(constants.Boltzmann),
|
||||
'Unit Charge': str(constants.elementary_charge)})
|
||||
|
||||
air = elmer.load_material('air', sim, 'coil_mag_materials.yml')
|
||||
fr4 = elmer.load_material('fr4', sim, 'coil_mag_materials.yml')
|
||||
copper = elmer.load_material('copper', sim, 'coil_mag_materials.yml')
|
||||
|
||||
solver_current = elmer.load_solver('StaticCurrent', sim, 'self_capacitance_solvers.yml')
|
||||
solver_estat = elmer.load_solver('Electrostatics', sim, 'self_capacitance_solvers.yml')
|
||||
|
||||
copper_eqn = elmer.Equation(sim, 'copperEqn', [solver_current, solver_estat])
|
||||
air_eqn = elmer.Equation(sim, 'airEqn', [solver_estat])
|
||||
|
||||
bdy_trace = elmer.Body(sim, 'trace', [physical['trace'][1]])
|
||||
bdy_trace.material = copper
|
||||
bdy_trace.equation = copper_eqn
|
||||
|
||||
bdy_sub = elmer.Body(sim, 'substrate', [physical['substrate'][1]])
|
||||
bdy_sub.material = fr4
|
||||
bdy_sub.equation = air_eqn
|
||||
|
||||
bdy_ab = elmer.Body(sim, 'airbox', [physical['airbox'][1]])
|
||||
bdy_ab.material = air
|
||||
bdy_ab.equation = air_eqn
|
||||
|
||||
bdy_if_top = elmer.Body(sim, 'interface_top', [physical['interface_top'][1]])
|
||||
bdy_if_top.material = copper
|
||||
bdy_if_top.equation = copper_eqn
|
||||
|
||||
bdy_if_bottom = elmer.Body(sim, 'interface_bottom', [physical['interface_bottom'][1]])
|
||||
bdy_if_bottom.material = copper
|
||||
bdy_if_bottom.equation = copper_eqn
|
||||
|
||||
potential_force = elmer.BodyForce(sim, 'electric_potential', {'Potential': 'Equals "PotentialStat"'})
|
||||
bdy_trace.body_force = potential_force
|
||||
|
||||
# boundaries
|
||||
boundary_airbox = elmer.Boundary(sim, 'FarField', [physical['airbox_surface'][1]])
|
||||
boundary_airbox.data['Electric Infinity BC'] = 'True'
|
||||
|
||||
boundary_vplus = elmer.Boundary(sim, 'Vplus', [physical['interface_top'][1]])
|
||||
boundary_vplus.data['PotentialStat'] = 'Real 1.0'
|
||||
boundary_vplus.data['Save Scalars'] = True
|
||||
|
||||
boundary_vminus = elmer.Boundary(sim, 'Vminus', [physical['interface_bottom'][1]])
|
||||
boundary_vminus.data['PotentialStat'] = 'Real 0.0'
|
||||
|
||||
with tempfile.TemporaryDirectory() as tmpdir:
|
||||
tmpdir = sim_dir if sim_dir else Path(tmpdir)
|
||||
|
||||
sim.write_startinfo(tmpdir)
|
||||
sim.write_sif(tmpdir)
|
||||
# Convert mesh from gmsh to elemer formats. Also scale it from 1 unit = 1 mm to 1 unit = 1 m (SI units)
|
||||
elmer_grid(mesh_file.absolute(), 'mesh', cwd=tmpdir, scale=[1e-3, 1e-3, 1e-3],
|
||||
stdout_log=(tmpdir / 'ElmerGrid_stdout.log'),
|
||||
stderr_log=(tmpdir / 'ElmerGrid_stderr.log'))
|
||||
solver_stdout, solver_stderr = (tmpdir / 'ElmerSolver_stdout.log'), (tmpdir / 'ElmerSolver_stderr.log')
|
||||
res = elmer_solver(tmpdir,
|
||||
stdout_log=solver_stdout,
|
||||
stderr_log=solver_stderr)
|
||||
|
||||
C, U_elec = None, None
|
||||
solver_error = False
|
||||
for l in res.stdout.splitlines():
|
||||
if (m := re.fullmatch(r'StatElecSolve:\s*Tot. Electric Energy\s*:\s*([0-9.+-Ee]+)\s*', l)):
|
||||
U_elec = float(m.group(1))
|
||||
elif (m := re.fullmatch(r'StatElecSolve:\s*Capacitance\s*:\s*([0-9.+-Ee]+)\s*', l)):
|
||||
C = float(m.group(1))
|
||||
elif re.fullmatch(r'IterSolve: Linear iteration did not converge to tolerance', l):
|
||||
solver_error = True
|
||||
|
||||
if solver_error:
|
||||
raise click.ClickException(f'Error: One of the solvers did not converge. See log files for details:\n{solver_stdout.absolute()}\n{solver_stderr.absolute()}')
|
||||
elif C is None or U_elec is None:
|
||||
raise click.ClickException(f'Error during solver execution. Electrical parameters could not be calculated. See log files for details:\n{solver_stdout.absolute()}\n{solver_stderr.absolute()}')
|
||||
|
||||
print(f'Total electric field energy: {format_si(U_elec, "J")}')
|
||||
print(f'Total parasitic capacitance: {format_si(C, "F")}')
|
||||
|
||||
@cli.command()
|
||||
@click.option('-d', '--sim-dir', type=click.Path(dir_okay=True, file_okay=False, path_type=Path))
|
||||
@click.option('--capacitance-matrix-file', type=click.Path(dir_okay=False, exists=True))
|
||||
@click.option('--total-inductance', type=float, required=True, help='Total inductance in Henry')
|
||||
@click.option('--total-resistance', type=float, required=True, help='Total resistance in Ohm')
|
||||
@click.option('--plot-out', type=click.Path(dir_okay=False, writable=True), help='Optional SVG plot output file')
|
||||
def resonance(sim_dir, capacitance_matrix_file, total_inductance, total_resistance, plot_out):
|
||||
import PySpice.Unit
|
||||
from PySpice.Spice.Library import SpiceLibrary
|
||||
from PySpice.Spice.Netlist import Circuit
|
||||
from PySpice.Plot.BodeDiagram import bode_diagram
|
||||
import scipy.signal
|
||||
from matplotlib import pyplot as plt
|
||||
|
||||
capacitance_matrix = np.loadtxt(capacitance_matrix_file)
|
||||
num_elements = capacitance_matrix.shape[0]
|
||||
|
||||
circ = Circuit('LC ladder parasitic sim')
|
||||
inputs = 'Vplus', circ.gnd
|
||||
coil_in = 'coil_in'
|
||||
|
||||
Rtest = circ.R('Rtest', inputs[0], coil_in, 50@PySpice.Unit.u_Ohm)
|
||||
|
||||
intermediate_nodes = [f'intermediate{i}' for i in range(num_elements-1)]
|
||||
inductor_nodes = [(a, b) for a, b in zip([coil_in, *intermediate_nodes], [*intermediate_nodes, inputs[1]])]
|
||||
inductor_midpoints = [f'midpoint{i}' for i in range(num_elements)]
|
||||
|
||||
circ.SinusoidalVoltageSource('input', inputs[0], inputs[1], amplitude=1@PySpice.Unit.u_V)
|
||||
|
||||
for i, ((a, b), m) in enumerate(zip(inductor_nodes, inductor_midpoints)):
|
||||
L = total_inductance / num_elements / 2
|
||||
R = total_resistance / num_elements / 2
|
||||
circ.L(f'L{i}A', a, f'R{i}A1', L@PySpice.Unit.u_H)
|
||||
circ.R(f'R{i}A', f'R{i}A1', m, R@PySpice.Unit.u_Ohm)
|
||||
circ.R(f'R{i}B', m, f'R{i}B1', R@PySpice.Unit.u_Ohm)
|
||||
circ.L(f'L{i}B', f'R{i}B1', b, L@PySpice.Unit.u_H)
|
||||
|
||||
for i in range(num_elements):
|
||||
for j in range(i):
|
||||
circ.C(f'C{i}_{j}', inductor_midpoints[i], inductor_midpoints[j], capacitance_matrix[i, j]@PySpice.Unit.u_F)
|
||||
|
||||
sim = circ.simulator(temperature=25, nominal_temperature=25)
|
||||
ana = sim.ac(start_frequency=10@PySpice.Unit.u_kHz, stop_frequency=1000@PySpice.Unit.u_MHz, number_of_points=1000, variation='dec')
|
||||
figure, axs = plt.subplots(2, figsize=(20, 10), sharex=True)
|
||||
|
||||
freq = ana.frequency
|
||||
gain = 20*np.log10(np.absolute(ana.coil_in))
|
||||
|
||||
peaks, peak_props = scipy.signal.find_peaks(-gain, height=20)
|
||||
for peak in peaks[:3]:
|
||||
print(f'Resonance at {float(freq[peak])/1e6:.3f} MHz')
|
||||
|
||||
if plot_out:
|
||||
plt.title("Bode Diagram of a Low-Pass RC Filter")
|
||||
bode_diagram(axes=axs,
|
||||
frequency=freq,
|
||||
gain=gain,
|
||||
phase=np.angle(ana.coil_in, deg=False),
|
||||
linestyle='-',
|
||||
)
|
||||
|
||||
for peak in peaks[:3]:
|
||||
for ax in axs:
|
||||
ax.axvline(float(freq[peak]), color='red', alpha=0.5)
|
||||
|
||||
plt.tight_layout()
|
||||
plt.savefig(plot_out)
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
cli()
|
||||
|
||||
|
|
@ -1,38 +0,0 @@
|
|||
air:
|
||||
Density: 1.1885 # 20°C
|
||||
Electric Conductivity: 0.0
|
||||
Heat Capacity: 1006.4 # 20°C
|
||||
Heat Conductivity: 0.025873 # 20°C
|
||||
Relative Permeability: 1
|
||||
Relative Permittivity: 1
|
||||
ro4003c:
|
||||
Density: 1790 # 23°C
|
||||
Relative Permeability: 1
|
||||
Relative Permittivity: 3.55
|
||||
fr4:
|
||||
Density: 1850 # 23°C
|
||||
Relative Permeability: 1
|
||||
Relative Permittivity: 4.4
|
||||
Heat Conductivity: 0.81 # in-plane
|
||||
ideal:
|
||||
Relative Permittivity: 1
|
||||
copper:
|
||||
Density: 8960.0 # 0°C
|
||||
Electric Conductivity: 59600000 # 20°C
|
||||
Emissivity: 0.012 # 327°C
|
||||
Heat Capacity: 415.0 # 200°C
|
||||
Heat Conductivity: 401.0 # 0°C
|
||||
Relative Permeability: 1
|
||||
Relative Permittivity: 1
|
||||
steel_1.4541:
|
||||
Density: 7900.0 # 20°C
|
||||
Electric Conductivity: 1370
|
||||
Emissivity: 0.111 # 200°C
|
||||
Heat Capacity: 470.0 # 20°C
|
||||
Heat Conductivity: 15.0 # 20°C
|
||||
Relative Permeability: 1
|
||||
Relative Permittivity: 1
|
||||
water:
|
||||
Density: 1000.0
|
||||
Heat Capacity: 4182.0
|
||||
Heat Conductivity: 0.6
|
||||
|
|
@ -1,50 +0,0 @@
|
|||
2D_steady:
|
||||
Max Output Level: 4
|
||||
Coordinate System: Cartesian 2D
|
||||
Simulation Type: Steady state
|
||||
Steady State Max Iterations: 10
|
||||
|
||||
3D_steady:
|
||||
Max Output Level: 5
|
||||
Coordinate System: Cartesian
|
||||
Coordinate Mapping(3): 1 2 3
|
||||
Simulation Type: Steady state
|
||||
Steady State Max Iterations: 1
|
||||
Output Intervals: 1
|
||||
Timestepping Method: BDF
|
||||
BDF Order: 1
|
||||
Solver Input File: case.sif
|
||||
Post File: case.vtu
|
||||
Output File: case.result
|
||||
|
||||
2D_transient:
|
||||
Max Output Level: 4
|
||||
Coordinate System: Cartesian 2D
|
||||
Simulation Type: Transient
|
||||
Steady State Max Iterations: 10
|
||||
Output File: case.result
|
||||
Output Intervals: 10
|
||||
Timestep Sizes: 0.1
|
||||
Timestep Intervals: 100
|
||||
Timestepping Method: BDF
|
||||
BDF Order: 1
|
||||
|
||||
|
||||
axi-symmetric_steady:
|
||||
Max Output Level: 4
|
||||
Coordinate System: Axi Symmetric
|
||||
Simulation Type: Steady state
|
||||
Steady State Max Iterations: 10
|
||||
Output File: case.result
|
||||
|
||||
axi-symmetric_transient:
|
||||
Max Output Level: 4
|
||||
Coordinate System: Axi Symmetric
|
||||
Simulation Type: Transient
|
||||
Steady State Max Iterations: 10
|
||||
Output File: case.result
|
||||
Output Intervals: 10
|
||||
Timestep Sizes: 0.1
|
||||
Timestep Intervals: 100
|
||||
Timestepping Method: BDF
|
||||
BDF Order: 1
|
||||
|
|
@ -1,203 +0,0 @@
|
|||
Electrostatics_Capacitance:
|
||||
Equation: Electrostatics
|
||||
Calculate Electric Field: True
|
||||
Calculate Capacitance Matrix: True
|
||||
Capacitance Matrix Filename: capacitance.txt
|
||||
Procedure: '"StatElecSolve" "StatElecSolver"'
|
||||
Variable: Potential
|
||||
Calculate Electric Energy: True
|
||||
Exec Solver: Always
|
||||
Stabilize: True
|
||||
Bubbles: False
|
||||
Lumped Mass Matrix: False
|
||||
Optimize Bandwidth: True
|
||||
Steady State Convergence Tolerance: 1.0e-5
|
||||
Nonlinear System Convergence Tolerance: 1.0e-7
|
||||
Nonlinear System Max Iterations: 20
|
||||
Nonlinear System Newton After Iterations: 3
|
||||
Nonlinear System Newton After Tolerance: 1.0e-3
|
||||
Nonlinear System Relaxation Factor: 1
|
||||
Linear System Solver: Iterative
|
||||
Linear System Iterative Method: BiCGStab
|
||||
Linear System Max Iterations: 500
|
||||
Linear System Convergence Tolerance: 1.0e-10
|
||||
BiCGstabl polynomial degree: 2
|
||||
Linear System Preconditioning: ILU0
|
||||
Linear System ILUT Tolerance: 1.0e-3
|
||||
Linear System Abort Not Converged: False
|
||||
Linear System Residual Output: 10
|
||||
Linear System Precondition Recompute: 1
|
||||
Electrostatics:
|
||||
Equation: Electrostatics
|
||||
Calculate Electric Field: True
|
||||
Procedure: '"StatElecSolve" "StatElecSolver"'
|
||||
Variable: Potential
|
||||
Calculate Electric Energy: True
|
||||
Exec Solver: Always
|
||||
Stabilize: True
|
||||
Bubbles: False
|
||||
Lumped Mass Matrix: False
|
||||
Optimize Bandwidth: True
|
||||
Steady State Convergence Tolerance: 1.0e-5
|
||||
Nonlinear System Convergence Tolerance: 1.0e-7
|
||||
Nonlinear System Max Iterations: 20
|
||||
Nonlinear System Newton After Iterations: 3
|
||||
Nonlinear System Newton After Tolerance: 1.0e-3
|
||||
Nonlinear System Relaxation Factor: 1
|
||||
Linear System Solver: Iterative
|
||||
Linear System Iterative Method: BiCGStab
|
||||
Linear System Max Iterations: 500
|
||||
Linear System Convergence Tolerance: 1.0e-10
|
||||
BiCGstabl polynomial degree: 2
|
||||
Linear System Preconditioning: ILU0
|
||||
Linear System ILUT Tolerance: 1.0e-3
|
||||
Linear System Abort Not Converged: False
|
||||
Linear System Residual Output: 10
|
||||
Linear System Precondition Recompute: 1
|
||||
ThermoElectricSolver:
|
||||
Equation: ThermoElectric
|
||||
Procedure: '"ThermoElectricSolver" "ThermoElectricSolver"'
|
||||
Variable: POT[Temperature:1 Potential:1]
|
||||
Element: '"p:1"'
|
||||
Calculate Loads: True
|
||||
Exec Solver: Always
|
||||
Nonlinear System Convergence Tolerance: 1.0e-6
|
||||
Nonlinear System Max Iterations: 100
|
||||
Nonlinear System Newton After Iterations : 1
|
||||
Nonlinear System Newton After Tolerance: 1e-9
|
||||
Linear System Solver: '"Iterative"'
|
||||
Linear System Iterative Method: BicgstabL
|
||||
Bicgstabl Polynomial Degree: 2
|
||||
Linear System Max Iterations: 200
|
||||
Linear System Residual Output: 40
|
||||
Linear System Preconditioning: Ilu
|
||||
Linear System Convergence Tolerance: 1e-8
|
||||
Steady State Convergence Tolerance: 1e-6
|
||||
HeatSolver:
|
||||
Equation: HeatSolver
|
||||
Procedure: '"HeatSolve" "HeatSolver"'
|
||||
Variable: '"Temperature"'
|
||||
Variable Dofs: 1
|
||||
Calculate Loads: True
|
||||
Exec Solver: Always
|
||||
Nonlinear System Convergence Tolerance: 1.0e-6
|
||||
Nonlinear System Max Iterations: 1000
|
||||
Nonlinear System Relaxation Factor: 0.7
|
||||
Steady State Convergence Tolerance: 1.0e-6
|
||||
Stabilize: True # Necessary in convection-dominated systems
|
||||
Optimize Bandwidth: True
|
||||
Linear System Solver: Iterative
|
||||
Linear System Iterative Method: BiCGStab
|
||||
Linear System Max Iterations: 1000
|
||||
Linear System Preconditioning: ILU
|
||||
Linear System Precondition Recompute: 1
|
||||
Linear System Convergence Tolerance: 1.0e-8
|
||||
Linear System Abort Not Converged: True
|
||||
Linear System Residual Output: 1
|
||||
Smart Heater Control After Tolerance: 1.0e-4
|
||||
MagnetoDynamics2DHarmonic:
|
||||
Equation: MagnetoDynamics2DHarmonic
|
||||
Procedure: '"MagnetoDynamics2D" "MagnetoDynamics2DHarmonic"'
|
||||
Variable: 'Potential[Potential Re:1 Potential Im:1]'
|
||||
Variable Dofs: 2
|
||||
Exec Solver: Always
|
||||
Nonlinear System Convergence Tolerance: 1.0e-6
|
||||
Nonlinear System Max Iterations: 1000
|
||||
Nonlinear System Relaxation Factor: 0.7
|
||||
Steady State Convergence Tolerance: 1.0e-6
|
||||
Optimize Bandwidth: True
|
||||
Linear System Solver: Iterative
|
||||
Linear System Iterative Method: BiCGStab
|
||||
Linear System Max Iterations: 1000
|
||||
Linear System Preconditioning: ILU
|
||||
Linear System Precondition Recompute: 1
|
||||
Linear System Convergence Tolerance: 1.0e-8
|
||||
Linear System Abort Not Converged: True
|
||||
Linear System Residual Output: 1
|
||||
MagnetoDynamicsCalcFields:
|
||||
Equation: MagnetoDynamicsCalcFields
|
||||
Procedure: '"MagnetoDynamics" "MagnetoDynamicsCalcFields"'
|
||||
Potential Variable: Potential
|
||||
Angular Frequency: 8.48e4
|
||||
Calculate Joule Heating: True
|
||||
Calculate Magnetic Field Strength: True
|
||||
Calculate Electric Field: True
|
||||
Exec Solver: Always
|
||||
Calculate Nodal Fields: Logical False
|
||||
Calculate Elemental Fields: Logical True
|
||||
StatMagSolver:
|
||||
Equation: StatMagSolver
|
||||
Procedure: '"StatMagSolve" "StatMagSolver"'
|
||||
Variable: Potential
|
||||
Variable DOFs: 2
|
||||
Calculate Joule Heating: 'Logical True'
|
||||
Calculate Magnetic Flux: 'Logical True'
|
||||
Exec Solver: Always
|
||||
Nonlinear System Convergence Tolerance: 1.0e-6
|
||||
Nonlinear System Max Iterations: 1000
|
||||
Nonlinear System Relaxation Factor: 0.7
|
||||
Steady State Convergence Tolerance: 1.0e-6
|
||||
Optimize Bandwidth: True
|
||||
Linear System Solver: Iterative
|
||||
Linear System Iterative Method: BiCGStab
|
||||
Linear System Max Iterations: 1000
|
||||
Linear System Preconditioning: ILU
|
||||
Linear System Precondition Recompute: 1
|
||||
Linear System Convergence Tolerance: 1.0e-8
|
||||
Linear System Abort Not Converged: True
|
||||
Linear System Residual Output: 1
|
||||
SaveMaterials:
|
||||
Exec Solver: 'After timestep'
|
||||
Procedure: 'File "SaveData" "SaveMaterials"'
|
||||
Parameter 1: 'String "Heat Conductivity"'
|
||||
ResultOutputSolver:
|
||||
Exec Solver: 'After timestep'
|
||||
Equation: ResultOutputSolver
|
||||
Procedure: '"ResultOutputSolve" "ResultOutputSolver"'
|
||||
VTU Format: True
|
||||
Save Geometry Ids: 'Logical True'
|
||||
FluxSolver:
|
||||
Exec Solver: 'After timestep'
|
||||
Equation: 'Flux Solver'
|
||||
Procedure: '"FluxSolver" "FluxSolver"'
|
||||
Calculate Grad: 'Logical True'
|
||||
Calculate Flux: 'Logical True'
|
||||
Target Variable: 'String "Temperature"'
|
||||
Flux Coefficient: 'String "Heat Conductivity"'
|
||||
Exec Solver: Always
|
||||
Nonlinear System Convergence Tolerance: 1.0e-6
|
||||
Nonlinear System Max Iterations: 1000
|
||||
Nonlinear System Relaxation Factor: 0.7
|
||||
Steady State Convergence Tolerance: 1.0e-6
|
||||
Optimize Bandwidth: True
|
||||
Linear System Solver: Iterative
|
||||
Linear System Iterative Method: BiCGStab
|
||||
Linear System Max Iterations: 1000
|
||||
Linear System Preconditioning: ILU
|
||||
Linear System Precondition Recompute: 1
|
||||
Linear System Convergence Tolerance: 1.0e-8
|
||||
Linear System Abort Not Converged: True
|
||||
Linear System Residual Output: 1
|
||||
SaveScalars:
|
||||
Exec Solver: 'After timestep'
|
||||
Equation: SaveScalars
|
||||
Procedure: '"SaveData" "SaveScalars"'
|
||||
Filename: '"boundary_scalars.dat"'
|
||||
Output Directory: './results'
|
||||
Operator 1: 'boundary sum'
|
||||
Variable 1: 'Temperature Loads'
|
||||
Operator 2: 'diffusive flux'
|
||||
Variable 2: Temperature
|
||||
Coefficient 2: 'Heat Conductivity'
|
||||
SaveLine:
|
||||
Exec Solver: 'After timestep'
|
||||
Equation: '"SaveLine"'
|
||||
Procedure: '"SaveData" "SaveLine"'
|
||||
Filename: '"boundary_lines.dat"'
|
||||
Output Directory: './results'
|
||||
Variable 1: Temperature Loads
|
||||
SteadyPhaseChange:
|
||||
Equation: SteadyPhaseChange
|
||||
Variable: '"Phase Surface"'
|
||||
Procedure: '"SteadyPhaseChange" "SteadyPhaseChange"'
|
||||
Internal Mesh Movement: 'Logical True'
|
||||
|
|
@ -1,483 +0,0 @@
|
|||
#!/usr/bin/env python3
|
||||
|
||||
import math
|
||||
import hashlib
|
||||
import re
|
||||
import itertools
|
||||
import datetime
|
||||
import tempfile
|
||||
import subprocess
|
||||
import sqlite3
|
||||
import json
|
||||
from pathlib import Path
|
||||
|
||||
import tqdm
|
||||
|
||||
import gerbonara.cad.kicad.pcb as pcb
|
||||
import gerbonara.cad.kicad.footprints as fp
|
||||
import gerbonara.cad.primitives as cad_pr
|
||||
import gerbonara.cad.kicad.graphical_primitives as kc_gr
|
||||
|
||||
|
||||
cols = 5
|
||||
rows = 5
|
||||
|
||||
coil_specs = [
|
||||
{'n': 1, 's': True, 't': 1, 'c': 0.20, 'w': 5.00, 'd': 3.00, 'v': 5.00},
|
||||
{'n': 2, 's': True, 't': 1, 'c': 0.20, 'w': 3.00, 'd': 1.50, 'v': 3.00},
|
||||
{'n': 3, 's': True, 't': 1, 'c': 0.20, 'w': 1.50, 'd': 1.20, 'v': 2.00},
|
||||
{'n': 5, 's': True, 't': 1, 'c': 0.20, 'w': 0.80, 'd': 0.40, 'v': 0.80},
|
||||
{'n': 10, 's': True, 't': 1, 'c': 0.20, 'w': 0.50, 'd': 0.30, 'v': 0.60},
|
||||
{'n': 25, 's': True, 't': 1, 'c': 0.15, 'w': 0.25, 'd': 0.30, 'v': 0.60},
|
||||
|
||||
{'n': 1, 's': False, 't': 3, 'c': 0.20, 'w': 5.00, 'd': 3.00, 'v': 5.00},
|
||||
{'n': 2, 's': False, 't': 1, 'c': 0.20, 'w': 3.00, 'd': 1.50, 'v': 3.00},
|
||||
{'n': 3, 's': False, 't': 1, 'c': 0.20, 'w': 2.50, 'd': 1.20, 'v': 2.00},
|
||||
{'n': 5, 's': False, 't': 1, 'c': 0.20, 'w': 2.50, 'd': 1.20, 'v': 0.80},
|
||||
{'n': 10, 's': False, 't': 1, 'c': 0.20, 'w': 1.50, 'd': 0.80, 'v': 0.60},
|
||||
{'n': 25, 's': False, 't': 1, 'c': 0.15, 'w': 0.50, 'd': 0.30, 'v': 0.60},
|
||||
|
||||
{'n': 1, 's': False, 't': 4, 'c': 0.20, 'w': 5.00, 'd': 3.00, 'v': 5.00},
|
||||
{'n': 2, 's': False, 't': 3, 'c': 0.20, 'w': 3.00, 'd': 1.50, 'v': 3.00},
|
||||
{'n': 3, 's': False, 't': 4, 'c': 0.20, 'w': 2.50, 'd': 1.20, 'v': 2.00},
|
||||
{'n': 5, 's': False, 't': 3, 'c': 0.20, 'w': 2.50, 'd': 1.20, 'v': 0.80},
|
||||
{'n': 10, 's': False, 't': 3, 'c': 0.20, 'w': 1.50, 'd': 0.80, 'v': 0.60},
|
||||
{'n': 25, 's': False, 't': 3, 'c': 0.15, 'w': 0.50, 'd': 0.30, 'v': 0.60},
|
||||
|
||||
{'n': 1, 's': False, 't': 5, 'c': 0.20, 'w': 5.00, 'd': 3.00, 'v': 5.00},
|
||||
{'n': 2, 's': False, 't': 5, 'c': 0.20, 'w': 3.00, 'd': 1.50, 'v': 3.00},
|
||||
{'n': 3, 's': False, 't': 4, 'c': 0.20, 'w': 2.50, 'd': 1.20, 'v': 2.00},
|
||||
{'n': 5, 's': False, 't': 7, 'c': 0.20, 'w': 2.50, 'd': 1.20, 'v': 0.80},
|
||||
{'n': 10, 's': False, 't': 7, 'c': 0.20, 'w': 1.50, 'd': 0.80, 'v': 0.60},
|
||||
{'n': 25, 's': False, 't': 13, 'c': 0.15, 'w': 0.50, 'd': 0.30, 'v': 0.60},
|
||||
|
||||
{'n': 25, 's': False, 't': 37, 'c': 0.15, 'w': 0.50, 'd': 0.30, 'v': 0.60},
|
||||
]
|
||||
|
||||
cachedir = Path('/tmp/coil_test_cache')
|
||||
version_string = 'v1.0'
|
||||
coil_border = 7 # mm
|
||||
cut_gap = 8 # mm
|
||||
tooling_border = 10 # mm
|
||||
vscore_extra = 10 # mm
|
||||
mouse_bite_width = 8 # mm
|
||||
mouse_bite_yoff = 0.175
|
||||
mouse_bite_hole_dia = 0.7
|
||||
mouse_bite_hole_spacing = 0.7
|
||||
hole_offset = 5
|
||||
hole_dia = 3.2
|
||||
coil_dia = 35 # mm
|
||||
coil_inner_dia = 15 # mm
|
||||
board_thickness = 0.80 # mm
|
||||
pad_offset = 2 # mm
|
||||
pad_dia = 2.0 # mm
|
||||
pad_length = 3.5 # mm
|
||||
pad_drill = 1.1 # mm
|
||||
pad_pitch = 2.54 # mm
|
||||
join_trace_w = 0.150 # mm
|
||||
do_v_cuts = False
|
||||
do_mouse_bites = False
|
||||
do_cut_gaps = False
|
||||
|
||||
db = sqlite3.connect('coil_parameters.sqlite3')
|
||||
db.execute('CREATE TABLE IF NOT EXISTS runs (run_id INTEGER PRIMARY KEY, timestamp TEXT, version TEXT)')
|
||||
db.execute('CREATE TABLE IF NOT EXISTS coils (coil_id INTEGER PRIMARY KEY, run_id INTEGER, FOREIGN KEY (run_id) REFERENCES runs(run_id))')
|
||||
db.execute('CREATE TABLE IF NOT EXISTS results (result_id INTEGER PRIMARY KEY, coil_id INTEGER, key TEXT, value TEXT, FOREIGN KEY (coil_id) REFERENCES coils(coil_id))')
|
||||
cur = db.cursor()
|
||||
cur.execute('INSERT INTO runs(timestamp, version) VALUES (datetime("now"), ?)', (version_string,))
|
||||
run_id = cur.lastrowid
|
||||
db.commit()
|
||||
|
||||
tile_width = tile_height = coil_dia + 2*coil_border
|
||||
coil_pitch_v = tile_width + cut_gap
|
||||
coil_pitch_h = tile_height + cut_gap
|
||||
|
||||
total_width = coil_pitch_h*cols + 2*tooling_border + cut_gap
|
||||
total_height = coil_pitch_v*rows + 2*tooling_border + cut_gap
|
||||
|
||||
drawing_text_size = 2.0
|
||||
|
||||
print(f'Calculated board size: {total_width:.2f} * {total_height:.2f} mm')
|
||||
print(f'Tile size: {tile_height:.2f} * {tile_height:.2f} mm')
|
||||
|
||||
x0, y0 = 100, 100
|
||||
|
||||
xy = pcb.XYCoord
|
||||
b = pcb.Board.empty_board(page=pcb.PageSettings(page_format='A2'))
|
||||
|
||||
b.add(kc_gr.Rectangle(xy(x0, y0), xy(x0+total_width, y0+total_height), layer='Edge.Cuts', stroke=pcb.Stroke(width=0.15)))
|
||||
|
||||
def do_line(x0, y0, x1, y1, off_x=0, off_y=0):
|
||||
b.add(kc_gr.Line(xy(x0+off_x, y0+off_y),
|
||||
xy(x1+off_x, y1+off_y),
|
||||
layer='Edge.Cuts', stroke=pcb.Stroke(width=0.15)))
|
||||
|
||||
if do_v_cuts:
|
||||
for y in range(rows):
|
||||
for off_y in [0, tile_height]:
|
||||
y_pos = y0 + tooling_border + cut_gap + off_y + y*coil_pitch_v
|
||||
do_line(x0 - vscore_extra, y_pos, x0 + total_width + vscore_extra, y_pos)
|
||||
b.add(kc_gr.Text(text='V-score',
|
||||
at=pcb.AtPos(x0 + total_width + vscore_extra + drawing_text_size/2, y_pos, 0),
|
||||
layer=kc_gr.TextLayer('Edge.Cuts'),
|
||||
effects=pcb.TextEffect(
|
||||
font=pcb.FontSpec(size=xy(drawing_text_size, drawing_text_size),
|
||||
thickness=drawing_text_size/10),
|
||||
justify=pcb.Justify(h=pcb.Atom.left))))
|
||||
|
||||
|
||||
for x in range(cols):
|
||||
for off_x in [0, tile_width]:
|
||||
x_pos = x0 + tooling_border + cut_gap + off_x + x*coil_pitch_h
|
||||
do_line(x_pos, y0 - vscore_extra, x_pos, y0 + total_height + vscore_extra)
|
||||
b.add(kc_gr.Text(text='V-score',
|
||||
at=pcb.AtPos(x_pos, y0 + total_height + vscore_extra + drawing_text_size/2, 90),
|
||||
layer=kc_gr.TextLayer('Edge.Cuts'),
|
||||
effects=pcb.TextEffect(
|
||||
font=pcb.FontSpec(size=xy(drawing_text_size, drawing_text_size),
|
||||
thickness=drawing_text_size/10),
|
||||
justify=pcb.Justify(h=pcb.Atom.right))))
|
||||
|
||||
def draw_corner(x0, y0, spokes):
|
||||
right, top, left, bottom = [True if c.lower() in 'y1' else False for c in spokes]
|
||||
|
||||
l = (tile_width - mouse_bite_width)/2 - cut_gap/2
|
||||
|
||||
if right:
|
||||
do_line(cut_gap/2, -cut_gap/2, cut_gap/2 + l, -cut_gap/2, x0, y0)
|
||||
do_line(cut_gap/2, cut_gap/2, cut_gap/2 + l, cut_gap/2, x0, y0)
|
||||
b.add(kc_gr.Arc(start=xy(x0+cut_gap/2+l, y0-cut_gap/2),
|
||||
end=xy(x0+cut_gap/2+l, y0+cut_gap/2),
|
||||
center=xy(x0+cut_gap/2+l, y0),
|
||||
layer='Edge.Cuts',
|
||||
stroke=pcb.Stroke(width=0.15)))
|
||||
|
||||
else:
|
||||
do_line(cut_gap/2, -cut_gap/2, cut_gap/2, cut_gap/2, x0, y0)
|
||||
|
||||
if left:
|
||||
do_line(-cut_gap/2, -cut_gap/2, -cut_gap/2 - l, -cut_gap/2, x0, y0)
|
||||
do_line(-cut_gap/2, cut_gap/2, -cut_gap/2 - l, cut_gap/2, x0, y0)
|
||||
b.add(kc_gr.Arc(end=xy(x0-cut_gap/2-l, y0-cut_gap/2),
|
||||
start=xy(x0-cut_gap/2-l, y0+cut_gap/2),
|
||||
center=xy(x0-cut_gap/2-l, y0),
|
||||
layer='Edge.Cuts',
|
||||
stroke=pcb.Stroke(width=0.15)))
|
||||
|
||||
else:
|
||||
do_line(-cut_gap/2, -cut_gap/2, -cut_gap/2, cut_gap/2, x0, y0)
|
||||
|
||||
if bottom:
|
||||
do_line(-cut_gap/2, cut_gap/2, -cut_gap/2, cut_gap/2 + l, x0, y0)
|
||||
do_line(cut_gap/2, cut_gap/2, cut_gap/2, cut_gap/2 + l, x0, y0)
|
||||
b.add(kc_gr.Arc(end=xy(x0-cut_gap/2, y0+cut_gap/2+l),
|
||||
start=xy(x0+cut_gap/2, y0+cut_gap/2+l),
|
||||
center=xy(x0, y0+cut_gap/2+l),
|
||||
layer='Edge.Cuts',
|
||||
stroke=pcb.Stroke(width=0.15)))
|
||||
|
||||
else:
|
||||
do_line(-cut_gap/2, cut_gap/2, cut_gap/2, cut_gap/2, x0, y0)
|
||||
|
||||
if top:
|
||||
do_line(-cut_gap/2, -cut_gap/2, -cut_gap/2, -cut_gap/2 - l, x0, y0)
|
||||
do_line(cut_gap/2, -cut_gap/2, cut_gap/2, -cut_gap/2 - l, x0, y0)
|
||||
b.add(kc_gr.Arc(start=xy(x0-cut_gap/2, y0-cut_gap/2-l),
|
||||
end=xy(x0+cut_gap/2, y0-cut_gap/2-l),
|
||||
center=xy(x0, y0-cut_gap/2-l),
|
||||
layer='Edge.Cuts',
|
||||
stroke=pcb.Stroke(width=0.15)))
|
||||
else:
|
||||
|
||||
do_line(-cut_gap/2, -cut_gap/2, cut_gap/2, -cut_gap/2, x0, y0)
|
||||
|
||||
|
||||
def make_mouse_bite(x, y, rot=0, width=mouse_bite_width, hole_dia=mouse_bite_hole_dia, hole_spacing=mouse_bite_hole_spacing, **kwargs):
|
||||
|
||||
pitch = hole_dia + hole_spacing
|
||||
num_holes = int(math.floor((width - hole_dia) / pitch)) + 1
|
||||
|
||||
actual_spacing = (width - num_holes*hole_dia) / (num_holes - 1)
|
||||
pitch = hole_dia + actual_spacing
|
||||
|
||||
f = fp.Footprint(name='mouse_bite', _version=None, generator=None, at=fp.AtPos(x, y, rot), **kwargs)
|
||||
for i in range(num_holes):
|
||||
f.pads.append(fp.Pad(
|
||||
number='1',
|
||||
type=fp.Atom.np_thru_hole,
|
||||
shape=fp.Atom.circle,
|
||||
at=fp.AtPos(-width/2 + i*pitch + hole_dia/2, 0, 0),
|
||||
size=xy(hole_dia, hole_dia),
|
||||
drill=fp.Drill(diameter=hole_dia),
|
||||
footprint=f))
|
||||
return f
|
||||
|
||||
|
||||
def make_hole(x, y, dia, **kwargs):
|
||||
f = fp.Footprint(name='hole', _version=None, generator=None, at=fp.AtPos(x, y, 0), **kwargs)
|
||||
f.pads.append(fp.Pad(
|
||||
number='1',
|
||||
type=fp.Atom.np_thru_hole,
|
||||
shape=fp.Atom.circle,
|
||||
at=fp.AtPos(0, 0, 0),
|
||||
size=xy(dia, dia),
|
||||
drill=fp.Drill(diameter=dia),
|
||||
footprint=f))
|
||||
return f
|
||||
|
||||
|
||||
def make_pads(x, y, rot, n, pad_dia, pad_length, drill, pitch, **kwargs):
|
||||
f = fp.Footprint(name=f'conn_gen_01x{n}', _version=None, generator=None, at=fp.AtPos(x, y, rot), **kwargs)
|
||||
|
||||
for i in range(n):
|
||||
f.pads.append(fp.Pad(
|
||||
number=str(i+1),
|
||||
type=fp.Atom.thru_hole,
|
||||
shape=fp.Atom.oval,
|
||||
at=fp.AtPos(-pitch*(n-1)/2 + i*pitch, 0, rot),
|
||||
size=xy(pad_dia, pad_length),
|
||||
drill=fp.Drill(diameter=drill),
|
||||
footprint=f))
|
||||
|
||||
return f
|
||||
|
||||
|
||||
corner_x0 = x0 + tooling_border + cut_gap/2
|
||||
corner_y0 = y0 + tooling_border + cut_gap/2
|
||||
corner_x1 = x0 + total_width - tooling_border - cut_gap/2
|
||||
corner_y1 = y0 + total_height - tooling_border - cut_gap/2
|
||||
|
||||
if do_cut_gaps:
|
||||
# Corners
|
||||
draw_corner(corner_x0, corner_y0, 'YNNY')
|
||||
draw_corner(corner_x0, corner_y1, 'YYNN')
|
||||
draw_corner(corner_x1, corner_y0, 'NNYY')
|
||||
draw_corner(corner_x1, corner_y1, 'NYYN')
|
||||
|
||||
# Top / bottom T junctions
|
||||
for x in range(1, cols):
|
||||
draw_corner(corner_x0 + x*coil_pitch_h, corner_y0, 'YYNY')
|
||||
draw_corner(corner_x0 + x*coil_pitch_h, corner_y1, 'NYYY')
|
||||
|
||||
# Left / right T junctions
|
||||
for y in range(1, rows):
|
||||
draw_corner(corner_x0, corner_y0 + y*coil_pitch_v, 'YYNY')
|
||||
draw_corner(corner_x1, corner_y0 + y*coil_pitch_v, 'NYYY')
|
||||
|
||||
# Middle X junctions
|
||||
for y in range(1, rows):
|
||||
for x in range(1, cols):
|
||||
draw_corner(corner_x0 + x*coil_pitch_h, corner_y0 + y*coil_pitch_v, 'YYYY')
|
||||
|
||||
else:
|
||||
for layer in ('F.SilkS', 'B.SilkS'):
|
||||
for x in range(0, cols+1):
|
||||
cx = x0 + tooling_border + cut_gap/2 + x*coil_pitch_h
|
||||
b.add(kc_gr.Line(xy(cx, corner_y0),
|
||||
xy(cx, corner_y1),
|
||||
layer=layer, stroke=pcb.Stroke(width=0.15)))
|
||||
|
||||
for y in range(0, rows+1):
|
||||
cy = y0 + tooling_border + cut_gap/2 + y*coil_pitch_v
|
||||
b.add(kc_gr.Line(xy(corner_x0, cy),
|
||||
xy(corner_x1, cy),
|
||||
layer=layer, stroke=pcb.Stroke(width=0.15)))
|
||||
|
||||
|
||||
# Mouse bites
|
||||
if do_mouse_bites:
|
||||
for x in range(0, cols):
|
||||
for y in range(0, rows):
|
||||
tile_x0 = x0 + tooling_border + cut_gap + x*coil_pitch_h
|
||||
tile_y0 = y0 + tooling_border + cut_gap + y*coil_pitch_v
|
||||
|
||||
b.add(make_mouse_bite(tile_x0 + tile_width/2, tile_y0 - mouse_bite_hole_dia/2, 0))
|
||||
b.add(make_mouse_bite(tile_x0 + tile_width/2, tile_y0 + tile_height + mouse_bite_hole_dia/2, 0))
|
||||
b.add(make_mouse_bite(tile_x0 - mouse_bite_hole_dia/2, tile_y0 + tile_height/2, 90))
|
||||
b.add(make_mouse_bite(tile_x0 + tile_width + mouse_bite_hole_dia/2, tile_y0 + tile_height/2, 90))
|
||||
|
||||
# Mounting holes
|
||||
for x in range(0, cols):
|
||||
for y in range(0, rows):
|
||||
tile_x0 = x0 + tooling_border + cut_gap + x*coil_pitch_h + tile_width/2
|
||||
tile_y0 = y0 + tooling_border + cut_gap + y*coil_pitch_v + tile_height/2
|
||||
|
||||
dx = tile_width/2 - hole_offset
|
||||
dy = tile_height/2 - hole_offset
|
||||
b.add(make_hole(tile_x0 - dx, tile_y0 - dy, hole_dia))
|
||||
b.add(make_hole(tile_x0 - dx, tile_y0 + dy, hole_dia))
|
||||
b.add(make_hole(tile_x0 + dx, tile_y0 - dy, hole_dia))
|
||||
b.add(make_hole(tile_x0 + dx, tile_y0 + dy, hole_dia))
|
||||
|
||||
# border graphics
|
||||
c = 3
|
||||
for layer in ['F.SilkS', 'B.SilkS']:
|
||||
b.add(kc_gr.Rectangle(start=xy(x0, y0), end=xy(x0+c, y0+total_height), layer=layer, stroke=pcb.Stroke(width=0),
|
||||
fill=kc_gr.FillMode(pcb.Atom.solid)))
|
||||
b.add(kc_gr.Rectangle(start=xy(x0, y0), end=xy(x0+total_width, y0+c), layer=layer, stroke=pcb.Stroke(width=0),
|
||||
fill=kc_gr.FillMode(pcb.Atom.solid)))
|
||||
b.add(kc_gr.Rectangle(start=xy(x0+total_width-c, y0), end=xy(x0+total_width, y0+total_height), layer=layer, stroke=pcb.Stroke(width=0),
|
||||
fill=kc_gr.FillMode(pcb.Atom.solid)))
|
||||
b.add(kc_gr.Rectangle(start=xy(x0, y0+total_height-c), end=xy(x0+total_width, y0+total_height), layer=layer, stroke=pcb.Stroke(width=0),
|
||||
fill=kc_gr.FillMode(pcb.Atom.solid)))
|
||||
|
||||
a = 3
|
||||
timestamp = datetime.datetime.now().strftime('%Y-%m-%d')
|
||||
b.add(kc_gr.Text(text=f'Planar inductor test panel',
|
||||
at=pcb.AtPos(x0 + tooling_border + cut_gap/2, y0 + c + 2*a/3),
|
||||
layer=kc_gr.TextLayer('F.SilkS'),
|
||||
effects=pcb.TextEffect(
|
||||
font=pcb.FontSpec(face="Inter Semi Bold",
|
||||
size=xy(6*a/3, 6*a/3),
|
||||
thickness=a/5),
|
||||
justify=pcb.Justify(h=pcb.Atom.left, v=pcb.Atom.top))))
|
||||
|
||||
b.add(kc_gr.Text(text=f'{version_string} {timestamp} © 2023 Jan Götte, FG KOM, TU Darmstadt',
|
||||
at=pcb.AtPos(x0 + total_width - tooling_border - cut_gap/2, y0 + c + 4*a/3),
|
||||
layer=kc_gr.TextLayer('F.SilkS'),
|
||||
effects=pcb.TextEffect(
|
||||
font=pcb.FontSpec(face="Inter Light",
|
||||
size=xy(a, a),
|
||||
thickness=a/5),
|
||||
justify=pcb.Justify(h=pcb.Atom.right, v=pcb.Atom.top))))
|
||||
|
||||
for index, ((y, x), spec) in tqdm.tqdm(enumerate(zip(itertools.product(range(rows), range(cols)), coil_specs), start=1)):
|
||||
pass
|
||||
with tempfile.NamedTemporaryFile(suffix='.kicad_mod') as f:
|
||||
tile_x0 = x0 + tooling_border + cut_gap + x*coil_pitch_h + tile_width/2
|
||||
tile_y0 = y0 + tooling_border + cut_gap + y*coil_pitch_v + tile_height/2
|
||||
|
||||
for key, alias in {
|
||||
'gen.inner_diameter': 'id',
|
||||
'gen.outer_diameter': 'od',
|
||||
'gen.trace_width': 'w',
|
||||
'gen.turns': 'n',
|
||||
'gen.twists': 't',
|
||||
'gen.clearance': 'c',
|
||||
'gen.single_layer': 's',
|
||||
'gen.via_drill': 'd',
|
||||
'gen.via_diameter': 'v'}.items():
|
||||
if alias in spec:
|
||||
spec[key] = spec.pop(alias)
|
||||
|
||||
if 'gen.via_diameter' not in spec:
|
||||
spec['gen.via_diameter'] = spec['gen.trace_width']
|
||||
|
||||
if 'gen.inner_diameter' not in spec:
|
||||
spec['gen.inner_diameter'] = coil_inner_dia
|
||||
|
||||
if 'gen.outer_diameter' not in spec:
|
||||
spec['gen.outer_diameter'] = coil_dia
|
||||
|
||||
args = ['python', '-m', 'twisted_coil_gen_twolayer', '--no-keepout-zone']
|
||||
for k, v in spec.items():
|
||||
prefix, _, k = k.partition('.')
|
||||
if (not isinstance(v, bool) or v) and prefix == 'gen':
|
||||
args.append('--' + k.replace('_', '-'))
|
||||
if v is not True:
|
||||
args.append(str(v))
|
||||
|
||||
arg_digest = hashlib.sha3_256(' / '.join(map(str, args)).encode()).hexdigest()
|
||||
cachedir.mkdir(exist_ok=True)
|
||||
cache_file = cachedir / f'C-{arg_digest}.kicad_mod'
|
||||
log_file = cachedir / f'Q-{arg_digest}.kicad_mod'
|
||||
if not cache_file.is_file():
|
||||
args.append(cache_file)
|
||||
try:
|
||||
res = subprocess.run(args, check=True, capture_output=True, text=True)
|
||||
log_file.write_text(res.stdout + res.stderr)
|
||||
except subprocess.CalledProcessError as e:
|
||||
print(f'Error generating coil with command line {args}, rc={e.returncode}')
|
||||
print(e.stdout)
|
||||
print(e.stderr)
|
||||
|
||||
coil = fp.Footprint.open_mod(cache_file)
|
||||
coil.at = fp.AtPos(tile_x0, tile_y0, 0)
|
||||
b.add(coil)
|
||||
|
||||
t = [f'n={spec["gen.turns"]}',
|
||||
f'{spec["gen.twists"]} twists',
|
||||
f'w={spec["gen.trace_width"]:.2f}mm']
|
||||
if spec.get('gen.single_layer'):
|
||||
t.append('single layer')
|
||||
|
||||
spec['gen.board_thickness'] = board_thickness
|
||||
cur.execute('INSERT INTO coils(run_id) VALUES (?)', (run_id,))
|
||||
coil_id = cur.lastrowid
|
||||
|
||||
for key, value in spec.items():
|
||||
if isinstance(value, bool):
|
||||
value = str(value)
|
||||
db.execute('INSERT INTO results(coil_id, key, value) VALUES (?, ?, ?)', (coil_id, key, value))
|
||||
|
||||
for l in log_file.read_text().splitlines():
|
||||
if (m := re.fullmatch(r'Approximate inductance:\s*([-+.0-9eE]+)\s*µH', l.strip())):
|
||||
val = float(m.group(1)) * 1e-6
|
||||
db.execute('INSERT INTO results(coil_id, key, value) VALUES (?, "calculated_approximate_inductance", ?)', (coil_id, val))
|
||||
if (m := re.fullmatch(r'Approximate track length:\s*([-+.0-9eE]+)\s*mm', l.strip())):
|
||||
val = float(m.group(1)) * 1e-3
|
||||
db.execute('INSERT INTO results(coil_id, key, value) VALUES (?, "calculated_trace_length", ?)', (coil_id, val))
|
||||
if (m := re.fullmatch(r'Approximate resistance:\s*([-+.0-9eE]+)\s*Ω', l.strip())):
|
||||
val = float(m.group(1))
|
||||
db.execute('INSERT INTO results(coil_id, key, value) VALUES (?, "calculated_approximate_resistance", ?)', (coil_id, val))
|
||||
if (m := re.fullmatch(r'Fill factor:\s*([-+.0-9eE]+)', l.strip())):
|
||||
val = float(m.group(1))
|
||||
db.execute('INSERT INTO results(coil_id, key, value) VALUES (?, "calculated_fill_factor", ?)', (coil_id, val))
|
||||
db.commit()
|
||||
|
||||
sz = 2
|
||||
b.add(kc_gr.Text(text='\\n'.join(t),
|
||||
at=pcb.AtPos(tile_x0, tile_y0),
|
||||
layer=kc_gr.TextLayer('B.SilkS'),
|
||||
effects=pcb.TextEffect(
|
||||
font=pcb.FontSpec(face='Inter Medium',
|
||||
size=xy(sz, sz),
|
||||
thickness=sz/5),
|
||||
justify=pcb.Justify(h=None, v=None, mirror=True))))
|
||||
|
||||
b.add(kc_gr.Text(text=f'Tile {index}',
|
||||
at=pcb.AtPos(tile_x0, tile_y0 - tile_height/2 + sz),
|
||||
layer=kc_gr.TextLayer('B.SilkS'),
|
||||
effects=pcb.TextEffect(
|
||||
font=pcb.FontSpec(face='Inter Semi Bold',
|
||||
size=xy(sz, sz),
|
||||
thickness=sz/5),
|
||||
justify=pcb.Justify(h=None, v=pcb.Atom.top, mirror=True))))
|
||||
|
||||
b.add(kc_gr.Text(text=f'{version_string} {timestamp}',
|
||||
at=pcb.AtPos(tile_x0, tile_y0 - tile_height/2 + sz*2.4),
|
||||
layer=kc_gr.TextLayer('B.SilkS'),
|
||||
effects=pcb.TextEffect(
|
||||
font=pcb.FontSpec(face='Inter Light',
|
||||
size=xy(sz, sz),
|
||||
thickness=sz/5),
|
||||
justify=pcb.Justify(h=None, v=pcb.Atom.top, mirror=True))))
|
||||
|
||||
b.add(kc_gr.Text(text=f'{index}',
|
||||
at=pcb.AtPos(tile_x0, tile_y0 - tile_height/2 + sz),
|
||||
layer=kc_gr.TextLayer('F.SilkS'),
|
||||
effects=pcb.TextEffect(
|
||||
font=pcb.FontSpec(face='Inter Medium',
|
||||
size=xy(sz, sz),
|
||||
thickness=sz/5),
|
||||
justify=pcb.Justify(h=None, v=pcb.Atom.top, mirror=False))))
|
||||
|
||||
pads_x0 = tile_x0 + tile_width/2 - pad_offset
|
||||
pads = make_pads(pads_x0, tile_y0, 270, 2, pad_dia, pad_length, pad_drill, pad_pitch)
|
||||
b.add(pads)
|
||||
|
||||
w = min(spec.get('gen.trace_width', pad_dia), pad_dia)
|
||||
wx, wy, _r, _f = pads.pad(2).abs_pos
|
||||
w2 = (wx - pad_length/2, wy)
|
||||
wx, wy, _r, _f = pads.pad(1).abs_pos
|
||||
w1 = (wx - pad_length/2, wy)
|
||||
b.add(cad_pr.Trace(w, coil.pad(1), pads.pad(1), waypoints=[w1], orientation=['ccw'], side='top'))
|
||||
b.add(cad_pr.Trace(w, coil.pad(2), pads.pad(2), waypoints=[w2], orientation=['cw'], side='bottom'))
|
||||
|
||||
k = 3
|
||||
for layer in ['F.SilkS', 'B.SilkS']:
|
||||
b.add(kc_gr.Rectangle(start=xy(wx-k/2, wy-pad_pitch-k/2), end=xy(wx+k/2, wy-pad_pitch), layer=layer, stroke=pcb.Stroke(width=0),
|
||||
fill=kc_gr.FillMode(pcb.Atom.solid)))
|
||||
|
||||
b.write('coil_test_board.kicad_pcb')
|
||||
|
||||
|
|
@ -1,173 +0,0 @@
|
|||
Header
|
||||
CHECK KEYWORDS "Warn"
|
||||
Mesh DB "." "mesh"
|
||||
End
|
||||
|
||||
Simulation
|
||||
Mesh Levels = 1
|
||||
Max Output Level = 7
|
||||
Coordinate System = Cartesian
|
||||
Coordinate Mapping(3) = 1 2 3
|
||||
Simulation Type = Steady state
|
||||
Steady State Max Iterations = 1
|
||||
Output Intervals = 1
|
||||
Timestepping Method = BDF
|
||||
Simulation Timing = True
|
||||
BDF Order = 1
|
||||
Solver Input File = case.sif
|
||||
Post File = case.vtu
|
||||
Output File = case.result
|
||||
End
|
||||
|
||||
Constants
|
||||
Stefan Boltzmann = 5.6704e-08
|
||||
Permittivity of Vacuum = 8.8541878128e-12
|
||||
Gravity(4) = 0 -1 0 9.80665
|
||||
Boltzmann Constant = 1.380649e-23
|
||||
Unit Charge = 1.602176634e-19
|
||||
End
|
||||
|
||||
! airEqn
|
||||
Equation 1
|
||||
Active Solvers(1) = 2
|
||||
End
|
||||
|
||||
! copperEqn
|
||||
Equation 2
|
||||
Active Solvers(1) = 1 2
|
||||
End
|
||||
|
||||
Solver 1
|
||||
Equation = Static Current Conduction
|
||||
Variable = Potential
|
||||
Variable DOFs = 1
|
||||
Procedure = "StatCurrentSolve" "StatCurrentSolver"
|
||||
Calculate Volume Current = True
|
||||
Calculate Joule Heating = False
|
||||
Optimize Bandwidth = True
|
||||
Nonlinear System Max Iterations = 1
|
||||
Linear System Solver = Iterative
|
||||
Linear System Iterative Method = CG
|
||||
Linear System Max Iterations = 10000
|
||||
Linear System Convergence Tolerance = 1e-10
|
||||
Linear System Preconditioning = ILU3
|
||||
Linear System ILUT Tolerance = 0.001
|
||||
Linear System Abort Not Converged = False
|
||||
Linear System Residual Output = 20
|
||||
Linear System Precondition Recompute = 1
|
||||
End
|
||||
|
||||
Solver 2
|
||||
Equation = Electrostatics
|
||||
Procedure = "StatElecSolve" "StatElecSolver"
|
||||
Variable = PotentialStat
|
||||
Calculate Electric Field = True
|
||||
Calculate Electric Energy = True
|
||||
Steady State Convergence Tolerance = 1e-05
|
||||
Nonlinear System Convergence Tolerance = 1e-07
|
||||
Nonlinear System Max Iterations = 20
|
||||
Nonlinear System Newton After Iterations = 3
|
||||
Nonlinear System Newton After Tolerance = 0.001
|
||||
Nonlinear System Relaxation Factor = 1
|
||||
Linear System Solver = Iterative
|
||||
Linear System Iterative Method = BiCGStab
|
||||
Linear System Max Iterations = 500
|
||||
Linear System Convergence Tolerance = 1e-10
|
||||
BiCGstabl polynomial degree = 2
|
||||
Linear System Preconditioning = none
|
||||
Linear System ILUT Tolerance = 0.001
|
||||
Linear System Abort Not Converged = False
|
||||
Linear System Residual Output = 10
|
||||
End
|
||||
|
||||
|
||||
|
||||
! air
|
||||
Material 1
|
||||
Density = 1.1885
|
||||
Electric Conductivity = 0.0
|
||||
Heat Capacity = 1006.4
|
||||
Heat Conductivity = 0.025873
|
||||
Relative Permeability = 1
|
||||
Relative Permittivity = 1
|
||||
End
|
||||
|
||||
! ro4003c
|
||||
Material 2
|
||||
Density = 1790
|
||||
Electric Conductivity = 0.0
|
||||
Relative Permeability = 1
|
||||
Relative Permittivity = 3.55
|
||||
End
|
||||
|
||||
! copper
|
||||
Material 3
|
||||
Density = 8960.0
|
||||
Electric Conductivity = 59600000
|
||||
Emissivity = 0.012
|
||||
Heat Capacity = 415.0
|
||||
Heat Conductivity = 401.0
|
||||
Relative Permeability = 1
|
||||
End
|
||||
|
||||
! trace
|
||||
Body 1
|
||||
Target Bodies(1) = 3
|
||||
Equation = 2 ! copperEqn
|
||||
Material = 3 ! copper
|
||||
Body Force = 1 ! electric_potential
|
||||
End
|
||||
|
||||
! substrate
|
||||
Body 2
|
||||
Target Bodies(1) = 1
|
||||
Equation = 1 ! airEqn
|
||||
Material = 2 ! ro4003c
|
||||
End
|
||||
|
||||
! airbox
|
||||
Body 3
|
||||
Target Bodies(1) = 2
|
||||
Equation = 1 ! airEqn
|
||||
Material = 1 ! air
|
||||
Body Force = 1
|
||||
End
|
||||
|
||||
! interface_top
|
||||
Body 4
|
||||
Target Bodies(1) = 4
|
||||
Material = 3 ! copper
|
||||
End
|
||||
|
||||
! interface_bottom
|
||||
Body 5
|
||||
Target Bodies(1) = 5
|
||||
Material = 3 ! copper
|
||||
End
|
||||
|
||||
|
||||
! FarField
|
||||
Boundary Condition 1
|
||||
Target Boundaries(1) = 6
|
||||
Electric Infinity BC = True
|
||||
End
|
||||
|
||||
! Vplus
|
||||
Boundary Condition 2
|
||||
Target Boundaries(1) = 4
|
||||
Potential = 1.0
|
||||
Save Scalars = True
|
||||
End
|
||||
|
||||
! Vminus
|
||||
Boundary Condition 3
|
||||
Target Boundaries(1) = 5
|
||||
Potential = 0.0
|
||||
End
|
||||
|
||||
|
||||
! electric_potential
|
||||
Body Force 1
|
||||
PotentialStat = Equals "Potential"
|
||||
End
|
||||
|
||||
|
|
@ -1,154 +0,0 @@
|
|||
Header
|
||||
CHECK KEYWORDS "Warn"
|
||||
Mesh DB "." "mesh"
|
||||
End
|
||||
|
||||
Simulation
|
||||
Mesh Levels = 1
|
||||
Max Output Level = 7
|
||||
Coordinate System = Cartesian
|
||||
Coordinate Mapping(3) = 1 2 3
|
||||
Simulation Type = Steady state
|
||||
Steady State Max Iterations = 1
|
||||
Output Intervals = 1
|
||||
Timestepping Method = BDF
|
||||
Simulation Timing = True
|
||||
BDF Order = 1
|
||||
Solver Input File = case.sif
|
||||
Post File = case.vtu
|
||||
Output File = case.result
|
||||
End
|
||||
|
||||
Constants
|
||||
Stefan Boltzmann = 5.6704e-08
|
||||
Permittivity of Vacuum = 8.8541878128e-12
|
||||
Gravity(4) = 0 -1 0 9.80665
|
||||
Boltzmann Constant = 1.380649e-23
|
||||
Unit Charge = 1.602176634e-19
|
||||
End
|
||||
|
||||
! airEqn
|
||||
Equation 1
|
||||
Active Solvers(1) = 1 ! Static_Current_Conduction, Magneto_Dynamics, Magneto_Dynamics_Calculations,
|
||||
End
|
||||
|
||||
! Static_Current_Conduction
|
||||
Solver 1
|
||||
Equation = Electrostatics
|
||||
Procedure = "StatElecSolve" "StatElecSolver"
|
||||
Variable = Potential
|
||||
Calculate Electric Field = True
|
||||
Calculate Electric Energy = True
|
||||
Exec Solver = Always
|
||||
Stabilize = True
|
||||
Bubbles = False
|
||||
Lumped Mass Matrix = False
|
||||
Optimize Bandwidth = True
|
||||
Steady State Convergence Tolerance = 1e-05
|
||||
Nonlinear System Convergence Tolerance = 1e-07
|
||||
Nonlinear System Max Iterations = 20
|
||||
Nonlinear System Newton After Iterations = 3
|
||||
Nonlinear System Newton After Tolerance = 0.001
|
||||
Nonlinear System Relaxation Factor = 1
|
||||
Linear System Solver = Iterative
|
||||
Linear System Iterative Method = BiCGStab
|
||||
Linear System Max Iterations = 500
|
||||
Linear System Convergence Tolerance = 1e-10
|
||||
BiCGstabl polynomial degree = 2
|
||||
Linear System Preconditioning = ILU0
|
||||
Linear System ILUT Tolerance = 0.001
|
||||
Linear System Abort Not Converged = False
|
||||
Linear System Residual Output = 10
|
||||
Linear System Precondition Recompute = 1
|
||||
End
|
||||
|
||||
|
||||
! air
|
||||
Material 1
|
||||
Density = 1.1885
|
||||
Electric Conductivity = 0.0
|
||||
Heat Capacity = 1006.4
|
||||
Heat Conductivity = 0.025873
|
||||
Relative Permeability = 1
|
||||
Relative Permittivity = 1
|
||||
End
|
||||
|
||||
! ro4003c
|
||||
Material 2
|
||||
Density = 1790
|
||||
Electric Conductivity = 0.0
|
||||
Relative Permeability = 1
|
||||
Relative Permittivity = 3.55
|
||||
End
|
||||
|
||||
! copper
|
||||
Material 3
|
||||
Density = 8960.0
|
||||
Electric Conductivity = 59600000
|
||||
Emissivity = 0.012
|
||||
Heat Capacity = 415.0
|
||||
Heat Conductivity = 401.0
|
||||
Relative Permeability = 1
|
||||
End
|
||||
|
||||
! trace
|
||||
Body 1
|
||||
Target Bodies(1) = 3
|
||||
Material = 3 ! copper
|
||||
Body Force = 1 ! electric_potential
|
||||
End
|
||||
|
||||
! substrate
|
||||
Body 2
|
||||
Target Bodies(1) = 1
|
||||
Equation = 1 ! airEqn
|
||||
Material = 2 ! ro4003c
|
||||
End
|
||||
|
||||
! airbox
|
||||
Body 3
|
||||
Target Bodies(1) = 2
|
||||
Equation = 1 ! airEqn
|
||||
Material = 1 ! air
|
||||
End
|
||||
|
||||
! interface_top
|
||||
Body 4
|
||||
Target Bodies(1) = 4
|
||||
Material = 3 ! copper
|
||||
End
|
||||
|
||||
! interface_bottom
|
||||
Body 5
|
||||
Target Bodies(1) = 5
|
||||
Material = 3 ! copper
|
||||
End
|
||||
|
||||
|
||||
! FarField
|
||||
Boundary Condition 1
|
||||
Target Boundaries(1) = 6
|
||||
Electric Infinity BC = True
|
||||
End
|
||||
|
||||
! Vplus
|
||||
Boundary Condition 2
|
||||
Target Boundaries(1) = 4
|
||||
Potential = 1.0
|
||||
Save Scalars = True
|
||||
End
|
||||
|
||||
! Vminus
|
||||
Boundary Condition 3
|
||||
Target Boundaries(1) = 5
|
||||
Potential = 0.0
|
||||
End
|
||||
|
||||
|
||||
! electric_potential
|
||||
Body Force 1
|
||||
Electric Potential = Equals "Potential"
|
||||
End
|
||||
|
||||
|
||||
|
||||
282
sim_runner.py
282
sim_runner.py
|
|
@ -1,282 +0,0 @@
|
|||
#!/usr/bin/env python3
|
||||
|
||||
import threading
|
||||
import queue
|
||||
import itertools
|
||||
import pathlib
|
||||
import tempfile
|
||||
import sys
|
||||
import sqlite3
|
||||
import time
|
||||
import math
|
||||
import json
|
||||
import subprocess
|
||||
|
||||
import tqdm
|
||||
import click
|
||||
from tabulate import tabulate
|
||||
|
||||
|
||||
def mesh_args(db, coil_id, mesh_type, mesh_file, outfile):
|
||||
mesh_type = {'split': '--mesh-split-out', 'normal': '--mesh-out', 'mutual': '--mesh-mutual-out'}[mesh_type]
|
||||
rows = db.execute('SELECT key, value FROM results WHERE coil_id=?', (coil_id,)).fetchall()
|
||||
args = ['python', '-m', 'twisted_coil_gen_twolayer', mesh_type, mesh_file, '--pcb']
|
||||
for k, v in rows:
|
||||
prefix, _, k = k.partition('.')
|
||||
if v != 'False' and prefix == 'gen':
|
||||
args.append('--' + k.replace('_', '-'))
|
||||
if v != 'True':
|
||||
args.append(str(v))
|
||||
args.append(outfile)
|
||||
return args
|
||||
|
||||
|
||||
def get_mesh_file(db, mesh_dir, run_id, coil_id, mesh_type):
|
||||
db.execute('CREATE TABLE IF NOT EXISTS meshes(coil_id INTEGER, mesh_type TEXT, error INTEGER, filename TEXT, timestamp TEXT DEFAULT current_timestamp, FOREIGN KEY (coil_id) REFERENCES coils(coil_id))')
|
||||
|
||||
row = db.execute('SELECT * FROM meshes WHERE coil_id=? AND mesh_type=? ORDER BY timestamp DESC LIMIT 1', (coil_id, mesh_type)).fetchone()
|
||||
if row is not None:
|
||||
mesh_file = mesh_dir / row['filename']
|
||||
if mesh_file.is_file():
|
||||
return mesh_file
|
||||
|
||||
timestamp = time.strftime('%Y-%m-%d_%H-%M-%S')
|
||||
return mesh_dir / f'mesh-{run_id}-{coil_id}-{mesh_type}-{timestamp}.msh'
|
||||
|
||||
|
||||
def ensure_mesh(db, mesh_dir, log_dir, run_id, coil_id, mesh_type):
|
||||
mesh_file = get_mesh_file(db, mesh_dir, run_id, coil_id, mesh_type)
|
||||
|
||||
if mesh_file.is_file():
|
||||
return mesh_file
|
||||
|
||||
db.execute('INSERT INTO meshes(coil_id, mesh_type, error, filename) VALUES (?, ?, 0, ?)', (coil_id, mesh_type, mesh_file.name))
|
||||
db.commit()
|
||||
|
||||
mesh_file.parent.mkdir(exist_ok=True)
|
||||
with tempfile.NamedTemporaryFile(suffix='.kicad_pcb') as f:
|
||||
args = mesh_args(db, coil_id, mesh_type, mesh_file, f.name)
|
||||
tqdm.tqdm.write(' '.join(map(str, args)))
|
||||
logfile = log_dir / mesh_file.with_suffix('.log').name
|
||||
logfile.parent.mkdir(exist_ok=True)
|
||||
try:
|
||||
res = subprocess.run(args, check=True, capture_output=True, text=True)
|
||||
logfile.write_text(res.stdout + res.stderr)
|
||||
|
||||
except subprocess.CalledProcessError as e:
|
||||
print('Mesh generation failed with exit code {e.returncode}', file=sys.stderr)
|
||||
logfile.write_text(e.stdout + e.stderr)
|
||||
print(e.stdout + e.stderr)
|
||||
raise
|
||||
|
||||
return mesh_file
|
||||
|
||||
|
||||
@click.group()
|
||||
@click.option('-d', '--database', default='coil_parameters.sqlite3')
|
||||
@click.pass_context
|
||||
def cli(ctx, database):
|
||||
ctx.ensure_object(dict)
|
||||
def connect():
|
||||
db = sqlite3.connect(database)
|
||||
db.row_factory = sqlite3.Row
|
||||
return db
|
||||
ctx.obj['db_connect'] = connect
|
||||
|
||||
|
||||
@cli.command()
|
||||
@click.pass_context
|
||||
def list_runs(ctx):
|
||||
for row in ctx.obj['db_connect']().execute('SELECT * FROM runs ORDER BY timestamp').fetchall():
|
||||
print(row['run_id'], row['timestamp'], row['version'])
|
||||
|
||||
|
||||
@cli.command()
|
||||
@click.pass_context
|
||||
def list_runs(ctx):
|
||||
for row in ctx.obj['db_connect']().execute('SELECT * FROM runs ORDER BY timestamp').fetchall():
|
||||
print(row['run_id'], row['timestamp'], row['version'])
|
||||
|
||||
|
||||
@cli.command()
|
||||
@click.option('-r', '--run-id')
|
||||
@click.option('-m', '--mesh-dir', default='meshes')
|
||||
@click.pass_context
|
||||
def list_coils(ctx, run_id, mesh_dir):
|
||||
db = ctx.obj['db_connect']()
|
||||
if run_id is None:
|
||||
run_id, = db.execute('SELECT run_id FROM runs ORDER BY timestamp DESC LIMIT 1').fetchone()
|
||||
timestamp, = db.execute('SELECT timestamp FROM runs WHERE run_id=?', (run_id,)).fetchone()
|
||||
mesh_dir = pathlib.Path(mesh_dir)
|
||||
|
||||
print(f'Listing meshes for run {run_id} at {timestamp}')
|
||||
print()
|
||||
|
||||
keys = {'gen.turns': 'N',
|
||||
'gen.twists': 'T',
|
||||
'gen.single_layer': '1L',
|
||||
'gen.inner_diameter': 'ID[mm]',
|
||||
'gen.outer_diameter': 'OD[mm]',
|
||||
'calculated_fill_factor': 'Fill factor',
|
||||
'calculated_approximate_inductance': 'L [µH]',
|
||||
'calculated_trace_length': 'track len [mm]',
|
||||
'calculated_approximate_resistance': 'R [mΩ]'}
|
||||
out = []
|
||||
for row in db.execute('SELECT *, MAX(meshes.timestamp) FROM coils LEFT JOIN meshes ON coils.coil_id=meshes.coil_id WHERE run_id=? GROUP BY coils.coil_id, mesh_type ORDER BY meshes.timestamp', (run_id,)).fetchall():
|
||||
if row['timestamp']:
|
||||
if row['error']:
|
||||
state = 'ERROR'
|
||||
elif not (mesh_dir / row['filename']).is_file():
|
||||
state = 'NOT FOUND'
|
||||
else:
|
||||
state = 'SUCCESS'
|
||||
else:
|
||||
state = 'NOT RUN'
|
||||
|
||||
params = dict(db.execute('SELECT key, value FROM results WHERE coil_id=?', (row['coil_id'],)).fetchall())
|
||||
|
||||
if 'calculated_approximate_inductance' in params:
|
||||
params['calculated_approximate_inductance'] = f'{float(params["calculated_approximate_inductance"])*1e6:.02f}'
|
||||
|
||||
if 'calculated_trace_length' in params:
|
||||
params['calculated_trace_length'] = f'{float(params["calculated_trace_length"])*1e3:.03f}'
|
||||
|
||||
if 'calculated_approximate_resistance' in params:
|
||||
params['calculated_approximate_resistance'] = f'{float(params["calculated_approximate_resistance"])*1e3:.03f}'
|
||||
|
||||
if 'calculated_fill_factor' in params:
|
||||
params['calculated_fill_factor'] = f'{float(params["calculated_fill_factor"]):.03f}'
|
||||
|
||||
out.append([row['coil_id'], row['mesh_type'], state, row['timestamp']] + [params.get(key, '-') for key in keys])
|
||||
|
||||
print(tabulate(out, headers=['coil', 'mesh', 'state', 'time'] + list(keys.values()), disable_numparse=True, stralign='right'))
|
||||
|
||||
@cli.command()
|
||||
@click.argument('coil_id', type=int)
|
||||
@click.argument('mesh_type', type=click.Choice(['normal', 'split', 'mutual']))
|
||||
@click.option('--mesh-file', default='/tmp/test.msh')
|
||||
@click.option('--pcb-file', default='/tmp/test.kicad_pcb')
|
||||
@click.pass_context
|
||||
def cmdline(ctx, coil_id, mesh_type, mesh_file, pcb_file):
|
||||
print(' '.join(mesh_args(ctx.obj['db_connect'](), coil_id, mesh_type, mesh_file, pcb_file)))
|
||||
|
||||
@cli.group()
|
||||
@click.option('-r', '--run-id')
|
||||
@click.option('-l', '--log-dir', default='logs')
|
||||
@click.option('-m', '--mesh-dir', default='meshes')
|
||||
@click.pass_context
|
||||
def run(ctx, run_id, log_dir, mesh_dir):
|
||||
if run_id is None:
|
||||
run_id, = ctx.obj['db_connect']().execute('SELECT run_id FROM runs ORDER BY timestamp DESC LIMIT 1').fetchone()
|
||||
ctx.obj['run_id'] = run_id
|
||||
ctx.obj['log_dir'] = pathlib.Path(log_dir)
|
||||
ctx.obj['mesh_dir'] = pathlib.Path(mesh_dir)
|
||||
|
||||
|
||||
@run.command()
|
||||
@click.option('-j', '--num-jobs', type=int, default=1, help='Number of jobs to run in parallel')
|
||||
@click.pass_context
|
||||
def generate_meshes(ctx, num_jobs):
|
||||
db = ctx.obj['db_connect']()
|
||||
rows = [row['coil_id'] for row in db.execute('SELECT coil_id FROM coils WHERE run_id=?', (ctx.obj['run_id'],)).fetchall()]
|
||||
mesh_types = ['split', 'normal', 'mutual']
|
||||
|
||||
params = list(itertools.product(rows, mesh_types))
|
||||
all_files = {get_mesh_file(db, ctx.obj['mesh_dir'], ctx.obj['run_id'], coil_id, mesh_type): (coil_id, mesh_type) for coil_id, mesh_type in params}
|
||||
todo = [(coil_id, mesh_type) for f, (coil_id, mesh_type) in all_files.items() if not f.is_file()]
|
||||
|
||||
q = queue.Queue()
|
||||
for elem in todo:
|
||||
q.put(elem)
|
||||
|
||||
tq = tqdm.tqdm(total=len(todo))
|
||||
def queue_worker():
|
||||
try:
|
||||
while True:
|
||||
coil_id, mesh_type = q.get_nowait()
|
||||
try:
|
||||
ensure_mesh(ctx.obj['db_connect'](), ctx.obj['mesh_dir'], ctx.obj['log_dir'], ctx.obj['run_id'], coil_id, mesh_type)
|
||||
except subprocess.CalledProcessError:
|
||||
tqdm.tqdm.write(f'Error generating {mesh_type} mesh for {coil_id=}')
|
||||
tq.update(1)
|
||||
q.task_done()
|
||||
except queue.Empty:
|
||||
pass
|
||||
|
||||
tqdm.tqdm.write(f'Found {len(params)-len(todo)} meshes out of a total of {len(params)}.')
|
||||
tqdm.tqdm.write(f'Processing the remaining {len(todo)} meshes on {num_jobs} workers in parallel.')
|
||||
threads = []
|
||||
for i in range(num_jobs):
|
||||
t = threading.Thread(target=queue_worker, daemon=True)
|
||||
t.start()
|
||||
threads.append(t)
|
||||
q.join()
|
||||
|
||||
@run.command()
|
||||
@click.option('-j', '--num-jobs', type=int, default=1, help='Number of jobs to run in parallel')
|
||||
@click.pass_context
|
||||
def self_inductance(ctx, num_jobs):
|
||||
db = ctx.obj['db_connect']()
|
||||
|
||||
q = queue.Queue()
|
||||
|
||||
def queue_worker():
|
||||
try:
|
||||
while True:
|
||||
mesh_file, logfile = q.get_nowait()
|
||||
with tempfile.TemporaryDirectory() as tmpdir:
|
||||
try:
|
||||
tqdm.tqdm.write(f'Processing {mesh_file}')
|
||||
res = subprocess.run(['python', '-m', 'coil_parasitics', 'inductance', '--sim-dir', tmpdir, mesh_file], check=True, capture_output=True)
|
||||
logfile.write_text(res.stdout+res.stderr)
|
||||
except subprocess.CalledProcessError as e:
|
||||
print(f'Error running simulation, rc={e.returncode}')
|
||||
logfile.write_text(e.stdout+e.stderr)
|
||||
tq.update(1)
|
||||
q.task_done()
|
||||
except queue.Empty:
|
||||
pass
|
||||
|
||||
num_meshes, num_params, num_completed = 0, 0, 0
|
||||
for coil_id, in db.execute('SELECT coil_id FROM coils WHERE run_id=?', (ctx.obj['run_id'],)).fetchall():
|
||||
num_params += 1
|
||||
mesh_file = get_mesh_file(ctx.obj['db_connect'](), ctx.obj['mesh_dir'], ctx.obj['run_id'], coil_id, 'normal')
|
||||
if mesh_file.is_file():
|
||||
num_meshes += 1
|
||||
logfile = ctx.obj['log_dir'] / (mesh_file.stem + '_elmer_self_inductance.log')
|
||||
if logfile.is_file():
|
||||
num_completed += 1
|
||||
else:
|
||||
q.put((mesh_file, logfile))
|
||||
|
||||
tqdm.tqdm.write(f'Found {num_meshes} meshes out of a total of {num_params} with {num_completed} completed simulations.')
|
||||
tqdm.tqdm.write(f'Processing the remaining {num_meshes-num_completed} simulations on {num_jobs} workers in parallel.')
|
||||
|
||||
tq = tqdm.tqdm(total=num_meshes-num_completed)
|
||||
threads = []
|
||||
for i in range(num_jobs):
|
||||
t = threading.Thread(target=queue_worker, daemon=True)
|
||||
t.start()
|
||||
threads.append(t)
|
||||
q.join()
|
||||
|
||||
@run.command()
|
||||
@click.pass_context
|
||||
def self_capacitance(ctx):
|
||||
db = ctx.obj['db_connect']()
|
||||
for coil_id, in tqdm.tqdm(db.execute('SELECT coil_id FROM coils WHERE run_id=?', (ctx.obj['run_id'],)).fetchall()):
|
||||
mesh_file = get_mesh_file(ctx.obj['db_connect'](), ctx.obj['mesh_dir'], ctx.obj['run_id'], coil_id, 'normal')
|
||||
if mesh_file.is_file():
|
||||
logfile = ctx.obj['log_dir'] / (mesh_file.stem + '_elmer_self_capacitance.log')
|
||||
with tempfile.TemporaryDirectory() as tmpdir:
|
||||
try:
|
||||
res = subprocess.run(['python', '-m', 'coil_parasitics', 'self-capacitance', '--sim-dir', tmpdir, mesh_file], check=True, capture_output=True)
|
||||
logfile.write_text(res.stdout+res.stderr)
|
||||
except subprocess.CalledProcessError as e:
|
||||
print(f'Error running simulation, rc={e.returncode}')
|
||||
logfile.write_text(e.stdout+e.stderr)
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
cli()
|
||||
|
||||
|
|
@ -1,295 +0,0 @@
|
|||
#!/usr/bin/env python3
|
||||
|
||||
import subprocess
|
||||
import sys
|
||||
import os
|
||||
from math import *
|
||||
from pathlib import Path
|
||||
from itertools import cycle
|
||||
from scipy.constants import mu_0
|
||||
|
||||
from gerbonara.cad.kicad import pcb as kicad_pcb
|
||||
from gerbonara.cad.kicad import footprints as kicad_fp
|
||||
from gerbonara.cad.kicad import graphical_primitives as kicad_gr
|
||||
from gerbonara.cad.kicad import primitives as kicad_pr
|
||||
import click
|
||||
|
||||
|
||||
__version__ = '1.0.0'
|
||||
|
||||
|
||||
def point_line_distance(p, l1, l2):
|
||||
x0, y0 = p
|
||||
x1, y1 = l1
|
||||
x2, y2 = l2
|
||||
# https://en.wikipedia.org/wiki/Distance_from_a_point_to_a_line
|
||||
return abs((x2-x1)*(y1-y0) - (x1-x0)*(y2-y1)) / sqrt((x2-x1)**2 + (y2-y1)**2)
|
||||
|
||||
def line_line_intersection(l1, l2):
|
||||
p1, p2 = l1
|
||||
p3, p4 = l2
|
||||
x1, y1 = p1
|
||||
x2, y2 = p2
|
||||
x3, y3 = p3
|
||||
x4, y4 = p4
|
||||
|
||||
# https://en.wikipedia.org/wiki/Line%E2%80%93line_intersection
|
||||
px = ((x1*y2-y1*x2)*(x3-x4)-(x1-x2)*(x3*y4-y3*x4))/((x1-x2)*(y3-y4)-(y1-y2)*(x3-x4))
|
||||
py = ((x1*y2-y1*x2)*(y3-y4)-(y1-y2)*(x3*y4-y3*x4))/((x1-x2)*(y3-y4)-(y1-y2)*(x3-x4))
|
||||
return px, py
|
||||
|
||||
def angle_between_vectors(va, vb):
|
||||
angle = atan2(vb[1], vb[0]) - atan2(va[1], va[0])
|
||||
if angle < 0:
|
||||
angle += 2*pi
|
||||
return angle
|
||||
|
||||
class SVGPath:
|
||||
def __init__(self, **attrs):
|
||||
self.d = ''
|
||||
self.attrs = attrs
|
||||
|
||||
def line(self, x, y):
|
||||
self.d += f'L {x} {y} '
|
||||
|
||||
def move(self, x, y):
|
||||
self.d += f'M {x} {y} '
|
||||
|
||||
def arc(self, x, y, r, large, sweep):
|
||||
self.d += f'A {r} {r} 0 {int(large)} {int(sweep)} {x} {y} '
|
||||
|
||||
def close(self):
|
||||
self.d += 'Z '
|
||||
|
||||
def __str__(self):
|
||||
attrs = ' '.join(f'{key.replace("_", "-")}="{value}"' for key, value in self.attrs.items())
|
||||
return f'<path {attrs} d="{self.d.rstrip()}"/>'
|
||||
|
||||
class SVGCircle:
|
||||
def __init__(self, r, cx, cy, **attrs):
|
||||
self.r = r
|
||||
self.cx, self.cy = cx, cy
|
||||
self.attrs = attrs
|
||||
|
||||
def __str__(self):
|
||||
attrs = ' '.join(f'{key.replace("_", "-")}="{value}"' for key, value in self.attrs.items())
|
||||
return f'<circle {attrs} r="{self.r}" cx="{self.cx}" cy="{self.cy}"/>'
|
||||
|
||||
def svg_file(fn, stuff, vbw, vbh, vbx=0, vby=0):
|
||||
with open(fn, 'w') as f:
|
||||
f.write('<?xml version="1.0" standalone="no"?>\n')
|
||||
f.write('<!DOCTYPE svg PUBLIC "-//W3C//DTD SVG 1.1//EN" "http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd">\n')
|
||||
f.write(f'<svg version="1.1" width="{vbw}mm" height="{vbh}mm" viewBox="{vbx} {vby} {vbw} {vbh}" xmlns="http://www.w3.org/2000/svg" xmlns:xlink="http://www.w3.org/1999/xlink">>\n')
|
||||
|
||||
for foo in stuff:
|
||||
f.write(str(foo))
|
||||
|
||||
f.write('</svg>\n')
|
||||
|
||||
@click.command()
|
||||
@click.argument('outfile', required=False, type=click.Path(writable=True, dir_okay=False, path_type=Path))
|
||||
@click.option('--footprint-name', help="Name for the generated footprint. Default: Output file name sans extension.")
|
||||
@click.option('--target-layer', default='F.Cu', help="Target KiCad layer for the generated footprint. Default: F.Cu.")
|
||||
@click.option('--jumper-layer', default='B.Cu', help="KiCad layer for jumper connections. Default: B.Cu.")
|
||||
@click.option('--turns', type=int, default=5, help='Number of turns')
|
||||
@click.option('--diameter', type=float, default=50, help='Outer diameter [mm]')
|
||||
@click.option('--trace-width', type=float, default=0.15)
|
||||
@click.option('--via-diameter', type=float, default=0.6)
|
||||
@click.option('--via-drill', type=float, default=0.3)
|
||||
@click.option('--keepout-zone/--no-keepout-zone', default=True, help='Add a keepout are to the footprint (default: yes)')
|
||||
@click.option('--keepout-margin', type=float, default=5, help='Margin between outside of coil and keepout area (mm, default: 5)')
|
||||
@click.option('--twist-width', type=float, default=20, help='Width of twist versus straight coil in percent (0-100, default: 20)')
|
||||
@click.option('--num-twists', type=int, default=1, help='Number of twists per revolution (default: 1)')
|
||||
@click.option('--clearance', type=float, default=0.15)
|
||||
@click.option('--clipboard/--no-clipboard', help='Use clipboard integration (requires wl-clipboard)')
|
||||
@click.option('--counter-clockwise/--clockwise', help='Direction of generated spiral. Default: clockwise when wound from the inside.')
|
||||
def generate(outfile, turns, diameter, via_diameter, via_drill, trace_width, clearance, footprint_name, target_layer,
|
||||
jumper_layer, twist_width, num_twists, clipboard, counter_clockwise, keepout_zone, keepout_margin):
|
||||
if 'WAYLAND_DISPLAY' in os.environ:
|
||||
copy, paste, cliputil = ['wl-copy'], ['wl-paste'], 'xclip'
|
||||
else:
|
||||
copy, paste, cliputil = ['xclip', '-i', '-sel', 'clipboard'], ['xclip', '-o', '-sel' 'clipboard'], 'wl-clipboard'
|
||||
|
||||
|
||||
out_path = SVGPath(fill='none', stroke='black', stroke_width=trace_width, stroke_linejoin='round', stroke_linecap='round')
|
||||
jumper_path = SVGPath(fill='none', stroke='gray', stroke_width=trace_width, stroke_linejoin='round', stroke_linecap='round')
|
||||
svg_stuff = [jumper_path, out_path]
|
||||
|
||||
pitch = clearance + trace_width
|
||||
twist_angle = 2*pi / (turns * num_twists - 1)
|
||||
twist_width = twist_angle * twist_width/100
|
||||
|
||||
via_diameter = max(trace_width, via_diameter)
|
||||
|
||||
# See https://coil32.net/pcb-coil.html for details
|
||||
|
||||
d_inside = diameter - 2*(pitch*turns - clearance)
|
||||
d_avg = (diameter + d_inside)/2
|
||||
phi = (diameter - d_inside) / (diameter + d_inside)
|
||||
c1, c2, c3, c4 = 1.00, 2.46, 0.00, 0.20
|
||||
L = mu_0 * turns**2 * d_avg*1e3 * c1 / 2 * (log(c2/phi) + c3*phi + c4*phi**2)
|
||||
print(f'Outer diameter: {diameter:g} mm')
|
||||
print(f'Average diameter: {d_avg:g} mm')
|
||||
print(f'Inner diameter: {d_inside:g} mm')
|
||||
print(f'Fill factor: {phi:g}')
|
||||
print(f'Approximate inductance: {L:g} µH')
|
||||
|
||||
|
||||
make_pad = lambda num, x, y: kicad_fp.Pad(
|
||||
number=str(num),
|
||||
type=kicad_fp.Atom.smd,
|
||||
shape=kicad_fp.Atom.circle,
|
||||
at=kicad_fp.AtPos(x=x, y=y),
|
||||
size=kicad_fp.XYCoord(x=trace_width, y=trace_width),
|
||||
layers=[target_layer],
|
||||
clearance=clearance,
|
||||
zone_connect=0)
|
||||
|
||||
make_line = lambda x1, y1, x2, y2, layer=target_layer: kicad_fp.Line(
|
||||
start=kicad_fp.XYCoord(x=x1, y=y1),
|
||||
end=kicad_fp.XYCoord(x=x2, y=y2),
|
||||
layer=layer,
|
||||
stroke=kicad_fp.Stroke(width=trace_width))
|
||||
|
||||
make_arc = lambda x1, y1, x2, y2, xc, yc, layer=target_layer: kicad_fp.Arc(
|
||||
start=kicad_fp.XYCoord(x=x1, y=y1),
|
||||
mid=kicad_fp.XYCoord(x=xc, y=yc),
|
||||
end=kicad_fp.XYCoord(x=x2, y=y2),
|
||||
layer=layer,
|
||||
stroke=kicad_fp.Stroke(width=trace_width))
|
||||
|
||||
|
||||
make_via = lambda x, y: kicad_fp.Pad(number="NC",
|
||||
type=kicad_fp.Atom.thru_hole,
|
||||
shape=kicad_fp.Atom.circle,
|
||||
at=kicad_fp.AtPos(x=x, y=y),
|
||||
size=kicad_fp.XYCoord(x=via_diameter, y=via_diameter),
|
||||
drill=kicad_fp.Drill(diameter=via_drill),
|
||||
layers=[target_layer, jumper_layer],
|
||||
clearance=clearance,
|
||||
zone_connect=0)
|
||||
|
||||
pads = []
|
||||
lines = []
|
||||
arcs = []
|
||||
|
||||
for n in range(turns * num_twists - 1):
|
||||
for k in range(turns):
|
||||
r = diameter/2 - trace_width/2 - k*pitch
|
||||
a1 = n*twist_angle + twist_width/2
|
||||
a2 = a1 + twist_angle - twist_width
|
||||
x1, y1 = r*cos(a1), r*sin(a1)
|
||||
out_path.move(x1, y1)
|
||||
x2, y2 = r*cos(a2), r*sin(a2)
|
||||
out_path.line(x2, y2)
|
||||
a3 = (a1 + a2) / 2
|
||||
xm, ym = r*cos(a3), r*sin(a3)
|
||||
arcs.append(make_arc(x2, y2, x1, y1, xm, ym))
|
||||
|
||||
for k in range(turns-1):
|
||||
r1 = diameter/2 - trace_width/2 - (k+1)*pitch
|
||||
r2 = diameter/2 - trace_width/2 - k*pitch
|
||||
a1 = n*twist_angle - twist_width/2
|
||||
a2 = a1 + twist_width
|
||||
x1, y1 = r1*cos(a1), r1*sin(a1)
|
||||
out_path.move(x1, y1)
|
||||
x2, y2 = r2*cos(a2), r2*sin(a2)
|
||||
out_path.line(x2, y2)
|
||||
a3 = (a1 + a2) / 2
|
||||
r3 = (r1 + r2) / 2
|
||||
xm, ym = r3*cos(a3), r3*sin(a3)
|
||||
arcs.append(make_arc(x2, y2, x1, y1, xm, ym))
|
||||
|
||||
rs = diameter/2 - trace_width/2
|
||||
rv = rs - trace_width/2 + via_diameter/2
|
||||
a = n*twist_angle - twist_width/2
|
||||
|
||||
x1, y1 = rs*cos(a), rs*sin(a)
|
||||
out_path.move(x1, y1)
|
||||
xv1, yv1 = rv*cos(a), rv*sin(a)
|
||||
out_path.line(xv1, yv1)
|
||||
svg_stuff.append(SVGCircle(via_diameter/2, xv1, yv1, fill='red'))
|
||||
pads.append(make_via(xv1, yv1))
|
||||
jumper_path.move(xv1, yv1)
|
||||
lines.append(make_line(x1, y1, xv1, yv1))
|
||||
|
||||
a += twist_width
|
||||
rs = diameter/2 - trace_width/2 - (turns-1)*pitch
|
||||
rv = rs + trace_width/2 - via_diameter/2
|
||||
|
||||
x1, y1 = rs*cos(a), rs*sin(a)
|
||||
out_path.move(x1, y1)
|
||||
xv2, yv2 = rv*cos(a), rv*sin(a)
|
||||
out_path.line(xv2, yv2)
|
||||
svg_stuff.append(SVGCircle(via_diameter/2, xv2, yv2, fill='red'))
|
||||
pads.append(make_via(xv2, yv2))
|
||||
lines.append(make_line(x1, y1, xv2, yv2))
|
||||
|
||||
if n > 0:
|
||||
jumper_path.line(xv2, yv2)
|
||||
lines.append(make_line(xv1, yv1, xv2, yv2, jumper_layer))
|
||||
else:
|
||||
pads.append(make_pad(1, xv1, yv1))
|
||||
pads.append(make_pad(2, xv2, yv2))
|
||||
|
||||
svg_file('/tmp/test.svg', svg_stuff, 100, 100, -50, -50)
|
||||
|
||||
if counter_clockwise:
|
||||
for p in pads:
|
||||
p.at.y = -p.at.y
|
||||
|
||||
for l in lines:
|
||||
l.start.y = -l.start.y
|
||||
l.end.y = -l.end.y
|
||||
|
||||
for a in arcs:
|
||||
a.start.y = -a.start.y
|
||||
a.end.y = -a.end.y
|
||||
|
||||
if footprint_name:
|
||||
name = footprint_name
|
||||
elif outfile:
|
||||
name = outfile.stem,
|
||||
else:
|
||||
name = 'generated_coil'
|
||||
|
||||
if keepout_zone:
|
||||
r = diameter/2 + keepout_margin
|
||||
tol = 0.05 # mm
|
||||
n = ceil(pi / acos(1 - tol/r))
|
||||
pts = [(r*cos(a*2*pi/n), r*sin(a*2*pi/n)) for a in range(n)]
|
||||
zones = [kicad_pr.Zone(layers=['*.Cu'],
|
||||
hatch=kicad_pr.Hatch(),
|
||||
filled_areas_thickness=False,
|
||||
keepout=kicad_pr.ZoneKeepout(copperpour_allowed=False),
|
||||
polygon=kicad_pr.ZonePolygon(pts=kicad_pr.PointList(xy=[kicad_pr.XYCoord(x=x, y=y) for x, y in pts])))]
|
||||
else:
|
||||
zones = []
|
||||
|
||||
fp = kicad_fp.Footprint(
|
||||
name=name,
|
||||
generator=kicad_fp.Atom('GerbonaraTwistedCoilGenV1'),
|
||||
layer='F.Cu',
|
||||
descr=f"{turns} turn {diameter:.2f} mm diameter twisted coil footprint, inductance approximately {L:.6f} µH. Generated by gerbonara'c Twisted Coil generator, version {__version__}.",
|
||||
clearance=clearance,
|
||||
zone_connect=0,
|
||||
lines=lines,
|
||||
arcs=arcs,
|
||||
pads=pads,
|
||||
zones=zones,
|
||||
)
|
||||
|
||||
if clipboard:
|
||||
try:
|
||||
print(f'Running {copy[0]}.')
|
||||
proc = subprocess.Popen(copy, stdin=subprocess.PIPE, text=True)
|
||||
proc.communicate(fp.serialize())
|
||||
except FileNotFoundError:
|
||||
print(f'Error: --clipboard requires the {copy[0]} and {paste[0]} utilities from {cliputil} to be installed.', file=sys.stderr)
|
||||
elif not outfile:
|
||||
print(fp.serialize())
|
||||
else:
|
||||
fp.write(outfile)
|
||||
|
||||
if __name__ == '__main__':
|
||||
generate()
|
||||
Plik diff jest za duży
Load Diff
Ładowanie…
Reference in New Issue