Computational Electromagnetics with MEEP Part 1: Getting Started
Lesson 3: (Extreme) Basics of FDTD
Before we actually start writing our first program, there's a couple more things you should know about how simulators work.
Garbage in, Garbage Out
When simulating anything, you should always start with something you can verify by hand, and you should constantly be asking yourself the question - does this make sense? This makes sure you know how to use the simulator, and that you aren't feeding it nonsense. There's a saying in the simulation world - "Garbage in, Garbage out". If you tell the simulator to simulate something ridiculous, it will do it, but the results will be worse than worthless (garbage), and even if they seem to make sense, you have no way of knowing. Let's start with something we can definitely verify by hand - a plane wave propagating in free space.
But what about infinity?
Planewaves are defined over all space, and we would like space to be infinite. MEEP (and most other simulators), unfortunately despite the best efforts of hordes of clever engineers, only have a finite amount of memory. So we have to restrict ourselves to a finite amount of space - this is called the "Simulation domain". So we have to choose how big our simulation domain should be. This is something of an art, but as a rule of thumb it should be as small as possible for you to get the answer you are looking for. Bigger simulation domains mean more computational time and more memory, which translates into unnecessary waiting.
Suppose you are interested in watching a plane wave of wavelength \(1\mu m\) propagate by \(5 \mu m\) in free space. How long should your simulation domain be (in real units)?
Planewaves are also defined for all time, and so I would love to run our simulation from \(-\infty\) to \(\infty \). Unfortunately, we run into the finite memory problem we had before. Besides, we won't be alive long enough to watch the results even if we had enough memory. This means we have to choose a time when we stop our simulation. As with the simulation domain size, longer simulations take longer to run, so we would like to choose a time that is as short as possible, but long enough to capture the important effects we are looking for. Let's do another example to illustrate the point.
Suppose you are interested in watching a plane wave of wavelength \(1\mu m\) propagate by \(5 \mu m\) in free space. How long should your simulation time be?
Split Space into pieces
As if things weren't bad enough, with only being able to simulate a finite time and a finite space, we can only simulate a finite number of points within our simulation domain. As you might imagine, if you don't choose enough points, you won't have an accurate representation of what's going on. For example, say you are simulating a planewave in one dimension. We actually want it to look like a continuous sinewave:
But, if we don't have enough points (say we choose 5 points per wavelength), that sinewave might end up looking like this: Looks more like a triangle to me than a sinewave. There are other, more subtle and serious consequences if we don't choose enough simulation points per wavelength - one being that waves in free space will not actually propagate at the speed of light. To avoid this and other unpleasantness, the developers of MEEP recommend having at least 8 points per wavelength (more will sometimes be necessary). In MEEP, you set this by changing the 'resolution' (let's call it \(R\)) - the number of sample points per characteristic length \(a\). If we want the number of points per wavelength, we just need to multiply the resolution by \(\frac{\lambda}{a}\):
\( points / \lambda = R * \frac{\lambda}{a} \)
So, as a starting point, by we should have a resolution \( R \geq 8*\frac{a}{\lambda} \). Note that this is not the free space wavelength, but the wavelength in your highest-index material (where the wavelength is the smallest). You want to make sure that you have enough points everywhere in your simulation domain. Let's do an example to hammer this home.
As before, suppose you have a planewave of wavelength \(1 \mu m\), and you chose your characteristic length \(a\) to be \(500nm\). What should your resolution be (at least)?
And time, too!
Time is also split into discrete chunks so that we can fit a simulation in computer memory. Each such discrete point is called a 'timestep', because your simulator steps through time, one point after another (so we start on timestep 0, and then timestep 1, and so on...). In FDTD simulators, the time between timesteps (let's call it \(\Delta t\)) is actually coupled the resolution by this thing called the 'Courant Factor'. This prevents numerical instability in the simulator. The actual timestep \(\Delta t\), in seconds is given by
\( \Delta t = \frac{S}{R}*\frac{a}{c}\)
Where \(a\) is the characteristic length, \(S\) is the Courant Factor, \(c\) is the speed of light, and \(R\) is the resolution. In MEEP, the default value of S is 0.5, and this should work for you unless you are doing something really weird (with refractive indices significantly smaller than 1). So, once you choose your resolution and characteristic length, you know what the time between timesteps must be.
Suppose your characteristic length was chosen to be 300nm, and your resolution was selected to be 10. What is the amount of time between timesteps?
Simulation Setup
OK, now with our newfound knowledge let's figure out the simulation parameters we want. I'm going to choose a plane wave of a nice color, green, which has a free-space wavelength \(\lambda_0\) of about \(500nm\). Because we don't have any devices, let's just choose a characteristic length reasonably close to our wavelength but easy to remember - 1 micron. Now, let's figure out the frequency of the planewave in MEEP units (which we will have to enter into MEEP). Since we know the frequency (in MEEP units) is just the characteristic length divided by the wavelength, we get a frequency of 2. What about the size of our simulation domain? Since the wavelength is half a micron, I would like to see the wave propagate in space, so let's say we want to give enough space to see 10 periods of the wave. This means our simulation domain should be 5 microns long. Since our characteristic length is 1 micron, this means we need to enter 5 MEEP units for the length of the simulation domain (let's call it \(L_{MEEP}\)).
How long should we run the simulation? Well, I would personally like to watch the planewave get from one side of the simulation domain to the other, so how long should that take? Since this is free space, the planewave travels at speed c, and it has to travel 10 wavelengths (the length of the simulation domain). This means we should run the simulation for about 16.7 femtoseconds, or exactly 5 MEEP time units (in keeping with the trend, let's call this \(T_{MEEP}\)).
Finally, what should our resolution be? Well, we want at least 8 points per wavelength, and we know the characteristic wavelength and the wavelength in our highest index material (free space). So this means we want a resolution of at least 16.
\(a = 1\mu m\), \(\lambda_0 = 500 nm \), \(L=5\mu m\), \(T=16.6782fs\)
\(f_{MEEP} = \frac{a}{\lambda_0} = 2\)
\(L_{MEEP} = \frac{L}{a} = 5 \)
\(T_{MEEP} = T * \frac{c}{a} = 5 \)
\(R \geq 8 \frac{a}{\lambda} = 16 \)
Next, we will actually use this to set up our simulation.
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!