停止与推力使用的odeint集成

Stop integration with odeint used with thrust

本文关键字：odeint 集成      更新时间：2023-10-16

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

是否存在一种方法，如果集合中的某些点未能通过某个测试，则仅对这些点停止积分过程?

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

在串行计算中，我认为这可以通过使用stepper.try_step函数的自定义while循环来执行，如本问题背后的思想所示。

如何在带推力的并行计算中实现这一点?

谢谢。

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

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

  1. 报告ODE已经停止的整型向量。因此，如果此向量的第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_step分支和/或给我留言以获取进一步的说明。

我认为这是一个很难在GPU上实现的问题。

我曾经做过一个类似的模拟，我必须对同一个系统的许多初始条件进行积分，但每次积分在不同的步数后停止。不仅仅是微小的变化，而是数量级的变化，比如在1000到10^6步之间。我用OpenMP编写了一个并行化程序，运行在48核上。我所做的对于OpenMP并行化非常简单:每当一个初始条件完成时，下一个初始条件就会启动。只要轨迹的总数比并行线程大得多，这是相当有效的。原则上，你也可以在GPU上实现这种方式。一旦一个轨迹结束，你用一个新的初始条件替换它。当然，你必须做一些记录，特别是如果你的系统是时间依赖的。但总的来说，我认为这是可行的

相关文章：