Getting started with MEEP - simple electric field simulations in python

by Matthew Baas

A quick intro to simulating the electric field for simple current distributions in python using MEEP.

TL;DR: While looking for a nice open source piece of electromagnetics simulation software, particularly ones which could be scripted using python, I found MEEP. It seems to work pretty well, and has quite a reasonable python interface, and allows for a wide range of kinds of electromagnetic simulations. However, the documentation and introductions appear to be more geared toward people who know a lot more about it than me, and seems to be more focused on photonics than simulation of basic things like electric/magnetic field intensity. So here I present a quick intro for those who have only taken one or two undergraduate modules in electromagnetics to get started with using MEEP.

1. Installation

TL;DR: use the docs. Bad news: MEEP does not support windows as of 2019-12-20. You can get around this hurdle on windows if you use the Windows Subsystem for Linux, but it is a bit buggy. Hopefully they add support soon.

Quickstart on Google Colab

If you want to get started using google colab notebooks (which are run on linux backends), then you need to do some special path stuffs to get it to work nicely. In short, to get MEEP working 100%, at the start of the notebook run:

!wget -c https://repo.continuum.io/miniconda/Miniconda3-latest-Linux-x86_64.sh
!chmod +x Miniconda3-latest-Linux-x86_64.sh
!bash ./Miniconda3-latest-Linux-x86_64.sh -b -u -p /usr/local 
import os
if os.path.isdir('/content/anaconda'): root_path = '/content/anaconda'
else: root_path = '/usr/local/'
os.environ['PATH'] += f":{root_path}/bin"
!conda create -n mp -c conda-forge pymeep python=3.7 -y
print(">> ", root_path)
import sys
sys.path.append(f'{root_path}envs/mp/lib/python3.7/site-packages/')

After that, you can just run import meep as mp and things will work as expected. Remember to confirm yes when prompted while creating the new conda environment.

Note: this is quite hacky and quite certainly not good code, but it works well enough to just get started with MEEP in colab.

Update 2021-02-19: I updated the quickstart code above after the miniconda installer changed how it handles args for installation directory.

Update 2021-08-19: Changed install method again after some changes in how meep and miniconda interact on colab. It should work on colab again.

Update 2023-06-26: Google colab has changed how it handles conda and dependencies in the last two years, and MEEP has evolved quite a bit from when this post was made. Depending on your setup and colab machine, this meep installation method might no longer work.

2. MEEP basics

MEEP can solve Maxwell’s equations using either in the time-domain using finite-difference time-domain (FDTD) or by using a frequency-domain solver. In some cases the one is faster than the other. Often for simple geometries and currents, the frequency domain solver is a bit quicker.

2.1 Units in MEEP

Basically, the official documentation on units and some of the great online resources (like this fan wiki page or this notebook) don’t have sufficient examples in python, are not in english (and most translators don’t work well with highly niche technical documents), or expect me to know a bit more than I actually do. So I’ll give a more Baas style explanation and notation here, coming from an engineering perspective using the SI system.

Defining base quantities

TL;DR: in electromagnetics, the 4 SI base quantities that are important are time, mass, length, and current. Just like how the SI provides base definitions for “1 unit of length is a meter”, or “1 unit of current is an ampere,” MEEP redefines these base quantities to construct the speed of light in vacuum, \(c\) to be equal to 1 in the MEEP base units chosen.

Length

First, a base unit od length is chosen. To use the same notation as the documentation, let’s call this base unit \(a\). This \(a\) can be anything you like, and then all distances are measured in units of ‘a’. For example, in the SI system, \(a = 1\text{m}\), so specifying a length as a distance of 13 in this unit system corresponds to a distance of \(13a = 13\) meters. Similarly, if we let \(a = 0.5\text{μm}\), then specifying some length variable as 13 will indicate a distance in the SI system as \(13a = 6.5\text{μm}\). Or to convert the other way, if we want to convert 7μm into the chosen \(a = 0.5\text{μm}\) MEEP units, then we just note that

Which is the dimensionless quantity of length we would put as the value in the python program.

Time

Next, by defining \(c = 1\), where we know that in the SI units, \(c\) is measured in m/s, using the reasoning from the above paragraph to equate the two values, we must find the base unit for time \(b\) to satisfy:

But we know that a ratio of \(a/\)(some measurement in meters) is dimensionless, since \(a\) is just equal to some multiple of meters (like 0.94μm for example). Namely, it is just the numerical (unitless) value of \(a\), which I will call \(|a|\) Thus:

In other words, 1 unit of time \(b\) in MEEP units is the time it takes for light (in a vacuum) to travel \(1\) units of distance in MEEP units (which is \(1\times a\) in meters). So for example, if \(a = 0.5\text{μm}\) as before, and we want to know how long 6 MEEP units of time was in seconds, then we would go:

Similarly, to convert 3 ps (picoseconds) to MEEP units using the same \(a = 0.5\text{μm}\), then we would simply go:

So to represent a time of 3ps if your base distance is 0.5μm, then you would enter 1798.8 as the time variable in the python code.

Velocity and frequency

Another corollary of this is that since the speed of light \(c\) is normalized to be 1, it means that the units of speed are all expressed as multiples of the speed of light (you can check the ratio of \(a/b\) to confirm this) – that is, a velocity of 0.2 in MEEP units is 0.2 multiplied by \(2.998\times 10^8\) in m/s. Also, since the unit of frequency is just the reciprocal of time, in MEEP units, the frequency is expressed in terms of multiples of a base unit \(1/b\). Example: to convert 30THz to MEEP units when \(a = 210\text{μm}\),

And with the same length base, to convert a frequency of 0.9 in MEEP units back to SI Hertz, then simply \(0.9\times \frac{2.998\times 10^8}{ | a | } = 1.28\ \text{THz}\) .

Current and mass

These last two base quantities need to be defined. However, they are actually related by the setting of \(c=1\), so we are free to chose one and the other is then set. TL;DR: pick the base current \(I_\text{base}\) as whatever you like (in SI units, this base is 1 Ampere). And then the base unit for mass, \(m_\text{base}\) is:

Where \(\epsilon_0, c\) are the vacuum permittivity and speed of light in SI units, respectively. There is some nice math showing this here.

Base units of field quantities

Now that we have defined the base units, all field quantities of the electric/magnetic field and power stuffs can be defined in terms of their base quantities, which can then easily be converted to MEEP base values and back to SI base values.

The MEEP fandom wiki did some nice derivations of the scaling factors for a few of the basic field quantities (link here). For the table below, constants like \(\epsilon_0, c\) are their values in the base SI unit system, and the MEEP base quantity is the factor by which to multiply a dimensionless MEEP value to convert it to the SI unit system. And:

• \(a\) is the base unit of length value you choose in MEEP. (e.g 5m)
• \(I_\text{base}\) is the base current you choose in MEEP. You can kind of keep this as 1A if you are not comfortable dealing with it.
• This table is adapted from this link, thanks to the MEEP fandom wiki!

Transformation factors from MEEP to SI units of basic field quantities

• To convert from MEEP to SI: multiply by the MEEP base unit (expressed in SI units).
• To convert from SI to MEEP: divide by the MEEP base unit (expressed in SI units).

And on a another note: for complex \(\epsilon\) with lossy materials, where the permittivity must be expressed as a complex number (or even a Fourier series), MEEP does support it as well, but then the conversions might get slightly more intricate. For the purposes of this post, let’s assume that the most ugly materials we deal with can be represented by just a complex \(\epsilon, \mu\) (see this the docs for more details).

2.2 Building simulations

Before we can simulate field quantities, we need to define the geometry of the simulation. That is, we need to define the material media, sources, and boundaries of what we intend to simulate. The docs give a great in-depth show of this part, so I will give a ‘getting started’ short version:

Decide base unit and config variables

TL;DR:

• Decide on the base unit of length \(a\). This is never explicitly programmed into the simulation code, but keep note of it if you want to convert between SI values and MEEP dimensionless values.
• MEEP calculates the fields at each pixel in a (2D or 3D) grid. First, define the resolution of the grid in terms of pixels/\(a\). That is, how many pixels should MEEP fit into a distance of 1 unit (in MEEP \(a\) units).
• Next, define the 2D or 3D cell/volume in which you are interested in calculating the fields for. MEEP has a 3D vector object Vector3(x, y, z), wherein you can define a 3D region in which to calculate the fields. This definition of the simulation cell, together with the resolution, determine the shape of the output tensor you will get from the simulation if you request a particular field component.

So, let’s simulate a simple linear array of 4 Hertzian dipoles (with z-directed current) aligned along x-axis. Let’s decide to make \(a = 1\text{mm}\), with a resolution of 50 pixels per 1mm. Further, let’s place the dipoles as shown in the image below, with each dipole separated by 1mm. So we have dipoles at x=2mm, 1mm, 0mm, and -1mm. Next let’s make the simulation cell be a 20x20mm block centered around the origin in the H-plane of the dipoles. The situation thus looks like:

The following code sets up the resolution and cell in the H-plane:

import meep as mp
resolution = 50 # 50 pixels per unit a.
cell = mp.Vector3(20, 20, 0) # just work in 2D for this

Add sources

Now we add sources, namely the Hertzian dipoles (which have an infinitesimal dipole length). We set the dipoles to all be in phase and have a current magnitude of 1. Remember in MEEP units, this current then becomes the current base unit equivalent to 1 ampere. For this simulation, lets assume 1.0 in MEEP units is equivalent to 1mA (thus our base unit for current is in mA). Also, Hertzian dipoles have a constant magnitude AC current at a fixed frequency, thus we use MEEP’s ContinuousSource.

Lets decide the frequency as 200GHz. With our chosen base of \(a=1\text{mm}\), this frequency in MEEP units is:

In MEEP, the direction of current in a current source is specified by a component argument for sources, which can be either Ex, Ey, Ez (for electric current sources in the x, y, z directions) or Hx, Hy, Hz (for magnetic current sources in the x, y, z directions). So for z-directed electric current, the component is mp.Ez. Thus, to define the sources as in the image above, the code is:

freq = 0.66713

sources = [mp.Source(src=mp.ContinuousSource(freq),
                     center=mp.Vector3(x=xi, y=0, z=0),
                     component=mp.Ez, 
                     amplitude=1.0) for xi in (2.0, 1.0, 0.0, -1.0)] # 1mA source

Note: to change the phase of each AC source relative to one another, just make amplitude a complex number, which will define the magnitude phasor for that source. For example, making amplitude=1.0+1.0j will make that source have a +45° phase and magnitude of \(\sqrt{2}\).

Add media

By default MEEP assumes everything to exist in a vacuum, and to keep things easy let’s assume for the example that the dipoles above exist in a vacuum as well. To define pieces of geometry, we go:

geometry = [mp.Block(mp.Vector3(mp.inf, mp.inf, mp.inf), # define an infinite block
                     center=mp.Vector3(0, 0, 0), # centered at the origin
                     material=mp.Medium(epsilon=1))] # that is a vacuum

Add cell boundary

Basically we just want to know what the near-fields of such an antenna array will look like, so ideally we should not have boundary conditions and the waves will propagate forever without reflection or termination. But, MEEP needs some concrete way to handle such a thing (without simulating an infinite cell size!). So to do this they introduce perfectly matched layers, which just absorb all fields thrown at it without reflection, as if the medium continued on forever. The docs describe these very concisely:

As for boundary conditions, we want to add absorbing boundaries around our cell. Absorbing boundaries in Meep are handled by perfectly matched layers (PML) — which aren’t really a boundary condition at all, but rather a fictitious absorbing material added around the edges of the cell.

They just need to be reasonably thick to prevent numerical errors, so let’s define a 1mm thick layer around the edge:

pml_layers = [mp.PML(1.0)]

Note: this PML will overlap the outer 1mm border of the cell (not increase the size of the cell by a 1mm border). i.e it will overlap whatever is in a 1mm rim of the cell’s border.

Define simulation

TL;DR: combine previous steps with

sim = mp.Simulation(cell_size=cell,
                    boundary_layers=pml_layers,
                    geometry=geometry,
                    sources=sources,
                    resolution=resolution)

2.3 Run the solver

Now that the simulation is defined, we just need to start the simulation! As mentioned earlier, this can be done either in the time-domain or frequency domain. Put simply, the ways to do both are given below (you only need to pick one):

Solve in the time-domain

With a 200GHz frequency, each period is 5ps. Let’s simulate in the time-domain for roughly 120ps, which is:

So, lets run the simulation for 35 units for time (35 lengths of the base time \(b\)):

sim.run(until=35)
# If we want to save the E-field z-components every 0.1 units of time,
# then, instead of the above, we can run:
sim.run(mp.to_appended("ez", mp.at_every(0.1, mp.output_efield_z)),
        until=35)

Solve in the frequency-domain

Basically this just finds the steady-state behavior of these constant current sources, so no time conversion or anything is needed. But, to use the frequency-domain solver, you must set the simulation object to use complex fields. i.e, when we define sim, we must add the keyword argument force_complex_fields=True. After that, we just go:

sim.init_sim()
sim.solve_cw()

3. Interpreting results

We can now plot the resulting field quantities however we want! Use sim.get_array(...) (and its associated documentation) to get arbitrary array quantities at the end of simulation. And, if we saved some field quantities to a file while the simulation was running, we can get a tensor of all the field quantities for every point in the cell for all instants in time where we recorded it.

Continuing with the running antenna array example, say we used the time-domain solver and appended the results of the z-component of the E-field to a file every 0.1 units of time (as in the second sim.run command given above). Note that 0.1 units of MEEP time corresponds to

which is well under a 5ps period. Next to note is that MEEP saves the tensor of field data as a .h5 file. So to read the data, you will require the h5py package. Once that is installed, we can go:

import h5py
f = h5py.File('ez.h5', 'r')
efield_z = np.array(f.get('ez')) # a float tensor of shape (600, 600, 350)
f.close()

The tensor given efield_z has 3 dimensions, the first two are the x- and y-axes, while the last dimension is the time axis. You can plot these any way you feel like, below is one such way which I came up with:

In this animation, each frame is 60ms long and represents the 0.1 time interval between recordings of the field value in the simulation. This plot was made with matplotlib, using its FuncAnimation function and using the jet colormap. Note it plots the absolute value/magnitude of \(E_z\), and not the signed value (which oscillates between positive and negative). An interesting factor to note is how the colorbar indicates the scale of the E-field in SI units. To do this, we just note that an E-field magnitude of 1.0 in MEEP units is converted to SI units using the table earlier:

The colorbar is then normalized such that a magnitude of 1.0 in the floating point array shown in the animation (MEEP units) has a magnitude of 376.6V/m (SI units) on the colorbar(using the norm argument for matplotlib colorbars). This value is somewhat reasonable for within a few millimeters of Hertzian dipole antennas with 1mA current flowing in them.

Broadside array

From the theory we know that a linear Hertzian dipole array of constant phase will be a broadside array – the main lobe/maximum of its radiation pattern is in the direction perpendicular to the linear array structure lying on the x-axis. Although the animation shows the near-field behavior of the antenna array, we can still see that the field is strongest in the broadside (y) direction, agreeing with the far-field theory.

There is also some angles at which the field seems to completely cancel and the magnitude drops to near zero. These are the nulls of the radiation pattern. Along the x-axis there also seems to be side lobes of lower peaks of E-field magnitude. The direction of the main lobe, nulls, and side lobes can be programmatically changed by adjusting the phase of each Hertzian dipole source. For example, if we keep the exact same setup as before but at 10 MEEP time units change the phases of each source to increase from 0 to 0.6π from right to left, then the field animation looks like:

The direction of the main lobe has certainly rotated, and the pattern of the nulls and side lobes have changed as well, albeit it is a bit hard to see them clearly since this is a near-field simulation. Note that the colormap is capped at 376.6V/m, so any magnitudes above that appear as the darkest maroon in the animation.

Conclusion

Anyway, I hope this helped a bit with using MEEP in python, and in particular converting between MEEP units and SI units. Have a comment, suggestion or noticed something wrong? Please let me know under the About section.