A Loop-Tree Conundrum and a note on Scientific Integrity
For the last couple of weeks, I have been trying to understand the so called loop-tree and loop-star basis functions used in CEM. I had a very simple case implemented in bbmm, but that was only tested at higher frequencies. Last week, I started testing it at very low frequencies, down to about 3Hz, on simple geometries. It produced nonsensical results. Naturally, intense debugging followed.
I started by trying to reproduce the condition number, currents and RCS of simple geometries in the literature. In particular, I used the paper by Wu, Glisson and Kajfez (Applied Computational Electromagnetic Society, No.3, 1995) as the primary reference. The first example I tried was a square plate of unit length at varying frequencies and I tried to reproduce figures 7a and 7b in the paper, which reported the inverse condition number.
Obviously my initial attempts produced the same results for RWG basis functions. But not for loop-tree. After reading more, especially papers by Burton and Kashyap (ACES Journal 1995) and Eibert (IEEE AP Magazine, June 2004), I realized that a diagonal scaling normalizing the diagonal entries to unit magnitude would produce better condition numbers. That worked. I constructed a new matrix B = D*A*D, where A being the original impedance matrix, and D = sqrt( 1 ./ abs(diag(A))) [using MATLAB notations]. That produced the nearly constant condition number at low frequencies.
But there was still a problem. Wu et. al did not mention anything about a scaling in their paper. Nothing. Does it mean I had a bug in the code? Finally, I wrote to Dr. Allen Glisson explaining this conundrum. He wrote a detailed reply. Since it was such a long time ago (I was still playing cricket and riding a clunky 10-speed BSA around Bangalore at that time), naturally, he could not recall if they were doing any scaling. He was leaning to the possibility of not using any scaling. However, he did suggest a few things that I could check with my code. That was very helpful. I checked everything again, and this time it became clearer that it is almost impossible to get good condition number at low frequencies without some sort of scaling. This is because the magnitudes of the rows/columns corresponding loop bases are many orders smaller than those corresponding to standard RWG functions. So, numerically these rows and columns can be almost zero, especially in single precision computations, leading to very large (10^14 or more) condition numbers at extremely low frequencies.
I then contacted Dr. Thomas Eibert, mentioned above, asking his experience with this. This was on Friday (06/11/2010). During the weekend, I read the paper by Zhao and Chew (IEEE TAP, Oct 2000). This was illuminating. Jun-Sheng discussed a frequency scaling (Section IIa, equation (8)) precisely to handle the issue of small numerical values corresponding to the loop bases. That gave me the first hint that there may not be any bug in the code after all.
I quickly tested the scaling for an 8×8 system by hand in Octave. That produced a very good condition number. Then I proceeded to implement it in the actual code and everything worked like a charm. The final scaling I used was this:
This wasn’t the end of the story. On Tuesday (06/15), Dr. Eibert replied that in his opinion “scaling is mandatory in order to achieve a good condition number, even if we were able to achieve a perfect Hodge decomposition.”
Shortly afterwards, I got another mail from Dr. Glisson. He went back to the code and discovered that there was a comment that he had missed earlier referring to a scaling suggested in a paper by Mautz and Harrington (1984) and “so it appears that he did do the frequency scaling.” (The ‘he’ Dr. Glisson refers is the student who did the work.)
Thus solved the mystery of nice condition numbers.
But what amazed and humbled me at the same time was Dr. Glisson’s act of going out of his way and digging through a 15 year old code to find out what was going on. That is, in my view, a great example of scientific integrity and commitment towards one’s work that is to be emulated. It is no surprise, then, that it is people like Dr. Glisson who have consistently set and raised the standards of research in computational EM.
Dr. Eibert’s insights were also extremely valuable in helping me resolve this issue. Thank you both.
Tags: EFIE, loop/tree, low frequency