Computational Electromagnetics with MEEP Part 3: 2D MEEP

Lesson 1: MEEP in Two Dimensions

←Previous: Resonators
Outline
Next: Gaussian Beams →

Finally!

Ah, at long last - learning how to do electromagnetics in two dimensions. We got really far with one dimension, but it has a serious limitation - it can only simulate structures which are infinite in two dimensions (also known as 'slabs'), and can really only work with fields that are, in some way, pointing 'towards' the slab (not fields propagating along the axis of a slab). Lots of interesting stuff happens in two dimensions that we could never simulate in one dimension - things like waveguides, ring resonators, optical beams, and structures which are generally finite along more than one dimension.

Baby Steps, As usual - Simulating a Simple Source

Let's start off with the simplest thing we can think of - a continuous source at a single point, just to check that everything makes sense. As a helpful reminder, here are the things we need to add to our simulation to get everything to run (discussed in lesson 4):

We need to add our sources, our geometry (if we have any, in free space we won't), our PMLs, and the other options the simulation needs to run. These include (bare minimum) the resolution, the size of the simulation domain, and when to terminate the simulation. We also want to be able to actually see what is going on, and we know from lesson 5 how to create simple animations of our fields. The only added wrinkle here is that we are working with an extra dimension, so when we store fields, we are storing a bunch of matrices, and we have to be mindful of that when accessing the data. Fortunately, other than this pretty much everything is the same as in working with one dimension.

Propagation

What should we check first? Well, let's just check that the plane wave sources we have been using up until this point behave as we expect. Let's make sure that they propagate at the speed of light and everything looks correct. Let's add a continous source at the origin, with a simulation domain of 10 units by 10 units (MEEP units, as usual).

How long should we run our simulation if it is 10 units by 10 units and the source is at the center for the waves to barely touch the edges of the simulation domain by the end?




Ok, now let's run our simulation (full code posted below as usual). What do we see, and does it make sense?


Well, if we were expecting a plane wave, we have clearly been surprised. What we actually have is a point source (or technically a line source, since it's a point source in two dimensions). But at least the speed of the source makes sense - the outer reddish ring barely starts to touch the boundary of the simulation at the end of the simulation. Point sources are all fine and good, and we will use them plenty, but I want a plane wave, because I know how to deal with those. To create one of these, we just need to extend the source along the y-axis, by adding a single argument to the Source() object: 

sourceSize = mp.Vector3(0, 10, 0)
sources = [mp.Source(mp.ContinuousSource(frequency=frequency),
        component = mp.Ex,
        center=sourceLocation,
        size=sourceSize
        )]

Now what does that look like?


Hmmm... it seems to be moving at the right speed, but it looks a little funky near the edges. This will actually cause problems. Why is this happening? Well, there's a PML now at all boundaries, not just the edges along the z-direction! Fortunately, there's an easy fix, we just need to set the is_integrated option for the source to True: 

mp.ContinuousSource(frequency=frequency, is_integrated=True)

Now what does it look like?


Ahh, much better. That's more like what I expect to see. Let's have some more fun with this plane wave, shall we? Let's try to tilt it, just like we did back in lesson 13, by adding the k_point argument to our simulations. Let's start with 45 degrees. 

theta = np.radians(45)
kVector = mp.Vector3(np.sin(theta), 0, np.cos(theta)).scale(frequency)

simulation = mp.Simulation(cell_size=cellSize,
        resolution=resolution,
        sources=sources,
        boundary_layers = pmlLayers,
        k_point=kVector)

??? What gives here? We specified our k-vector, but the plane wave doesn't appear to have tilted at all! It looks almost exactly the same!

Boundary Conditions vs. Phase

The problem is that we've only specified the boundary conditions - what happens when the plane wave crosses from the bottom of the simulation back to the top. But because it starts off as a planewave propagating at normal incidence, it never has a chance to cross the boundaries and become the tilted plane wave it was meant to be. We need to explicitly add the tilt to the plane wave itself, and we do this by including another argument, amp_func, in the Source definition. 

def tiltFunction(vec):
    return np.exp(1j*2*np.pi*(kVector.x*vec.x + kVector.y*vec.y + kVector.z*vec.z))

sources = [mp.Source(mp.ContinuousSource(frequency=frequency, is_integrated=True),
        component = mp.Ey,
        center=sourceLocation,
        size=sourceSize,
        amp_func=tiltFunction
        )]

This might be a somewhat awkward way of doing things (if you aren't already aquainted with Fourier Optics). We're basically just telling MEEP to multiply our line source by \(e^{j \overrightarrow{k}\cdot \overrightarrow{r}\), which is a plane wave with k-vector \(\overrightarrow{k}\) in phasor notation.


Ah, finally! It's got some tilt to it. Huh. That's weird. It looks like near the bottom of the simulation domain this looks like a planewave, but near the top it gets weaker and isn't even propagating in the right direction! Why might this be? Well, it's as if our source is acting like a point source, with waves near the top spreading out sort-of-in-a-circle. This is because this source is finite. To fix this, we need to make the source infinite. We do this by applying bloch-periodic boundary conditions in the x (vertical) direction. But didn't we already do that?

That pesky PML

We did, but we also included a PML. By default, MEEP will put the PML around all the simulation boundaries, not just those along z. So we need to remove the PML from the \(+x\) and \(-x\) boundaries. How do we do that? Well, we specify explicitly that we only want the PML to exist on the z boundaries: 

pmlLayers = [mp.PML(thickness=pmlThickness, direction=mp.Z)]

Let's run things again and see how they look:


Finally! We finally got a tilted plane wave. This plane wave is infinite along the x- direction, rather than being finite as it was before. However, this also means that whatever devices we simulate will also be infinite, or more precisely, periodic along the x-direction. Often this is desired, but if you want to do a 2D simulation this way of a finite device, you have to resort to more complex ways of creating a plane wave, which I might make an article on if someone harasses me :). In the next lesson, we're going to use our newfound friend the amplitude function to start shaping waves. In particular, we're going to revisit our best friend the Gaussian.

Full Code

You can find the full code here. 

import numpy as np
import meep as mp
import matplotlib.pyplot as plt
from matplotlib import animation

frequency = 1
resolution = 32
pmlThickness = 1.0
cellSize = mp.Vector3(10,0,10)
sourceLocation = mp.Vector3(0,0,0)
sourceSize = mp.Vector3(10, 0, 0)
animationTimestepDuration = 0.05
endTime = 8
theta = np.radians(45)
kVector = mp.Vector3(np.sin(theta), 0, np.cos(theta)).scale(frequency)

def tiltFunction(vec):
    return np.exp(1j*2*np.pi*(kVector.x*vec.x + kVector.y*vec.y + kVector.z*vec.z))

sources = [mp.Source(mp.ContinuousSource(frequency=frequency, is_integrated=True),
        component = mp.Ey,
        center=sourceLocation,
        size=sourceSize,
        amp_func=tiltFunction
        )]

pmlLayers = [mp.PML(thickness=pmlThickness, direction=mp.Z)]

# First, run a simulation to get the shape of the field array
simulation = mp.Simulation(cell_size=cellSize,
        resolution=resolution,
        sources=sources,
        boundary_layers = pmlLayers,
        k_point=kVector)

simulation.run(until=0)
fieldEy = np.array(simulation.get_array(center=sourceLocation,size=cellSize, component=mp.Ey))
fieldData = np.array([np.zeros(fieldEy.shape)])

def updateField(sim):
    global fieldData
    fieldEy = np.array([sim.get_array(center=sourceLocation,size=cellSize,component=mp.Ey)])
    fieldData = np.vstack((fieldData, fieldEy))

simulation = mp.Simulation(cell_size=cellSize,
        resolution=resolution,
        sources=sources,
        boundary_layers=pmlLayers,
        k_point=kVector)

def animate(i):
    im.set_data(fieldData[i])
    return im,

simulation.run(mp.at_every(animationTimestepDuration, updateField), until=endTime)
fieldData = np.real(fieldData)
print(fieldData.shape)

fig, ax = plt.subplots()
im = ax.imshow(fieldData[0], cmap='RdBu', vmin=-1, vmax=1, extent=[cellSize.x/2, -cellSize.x/2, cellSize.z/2, -cellSize.z/2])
fig.colorbar(im)
numberAnimationTimesteps = len(fieldData)
fieldAnimation = animation.FuncAnimation(fig, animate, frames=numberAnimationTimesteps, interval=20, blit=True)
fieldAnimation.save('basic_animation.mp4', fps=30, extra_args=['-vcodec', 'libx264'])
plt.show()
←Previous: Resonators
Outline
Next: Gaussian Beams →

If you found this content helpful, it would mean a lot to me if you would support me on Patreon. Help keep this content ad-free, get access to my Discord server, exclusive content, and receive my personal thanks for as little as $2. :)

Become a Patron!