reduce round-off errors in field integrals and fluxes with pairwise summation

[oskooi]

One approach, not currently implemented, to reducing round-off errors in the fields which are accumulated during the time stepping is to use pairwise summation. This would mainly involve rewriting parts of fields::loop_in_chunks in src/loop_in_chunks.cpp. It is possible that the recursion involved in the implementation would lead to an overall performance penalty (i.e., longer run times).

Since the implementation would likely be a major undertaking, perhaps it would be useful as a first step to investigate this accuracy-performance tradeoff using some simple benchmarking.

[stevengj]

The only sums we could implement this way are things like field integrals and fluxes. You can't implement timestepping using pairwise summation because timestepping is inherently in-order; you also can't use pairwise summation for Fourier-transform accumulation without incurring the storage cost of saving all of the fields, for the same reason — pairwise summation requires a random-access array.

Are roundoff errors in loop_in_chunks sums typically a limiting factor? Aren't they usually dwarfed by the discretization error? Even if you are doing a summation with 10^10 points, roundoff errors will usually only lose you about 5 significant digits compared to machine precision. And almost all loop_in_chunks sums will be much smaller than this, summing over a small subset of the domain like a flux region.

(Pairwise summation can, in principle, be implemented with negligible performance penalty by enlarging the base case. We did this successfully in the Julia implementation of sum, for example (JuliaLang/julia#4039). I'm more worried about the complexity of replacing LOOP_OVER_IVECS with a recursive version, for limited benefit in our use case.)

[stevengj]

(Fun fact: I wrote most of the Pairwise Summation Wikipedia article that you linked.)

[oskooi]

In cases where Meep is compiled with single precision floating point (i.e., on machines with limited memory, processing, and power capabilities), it's possible that roundoff errors in the time stepping could become comparable to the discretization error and thus affect overall accuracy.

For the near term, given the challenges of replacing LOOP_VER_IVECS, rewriting just the field integrals and fluxes would be a useful start.

[stevengj]

Again, you can't use pairwise summation for timestepping.

Even in single precision, you need to be computing field integrals over a lot of points (> 10^8, probably), before roundoff errors become comparable to discretization error...