FDTD: Total Field Scattered Field Hard Source and PML

Currently, I have a working 3D FDTD code with the Total Field/Scattered Field formulation. I am using a 3D plane wave constrained to propagate in one direction. However, I get severe leakage/reflections at the TFSF box edges (most likely a bounds issue, but I cannot figure out where I'm going wrong). Therefore, I am attempting to implement a 1D auxiliary field as a hard source (method presented in Taflove). However, I cannot figure out how to implement this 1D auxiliary field into my code.

Also, if anyone has any knowledge with a 3D Perfectly Matched Layer Absorbing Boundary Conditions, it would be much appreciated!


P.S. I wrote my code in Fortran90.

Attached screenshots of the total electric field (incident) propagating in the z-direction.

https://obrazki.elektroda.pl/7291692700_1416344333.png
https://obrazki.elektroda.pl/3198694900_1416344334.png
https://obrazki.elektroda.pl/2081102900_1416344335.png
https://obrazki.elektroda.pl/5665320500_1416344336.png

I lack the knowledge to help, but your screenshots remind me of Paul Falstad's simulations. He makes several which have to do with electrodynamics (as well as electronics and other physics areas).

https://www.falstad.com/mathphysics.html

He programs in Java. He puts comments into his code. You might pick up something even without knowing Java.

If you get immediate "leakage" its most likely due to not having the time delay correct between the E and H fields in your source. If the "leakage" gets more severe at the opposite side of the grid than you launch your source, it is probably a numerical dispersion problem. This can be helped with better grid resolution or by compensating for it through the constitutive values. Otherwise, maybe there is a mistake in how you implemented your TF/SF.

Try using finer resolution and see if the problem gets better. That will tell you a lot about what is happening.

I'm checking with finer grid now. I guess the term "leakage" isn't exactly correct. It's almost like the plane wave is scattering slightly off of the TF/SF boundary. Is this normal?

Also, it is important to mention that I am not using the 1D auxiliary field currently. I set up a 3D plane wave but constrained it to only travel in the z direction (x polarized).

As the plane wave enters the TFSF (k=k0 plane), i get slight scattering off of the boundary which continues until the plane wave approaches the opposite tfsf boundary. Once it reaches the other side (k=k1 plane), it appears that it does not fully exit and scatters quite a bit in each direction. This becomes a problem when I try set up monitors to collect the scattering and transmission of a gold nanosphere.

If your grid is filled entirely with air, you should never see anything in the SF regions, which I assume are running around the entire outside perimeter of your grid. Is that correct? Don't let your TF/SF slice through a PML or you will have trouble.

I really think you are experiencing a numerical dispersion problem. The wave in the FDTD grid is travelling at a slightly slower speed than a physical wave would due to the interaction with the grid. So at the edges of your TF/SF, you are adding/subtracting a plane wave with a different phase than the actual wave in the grid. This phase difference causes what you are seeing. The sign this problem is happening is that it gets worse as the plane wave propagates and using finer grid resolution will help it. The fix is to compensate for grid dispersion by adjusting the values for permittivity and permeability across the grid. This "boosts" the speed of a wave to compensate.

Not sure if you are aware, but I have a lecture series on FDTD here: **broken link removed**. Sounds like you have advanced passed what is in this course, but it may help you. Lecture 10 contains the grid dispersion and compensation topic. See slides 16-18.

Turns out, all of my derivatives were backwards from what is presented in Taflove (which is still valid); however, I was using the indices found in Taflove for my TF/SF boundaries so they were shifted by one grid point in the wrong direction. I corrected the code to match the forward derivatives in Taflove and the corresponding TF/SF indices.

Although, I still see these cassette tape looking patterns arising after the plane wave enters the TF/SF (see image below). I completely understand the dispersion issue with the other side, but why does this occur upon entry into the TF/SF? Is it because I'm not using a 1D auxiliary field for my plane wave? I used a 3D plane wave but set it up to only propagate in one direction (i.e. a pseudo 1D plane wave). Even still I wouldn't think I much dispersion since it's directed along the main grid (no angles).

I guess long story short, how do you correct for the phase difference where there is no delay from propagating your plane wave at an angle?

Also, very nice lecture series. I have watched a few of them within the last month or so.


Cassette Tape Scatter Pattern from PW hitting the TFSF boundary
https://obrazki.elektroda.pl/4109325200_1416552539.png

Dispersion pattern upon entry out of TFSF boundary
https://obrazki.elektroda.pl/9287039900_1416552539.png



As an aside, how hard is it to implement a PML boundary? Currently I'm using a damping function I cooked up which seems to work okay but there is a lot of wasted real estate. I have been through the formalism in Taflove but it looks very daunting to code up.

It is hard to see what is happening with still pictures. You may want to make movies. I can very often figure out what is wrong with that. My best guess with the scattering at the beginning is not having the time exactly right between the sources you use on either side of the interface.

I think the PML is the hardest thing in FDTD. They are a pain, but worth your time. Depending on what you are modeling you may be able to get away with a Mur boundary condition that is a lot easier to implement. If you are in 3D, go immediately for the stretched coordinate convolutional PML. Don't both with the UPML or anything else. You may want to check out ToyFDTD. I think they use the convolutional PML, but I could be wrong. I know I have seen code with it implemented in 3D.

Here is a .gif of the plane wave entering the TF/SF and exiting. There is a slight reflection upon entry and quite a bit of reflection/scattering upon exiting the boundary.

Still not sure what is going on here.

What happened when you used finer grid resolution? Be sure to do that representing the same physical size of the grid, just with smaller cell size so there will be more points. Did that not fix the problem? Otherwise, I suspect it is the delays in your calculations for the TF/SF interface.

The finer grid resolution decreases the amount of reflection at the left TF/SF boundary. However, it doesn't make much of a difference for the right TF/SF boundary. How do you correct for the delays in the TF/SF interface? I read over Taflove's discussion for an arbitrary 1D auxiliary plane wave as a source. However, his is appears to be generalized for any angle of incidence.

Use super find grid resolution for this test to eliminate that as a possibility. Set the longitudinal size of your box to something very short, maybe on 4 cells or so. Does that work much better? If so, I think you have a numerical dispersion problem. It not, a delay problem. By delay problem, I am not talking having to incorporate something artificial. I am talking about you possibly having a mistake where you forgot to account for a half time-step or delay through half a grid cell or something like that.

BTW, the auxiliary method you keep mentioning is very good for eliminating the numerical dispersion problem, but I don't think that is what is happening to you. Are you using this? If so, maybe the problem is there.

Of course, this is my best guess as to what is going wrong.

Attached are two animations illustrating the effect of making the grid finer. In this example, I went from a dx=dy=dz=16 down to a dx=dy=dz=8. As you can see, the negative wave that results as the plane wave travels through the Total Field Region is much smaller for the finer grid spacing. Is there any way to get rid of this completely without making the grid spacing infinitesimally smaller? Because of my home-cooked Absorbing Boundary Condition, I have a lot of wasted real estate and I'm overloading the memory on my server. I can try it with a MUCH finer grid spacing on our large memory node of the supercomputer, but I think this will take too long to diagnose.

Also, why is the plane wave scattering off of the TF/SF boundary as it enters the Total Field region? At this point, there shouldn't be much dispersion. I went back to check my indices for the TF/SF boundary, but I couldn't find any mistakes. However, it could be possible that the incident is not on the same time grid (Off by half-time step) as I had to implement a multi-lorentz material for gold and silver. This calls for an Auxiliary Differential Equation for the Electric Field. Although I'm not using an material in this test, the ADE method requires E, E(n+1), and E(n-1). I am assuming that I have chosen the correct E (in time) to apply the TF/SF boundary to.

It is important to note that I am simulating gold and silver nanoparticles, as well as, gold-silver core-shell nanoparticles. My code works very well for the gold or silver nanoparticle alone, however, I run into issues when I create a core (gold) encapsulated by a shell (silver). The optical response is not even close to what I see experimentally. I just wanted to rule out this TF/SF problem before I proceeded.

I am getting over stomach bugs so my brain is not working very well, but I will try to keep helping here. Read a bit about numerical dispersion in the Yee grid. You can calculate a phase velocity on the grid and base your analytical wave on that phase velocity rather than the physical phase velocity I think you are using. You could also adjust your mu and eps parameters to compensate for the dispersion. This is probably the better approach since your simulation will be more accurate as well.

A whole different approach is to just have one TF/SF interface across the grid and not have to worry about three other boundaries. Is this possible or are you doing something where this would not work for you?

I will read more in Taflove. What is the best approach to adjusting the mu and eps parameters? Is there some systematic approach?

Unfortunately, I have to stick with the 4 TF/SF boundaries as I am measuring the extinction of the metallic nanoparticles. Later on, I will be focusing on scattering and absorbance individually.

To adjust mu and eps, you end up calculating a "fudge factor" f. Then you simply multiple all your mu's and eps's by the fudge factor. f will be on the order of 0.99. This will speed up the wave in your grid to match the physical wave.

Finally fixed the issue. Turns out it was a combination of few problems. First and foremost, you were correct about being off by a timestep. I have to reorder my curl updates and break my incident subroutine into two separate parts (E and H incident). Once I did that, I quit having reflections or scattering at the box edge. My other problem was that the gaussian was too narrow in time which lead to dispersion within the grid. Thank you for all of your help!

Great! Congrats!