4

ODEのシステムをodeintライブラリと統合し、一連のポイントで並列に推力をかけようとしています（これは、多くの異なる初期条件を持つ同じODEを意味します）。特に、私は適応ステップサイズアルゴリズムrunge_kutta_dopri5を使用しています。アルゴリズムが失敗し、ステップサイズが小さくなり、統合プロセス全体が非常に遅くなるポイントがいくつかあります。

特定のテストに失敗した場合に、セットの一部のポイントでのみ統合プロセスを停止する方法はありますか？

私の特定のケースでは、重力の問題を積分しているので、ポイントがアトラクタに近づいたときに積分を停止したいので、距離は一定の制限未満です。

シリアルコンピューティングでは、この質問stepper.try_stepの背後にある考え方に多かれ少なかれ示されているように、これは関数を使用したカスタムwhileループによって実行できると思います。

推力を使用した並列コンピューティングでこれをどのように実行できますか？

ありがとう。

4

2 に答える 2

2

すでに述べたように、この問題に対する簡単な解決策はありません。

考えられる解決策の1つは、

  1. ODEがすでに停止していることを報告するintのベクトル。したがって、このベクトルの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;
};

このようなチェッカーは、制御されたルンゲクッタ法で使用できます。

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_driven_stepperブランチを確認するか、メッセージをドロップして詳細な手順を確認してください。

于 2013-03-27T21:07:41.583 に答える
2

これは、推力のあるGPUに実装するのは非常に難しい問題だと思います。

私はかつて同様のシミュレーションを行いました。同じシステムの多くの初期条件を統合する必要がありましたが、各統合は異なるステップ数の後に停止しました。わずかな変動だけでなく、1000〜10^6ステップのような桁違いの変動もあります。48コアで実行されるOpenMPを使用して、その並列化を作成しました。私が行ったことは、OpenMP並列化では非常に単純でした。ある初期条件が終了するたびに、次の条件が開始されます。軌道の総数が並列スレッドよりもはるかに多い限り、これはかなり効率的です。原則として、GPUにもこの方法で実装できます。1つの軌道が終了するとすぐに、それを新しい初期条件に置き換えます。もちろん、特にシステムが時間に依存している場合は、簿記を行う必要があります。しかし、一般的にこれはうまくいくかもしれないと私は思います。

于 2013-03-30T19:33:08.003 に答える