4

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

如果某些测试失败，有没有办法仅针对集合中的某些点停止集成过程？

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

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

如何在带推力的并行计算中执行此操作？

谢谢。

4

2 回答 2

2

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

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

  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_stepper 分支和/或给我留言以获取更多说明。

于 2013-03-27T21:07:41.583 回答
2

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

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

于 2013-03-30T19:33:08.003 回答