Merge pull request #470 from bluescarni/pr/erfa_fix

Fix for thread safety issue in erfa

CMakeLists.txt

@@ -11,7 +11,7 @@ if(NOT CMAKE_BUILD_TYPE)
 	FORCE)
 endif()
 
-project(heyoka VERSION 7.2.0 LANGUAGES CXX C)
+project(heyoka VERSION 7.2.1 LANGUAGES CXX C)
 
 list(APPEND CMAKE_MODULE_PATH "${CMAKE_CURRENT_SOURCE_DIR}/cmake" "${CMAKE_CURRENT_SOURCE_DIR}/cmake/yacma")
 

doc/changelog.rst

@@ -1,6 +1,16 @@
 Changelog
 =========
 
+7.2.1 (2025-01-07)
+------------------
+
+Fix
+~~~
+
+- Work around a thread-safety issue in the erfa routines used
+  to convert between UTC and TAI Julian dates
+  (`#470 <https://github.com/bluescarni/heyoka/pull/470>`__).
+
 7.2.0 (2025-01-02)
 ------------------
 

src/detail/erfa/.clang-format

@@ -0,0 +1 @@
+DisableFormat: true

src/detail/erfa/dat.c

@@ -199,7 +199,26 @@ int eraDat(int iy, int im, int id, double fd, double *deltat)
    int NDAT;
    eraLEAPSECOND *changes;
 
-   NDAT = eraDatini(_changes, _NDAT, &changes);
+   // NOTE: erfa supports changing at runtime the leap seconds table.
+   // This is accomplished with the help of a couple of global variables
+   // that can be set by the user via eraSetLeapSeconds(). If the user has not
+   // set a custom leap seconds table, the line
+   //
+   // eraDatini(_changes, _NDAT, &changes);
+   //
+   // will set these global variables to _NDAT and _changes (i.e.,
+   // the builtin leap seconds table defined a few lines earlier).
+   //
+   // The problem with this is that access to these global variables is not
+   // thread safe, and invoking the current function from multiple threads
+   // at the same time will result in data races.
+   //
+   // Hence, we remove the invocation of eraDatini() in favour of always using
+   // the builtin leap seconds table.
+   //
+   //NDAT = eraDatini(_changes, _NDAT, &changes);
+   NDAT = _NDAT;
+   changes = (eraLEAPSECOND *)_changes;
 
 /* Miscellaneous local variables */
    int j, i, m;

src/model/sgp4.cpp

@@ -915,10 +915,10 @@ void sgp4_propagator<T>::operator()(out_2d out, in_1d<date> dates)
     using tms_vec_size_t = decltype(tms_vec.size());
     tms_vec.resize(boost::numeric_cast<tms_vec_size_t>(dates.extent(0)));
 
-    // NOTE: we have a rough estimate of ~200 flops for a single execution of sgp4_date_to_tdelta().
+    // NOTE: we have a rough estimate of ~500 flops for a single execution of sgp4_date_to_tdelta().
     // Taking the usual figure of ~10'000 clock cycles as minimum threshold to parallelise, we enable
-    // parallelisation if we have 50 satellites or more.
-    if (n_sats >= 50u) {
+    // parallelisation if we have 20 satellites or more.
+    if (n_sats >= 20u) {
         oneapi::tbb::parallel_for(oneapi::tbb::blocked_range<tms_vec_size_t>(0, tms_vec.size()),
                                   [&tms_vec, dates, this, n_sats](const auto &range) {
                                       for (auto i = range.begin(); i != range.end(); ++i) {
@@ -1033,10 +1033,10 @@ void sgp4_propagator<T>::operator()(out_3d out, in_2d<date> dates)
     using tms_vec_size_t = decltype(tms_vec.size());
     tms_vec.resize(boost::safe_numerics::safe<tms_vec_size_t>(n_evals) * dates.extent(1));
 
-    // NOTE: we have a rough estimate of ~200 flops for a single execution of sgp4_date_to_tdelta().
+    // NOTE: we have a rough estimate of ~500 flops for a single execution of sgp4_date_to_tdelta().
     // Taking the usual figure of ~10'000 clock cycles as minimum threshold to parallelise, we enable
-    // parallelisation if tms_vec.size() (i.e., n_evals * n_sats) is 50 or more.
-    if (tms_vec.size() >= 50u) {
+    // parallelisation if tms_vec.size() (i.e., n_evals * n_sats) is 20 or more.
+    if (tms_vec.size() >= 20u) {
         oneapi::tbb::parallel_for(
             oneapi::tbb::blocked_range2d<std::size_t>(0, dates.extent(0), 0, dates.extent(1)),
             [&tms_vec, dates, this, n_sats](const auto &range) {