EinsteinInitialData/TwoPunctures segfaults and miscalcs with SP, maybe DP too.

Issue #1894 closed
David Koppelman created an issue

When using EinsteinInitialData/TwoPunctures with CCTK_REAL set to single-precision memory-access errors and incorrect results were encountered. The incorrect results could potentially affect double-precision code too.

The memory access errors occur due to double / CCTK_REAL confusion.

The incorrect results are due to rounding errors in PunctTaylorExpandAtArbitPosition that push an acos argument outside of [-1,1]. Though this was encountered with single-precision code, perhaps it could occur with DP code too.

The attached patch fixes the type mismatch and clamps the acos argument to [-1,1].

Keyword:

Comments (8)

  1. Roland Haas
    • changed status to open
    • removed comment

    Looks mostly ok to me. Minor issues that should be addressed:

    • layout of whitespace in the clamp routine should mimic what TwoPunctures' other code uses
    • I would rename the routine clamptoone to make it clear to which range it clamps
    • I would add a comment to the clamp on tan() stating that yes indeed the result of this compuation (ie B) is expected to be in the range of [-1,1] which tan is not normally guaranteed to be
  2. David Koppelman reporter
    • removed comment

    I've addressed the comments and updated the patch.

    clamp -> clamp_pm_one. I've also harmonized the formatting with surrounding code.

    I've added comments explaining why B should be in [-1,1].

    And, just to be safe I added assertions checking for NaNs to PunctTaylorExpandAtArbitPosition and also to PunctEvalAtArbitPositionFast since there are several functions, including square roots, that might cause trouble.

  3. David Koppelman reporter
    • changed status to open
    • removed comment

    I'm not sure whether I'm the one that should set the review flag.

  4. Roland Haas
    • changed status to open
    • removed comment

    Yes :-). Review basically means that you think it is now ready for someone else to look at. Looks good to me. If you can: please apply. Otherwise I will apply it for you (and you'll be listed as author in git).

  5. David Koppelman reporter
    • removed comment

    AFAIK I don't have check-in privileges on the ET repo, so please check it in for me.

  6. Log in to comment