Merge pull request #124 from bluescarni/pr/ts_event_fix

Don't ignore the event equations when determining the timestep

doc/changelog.rst

@@ -4,9 +4,18 @@ Changelog
 0.6.1 (unreleased)
 ------------------
 
+Changes
+~~~~~~~
+
+- The event equations are now taken into account in the
+  determination of the adaptive timestep
+  (`#124 <https://github.com/bluescarni/heyoka/pull/124>`__).
+
 Fix
 ~~~
 
+- Fix an initialisation order issue in the event detection code
+  (`#124 <https://github.com/bluescarni/heyoka/pull/124>`__).
 - Fix an assertion misfiring in the event detection function
   (`#123 <https://github.com/bluescarni/heyoka/pull/123>`__).
 

doc/tut_adaptive.rst

@@ -128,6 +128,21 @@ then the timestep will be clamped to ``max_delta_t``:
 
 Note that the integration outcome is now ``taylor_outcome::time_limit``, instead of ``taylor_outcome::success``.
 
+Before moving on, we need to point out an important caveat when using the single
+step functions:
+
+.. warning::
+
+   If the exact solution of the ODE system can be expressed as a polynomial function
+   of time, the automatic timestep deduction algorithm may return a timestep of infinity.
+   This is the case, for instance, when integrating the rectilinear motion of a free
+   particle or the constant-gravity free-fall equation. In such cases, the step functions
+   should be called with a finite ``max_delta_t`` argument, in order to clamp the timestep
+   to a finite value.
+
+   Note that the ``propagate_*()`` functions (described :ref:`below <tlimited_prop>`)
+   are not affected by this issue.
+
 Accessing state and time
 ------------------------
 
@@ -144,6 +159,8 @@ Note that ``get_state()`` returns a const reference to the ``std::vector``
 holding the integrator state, while ``get_state_data()`` returns a naked pointer
 to the state data. Only ``get_state_data()`` can be used to mutate the state.
 
+.. _tlimited_prop:
+
 Time-limited propagation
 ------------------------
 

src/detail/event_detection.cpp

@@ -543,6 +543,16 @@ void taylor_detect_events_impl(std::vector<std::tuple<std::uint32_t, T, bool>> &
 
     assert(order >= 2u);
 
+    // NOTE: trigger the creation of the poly cache
+    // *before* triggering the creation of the other
+    // thread-local variables. This is important because
+    // the working list contains pwrap objects, which in turn
+    // interact with the poly cache. If we don't do this,
+    // the poly cache will be destroyed before the working
+    // list, and the destructors in the working list will
+    // interact with invalid memory.
+    [[maybe_unused]] auto *pc = &get_poly_cache<T>();
+
     // Fetch a reference to the list of isolating intervals.
     auto &isol = get_isol<T>();
 

src/taylor.cpp

@@ -1492,16 +1492,54 @@ llvm::Value *taylor_step_abs(llvm_state &s, llvm::Value *x_v)
 #endif
 }
 
+// Create a global read-only array containing the values in sv_funcs_dc, if any
+// (otherwise, the return value will be null). This is for use in the adaptive steppers
+// when employing compact mode.
+llvm::Value *taylor_c_make_sv_funcs_arr(llvm_state &s, const std::vector<std::uint32_t> &sv_funcs_dc)
+{
+    auto &builder = s.builder();
+
+    if (sv_funcs_dc.empty()) {
+        return nullptr;
+    } else {
+        auto *arr_type = llvm::ArrayType::get(llvm::Type::getInt32Ty(s.context()),
+                                              boost::numeric_cast<std::uint64_t>(sv_funcs_dc.size()));
+        std::vector<llvm::Constant *> sv_funcs_dc_const;
+        sv_funcs_dc_const.reserve(sv_funcs_dc.size());
+        for (auto idx : sv_funcs_dc) {
+            sv_funcs_dc_const.emplace_back(builder.getInt32(idx));
+        }
+        auto *sv_funcs_dc_arr = llvm::ConstantArray::get(arr_type, sv_funcs_dc_const);
+        auto *g_sv_funcs_dc = new llvm::GlobalVariable(s.module(), sv_funcs_dc_arr->getType(), true,
+                                                       llvm::GlobalVariable::InternalLinkage, sv_funcs_dc_arr);
+
+        // Get out a pointer to the beginning of the array.
+        return builder.CreateInBoundsGEP(g_sv_funcs_dc, {builder.getInt32(0), builder.getInt32(0)});
+    }
+}
+
 // Helper to generate the LLVM code to determine the timestep in an adaptive Taylor integrator,
 // following Jorba's prescription. diff_variant is the output of taylor_compute_jet(), and it contains
 // the jet of derivatives for the state variables and the sv_funcs. h_ptr is a pointer containing
-// the clamping values for the timesteps.
+// the clamping values for the timesteps. svf_ptr is a pointer to an LLVM array containing the
+// values in sv_funcs_dc.
 template <typename T>
-llvm::Value *taylor_determine_h(llvm_state &s,
-                                const std::variant<llvm::Value *, std::vector<llvm::Value *>> &diff_variant,
-                                const std::vector<std::uint32_t> &sv_funcs_dc, llvm::Value *h_ptr, std::uint32_t n_eq,
-                                std::uint32_t n_uvars, std::uint32_t order, std::uint32_t batch_size)
+llvm::Value *
+taylor_determine_h(llvm_state &s, const std::variant<llvm::Value *, std::vector<llvm::Value *>> &diff_variant,
+                   const std::vector<std::uint32_t> &sv_funcs_dc, llvm::Value *svf_ptr, llvm::Value *h_ptr,
+                   std::uint32_t n_eq, std::uint32_t n_uvars, std::uint32_t order, std::uint32_t batch_size)
 {
+    assert(batch_size != 0u);
+#if !defined(NDEBUG)
+    if (diff_variant.index() == 0u) {
+        // Compact mode.
+        assert(sv_funcs_dc.empty() == !svf_ptr);
+    } else {
+        // Non-compact mode.
+        assert(svf_ptr == nullptr);
+    }
+#endif
+
     using std::exp;
 
     auto &builder = s.builder();
@@ -1514,7 +1552,7 @@ llvm::Value *taylor_determine_h(llvm_state &s,
 
         auto *diff_arr = std::get<llvm::Value *>(diff_variant);
 
-        // These will end up containing the norm infinity of the state vector and the
+        // These will end up containing the norm infinity of the state vector + sv_funcs and the
         // norm infinity of the derivatives at orders order and order - 1.
         auto vec_t = to_llvm_vector_type<T>(context, batch_size);
         max_abs_state = builder.CreateAlloca(vec_t);
@@ -1548,6 +1586,30 @@ llvm::Value *taylor_determine_h(llvm_state &s,
                 max_abs_diff_om1);
         });
 
+        if (svf_ptr != nullptr) {
+            // Consider also the functions of state variables for
+            // the computation of the timestep.
+            llvm_loop_u32(
+                s, builder.getInt32(0), builder.getInt32(boost::numeric_cast<std::uint32_t>(sv_funcs_dc.size())),
+                [&](llvm::Value *arr_idx) {
+                    // Fetch the index value from the array.
+                    auto cur_idx = builder.CreateLoad(builder.CreateInBoundsGEP(svf_ptr, {arr_idx}));
+
+                    builder.CreateStore(
+                        taylor_step_maxabs(s, builder.CreateLoad(max_abs_state),
+                                           taylor_c_load_diff(s, diff_arr, n_uvars, builder.getInt32(0), cur_idx)),
+                        max_abs_state);
+                    builder.CreateStore(
+                        taylor_step_maxabs(s, builder.CreateLoad(max_abs_diff_o),
+                                           taylor_c_load_diff(s, diff_arr, n_uvars, builder.getInt32(order), cur_idx)),
+                        max_abs_diff_o);
+                    builder.CreateStore(taylor_step_maxabs(s, builder.CreateLoad(max_abs_diff_om1),
+                                                           taylor_c_load_diff(s, diff_arr, n_uvars,
+                                                                              builder.getInt32(order - 1u), cur_idx)),
+                                        max_abs_diff_om1);
+                });
+        }
+
         // Load the values for later use.
         max_abs_state = builder.CreateLoad(max_abs_state);
         max_abs_diff_o = builder.CreateLoad(max_abs_diff_o);
@@ -1567,7 +1629,10 @@ llvm::Value *taylor_determine_h(llvm_state &s,
         // order * (n_eq + n_sv_funcs).
         max_abs_diff_o = taylor_step_abs(s, diff_arr[order * (n_eq + n_sv_funcs)]);
         max_abs_diff_om1 = taylor_step_abs(s, diff_arr[(order - 1u) * (n_eq + n_sv_funcs)]);
-        for (std::uint32_t i = 1; i < n_eq; ++i) {
+        // NOTE: iterate up to n_eq + n_sv_funcs in order to
+        // consider also the functions of state variables for
+        // the computation of the timestep.
+        for (std::uint32_t i = 1; i < n_eq + n_sv_funcs; ++i) {
             max_abs_state = taylor_step_maxabs(s, max_abs_state, diff_arr[i]);
             max_abs_diff_o = taylor_step_maxabs(s, max_abs_diff_o, diff_arr[order * (n_eq + n_sv_funcs) + i]);
             max_abs_diff_om1
@@ -2637,42 +2702,25 @@ taylor_compute_jet(llvm_state &s, llvm::Value *order0, llvm::Value *par_ptr, llv
     }
 }
 
-// Create a global read-only array containing the values in sv_funcs_dc, if any
-// (otherwise, the return value will be null). This is for use in the adaptive steppers
-// when employing compact mode.
-llvm::Value *taylor_c_make_sv_funcs_arr(llvm_state &s, const std::vector<std::uint32_t> &sv_funcs_dc)
-{
-    auto &builder = s.builder();
-
-    if (sv_funcs_dc.empty()) {
-        return nullptr;
-    } else {
-        auto *arr_type = llvm::ArrayType::get(llvm::Type::getInt32Ty(s.context()),
-                                              boost::numeric_cast<std::uint64_t>(sv_funcs_dc.size()));
-        std::vector<llvm::Constant *> sv_funcs_dc_const;
-        sv_funcs_dc_const.reserve(sv_funcs_dc.size());
-        for (auto idx : sv_funcs_dc) {
-            sv_funcs_dc_const.emplace_back(builder.getInt32(idx));
-        }
-        auto *sv_funcs_dc_arr = llvm::ConstantArray::get(arr_type, sv_funcs_dc_const);
-        auto *g_sv_funcs_dc = new llvm::GlobalVariable(s.module(), sv_funcs_dc_arr->getType(), true,
-                                                       llvm::GlobalVariable::InternalLinkage, sv_funcs_dc_arr);
-
-        // Get out a pointer to the beginning of the array.
-        return builder.CreateInBoundsGEP(g_sv_funcs_dc, {builder.getInt32(0), builder.getInt32(0)});
-    }
-}
-
 // Helper to generate the LLVM code to store the Taylor coefficients of the state variables and
 // the sv funcs into an external array. The Taylor polynomials are stored in row-major order,
 // first the state variables and after the sv funcs. For use in the adaptive timestepper implementations.
 void taylor_write_tc(llvm_state &s, const std::variant<llvm::Value *, std::vector<llvm::Value *>> &diff_variant,
-                     const std::vector<std::uint32_t> &sv_funcs_dc, llvm::Value *tc_ptr, std::uint32_t n_eq,
-                     std::uint32_t n_uvars, std::uint32_t order, std::uint32_t batch_size)
+                     const std::vector<std::uint32_t> &sv_funcs_dc, llvm::Value *svf_ptr, llvm::Value *tc_ptr,
+                     std::uint32_t n_eq, std::uint32_t n_uvars, std::uint32_t order, std::uint32_t batch_size)
 {
-    auto &builder = s.builder();
-
     assert(batch_size != 0u);
+#if !defined(NDEBUG)
+    if (diff_variant.index() == 0u) {
+        // Compact mode.
+        assert(sv_funcs_dc.empty() == !svf_ptr);
+    } else {
+        // Non-compact mode.
+        assert(svf_ptr == nullptr);
+    }
+#endif
+
+    auto &builder = s.builder();
 
     // Convert to std::uint32_t for overflow checking and use below.
     const auto n_sv_funcs = boost::numeric_cast<std::uint32_t>(sv_funcs_dc.size());
@@ -2693,10 +2741,6 @@ void taylor_write_tc(llvm_state &s, const std::variant<llvm::Value *, std::vecto
 
         auto *diff_arr = std::get<llvm::Value *>(diff_variant);
 
-        // Create a global read-only array containing the values in sv_funcs_dc, if any
-        // (otherwise, svf_ptr will be null).
-        auto *svf_ptr = taylor_c_make_sv_funcs_arr(s, sv_funcs_dc);
-
         // Write out the Taylor coefficients for the state variables.
         llvm_loop_u32(s, builder.getInt32(0), builder.getInt32(n_eq), [&](llvm::Value *cur_var) {
             llvm_loop_u32(s, builder.getInt32(0), builder.CreateAdd(builder.getInt32(order), builder.getInt32(1)),
@@ -2864,18 +2908,22 @@ auto taylor_add_adaptive_step_with_events(llvm_state &s, const std::string &name
     assert(bb != nullptr);
     builder.SetInsertPoint(bb);
 
+    // Create a global read-only array containing the values in ev_dc, if there
+    // are any and we are in compact mode (otherwise, svf_ptr will be null).
+    auto *svf_ptr = compact_mode ? taylor_c_make_sv_funcs_arr(s, ev_dc) : nullptr;
+
     // Compute the jet of derivatives at the given order.
     auto diff_variant = taylor_compute_jet<T>(s, state_ptr, par_ptr, time_ptr, dc, ev_dc, n_eq, n_uvars, order,
                                               batch_size, compact_mode);
 
     // Determine the integration timestep.
-    auto h = taylor_determine_h<T>(s, diff_variant, ev_dc, h_ptr, n_eq, n_uvars, order, batch_size);
+    auto h = taylor_determine_h<T>(s, diff_variant, ev_dc, svf_ptr, h_ptr, n_eq, n_uvars, order, batch_size);
 
     // Store h to memory.
     store_vector_to_memory(builder, h_ptr, h);
 
     // Copy the jet of derivatives to jet_ptr.
-    taylor_write_tc(s, diff_variant, ev_dc, jet_ptr, n_eq, n_uvars, order, batch_size);
+    taylor_write_tc(s, diff_variant, ev_dc, svf_ptr, jet_ptr, n_eq, n_uvars, order, batch_size);
 
     // Create the return value.
     builder.CreateRetVoid();
@@ -3390,7 +3438,7 @@ auto taylor_add_adaptive_step(llvm_state &s, const std::string &name, U sys, T t
                                               compact_mode);
 
     // Determine the integration timestep.
-    auto h = taylor_determine_h<T>(s, diff_variant, sv_funcs_dc, h_ptr, n_eq, n_uvars, order, batch_size);
+    auto h = taylor_determine_h<T>(s, diff_variant, sv_funcs_dc, nullptr, h_ptr, n_eq, n_uvars, order, batch_size);
 
     // Evaluate the Taylor polynomials, producing the updated state of the system.
     auto new_state_var
@@ -3433,7 +3481,7 @@ auto taylor_add_adaptive_step(llvm_state &s, const std::string &name, U sys, T t
         [&]() {
             // tc_ptr is not null: copy the Taylor coefficients
             // for the state variables.
-            taylor_write_tc(s, diff_variant, {}, tc_ptr, n_eq, n_uvars, order, batch_size);
+            taylor_write_tc(s, diff_variant, {}, nullptr, tc_ptr, n_eq, n_uvars, order, batch_size);
         },
         [&]() {
             // Taylor coefficients were not requested,

test/CMakeLists.txt

@@ -91,3 +91,4 @@ ADD_HEYOKA_TESTCASE(cosh)
 ADD_HEYOKA_TESTCASE(sinh)
 ADD_HEYOKA_TESTCASE(tan)
 ADD_HEYOKA_TESTCASE(tanh)
+ADD_HEYOKA_TESTCASE(wavy_ramp)

test/taylor_adaptive.cpp

@@ -904,3 +904,16 @@ TEST_CASE("taylor scalar move")
     REQUIRE(tes_data == ta.get_t_events().data());
     REQUIRE(ntes_data == ta.get_nt_events().data());
 }
+
+// A test to make sure the propagate functions deal correctly
+// with trivial dynamics.
+TEST_CASE("propagate trivial")
+{
+    auto [x, v] = make_vars("x", "v");
+
+    auto ta = taylor_adaptive<double>{{prime(x) = v, prime(v) = 1_dbl}, {0, 0}};
+
+    REQUIRE(std::get<0>(ta.propagate_for(1.2)) == taylor_outcome::time_limit);
+    REQUIRE(std::get<0>(ta.propagate_until(2.3)) == taylor_outcome::time_limit);
+    REQUIRE(std::get<0>(ta.propagate_grid({3, 4, 5, 6, 7.})) == taylor_outcome::time_limit);
+}

test/taylor_adaptive_batch.cpp

@@ -275,3 +275,24 @@ TEST_CASE("propagate grid")
         }
     }
 }
+
+// A test to make sure the propagate functions deal correctly
+// with trivial dynamics.
+TEST_CASE("propagate trivial")
+{
+    auto [x, v] = make_vars("x", "v");
+
+    auto ta = taylor_adaptive_batch<double>{{prime(x) = v, prime(v) = 1_dbl}, {0, 0, 0.1, 0.1}, 2};
+
+    ta.propagate_for({1.2, 1.3});
+    REQUIRE(std::all_of(ta.get_propagate_res().begin(), ta.get_propagate_res().end(),
+                        [](const auto &t) { return std::get<0>(t) == taylor_outcome::time_limit; }));
+
+    ta.propagate_until({2.3, 4.5});
+    REQUIRE(std::all_of(ta.get_propagate_res().begin(), ta.get_propagate_res().end(),
+                        [](const auto &t) { return std::get<0>(t) == taylor_outcome::time_limit; }));
+
+    ta.propagate_grid({5, 6, 7, 8.});
+    REQUIRE(std::all_of(ta.get_propagate_res().begin(), ta.get_propagate_res().end(),
+                        [](const auto &t) { return std::get<0>(t) == taylor_outcome::time_limit; }));
+}

test/taylor_nt_event.cpp

@@ -278,9 +278,7 @@ TEST_CASE("taylor nte multizero")
                                                cur_time = t;
                                            })}};
 
-        for (auto i = 0; i < 20; ++i) {
-            REQUIRE(std::get<0>(ta.step()) == taylor_outcome::success);
-        }
+        REQUIRE(std::get<0>(ta.propagate_until(fp_t(4))) == taylor_outcome::time_limit);
 
         REQUIRE(counter == 12u);
 
@@ -350,9 +348,7 @@ TEST_CASE("taylor nte multizero")
                                           },
                                           event_direction::negative)}};
 
-        for (auto i = 0; i < 20; ++i) {
-            REQUIRE(std::get<0>(ta.step()) == taylor_outcome::success);
-        }
+        REQUIRE(std::get<0>(ta.propagate_until(fp_t(4))) == taylor_outcome::time_limit);
 
         REQUIRE(counter == 10u);
     };
@@ -440,9 +436,7 @@ TEST_CASE("taylor nte multizero negative timestep")
                                                cur_time = t;
                                            })}};
 
-        for (auto i = 0; i < 20; ++i) {
-            REQUIRE(std::get<0>(ta.step_backward()) == taylor_outcome::success);
-        }
+        REQUIRE(std::get<0>(ta.propagate_until(fp_t(-4))) == taylor_outcome::time_limit);
 
         REQUIRE(counter == 12u);
     };