Week 5: Making a (Somewhat) Complete Implementation

April 12, 2024

Hello everyone, and welcome back to my blog. Today I have some very exciting news: most of the core aspects of my simulation program have been implemented! My blog today will focus on the process of implementation and the challenges I faced along the way.

Implementing the Grid class

The most important part of the simulation program is the Grid class, whose task is to calculate and update the fields using the discretized Ampere and Faraday’s Laws (essentially a 3D implementation of Yee’s algorithm). I’ll refrain from discussing any theoretical details of Yee’s algorithm here, but here is a good source if any readers are feeling curious.

In 3D, there are 6 field components to be calculated. Thus, every Yee cell now has six different field values to be stored (shown above). In theory, all of the components are offset by half spatial steps, but the halves are suppressed in practice so that all of the members of the Yee cell actually have the same indices.

Another aspect of indexing that I had to be careful of was setting the proper 3D array sizes for the fields. Notice that the E fields (black circles in the image above) “stick out” in the direction they face while the H fields (white squares) “stick out” in the two directions the field doesn’t face. For example, the Ex field is offset a half-step in the x-direction, while the Hx field is offset a half-step in the y- and z-directions. This means that if the grid is M blocks wide in the x-direction, the grid of Ex field values would only be M – 1 blocks wide. Applying this logic to all the fields, we have the below code for initializing the grid:

(Here, _init_3d_array is a helper function to initialize a 3D array with the given x, y, and z dimensions)

Moving on, I then implemented the update equations for the fields. The code below is the update equation for the x-component of the E field. ExE and ExH are grids storing scalar coefficients, while Ex, Hy, Hz are grids for storing the values of the respective field components. The rest is pretty self-explanatory, and particularly curious/motivated/masochistic readers can compare the code to the theoretical update equation for some more clarity.

Adding Boundary Conditions and Sources

After implementing the Grid class, I was ready to add some basic field sources and boundary conditions to test the simulation. As promised last week, I implemented a simple ABC (absorbing boundary condition) based on the advection equations. To make coding easier in the future, I first coded a BoundaryCond class that all boundary conditions have to inherit from. Notice that each boundary condition has to bind itself to a Grid in order to be usable. Moreover, the update() method allows for each unique boundary condition to have its own custom update equations.

The FirstOrderABC class inherits from BoundaryCond (nothing special there). The hardest part of implementing the ABC was making sure the implementation of the update equations was bug-free. To show a bit of what I mean, here’s the code for updating the fields along the YZ oriented boundaries of the simulation space:

For the more insane readers out there, here is a good explanation of how ABCs work in full detail. Although it may be hard to tell at first, the code above resembles the advection equations I showed last week, but this time discretized and rearranged to solve for a particular term (it might be helpful to note that EyX0 stores the Ey fields at the X = 0 boundary of the simulation space at the previous timestep, and so on).

I also implemented (with much less frustration, to my relief) some EM wave sources. Right now, all the sources I use are additive sources, meaning that a value is simply being added onto certain points in space at given times (as opposed to using a TFSF boundary).

The code is pretty self-explanatory; I’ve implemented one sinusoidal waveform and one Ricker wavelet:

Visualizing the Simulation

Visualization is one of the most important tasks of the simulation program. Without proper visualization, it would be hard for users to decipher the results of the simulation.

To effectively visualize the simulation, I wrote a program to take “snapshots” of 2D slices in the 3D simulation grid at given instants in time. The snapshotting is performed by the snapshot method of the GridVisualizer class, as shown below:

The argument “dim” refers to the dimension along which the 2D snapshot slice is taken and “slice_index” refers to the coordinate along the chosen dimension to slice the grid. So for instance, if dim = 0 and slice_index = 10, the snapshot would be of the YZ plane at the coordinate x = 10. An example of what a snapshot would look like is shown below.

By compiling many snapshots together, a GIF of the field’s evolution over time is created. Shown below is the visualization of the Ex field of an electric dipole:

Here’s another one, but it’s the Hy field of the electric dipole:

Conclusion

Thanks for reading my blog all the way through! If you want more insight into the details of all the theory behind my code, you can take a look at the sources below (or look at my GitHub repo https://github.com/switchpiggy/yee-algorithm). Have a great day! (Make sure you get some sunlight, especially after experiencing the ordeal that is reading my blog).