On Mon, Nov 15, 2010 at 11:46 AM, Roland Haas roland.haas@physics.gatech.edu wrote:
Meeting minutes:
Present were: Frank, Roland, Bruno, Christian, Ian, Erik
Testsuites
- AHFinderDirect/AEILocalInterp: Carpet/mercurial with full grid gives
identical results on 1 or 2 processors, PUGH with bitant grid does not. This points to PUGHInterp to be th culprit. However Ian sees differences when using Multipole and Carpet as well (for 32 processors). This points to the symmetry handling in AEILocalInterp being the culprit.
I have looked at the interpolation coordinates that AHFinderDirect requests from the interpolator, and how they change during the nonlinear elliptic solve. I have used PUGH, and compared runs on 1 and 2 processes.
The first three times the interpolator is called, the coordinates are almost identical, differing by round-off only (1e-16). This round-off error exists only in the domain of the second process. It seems to be cause by the floating-point representation of the lower domain boundary -0.4, which cannot be represented exactly, thus many coordinate values differ by round-off.
When the interpolator is called the fourth time, there are suddenly substantial differences of the order of 1e-7 everywhere in the domain. My best guess this is that the linear solver, which uses an ILU decomposition as preconditioner (I believe) and thus makes some discrete choices depending on its input, yields a different approximation because of the abovementioned floating-point differences.
This is also evident in the output of AHFinderDirect, if one increases the precision. The trial surfaces' Theta value (which depends on the interpolation result) begins to differ significantly in the same iteration in which the average surface radius (which depends on the interpolator input) begins to differ as well. That is, the interpolation result only differs if the interpolator input differs as well. (Note that there are two interpolator calls per iteration; this is used to evaluate the radial dependence of the Jacobian.)
grep 'proc 0/horizon 1' np?/Kerr*.out | awk '{ print $6,$7,$8; }'
1 process:
1 r_grid=1.8888370318372081 ||Theta||=0.58245864115033152 2 r_grid=1.8382954223001324 ||Theta||=0.41707953161957556 3 r_grid=1.7988326341366765 ||Theta||=0.19206874617879019 4 r_grid=1.7970186650813291 ||Theta||=0.020839031289570303 5 r_grid=1.8003073662257811 ||Theta||=0.00049321965075897409 6 r_grid=1.8003202072919120 ||Theta||=5.0676410161959451e-08 7 r_grid=1.8003202078276501 ||Theta||=3.8071282237872595e-13
2 processes:
1 r_grid=1.8888370318372081 ||Theta||=0.5824586411503323 2 r_grid=1.8382954098754503 ||Theta||=0.41707944199389968 3 r_grid=1.7988327249094958 ||Theta||=0.19206838061444106 4 r_grid=1.7970189739199420 ||Theta||=0.020838947150236913 5 r_grid=1.8003076389626189 ||Theta||=0.00049321768953301467 6 r_grid=1.8003204797819821 ||Theta||=5.0679888046606392e-08 7 r_grid=1.8003204803171862 ||Theta||=3.8437655835998896e-13
In conclusion, I could not find any problem with AHFinderDirect, AEILocalInterp, CarpetInterp, or PUGHInterp. It may in principle be good to have more thorough test suites for these, but this would be independent of the horizon finding problem discussed here.
-erik