PML behaviour and how to reduce reflections

[galcerte]

I've been experimenting with PMLs in FDTDX and I still haven't quite gotten to reduce reflections significantly. I was looking for a pretty significant drop in the reflected radiation's power flux, maybe something like -60 dB or even less. However, I've only managed about -35 dB.

Now, I am new to computational electromagnetism, so bear with me when I ask this but; if I understand correctly, the kappa parameters are for controlling how much grazing incidence radiation gets absorbed right? And alpha is more of the same? If that is the case, how should I adjust the PML such that the reflected radiation gets weaker (PML absorbs more)? These two don't really serve that much of a purpose to me as my radiation incides normally. I read that sigma should be scaled up, as the negative natural log of reflectivity is approx the sum of sigma/kappa over all Yee cells. I tried this, and didn't get any decent results. I also read that increasing thickness works; I tried massively increasing it (I was using 24-48 Yee cells before; went from that to 150) and that did give me results, but still not enough (-40, -45 dB). However, that's not feasible, as I need to keep the PML thickness down, since that adds a lot of computation time.

On that note, I came across several conflicting suggestions on the PML thickness. ChatGPT says 18-36 Yee cells; my colleagues say 2-4 lambda (which in my case would equate 150-600 cells so completely unworkable with a uniform mesh), and Meep's docs say 0.5 lambda. I couldn't even simulate using more than 1 lambda of thickness, as any more would induce an integer overflow. Maybe the reason I am having so much trouble is because of how small my Yee cell resolution is with respect to the wavelength (resolution is 0.1 mm and lambda ranges from 14.9 mm to 149.0 mm), but I'm afraid I need to keep resolution that way.

Pictured below is a plot of the absolute value of the Poynting flux in an empty simulation volume with just PMLs in the top and bottom, periodic boundaries elsewhere. I have a plane wave source with Gaussian temporal profile up against a PML which launches a pulse towards a detector first and then the opposite PML behind it. You can clearly observe the incident pulse, the reflection on the PML, and the twice reflected pulse. No idea where that background flux comes from, that's something I'd also like to know about. This is with 40 GHz, detector is about 12.2 mm away from source.

[ymahlau]

I am afraid I cannot give you good advice on this issue. Some ideas are:

• are you using float32? I think for such small numbers you will need float64. We have not tested a lot with float64 data type, maybe there is an issue there?
• Maybe other types of PML are better suited for this problem? Or absorbing boundary conditions? We just chose to implement the convolutional PML because it seemed decent for most use cases, but probably there are better variants for your use-case in literature.
• Regarding PML thickness, the thicker the better but at some point you get diminishing returns and just waste computational power. I think it always depends on the use-case.

I have never tested the PML for such fine resolution, but it is good to know the limitations. Currently, I am doing an internship and not very active for fdtdx. From February next year you can expect to see more activity here again. But, if you need advice or help implementing something feel free to ask!

[galcerte]

So I'm studying the CPML itself and its implementation in FDTDX, and I'm a bit confused at times. Do you happen to have handy references for the papers or books that you used to write the CPML? I'm kinda lost with some aspects, especially the dimensions of each quantity. They seem to be normalized, but I am not quite sure what constants they are divided by. So far I have located the Schneider book https://eecs.wsu.edu/~schneidj/ufdtd/, mentioned in update.py and it looks like in some places you follow the sign convention of Taflove and Hagness' Computational Electromagnetics: The Finite-Difference Time-domain Method. I guess I should also look at the Roden and Gedney paper I mentioned in #164 which I found in the FDTDX whitepaper; is there any other reference that could be of use?

[ymahlau]

Like I said, I am definitely not an expert on the units myself, mainly because a colleague implemented the PML in fdtdx. With the upcoming unit system changes (see here, we will add units to the pml parameter values as well, though I am currently not sure what the correct units for that would be. If you reach any insights on that please let me know! Since I am currently a bit involved with my internship, you may expect the unit system to come at the beginning of next year.

Regarding the PML implementation, we have gotten most of our implementation from the FDTD software of Floris Laporte: https://github.com/flaport/fdtd?tab=readme-ov-file#perfectly-matched-layer

There is a pretty good explanation in the readme on how the PML works.

[galcerte]

though I am currently not sure what the correct units for that would be

alpha and kappa are both dimensionless in FDTDX, but alpha usually has units of conductivity (as per Schneider and other papers), so FDTDX's alpha must be implied to have been multiplied by a constant with units of resistivity. kappa is dimensionless either way, as it has the same dimensions as the stretching coefficients of the coordinate system. If you look at the definition of the bE and bH coeffs in perfectly_matched_layer.py, there's an exponent, its argument must be dimensionless always. Assuming the Courant number is also dimensionless in FDTDX's unit system, then sigma_H and sigma_E must also have the same dimensions as kappa, that is, none. But if you took a look at how these sigmas are defined, you end up in sigma_fn() in utils.py, where the actual sigma profile is defined as 40 * x**3 / (thickness + 1) ** 4, and this has dimensions of inverse length, not dimensionless, if 40 has no dimensions. However, due to this contradiction, 40 certainly should have dimensions... well, now you can see where part of my confusion comes from.

If you reach any insights on that please let me know!

Of course, as a matter of fact, I think I will open an issue, because as a result of my digging into this, I came to the conclusion that there are a few things in the PML implementation that don't quite match the literature. I have a lot of things to say on that regard, but I hope I get some scrutiny from others to make sure I am on the right track.

Regarding the PML implementation, we have gotten most of our implementation from the FDTD software of Floris Laporte: https://github.com/flaport/fdtd?tab=readme-ov-file#perfectly-matched-layer

I will take a look, thank you.

[galcerte]

• I have not tested with float64. I'd rather not use this data type, but I will try to use it anyways and report back, maybe there is some issue with numerical errors propagating due to insufficient precision? No idea.
• Yes, I did see that there are other kinds of boundary conditions out there, like adiabatic absorbers. I will tinker around further with the CPML a bit, and if I don't see any results maybe I will try to implement others, though that will take a good while...
• Yeah, I noticed that... I also tried lowering the wavelength and that led to higher absorption too, but again, I can't lower wavelength that much. I don't understand why this doesn't work well enough, while examples such as the vertical coupler in speed_benchmark.py have a set thickness much smaller than the wavelength (10 cells, one seventh of the wavelength). Well, the comments in the script did say that it doesn't reproduce the result on the listed paper. Maybe there's too much reflectivity there in the PMLs as well?

I have never tested the PML for such fine resolution, but it is good to know the limitations. Currently, I am doing an internship and not very active for fdtdx. From February next year you can expect to see more activity here again. But, if you need advice or help implementing something feel free to ask!

Okay, thank you very much, I'll only bother you regarding implementations probably. I will keep testing and investigating these issues and share any solutions or at least useful observations.

For now, I do have a couple; raising alpha eliminates this weird background flux that is observed in the plot above. I raised it quite a bit, from the default 1e-8 to 0.1. Lowering kappa also increases absorption. I tried increasing it to 2.5, but that worsened things, so I used 1.5 again (it's also used above), which is the default. However, I would not lower it beyond that value, as the reflections get distorted again. No idea what that would entail in a simulation with objects in the volume, but that is probably not entirely good. Here is a plot having simulated using these values.

[galcerte]

I just tried out simulating using float64 and apparently it does work. However, I had to set the environment variable JAX_ENABLE_X64 to True, and even then I get a warning, /home/user/project_folder/.venv/lib/python3.12/site-packages/jax/_src/ops/scatter.py:92: FutureWarning: scatter inputs have incompatible types: cannot safely cast value from dtype=float64 to dtype=float32 with jax_numpy_dtype_promotion='standard'. In future JAX releases this will result in an error. warnings.warn(, but it does seem to simulate using 64-bit floats, as the runtime is noticeably higher, and the result is different. Specifically, it seems that the noisy flux readings are no more, they are very smooth even as the value decreases by several orders of magnitude. Sadly however, the PML still reflects too much power.

[ymahlau]

The warning might also indicate, that some places in the simulation do not handle float64 correctly. This is what I was talking about above, you can try it out both data types. In the end it probably does not make a big difference.