search

adda-team / adda

ADDA - light scattering simulator based on the discrete dipole approximation
GNU General Public License v3.0
93 stars 56 forks source link

Optimization of integration of Green's tensor #252

Open myurkin opened 5 years ago

myurkin commented 5 years ago

There are two problems with current implementation of IGT. It is slow in general, and unstable for distances much smaller than the wavelength (static limit, #36 ). The latter causes errors IFAIL, which indicate either the imperfection of the convergence test or a true error (inaccuracy), not clear which one. Anyway, it is caused by strong variation (and large absolute values) of the Green's tensor in the static limit.

IGT_SO (#37 ) partly addresses those issues, but it is inherently limited to the cubical grid, while general rectangular dipoles (#196 ) makes tabulation impractical.

The new idea stems from the fact that integral of G can be decomposed as integral of G-Gst + integral of Gst. The former is similar to M-term and decreases with kd->0 (the order of smallness need to be determined). And the latter can be brought down to surface integrals (as in L term), where the integrand is bounded (solid-angle type). Maybe it can even be further simplified (e.g., to contour integrals). So we may build IGT_SO formulation without tables - numerical integration will remain only for Gst and in much simpler form, while G-Gst would be considered accurately up to (kd)^2.

Alternative (and may be equivalent) approach is to start with current IGT_SO formulae and to numerically calculate remaining (non-oscillating) integrals, by first simplifying them as much as possible (reducing dimensionality, etc.).

Or we may reduce dimensionality in the original integral of G - see Polimeridis A.G., Fernandez Villena J., Daniel L., and White J.K. Robust J-EFVIE solvers based on purely surface integrals, 2013 International Conference on Electromagnetics in Advanced Applications (ICEAA), pp. 379–381.

myurkin commented 5 years ago

In terms of unstable behavior, trivial renormalization in file propaesplibreintadda.f can help. Now it computes very small value and then divides by the dipole volume.

myurkin commented 2 years ago

299 constitutes first working implementation, which should be sufficient for most purposes in ADDA. However, we still want to improve it further, including making it more robust for arbitrary distance (not only between two voxels).