Assignment II: Finite Differences, FDTD & MEEP

Prepare a .pdf file presenting your work in this assignment and upload it to Canvas.

1. The finite difference for Poisson Equation and Capacitance of Square Coaxial Cable

Objective: perform the following tasks and answer the questions:

• 1.0 Plot the potential on the grid obtained using the Gauss-Seidel solver we implemented in Matlab in the lecture (code here).
• 1.1 Perform a convergence test for the square coaxial cable. Explain how you perform the convergence tests and your expectation.
• 1.2 Determine the accuracy of the numerical solution experimentally: plot the error (wrt to a high-resolution reference run) vs. grid spacing and superimpose the plot of errors decreasing with h, h^2, or h^3. Which truncation error are you expecting?
• 1.3 Formulate the Poisson equation for the coaxial problem as a matrix problem. Your task is to create the A matrix and g vector to solve it using f = A\g Matlab command. No implementation is needed, only the formulation in matrix form.
  • What is the structure of the matrix, e.g., diagonal, tridiagonal,  block, ...?
  • What is the equivalent of symmetry BC in matrix form?

2. Finite-Difference / Complex Exponential / High-Order Differentiation / Stability

In the lecture Lecture: Finite Difference of Complex Exponentials, we saw that there are two possible ways of discretizing first-order equations with centered differencing, e.g., using one cell or two cells. Although they look pretty similar, they have different numerical properties.

Answer the following questions:

• 2.0 Which scheme should we use for the discretization of first-order equations with centered differencing? What is the reason? Motivate your answer.
• 2.1 Derive the numerical wavenumber for second-order derivative using the complex exponential technique (in the lecture, we used it to analyze the numerical property and numerical wave number for the first-order equation). For this exercise, check textbook section 3.2 and, in particular, section 3.2.3. 

Regarding the high-order (>2nd order) discretization:

• 2.3 What are the advantages of high-order finite difference schemes? What are the disadvantages of high-order finite difference schemes?

Regarding solving eigenvalue problems with the finite difference technique:

• 2.4 Use the 1D Maxwell solver code to calculate eigenfrequencies in the time-domain for a simulation 30*pi long and 31*pi long in time. What are the differences between the results? What is the reason? Check the textbook Section 4.4.3.
• 2.5 Use the code that we using a different two different time steps dt: d2tmax = 1.9*dx and d2tmax = 3.0*dx? Plot the eigenfrequencies s1 and s2, if possible. What is the stability condition for this numerical scheme? Check the lecture and textbook sections 4.4.1 and 4.4.2.

3. FDTD

In the lecture, we studied the propagation of a square pulse modeled with FDTD. 

Objective: Complete the following tasks and answer the questions:

3.0 Implement in Matlab, Python, or Julia the propagation of a Gaussian pulse. For this exercise, you can refer to the Matlab implementation presented in the lecture and change the initial condition to a Gaussian pulse from the square pulse.

3.1 What is the magic time step, and how is it defined?

3.2 Plot the pulse propagation at different time steps in simulations using

- the magic time-step

- 0.9*the magic time

What is the impact of not using the magic time-step, and why is the simulation affected?

3.3 Determine the numerical dispersion relation experimentally using a 2D FFT (space and time) of the history of Ez. For doing that, you can refer to the following snippet code in Matlab:

a) Record the values of the electric field Er at each grid point at each time step (we call this variable ErHist). Matlab snippet:

ErHist = [];
% Start loop
for t=1:nsteps
  ...
   ErHist = [ErHist; Er];
end

2) perform a 2D FFT on the history of the electric field (ErHist) and plot the logarithm of the absolute value of the electric field history. Matlab snippet:

spectrum = fft2(ErHist);
plot_spectrum = log(abs(spectrum));
pcolor(plot_spectrum(1:ke/2,1:nsteps/2))
shading interp

Question: How do you calculate the x and y-axis values?

Hint: You need to check the fft2 function online or check the textbook where it presents FFT with Matlab (Sections 4.4.2 and 4.4.3)

3.4 Compare the numerical dispersion relations obtained experimentally with the two methods and time steps (1-magic number) and (2-0.9*magic number). Which differences do you observe?

3.5 What is the Yee lattice, and what is its main characteristic?

4. Use MEEP to model a bent waveguide

Objective: The overall goal of this task is to expose you to MEEP and a step-by-step solution using an online tutorial on using MEEP. As part of this task, you need to install the MEEP Python code, familiarize yourself with it and implement the bent waveguide problem.

Complete the following tasks:

4.1: Install the MEEP Python using Conda (Official Releases) and the instructions at https://meep.readthedocs.io/en/latest/Installation/

Hint: You might get an error when running the code because you don't have six installed. To install six, use the command: 

 pip install six

4.2: Follow the step-by-step instructions at https://github.com/NanoComp/meep/blob/master/python/examples/bent-waveguide.ipynb (Jupyter notebook) to implement a simulation of the 90-degree bent waveguide using MEEP. Use the same geometry and material properties as the example. In this case, you can create a standalone Python code from the Jupyter notebook instructions and run it.

4.3: Plot the Ez value after 2,000 steps (time=200) in MEEP units.

4.3: Answer the questions: 1) what are the PML boundary conditions and how do you implement them in MEEP? 2) How is the time defined in MEEP?

Bonus Exercise (+1) -  Develop a 3D Cubical Cavity Solver and calculate the eigenfrequencies

B.1 Develop a 3D FDTD Maxwell solver using the Yee Lattice in Matlab, Python, or Julia with PEC boundary conditions and random initial values within the 3D cavity. Make sure to read Sections 5.2.3 of the textbook. The textbook guides the different steps in developing the 3D Maxwell's solver.

For the initial conditions and parameters, you can use the following setup (using Matlab):

% Physical constants
eps0 = 8.8541878e-12;         % Permittivity of vacuum
mu0  = 4e-7 * pi;             % Permeability of vacuum
c0   = 299792458;             % Speed of light in vacuum

% Parameter initiation
Lx = .05; Ly = .04; Lz = .03; % Cavity dimensions in meters
Nx =  25; Ny =  20; Nz =  15; % Number of cells along each axis
Nt = 8192;                    % Number of time steps

You will need to calculate the time step.

The electric field at the grid points can be initialized randomly (in Matlab):

% Initiate fields with noise (except on the boundary)
Ex( :  , 2:Ny, 2:Nz) = rand(Nx  , Ny-1, Nz-1) - 0.5;
...

B.2 Select a point within the domain and record the evolution of the electric field. Perform a spectral analysis (FFT) of the electric field evolution.

Plot the eigenfrequencies and answer the following questions:

    B.2.1 How do you initialize the fields?

    B.2.2 How do you choose the time step?

    B.2.3 How do you choose the sampling points?

    B.2.3 How do you calculate the eigenfrequencies?

B.3 Implement the 3D FDTD Maxwell solver using MEEP.  The objective of the exercise is to get some experience with MEEP, understand what the challenges are, and try to solve them using resources online. MEEP is a difficult code itself and it is not easy to navigate it is not expected that you to become an expert. It is not expected that you compare the results with Matlab as it would require a lot of work on making the configurations exactly the same.

For this task, you will need to consult the concepts you learned in exercise 4 and MEEP online resources.

There is no need to perform the spectral analysis using MEEP. You can use MEEP to calculate the evolution of the fields and use Python to perform the spectral analysis. 

Hints for the MEEP implementation

a. By default MEEP uses PEC or "metallic" boundary conditions so you don't need to set up them.

b. To initialize the electric field to random values, you can use the MEEP function

sim = mp.Simulation(cell_size=cell, resolution=resolution)

from random import random
from random import seed
seed(1)
def field_func(p):
    rr = random() - 0.5
    return rr

sim.initialize_field(mp.Ex, field_func)
sim.initialize_field(mp.Ey, field_func)
sim.initialize_field(mp.Ez, field_func)

c. To extract 3D data, e,g, the E_z

eps_data = sim.get_efield_z()

d. To print the value of a component (Ez) at each timestep (after 100 time units until 200)

def print_field(sim):
    fp = sim.get_field_point(mp.Ez, mp.Vector3(0,0,0)).real
    print("field:, {}, {}".format(sim.meep_time(), fp))
sim.run(mp.after_time(100, print_field),until=200)

This is not necessary and requires the usage of extra packages.

e. To visualize 3D data with a contour plot, you can use the mayavi module (you need to install it first with pip install mayavi and pip install pyqt5)

eps_data = sim.get_efield_z()
from mayavi import mlab
s = mlab.contour3d(eps_data,contours=10)
mlab.colorbar()
mlab.show()