Pyems is a Python interface to the electromagnetic field solver, OpenEMS. It uses OpenEMS’s own Python interface to hook into CSXCAD and OpenEMS. However, unlike that interface, whose primary purpose is to expose the underlying C++ interface as a Python API, pyems provides high-level functionality intended to facilitate and accelerate the simulation process. Moreover, OpenEMS contains a number of subtle usage traps that can be confusing to new users, and requires that those users understand certain limitations of the FDTD algorithm. Pyems attempts (when possible) to enforce correct usage in these situations and in other cases to make the usage tradeoffs more visible.
To accomplish this goal, pyems provides:
- Configurable, automatic mesh generation.
- KiCad footprint creation.
- Simple port impedance and S-parameter calculations.
- A collection of cooperative classes for building commonly-needed microwave structures (many of which are PCB-based).
- Functions performing frequently-needed calculations involved with microwave design.
- Methods to optimize a simulation structure to achieve a desired result for any arbitrary parameter.
- A simple, expressive interface intended to make the simulation process more intuitive.
Although ease of use is one of pyems’s primary goals, another one of its major goals is to not restrict the power of OpenEMS. To this effect, the underlying OpenEMS Python interface is still directly accessible to the user. Indeed, there are cases in which it is necessary to use the OpenEMS interface (e.g. applying transformations and constructing some simple shapes). Additionally, there are some cases where the interface has been designed to allow feature accessibility in a way that can be misused by the casual user. Pyems will always prioritize expressivity over protection against misuse. These cases should be properly documented.
Pypi installation is planned but not yet available. Pyems must currently be installed manually. It requires a working installation of OpenEMS, including its Python interfaces for CSXCAD and OpenEMS. Pyems also requires the NumPy and SciPy Python libraries. If you’d like to view simulation field dumps (recommended) you must also have a copy of ParaView. All of this software is (of course) open source and free to download and use.
I have not yet written a proper reference documentation for pyems. To help new users get acquainted with it, I’ve written a tutorial presented in this section. This should help users understand the API and what can be done with pyems. Additionally, there are numerous examples present in the examples directory that can be referenced. Finally, I’ve made a significant effort to thoroughly document class initialization and function parameters in docstrings. Please reference them and file bugs when ambiguities exist.
This tutorial will describe how to simulate a directional coupler fabricated on the top layer of a PCB fabricated with OSHPark’s 4-layer process. The simulation file is available in the examples directory. The user will be able to view the simulation structure as well as a movie of the current density. The post-simulation analysis will include the calculated scattering parameters for all 4 ports. Finally, the simulation will generate a PCB footprint that can be imported directly into KiCAD.
import numpy as np
from pyems.structure import Microstrip, PCB, MicrostripCoupler
from pyems.simulation import Simulation
from pyems.pcb import common_pcbs
from pyems.calc import phase_shift_length, microstrip_effective_dielectric
from pyems.utilities import print_table
from pyems.coordinate import Coordinate2, Axis, Box3, Coordinate3
from pyems.mesh import Mesh
from pyems.field_dump import FieldDump, DumpType
from pyems.kicad import write_footprint
freq = np.linspace(0e9, 18e9, 501)
ref_freq = 5.6e9
unit = 1e-3
sim = Simulation(freq=freq, unit=unit, reference_frequency=ref_freq)This first code block instantiates a Simulation()
object. Simulation initializes the underlying OpenEMS objects for
the simulation structure and the FDTD algorithm. freq specifies the
frequency values to simulate. This simulation uses a granularity of
501 values from DC to 18GHz. The post-simulation analysis results
(such as the S-parameters) will be given for these frequency
values. Bear in mind that using a smaller granularity results in a
longer simulation time. unit specifies the length unit to use in
meters. By using 1e-3, we’ve told pyems that all of our subsequent
dimensions will be given in mm.
Certain dielectric properties are frequency-dependent. However,
OpenEMS requires that all dielectric properties be constant in a
simulation. The reference_frequency parameter specifies the
frequency used to determine these values. It can be ommitted in which
case the center frequency of freq is used.
Simulation also takes a number of other parameters. They have been
given reasonable defaults but will sometimes need to be changed. See
the code documentation for these parameters. In particular, it will
frequently be necessary to specify the boundary_conditions
parameter.
pcb_prop = common_pcbs["oshpark4"]
trace_width = 0.38
trace_spacing = 0.2
eeff = microstrip_effective_dielectric(
pcb_prop.substrate.epsr_at_freq(ref_freq),
substrate_height=pcb_prop.substrate_thickness(0, unit=unit),
trace_width=trace_width,
)
coupler_length = phase_shift_length(90, eeff, ref_freq)
print("coupler_length: {:.4f}".format(coupler_length))pcb_prop is an object that contains details about the OSHPark
4-layer process. It knows the thickness of all metal and dielectric
layers as well as the dielectric frequency-dependent electrical
properties. Only a few PCB processes are supported at the moment, but
more will be added in the future.
eeff is the effective dielectric of the top PCB layer. It correctly
accounts for the fact that the microstrip is bounded below by the
substrate and above by air.
coupler_length is the length (in our chosen unit, which is mm)
required for a signal (specified by the reference frequency) to
undergo a quarter-wavelength phase shift. Since this coupler is a
backward-wave directional coupler, the quarter wave maximizes the
coupling coefficient and bandwidth at our reference frequency.
The effective dielectric equation (and by virtue the coupler length) is approximate, not based on a proper simulation. Although the approximation should be more than adequate for most cases, we could optimize the length later (and calculate a more precise effective dielectric) with the OpenEMS simulation if we wanted.
pcb_len = 2 * coupler_length
pcb_width = 0.5 * pcb_len
pcb = PCB(
sim=sim,
pcb_prop=pcb_prop,
length=pcb_len,
width=pcb_width,
layers=range(3),
omit_copper=[0],
)PCB creates a PCB object as part of the simulation structure. PCB
is our first example of what pyems refers to as a structure, which is
a collection of primitives (the OpenEMS terminology for simple shapes
with associated electrical properties) and other pyems structures that
present a useful abstraction as a single object. In practice,
structures allow you to quickly instantiate frequently-needed physical
objects while using OpenEMS best-practices. They also make it easy to
apply transformations (physical rotations and translations) to these
objects.
Structures play well together. For instance, there is a via structure which requires an associated PCB structure. Instead of having to worry about the 3-dimensional position and orientation of the via, you can simply specify its 2-dimensional coordinates on the PCB. The via will then be automatically oriented correctly on the PCB.
The via also serves to illustrate the benefits of structures over the
underlying OpenEMS primitives. Instead of having to instantiate a
cylinder for the via drill, another cylinder or cylindrical shell for
the via plating and then flat cylindrical shells for the each of the
pads and antipads, we can simply instantiate a Via object with the
desired attributes. Pyems fully supports blind and buried vias too, as
well as physically-inaccurate approximations of vias that shorten
simulation time.
Let’s return to the PCB object we instantiated above. This is a core
structure of many simulations, since many simulations instantiate
microwave structures on a PCB. We must tell the PCB object what
process we are using (so that it can automatically determine certain
dimensional and electrical properties) as well as the simulation
object we instantiated at the beginning. Additionally, we must specify
the x-dimensional length and y-dimensional width of the PCB. Although
our PCB process is a 4-layer process, by building a microstrip
directional coupler, we really only care about the first and second
metal layers and the substrate layer in-between. This is what the
layers parameter does. range(3) specifies that we only want to
include layers 0, 1, and 2, where 0 and 2 correspond to the first and
second metal layers and 1 corresponds to the top substrate layer. This
is an important feature since it leads to shorter simulation times
with virtually zero accuracy cost. By default all layers are
included. Pyems does not presently support layers other than
dielectric and metal layers (such as soldermask or silkscreen
layers). These may be added later if desired.
Finally, PCB by default fills all metal layers with a copper
pour. This is often useful and obviates the need for the user to do
this manually. We can use the omit_copper parameter to specify metal
layers where all the metal should be etched away. Although the
layers and omit_copper parameters may seem similar, there are a
few subtle differences. Firstly, layers requires a Python range
object wherease omit_copper requires a list. While it is reasonable
for us to include/disclude a copper pour on any metal layer, it
doesn’t make sense for us to use construct our PCB from the first and
second metal layers and the second substrate layer (omitting the first
substrate layer). Secondly, layers considers all layers (metal and
dielectric) when considering indices for the layers. By contrast,
omit_copper only cares about the metal layers and thus ignores
dielectric layers. As a result, the first and second metal layers are
indicated by 0 and 2 when passed to layers and by 0 and 1 when
passed to omit_copper.
coupler = MicrostripCoupler(
pcb=pcb,
position=Coordinate2(0, 0),
trace_layer=0,
gnd_layer=1,
trace_width=trace_width,
trace_gap=trace_spacing,
length=coupler_length,
miter=None,
)MicrostripCoupler instantiates coupled microstrip lines. It is
another example of a pyems structure. It acquires information about
the PCB object and simulation via the pcb parameter, since
microstrip couplers will always be instantiated on a PCB. position
specifies its center position. trace_layer and gnd_layer specify
the PCB metal layers of the trace and backing ground
plane. trace_width is the width of each microstrip and trace_gap
is the perpendicular distance between the inside of each
trace. length is the x-dimensional length, which we set to the
desired coupler length. The last parameter, miter specifies the
amount to miter the corners of ports 3 and 4. By specifying None
we’ve chosen an approximate, optimial miter (see the Miter structure
for more information). The use of miter here may be changed in the
future for something more general, since it is conceivable that a user
might not want to miter these corners, or do something else to them
like rounding. It is worth mentioning that MicrostripCoupler also
takes a transform parameter that we could use to rotate it.
coupler_port_positions = coupler.port_positions()
port0_x = coupler_port_positions[0].x
port0_y = coupler_port_positions[0].y
Microstrip(
pcb=pcb,
position=Coordinate2(np.average([port0_x, -pcb_len / 2]), port0_y),
length=port0_x + pcb_len / 2,
width=trace_width,
propagation_axis=Axis("x"),
port_number=1,
excite=True,
ref_impedance=50,
feed_shift=0.3,
)Microstrip creates a microstrip port. Microstrip is another
structure, but it is also an example of another important concept in
pyems: a port. Ports are conceptually identical to the OpenEMS concept
(and there is a significant degree of overlap in the implementation)
except that they integrate better with the rest of pyems. A port is
essentially a point of interface to the outside world. Ports are
locations where signal excitations are created and where voltages and
currents are measured.
The notion of ports used here is analogous to the notion of ports used by a VNA. For instance (although it is not the case in this simulation) we might have added SMA connectors at each port (pyems provides a structure for this too). Then, if we wanted to measure S₂₁ we’d terminate ports 3 and 4 with matching loads, attach the transmission port of the VNA to port 1 via an SMA cable and the other port of the VNA (assuming a 2-port VNA) to port 2. If the VNA is properly calibrated for the SMA cables, it will measure the signal as “starting” at the SMA connector of port 1 and “ending” at the SMA connector of port 2. Pyems will do exactly the same thing and should yield the same results.
There are a few aspects to the instantiation of Microstrip that
indicate this is used as a port. The first (and most obvious) is
port_number. As should be evident, this tells the simulation that
this microstrip structure acts as port 1. The numbering will be
important in the post-simulation analysis when calculating our
S-parameters. Next, the excite parameter tells the simulation that
we’d like to perform a signal excitation at this port. The excitation
is a Gaussian excitation whose frequency range is determined by the
Simulation freq parameter used at the beginning of this
tutorial. ref_impedance specifies the impedance value to use when
calculating the port’s voltage and current values. We could also have
omitted this parameter in which case the calculated value of the
microstrip’s characteristic impedance would have been used. Typically,
this should be set to the desired characteristic impedance as is done
here. feed_shift specifies the position of the signal excitation
along the port as a fraction of the port length. The feed needs to be
placed far enough along the port such that it is not contained within
a boundary (see the OpenEMS documentation for boundary
conditions). Pyems will notify you if the excitation is placed in a
boundary.
The propagation_axis parameter specifies the direction the port
faces. Because of the way the FDTD rectilinear grid works, we cannot
place the port in any arbitrary orientation. Finally, we can see that
the position and length parameters were used to place the port as
extending from the lowermost x-position of the PCB to the edge of the
MicrostripCoupler structure.
port1_x = coupler_port_positions[1].x
Microstrip(
pcb=pcb,
position=Coordinate2(np.average([port1_x, pcb_len / 2]), port0_y),
length=pcb_len / 2 - port1_x,
width=trace_width,
propagation_axis=Axis("x", direction=-1),
port_number=2,
excite=False,
ref_impedance=50,
)
port2_x = coupler_port_positions[2].x
port2_y = coupler_port_positions[2].y
Microstrip(
pcb=pcb,
position=Coordinate2(port2_x, np.average([port2_y, -pcb_width / 2])),
length=port2_y + pcb_width / 2,
width=trace_width,
propagation_axis=Axis("y"),
ref_impedance=50,
port_number=3,
)
port3_x = coupler_port_positions[3].x
Microstrip(
pcb=pcb,
position=Coordinate2(port3_x, np.average([port2_y, -pcb_width / 2])),
length=port2_y + pcb_width / 2,
width=trace_width,
propagation_axis=Axis("y"),
ref_impedance=50,
port_number=4,
)Ports 2, 3 and 4 are instantiated in much the same way as port 1. There are two main differences, however. The first is that ports 3 and 4 face in the y-direction. This rotates the structure and measurement probes by 90 degrees relative to an x-orientation. The other difference is that port 2 faces in the negative x-direction. This ensures that the voltage and current calculations are performed correctly for its orientation.
Mesh(
sim=sim,
metal_res=1 / 80,
nonmetal_res=1 / 40,
smooth=(1.1, 1.5, 1.5),
min_lines=5,
expand_bounds=((0, 0), (0, 0), (10, 40)),
)At this point we’ve finished the entire physical structure used in the simulation. In other words if we viewed the structure with AppCSXCAD (which we’ll do shortly), it would look like it would if you were holding the PCB in front of you. Additionally, we’ve imbued that structure with all the electrical properties it needs for simulation.
However, OpenEMS’s FDTD algorithm needs to know where in that structure it should be calculating the solutions to Maxwell’s equations at each timestep. This is where the simulation mesh comes in and is, in my opinion, one of the greatest advantages of pyems over OpenEMS’s default Python interface. Traditionally, creating the mesh has been one of the hardest and most cumbersome parts of the OpenEMS simulation process. There are a number of implementation-specific reasons for this. For instance, the FDTD algorithm performs badly when a mesh line is placed at the boundary of a conductor and insulator. Instead, something called the thirds rule should be applied to achieve a more accurate simulation result without simply adding more mesh lines (which would increase the simulation time). Pyems takes care of this and a bunch of other implementation-specific details for you. For instance it ensures a proper smoothness between adjacent mesh line spacings and makes sure that mesh lines work well with voltage and current probes (there are a number of important considerations in this regard that I won’t go into now).
metal_res specifies the maximum spacing between mesh lines inside a
metal. It is specified as a fraction of the minimum simulation
wavelength, which in turn is determined by the maximum frequency of
freq from the beginning of this tutorial. nonmetal_res does the
same thing but for non-metal areas such as the substrate and
surrounding air. smooth ensures that adjacent spacings are within a
multiplicative factor of one another. Each dimension abides by its own
smoothness factor, which is why we pass a tuple of 3 elements
corresponding to (x, y, z). In this example, we’ve kept the x lines
“smoother” than the y or z lines since the signal propagates primarily
in the x-direction. The min_lines parameter specifies the minimum
number of mesh lines that must be present in one dimension of a
primitive. For instance, the width of a microstrip trace (given the
resolution we’ve provided) would normally contain fewer than 5 mesh
lines. However, if there are too few mesh lines the simulation will
give incorrect results, believing that the microstrip structure is a
different width than it actually is. Finally, expand_bounds
specifies the number of additional lines we’d like outside our
simulation structure. This creates an air layer between the structure
and the boundary. The parameter is passes as a tuple of 3 tuples each
of 2 elements. It signifies
((xmin, xmax), (ymin, ymax), (zmin, zmax))
We can see from our example that we’ve only added an air layer in the z-dimension. We haven’t done this in the x-, or y-dimensions because the ports must terminate in a perfectly-matched layer (PML). This ensures that we don’t get signal reflections at the ports, making our post-simulation analysis more accurate.
FieldDump(
sim=sim,
box=Box3(
Coordinate3(-pcb_len / 2, -pcb_width / 2, 0),
Coordinate3(pcb_len / 2, pcb_width / 2, 0),
),
dump_type=DumpType.current_density_time,
)FieldDump adds a non-physical structure to our simulation, which
will record and allow us to view the current density at the top PCB
metal layer. box specifies the region to record. We have made it
2-dimensional though we could have made it 3-dimensional. dump_type
specifies the type of field to record, for which there are a number of
possibilities. See DumpType for other options.
write_footprint(coupler, "coupler_20db", "coupler_20db.kicad_mod")write_footprint writes a KiCAD-compatible footprint relative to the
current directory.
sim.run()Calling the run method of our Simulation object first displays our
CSXCAD object with AppCSXCAD (this can be turned off for usage in
scripts) and then asks us if we’d like to proceed with the OpenEMS
simulation.
At this point you should have an AppCSXCAD window open with the following structure
sim.view_field()view_field() runs ParaView on the recorded field dump. Here’s a GIF
of the result
print_table(
data=[
sim.freq / 1e9,
np.abs(sim.ports[0].impedance()),
sim.s_param(1, 1),
sim.s_param(2, 1),
sim.s_param(3, 1),
sim.s_param(4, 1),
],
col_names=["freq", "z0", "s11", "s21", "s31", "s41"],
prec=[4, 4, 4, 4, 4, 4],
)print_table is a convenience method to print tabular data in
nicely-spaced columns. This displays the calculated port 1 impedance
and all S-parameters for each frequency value of the simulation.
If we had plotted this and additionally computed the directivity, we would see
This section is very incomplete.
Many structure objects accept optional transformation parameters. They also generally accept position coordinates. The object is first created at the origin, then the transform is applied, finally followed by a translation of the center of the structure to the supplied position. As a result translation transformations should not be needed, although pyems will accept them.
This section describes how the automatic mesh generation algorithm works. Although I intend to keep it up to date, since the mesher is still evolving this description may at times lag behind development. If you find an inconsistency, please submit an issue.
In order to generate a mesh from a CSXCAD geometry, pyems starts by
getting a list of all physical CSXCAD primitives (i.e., CSXCAD
primitives that have an effect on the simulation). Then, for each
dimension it extracts a list of locations for mesh lines that must be
fixed at those locations. These correspond to zero-dimension
structures (e.g., a planar structure created by
AddConductingSheet). Next, pyems iterates through the full list of
physical primitives and extracts 3 lists of boundary positions, one
for each dimension. For example, a boundary position in the
x-dimension would correspond to a location where the physical
properties of the structure change anywhere at that location. This
change could occur in the y-dimension or the z-dimension. Boundary
positions in the y-dimension and z-dimension lists are analogous.
Pyems then converts each element of these lists of boundary positions
in each dimension to a type that associates a boundary region
(consisting of lower and upper bounds) with a CSXCAD property type,
which it classifies as metal, nonmetal or air. This is called a
BoundedType. To associate the property type it finds the type of the
primitive corresponding to the smallest length in that dimension
encompassing the bounded region. Where ties exist, the metal type gets
priority. There is still at least one issue with this part of the
process, which I will fix (e.g., see issue #2).
Pyems then adds peripheral BoundedType’s for simulation air space,
records the location of boundaries between a metal, and nonmetal and
orders the BoundedType’s by size (smallest first). Then, it iterates
through this list of BoundedType’s and generates a list of mesh line
locations inside each.
Generating this list of mesh lines in each BoundedType is, of
course, where most of the work happens. Pyems starts by finding the
mesh line spacing at the lower and upper boundaries of the boundary
region. It also determines the maximum spacing in this region
according to metal_res and nonmetal_res specified by the user and
whether this BoundedType is a metal or not. Then it adjusts the
upper and lower line positions if they correspond to metal-nonmetal
boundaries, which must satisfy the thirds rule. For instance, if the
upper boundary position corresponds to a metal-nonmetal boundary, we
would move the upper position 1/3 the mesh spacing inside the
boundary. Then pyems computes a geometric series whose constant factor
is between 1 and the smoothness factor specified by the user for that
dimension and whose distance is equal to the distance of the bounded
region. Computing the geometric series uses a Scipy optimization
routine and accounts for most of the time spent generating the mesh.
Finally, pyems trims the air mesh to the desired number of cells, smooths the mesh so that the mesh spacing inside the PML is uniform, calls hooks to other parts of pyems that need to know the final mesh location (e.g., probes need to align to the mesh), and generates the actual mesh in the CSXCAD structure. It then performs a number of checks for correctness.
The following set of features is planned, but not currently implemented.
- A tolerance analysis that incorporates variation in the input simulation parameters (e.g. prepreg thickness, etching precision, etc.).
- Support for independent dielectric properties for each substrate layer. Many PCB processes (especially in microwave contexts) require this. This is not difficult to implement. Please raise an issue if you’d like this.
A number of equations in this code base come from microwave design and theory textbooks. I’ve made an effort to make a comment in the code whenever an equation is taken from one of these textbooks so that users can look up the corresponding theory and to make it easier to find bugs in the code.
Here’s a list of the textbooks referenced:
- Pozar refers to “Microwave Engineering” by David Pozar, Fourth Edition.
- Wadell refers to “Transmission Line Design Handbook” by Brian Wadell, published 1991.
If you find a reference to a text not mentioned here, please submit a bug report or pull request.
The via wall otherwise often gets ignored. I believe this is a result of the floating point precision errors.
These methods are poorly named. freq_data and time_data are better names. Additionally, they shouldn’t pass back frequency and time values. This should be retreived with other methods. Note that this will require adjustments to port.py too.
This should use the Axis object.
self._propagation_axis is not currently required for the port base class. The interface must be changed in some way that is also compatible with the derived classes.