For ODE integration, we suggest using the boost libraries whenever possible. More...

Detailed Description

For ODE integration, we suggest using the boost libraries whenever possible.

Details

Here we describe briefly the suggested setup to use a boost ODE integrator in SpECTRE. The boost utilities have a number of more elaborate features, including features that may result in simpler function calls than the specific recipes below. For full documentation, see the Boost odeint documentation: https://www.boost.org/doc/libs/1_72_0/libs/numeric/odeint/doc/html/boost_numeric_odeint/

The integration methods can be used largely as in the boost documentation. We summarize the salient points necessary to get an example running. The main steps in using a boost integrator are:

define the system
construct the stepper
initialize (might be performed implicitly, depending on the stepper)
perform steps
(compute dense output, for dense-output steppers)

For most cases the Dormand-Prince fifth order controlled, dense-output stepper is recommended. That stepper is used in the section "SpECTRE vectors or `std::array` thereof in boost integrators" below.

Fundamental types or <tt>std::array</tt> thereof in fixed-step integrators

Let us consider a simple oscillator system, which we'll declare as a lambda:

  const auto oscillatory_array_system = [](const std::array<double, 2>& state,
                                           std::array<double, 2>& dt_state,
                                           const double /*time*/) {
    dt_state[0] = state[1];
    dt_state[1] = -state[0];
  };

Note that the first argument is a const lvalue reference to the state type, std::array<double, 2>, the second argument is an lvalue reference for the time derivatives that is written to by the system function, and the final argument is the current time.

The construction and initialization of the stepper is simple:

boost::numeric::odeint::runge_kutta4<std::array<double, 2>>

fixed_array_stepper;

Finally, we can perform the steps and examine the output,

  const double fixed_step = 0.001;
  std::array<double, 2> fixed_step_array_state{{1.0, 0.0}};
  for (size_t i = 0; i < 1000; ++i) {
    CHECK(approx(fixed_step_array_state[0]) == cos(i * fixed_step));
    CHECK(approx(fixed_step_array_state[1]) == -sin(i * fixed_step));
    fixed_array_stepper.do_step(oscillatory_array_system,
                                fixed_step_array_state, i * fixed_step,
                                fixed_step);
  }

Fundamental types or <tt>std::array</tt> thereof in Dense-output integrators

The dense-output and controlled-step-size ODE integrators in boost comply with a somewhat different interface, as significantly more state information must be managed. The result is a somewhat simpler, but more opaque, interface to the user code.

Once again, we start by constructing the system we'd like to integrate,

  const auto quadratic_system = [](const double /*state*/, double& dt_state,
                                   const double time) {
    dt_state = 2.0 * time;
  };

The constructor for dense steppers takes optional arguments for tolerance, and dense steppers need to be initialized.

  boost::numeric::odeint::dense_output_runge_kutta<
      boost::numeric::odeint::controlled_runge_kutta<
          boost::numeric::odeint::runge_kutta_dopri5<double>>>
      dense_fundamental_stepper = boost::numeric::odeint::make_dense_output(
          1.0e-12, 1.0e-12,
          boost::numeric::odeint::runge_kutta_dopri5<double>{});
  dense_fundamental_stepper.initialize(1.0, 1.0, 0.01);

It is important to take note of the tolerance arguments, as the defaults are 1.0e-6, which is often looser than we want for calculations in SpECTRE.

We then perform the step supplying the system function, and can retrieve the dense output state with calc_state, which returns by reference by the second argument.

  std::pair<double, double> step_range =
      dense_fundamental_stepper.do_step(quadratic_system);
  double state_output;
  for(size_t i = 0; i < 50; ++i) {
    while(step_range.second < 1.0 + 0.01 * square(i) ) {
      step_range = dense_fundamental_stepper.do_step(quadratic_system);
    }
    dense_fundamental_stepper.calc_state(1.0 + 0.01 * square(i), state_output);
    CHECK(approx(state_output) == square(1.0 + 0.01 * square(i)));
  }

SpECTRE vectors or <tt>std::array</tt> thereof in boost integrators

The additional template specializations present in the OdeIntegration.hpp file ensure that the boost ODE methods work transparently for SpECTRE vectors as well. We'll run through a brief example to emphasize the functionality.

  const auto oscillatory_vector_system =
      [](const std::array<DataVector, 2>& state,
         std::array<DataVector, 2>& dt_state, const double /*time*/) {
        dt_state[0] = state[1];
        dt_state[1] = -state[0];
      };
  boost::numeric::odeint::dense_output_runge_kutta<
      boost::numeric::odeint::controlled_runge_kutta<
          boost::numeric::odeint::runge_kutta_dopri5<
              std::array<DataVector, 2>>>>
      dense_stepper =
      boost::numeric::odeint::make_dense_output(1.0e-14, 1.0e-14,
                            boost::numeric::odeint::runge_kutta_dopri5<
                                std::array<DataVector, 2>>{});
  const DataVector initial_vector = {{0.1, 0.2, 0.3, 0.4, 0.5}};
  dense_stepper.initialize(
      std::array<DataVector, 2>{{initial_vector, DataVector{5, 0.0}}}, 0.0,
      0.01);
  std::array<DataVector, 2> dense_output_vector_state{
      {DataVector{5}, DataVector{5}}};
  step_range = dense_stepper.do_step(oscillatory_vector_system);
  for(size_t i = 0; i < 100; ++i) {
    while (step_range.second < 0.01 * i) {
      step_range = dense_stepper.do_step(oscillatory_vector_system);
    }
    dense_stepper.calc_state(0.01 * i, dense_output_vector_state);
    for(size_t j = 0; j < 5; ++j) {
      CHECK(approx(dense_output_vector_state[0][j]) ==
            (0.1 * j + 0.1) * cos(0.01 * i));
      CHECK(approx(dense_output_vector_state[1][j]) ==
            -(0.1 * j + 0.1) * sin(0.01 * i));
    }
  }