Intrepid2: Tolerance for mapToReferenceFrameInitGuess

[dridzal]

Question

@trilinos/intrepid2

I'm confused about the tolerance for the Newton solve in the function CellTools::mapToReferenceFrameInitGuess(). It looks like the tolerance is 100 time epsilon(), where the double version of epsilon() is defined below. I believe that the output of this function is zero. Is that correct? If so, then the Newton solve would proceed for 15 iterations, even if the mapping is linear. Is that the intent? Also, even if the mapping is nonlinear, using 15 iterations could be quite inefficient.

  template<>
  KOKKOS_FORCEINLINE_FUNCTION
  double epsilon<double>() {
    typedef union {
      long long i64;
      double d64;
    } dbl_64;

    dbl_64 s;
    s.d64 = 1;
    s.i64++;
    return (s.i64 < 0 ? 1 - s.d64 : s.d64 - 1);
  }

[CamelliaDPG]

@dridzal, I just tried running CellTools test_04, which invokes this method, in the debugger on my local machine. I see tol being set to 2.2204460492503131E-14.

[mperego]

I guess we could allow the user to specify a different tolerance by adding an extra parameter (defaulted to tolerance), if you think it's important.

[dridzal]

@CamelliaDPG , @mperego : A reasonable non-zero tolerance would be fine, as the above experiment indicates. I'm just trying to understand how 1 - s.d64 results in 2.2204460492503131E-16.

[CamelliaDPG]

So, today I learned about unions in C++. What's happening here appears to be the exploitation of technically undefined (but commonly implemented) behavior; we write 1.0 to the double member. Then we increment the integer member, which uses the same memory, so you're basically adding the smallest increment you can to the 1.0 double value. Then you subtract that from 1.0, or 1.0 from that, whichever yields a positive result. And thus you get machine epsilon, the smallest increment you can add to 1.0.

Does that make sense? It's not my code; I'm not sure whose code it is.

[dridzal]

@CamelliaDPG : Thanks! That makes sense. I also didn't know this about unions. It may be worth adding a comment to this function, with your explanation.

[jhux2]

So, today I learned about unions in C++. What's happening here appears to be the exploitation of technically undefined (but commonly implemented) behavior; we write 1.0 to the double member. Then we increment the integer member, which uses the same memory, so you're basically adding the smallest increment you can to the 1.0 double value. Then you subtract that from 1.0, or 1.0 from that, whichever yields a positive result. And thus you get machine epsilon, the smallest increment you can add to 1.0.

That's clever.