Table of Contents
About Pyems is a Python interface to the electromagnetic field solver, [[https://openems.de/start/index.php][OpenEMS]]. It uses OpenEMS's own Python interface to hook into [[https://github.com/thliebig/CSXCAD][CSXCAD]] and [[https://github.com/thliebig/openEMS][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:
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.
Installation 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 [[https://numpy.org/][NumPy]] and [[https://scipy.org/][SciPy]] Python libraries. If you'd like to view simulation field dumps (recommended) you must also have a copy of [[https://www.paraview.org/][ParaView]]. All of this software is (of course) open source and free to download and use.
Tutorial 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 [[https://github.com/matthuszagh/pyems/tree/master/examples][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 [[https://en.wikipedia.org/wiki/Power_dividers_and_directional_couplers#Directional_couplers][directional coupler]] fabricated on the top layer of a PCB fabricated with [[https://docs.oshpark.com/services/four-layer/][OSHPark's 4-layer process]]. The simulation file is available [[https://github.com/matthuszagh/pyems/blob/master/examples/coupler.py][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 [[http://openems.de/index.php/FDTD_Boundary_Conditions][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 [[https://en.wikipedia.org/wiki/Regular_grid][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 [[https://openems.de/index.php/FDTD_Mesh.html][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
[[file:.img/coupler_csxcad.png]]
sim.view_field()
~view_field()~ runs ParaView on the recorded field dump. Here's a GIF of the result
[[file:.img/coupler_current_time.gif]]
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
[[file:.img/coupler_plot.svg]]
Structures :PROPERTIES: :ID: 556e1040-a5bd-4175-8b80-d7d613fea8ba :END: Transformations Many [[id:556e1040-a5bd-4175-8b80-d7d613fea8ba][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.
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 [[https://github.com/matthuszagh/pyems/issues/2][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.
Here's a list of the textbooks referenced:
If you find a reference to a text not mentioned here, please submit a bug report or pull request.
TODO probe should not hold onto freq TODO probe get_freq_data and get_time_data 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.
** TODO rectwaveguideport propagation axis This should use the Axis object.
** TODO port calc requires self._propagation_axis set 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.
** HOLD mesh should support primitive priorities