 removed comment
If I set REAL_PRECISION=4 in my option list (single precision), Carpet will not run. This is because it performs selftests on startup of the prolongation operators, and my guess is that single precision is insufficient for the 9th order operators.
WARNING[L1,P0] (CarpetLib): Error in prolongate_3d_rf2::coeffs_3d_rf2 RT=CCTK_REAL4 ORDER=9 n=6 y0=0, res=9.53674e07 WARNING[L1,P0] (CarpetLib): Error in prolongate_3d_rf2::coeffs_3d_rf2 RT=CCTK_REAL4 ORDER=9 n=6 y0=0, res=9.53674e07 WARNING[L1,P0] (CarpetLib): Error in prolongate_3d_rf2::coeffs_3d_rf2 RT=CCTK_REAL4 ORDER=9 n=8 y0=0, res=1.52588e05 WARNING[L1,P0] (CarpetLib): Error in prolongate_3d_rf2::coeffs_3d_rf2 RT=CCTK_REAL4 ORDER=9 n=8 y0=0, res=1.52588e05 WARNING level 0 in thorn CarpetLib processor 0 host MacBook2.local (line 98 of /Users/ian/Cactus/EinsteinToolkit/arrangements/Carpet/CarpetLib/src/prolongate_3d_rf2.cc): > Aborting. WARNING level 0 in thorn CarpetLib processor 0 host MacBook2.local (line 98 of /Users/ian/Cactus/EinsteinToolkit/arrangements/Carpet/CarpetLib/src/prolongate_3d_rf2.cc): > Aborting.
I propose that the final warning should be changed to level 1 unless that order of prolongation has actually be requested for the run. In my case, I'm doing unigrid and don't care about prolongation. I don't know how to test the type (what type of object is "RT"?) in this template function.
Keyword:
Comments (12)


reporter  removed comment
I looked up numeric_limits, and according to http://www.cplusplus.com/reference/std/limits/numeric_limits, numeric_limits<T>::min() is of type T and is the minimum positive value that T can take. i.e. it would be 1.17549435e38F for single precision and 2.2250738585072014e308 for double precision. Do you mean numeric_limits<T>::epsilon() instead? This would be 1.19209290e07F for single precision and 2.2204460492503131e16 for double precision. eps^3/4^ would be correspondingly 6.4e6 and 1.8e12, which sound more like what we want. In the code, with single precision, I see:
INFO (CarpetLib): numeric_limits<RT>::min() = 1.1920928955078125e07 INFO (CarpetLib): good::abs (res  y0) = 1.52587890625e05
Is this what you mean by borderline? Should we use eps^4/7^? This gives (1.19209290e07)^(4/7)^ = 1.1e4 as the tolerance for single precision and 1.13e9 for double precision.

 removed comment
Right, it should be epsilon, not min.
4/7 sound like a very finetuned choice. What about 100*eps^(3/4)? Or maybe a lookup table, manually choosing 1e12 and 1e4? Or maybe we want 1e+3*eps, acknowledging that we lose 3 digits in our calculations.

reporter  removed comment
I agree it is finetuned. I like your last suggestion the best. That will make the test stricter than it currently is for double precision (2.2e13) and allows single precision to work (1.2e4). 1e13 seems a little small, but then we are getting to finetuning again...

 changed status to resolved
 removed comment
Applied.

 removed comment
This fix seems to have only been applied to the cellcentered version of the code. The original problem reported by Ian still exists. Can we apply the same fix to CarpetLib/src/prolongate_3d_rf2.cc (and possibly any other affected versions)?

 removed comment
Replying to [comment:6 barry.wardell]:
This fix seems to have only been applied to the cellcentered version of the code. The original problem reported by Ian still exists. Can we apply the same fix to CarpetLib/src/prolongate_3d_rf2.cc (and possibly any other affected versions)?
To answer my own question: no. The fix gets rid of the warning for 9th order, but a similar one appears for the 11th order operators:
WARNING[L1,P0] (CarpetLib): Error in prolongate_3d_rf2::coeffs_3d_rf2 RT=CCTK_REAL4 ORDER=11 n=10 y0=0, res=0.000732422 WARNING[L1,P0] (CarpetLib): Error in prolongate_3d_rf2::coeffs_3d_rf2 RT=CCTK_REAL4 ORDER=11 n=10 y0=0, res=0.000732422

 removed comment
At some point we just have to accept that higher order numerical methods don't make sense with single precision. I suggest to disable this operator for single precision; patches welcome.

 removed comment
How about this?

 removed comment
This seems more invasive than necessary. I think it suffices to disable the instantiation and the testing of the 11th order operators; the coefficient arrays can still be defined. However, your patch is also fine.
Please commit.

 removed comment
I agree. I have committed a slightly different patch which is less invasive.

 edited description
 changed status to closed
 Log in to comment
I don't think changing things to an L1 warning is a good solution. Either 9th order prolongation works with single precision (then we fix it), or it doesn't because of lack of precision (then we disable it). Or maybe we just have to lower our expectations with single precision.
For example, the statement
if (not (good::abs (res  y0) < 1.0e12))
may need adjusting. Replacing 1.0e12 by eps^(3/4) where eps is the precision of CCTK_REAL may be a way.
RT is a type; it is the real type corresponding to T, which may be a complex type. In C++, eps = numeric_limits<RT>::min(). Can you check whether eps^(3/4) works for single precision? 1.5e5 seems borderline for this.