From 064d970c11b5587d8041038d655bfd755881a1e2 Mon Sep 17 00:00:00 2001 From: Francesco Biscani Date: Wed, 7 Apr 2021 23:32:34 +0200 Subject: [PATCH 1/9] Initial code to consider event equations when determining the timestep. --- src/taylor.cpp | 88 +++++++++++++++++++++++++++++++++----------------- 1 file changed, 59 insertions(+), 29 deletions(-) diff --git a/src/taylor.cpp b/src/taylor.cpp index cb7c8aa47..3200178f1 100644 --- a/src/taylor.cpp +++ b/src/taylor.cpp @@ -1492,6 +1492,32 @@ 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 &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(sv_funcs_dc.size())); + std::vector 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 @@ -1548,6 +1574,36 @@ llvm::Value *taylor_determine_h(llvm_state &s, max_abs_diff_om1); }); + // Create a global read-only array containing the values in sv_funcs_dc, if any + // (otherwise, svf_ptr will be null). + // NOTE: in some cases we are creating another copy of this global + // array. See if we can re-arrange the API to avoid this duplication. + auto *svf_ptr = taylor_c_make_sv_funcs_arr(s, sv_funcs_dc); + + 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(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 +1623,7 @@ 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) { + 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,32 +2693,6 @@ 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 &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(sv_funcs_dc.size())); - std::vector 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. @@ -2670,10 +2700,10 @@ void taylor_write_tc(llvm_state &s, const std::variant &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) { - auto &builder = s.builder(); - assert(batch_size != 0u); + auto &builder = s.builder(); + // Convert to std::uint32_t for overflow checking and use below. const auto n_sv_funcs = boost::numeric_cast(sv_funcs_dc.size()); From 95e0a90a07a8975ef0bfd732e742c7bb0214da85 Mon Sep 17 00:00:00 2001 From: Francesco Biscani Date: Thu, 8 Apr 2021 10:10:32 +0200 Subject: [PATCH 2/9] Update internal docs. --- src/taylor.cpp | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/taylor.cpp b/src/taylor.cpp index 3200178f1..bdc8097c7 100644 --- a/src/taylor.cpp +++ b/src/taylor.cpp @@ -1540,7 +1540,7 @@ llvm::Value *taylor_determine_h(llvm_state &s, auto *diff_arr = std::get(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(context, batch_size); max_abs_state = builder.CreateAlloca(vec_t); @@ -1623,6 +1623,9 @@ 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)]); + // 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]); From b0a74ea42ae99b31548b4203140bd922a1da20d0 Mon Sep 17 00:00:00 2001 From: Francesco Biscani Date: Thu, 8 Apr 2021 11:25:22 +0200 Subject: [PATCH 3/9] Avoid creating the same array twice, and adapt test for changes. --- src/taylor.cpp | 54 +++++++++++++++++++++++++--------------- test/taylor_nt_event.cpp | 12 +++------ 2 files changed, 37 insertions(+), 29 deletions(-) diff --git a/src/taylor.cpp b/src/taylor.cpp index bdc8097c7..8c9771da2 100644 --- a/src/taylor.cpp +++ b/src/taylor.cpp @@ -1523,11 +1523,22 @@ llvm::Value *taylor_c_make_sv_funcs_arr(llvm_state &s, const std::vector -llvm::Value *taylor_determine_h(llvm_state &s, - const std::variant> &diff_variant, - const std::vector &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> &diff_variant, + const std::vector &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(); @@ -1574,12 +1585,6 @@ llvm::Value *taylor_determine_h(llvm_state &s, max_abs_diff_om1); }); - // Create a global read-only array containing the values in sv_funcs_dc, if any - // (otherwise, svf_ptr will be null). - // NOTE: in some cases we are creating another copy of this global - // array. See if we can re-arrange the API to avoid this duplication. - auto *svf_ptr = taylor_c_make_sv_funcs_arr(s, sv_funcs_dc); - if (svf_ptr != nullptr) { // Consider also the functions of state variables for // the computation of the timestep. @@ -2700,10 +2705,19 @@ taylor_compute_jet(llvm_state &s, llvm::Value *order0, llvm::Value *par_ptr, llv // 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> &diff_variant, - const std::vector &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 &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) { 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(); @@ -2726,10 +2740,6 @@ void taylor_write_tc(llvm_state &s, const std::variant(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)), @@ -2897,18 +2907,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(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(s, diff_variant, ev_dc, h_ptr, n_eq, n_uvars, order, batch_size); + auto h = taylor_determine_h(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(); @@ -3423,7 +3437,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(s, diff_variant, sv_funcs_dc, h_ptr, n_eq, n_uvars, order, batch_size); + auto h = taylor_determine_h(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 @@ -3466,7 +3480,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, diff --git a/test/taylor_nt_event.cpp b/test/taylor_nt_event.cpp index 232eec053..77a8aefc2 100644 --- a/test/taylor_nt_event.cpp +++ b/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); }; From da10b978371e58e9aea0157e7eff4e9a3abb847d Mon Sep 17 00:00:00 2001 From: Francesco Biscani Date: Thu, 8 Apr 2021 12:30:45 +0200 Subject: [PATCH 4/9] Fix order of initialisation issue in the event detection code. --- src/detail/event_detection.cpp | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/src/detail/event_detection.cpp b/src/detail/event_detection.cpp index 47b633614..2cfb7af2d 100644 --- a/src/detail/event_detection.cpp +++ b/src/detail/event_detection.cpp @@ -543,6 +543,16 @@ void taylor_detect_events_impl(std::vector> & 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(); + // Fetch a reference to the list of isolating intervals. auto &isol = get_isol(); From 69d9a785c15751c884217644dfd5b05958aeb80e Mon Sep 17 00:00:00 2001 From: Francesco Biscani Date: Thu, 8 Apr 2021 12:31:20 +0200 Subject: [PATCH 5/9] Add tests inspired by the wavy ramp. --- test/CMakeLists.txt | 1 + test/wavy_ramp.cpp | 98 +++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 99 insertions(+) create mode 100644 test/wavy_ramp.cpp diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index d215d804d..347d82a78 100644 --- a/test/CMakeLists.txt +++ b/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) diff --git a/test/wavy_ramp.cpp b/test/wavy_ramp.cpp new file mode 100644 index 000000000..fb06f2e28 --- /dev/null +++ b/test/wavy_ramp.cpp @@ -0,0 +1,98 @@ +// Copyright 2020, 2021 Francesco Biscani (bluescarni@gmail.com), Dario Izzo (dario.izzo@gmail.com) +// +// This file is part of the heyoka library. +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +// NOTE: a series of tests for bugs detected when implementing the wavy ramp example. + +#include +#include +#include + +#include + +#include +#include +#include + +#include "catch.hpp" +#include "test_utils.hpp" + +using namespace heyoka_test; +using namespace heyoka; + +const auto pi_const = boost::math::constants::pi(); + +auto make_wavy_ramp(bool check_event) +{ + auto cb_curve = [check_event](taylor_adaptive &ta, bool) { + auto x = ta.get_state()[0]; + auto y = ta.get_state()[1]; + + if (check_event) { + REQUIRE(y - (1. - x + 0.05 * std::cos(11 * pi_const * x)) == approximately(0.)); + } + + auto grad_x = 1 + 0.05 * 11 * pi_const * std::sin(11 * pi_const * x); + auto grad_y = 1.; + auto grad_norm = std::sqrt(grad_x * grad_x + grad_y * grad_y); + grad_x /= grad_norm; + grad_y /= grad_norm; + + auto vx = ta.get_state()[2]; + auto vy = ta.get_state()[3]; + + auto vproj = vx * grad_x + vy * grad_y; + + ta.get_state_data()[2] -= 1.8 * vproj * grad_x; + ta.get_state_data()[3] -= 1.8 * vproj * grad_y; + + return true; + }; + + auto cb_bottom = [check_event](taylor_adaptive &ta, bool) { + if (check_event) { + auto y = ta.get_state()[1]; + + REQUIRE(y == approximately(0.)); + } + + ta.get_state_data()[3] = -0.8 * ta.get_state_data()[3]; + + return true; + }; + + auto [x, y, vx, vy] = make_vars("x", "y", "vx", "vy"); + + auto eq_w_curve = y - (1. - x + 0.05 * cos(11 * pi_const * x)); + auto eq_bottom = y; + + auto ta + = taylor_adaptive({prime(x) = vx, prime(y) = vy, prime(vx) = 0_dbl, prime(vy) = -1_dbl}, {0, 1.2, 0, 0}, + kw::t_events = {t_event(eq_w_curve, kw::callback = cb_curve), + t_event(eq_bottom, kw::callback = cb_bottom)}); + + return ta; +} + +// This test case would trigger an assertion misfire +// in the event detection code. +TEST_CASE("assertion misfire") +{ + auto ta = make_wavy_ramp(false); + + REQUIRE(std::get<0>(ta.step()) == taylor_outcome::err_nf_state); + REQUIRE(std::get<0>(ta.step(10.)) == taylor_outcome::err_nf_state); +} + +// Make sure that the nonlinear event equation is correctly propagated +// at the desired precision. +TEST_CASE("accurate event propagation") +{ + auto ta = make_wavy_ramp(true); + + REQUIRE(std::get<0>(ta.propagate_until(10.)) == taylor_outcome::time_limit); +} From 16c74acc48b69b94983bbf42e6ac7af604d6b4e9 Mon Sep 17 00:00:00 2001 From: Francesco Biscani Date: Thu, 8 Apr 2021 13:49:06 +0200 Subject: [PATCH 6/9] A couple of test additions. --- test/taylor_adaptive.cpp | 13 +++++++++++++ test/taylor_adaptive_batch.cpp | 21 +++++++++++++++++++++ 2 files changed, 34 insertions(+) diff --git a/test/taylor_adaptive.cpp b/test/taylor_adaptive.cpp index 89f619ade..3e3ed6b81 100644 --- a/test/taylor_adaptive.cpp +++ b/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{{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); +} diff --git a/test/taylor_adaptive_batch.cpp b/test/taylor_adaptive_batch.cpp index d5626acf8..a6980fe43 100644 --- a/test/taylor_adaptive_batch.cpp +++ b/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{{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; })); +} From 4756ff4b6baa467551d3bcffbfc63478517a5644 Mon Sep 17 00:00:00 2001 From: Francesco Biscani Date: Thu, 8 Apr 2021 13:50:04 +0200 Subject: [PATCH 7/9] Update changelog. --- doc/changelog.rst | 2 ++ 1 file changed, 2 insertions(+) diff --git a/doc/changelog.rst b/doc/changelog.rst index f77ec463d..3f22f867a 100644 --- a/doc/changelog.rst +++ b/doc/changelog.rst @@ -7,6 +7,8 @@ Changelog Fix ~~~ +- Fix an initialisation order issue in the event detection code + (`#124 `__). - Fix an assertion misfiring in the event detection function (`#123 `__). From c09cbac6165ba20688209b01b4da2462f66519e5 Mon Sep 17 00:00:00 2001 From: Francesco Biscani Date: Thu, 8 Apr 2021 14:59:08 +0200 Subject: [PATCH 8/9] Minor doc addition. --- src/taylor.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/taylor.cpp b/src/taylor.cpp index 8c9771da2..ee372c484 100644 --- a/src/taylor.cpp +++ b/src/taylor.cpp @@ -1521,7 +1521,8 @@ llvm::Value *taylor_c_make_sv_funcs_arr(llvm_state &s, const std::vector llvm::Value * taylor_determine_h(llvm_state &s, const std::variant> &diff_variant, From 4137faf8528570733b815a3d5293dd2da913db2e Mon Sep 17 00:00:00 2001 From: Francesco Biscani Date: Thu, 8 Apr 2021 15:09:39 +0200 Subject: [PATCH 9/9] More docs, changelog update. --- doc/changelog.rst | 7 +++++++ doc/tut_adaptive.rst | 17 +++++++++++++++++ 2 files changed, 24 insertions(+) diff --git a/doc/changelog.rst b/doc/changelog.rst index 3f22f867a..25b2b0df5 100644 --- a/doc/changelog.rst +++ b/doc/changelog.rst @@ -4,6 +4,13 @@ Changelog 0.6.1 (unreleased) ------------------ +Changes +~~~~~~~ + +- The event equations are now taken into account in the + determination of the adaptive timestep + (`#124 `__). + Fix ~~~ diff --git a/doc/tut_adaptive.rst b/doc/tut_adaptive.rst index c3ee8a222..508db0861 100644 --- a/doc/tut_adaptive.rst +++ b/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 `) + 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 ------------------------