c++ - 停止与推力一起使用的 odeint 集成

Question

我正在尝试将 ODE 系统与 odeint 库集成，并在一组点上并行推进(这意味着具有许多不同初始条件的相同 ODE)。特别是我正在使用自适应步长算法 runge_kutta_dopri5。 算法在某些点上失败，减少了步长并极大地减慢了整个集成过程。

是否有一种方法可以仅针对未通过特定测试的集合中的某些点停止集成过程？

在我的特殊情况下，因为我正在整合重力问题，所以我想在点接近吸引子时停止积分，因此距离小于某个限制。

在串行计算中，我认为这可以通过使用 stepper.try_step 函数的自定义 while 循环来执行，正如 this question 背后的想法或多或少所说明的那样.

这如何在具有推力的并行计算中执行？

谢谢。

Answer 1

如前所述，这个问题没有直接的解决方案。

一个可能的解决方案是提供一个

1. 报告 ODE 已经停止的整数 vector 。因此，如果该 vector 的第 i 个元素为零，则意味着第 i 个 ODE 仍在运行。
2. runge_kutta_dopri5 中的自定义错误检查器，如果当前系统已经停止，它将错误设置为 0。从而避免此错误影响步长控制机制，至少它不会对步长控制产生负面影响。请参阅下面的草图，了解如何实现这一点
3. 检查集成是否停止的函数。这原则上可以放在观察者中。

自定义错误检查器的草图可以是(它是从 default_error_checker 复制的):

class custom_error_checker
{
public:

    typedef double value_type;
    typedef thrust_algebra algebra_type;
    typedef thrust_operations operations_type;

    default_error_checker(
            const thrust::device_vector< int > &is_stopped ,
            value_type eps_abs = static_cast< value_type >( 1.0e-6 ) ,
            value_type eps_rel = static_cast< value_type >( 1.0e-6 ) ,
            value_type a_x = static_cast< value_type >( 1 ) ,
            value_type a_dxdt = static_cast< value_type >( 1 ) )
    : m_is_stopped( is_stopped ) ,
      m_eps_abs( eps_abs ) , m_eps_rel( eps_rel ) ,
      m_a_x( a_x ) , m_a_dxdt( a_dxdt )
    { }


    template< class State , class Deriv , class Err , class Time >
    value_type error( const State &x_old ,
                      const Deriv &dxdt_old ,
                      Err &x_err , Time dt ) const
    {
        return error( algebra_type() , x_old , dxdt_old , x_err , dt );
    }

    template< class State , class Deriv , class Err , class Time >
    value_type error( algebra_type &algebra ,
                      const State &x_old ,
                      const Deriv &dxdt_old ,
                      Err &x_err , Time dt ) const
    {
        // this overwrites x_err !
        algebra.for_each3( x_err , x_old , dxdt_old ,
            typename operations_type::template rel_error< value_type >(
                m_eps_abs , m_eps_rel , m_a_x ,
                m_a_dxdt * get_unit_value( dt ) ) );

        thrust::replace_if( x_err.begin() ,
                            x_err.end() ,
                            m_is_stopped.begin() ,
                            thrust::identity< double >()
                            0.0 );

        value_type res = algebra.reduce( x_err ,
                typename operations_type::template maximum< value_type >() ,
                    static_cast< value_type >( 0 ) );
        return res;
    }

private:

    thrust::device_vector< int > m_is_stopped;
    value_type m_eps_abs;
    value_type m_eps_rel;
    value_type m_a_x;
    value_type m_a_dxdt;
};

您可以通过以下方式在受控的 runge kutta 中使用这样的检查器

typedef runge_kutta_dopri5< double , state_type , double , state_type ,
    thrust_algebra , thrust_operation > stepper;
typedef controlled_runge_kutta< stepper ,
    custom_error_checker > controlled_stepper ;
typedef dense_output_runge_kutta< controlled_stepper > dense_output_stepper;

thrust::device_vector< int > is_stopped;
// initialize an is_finished
dense_output_stepper s(
    controlled_stepper(
        custom_controller( is_stopped , ... )));

最后，您需要一个函数来检查哪个 ODE 已经停止。让我们调用这个函数check_finish_status:

void check_finish_status(
    const state_type &x ,
    thrust::device_vector< int > &is_stopped );

可以在观察者内部调用该函数，需要将is_stopped传递给观察者。

我们还有一个实验性的脏分支，其中步长控制直接在 GPU 上运行并分别控制每个 ODE。这真的很强大而且性能很高。不幸的是，此功能无法轻松集成到主分支中，因为需要将许多 __device__ __host__ 说明符添加到 odeint 的方法中。如果愿意，您可以查看 odeint 存储库中的 cuda_controlled_stepper 分支和/或给我发消息以获取更多说明。