Остановить интеграцию с odeint, используемым с тягой

Question

Я пытаюсь интегрировать систему ODE с библиотекой odeint и параллельно параллельно на множестве точек (это означает то же ODE со многими разными начальными условиями). В частности, я использую алгоритм адаптивного шага размера runge_kutta_dopri5. Есть некоторые моменты, в которых алгоритм терпит неудачу, уменьшая размер шага и чрезвычайно замедляет весь процесс интеграции.

Есть ли способ остановить процесс интеграции только для некоторых точек набора, если они не проходят определенный тест?

В моем конкретном случае, поскольку я интегрирую проблему гравитации, я хотел бы остановить интеграцию, когда точка приближается к аттрактору, поэтому расстояние меньше определенного предела.

В последовательных вычислениях я думаю, что это может быть выполнено с помощью настраиваемого цикла while с функцией stepper.try_step, как показано более или менее в идее за this question.

Как это можно выполнить при параллельных вычислениях с тягой?

Спасибо.

Answer 1

Как уже упоминалось, прямого решения этой проблемы нет.

Одним из возможных решений является создание

1. вектора Интс отчетности, который ОДА уже остановился. Итак, если 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 к наблюдателю.

У нас также есть экспериментальная и грязная ветка, где управление размером шага выполняется непосредственно на графическом процессоре и контролирует каждый ОДУ отдельно. Это действительно мощный и высокопроизводительный. К сожалению, эта функция не может быть легко интегрирована в основную ветвь, так как многие из __device__ __host__ спецификаторов должны быть добавлены к способам odeint в. Если хотите, вы можете проверить ветку cuda_controlled_stepper в репозитории odeint и/или отправить мне сообщение для получения дополнительных инструкций.

Answer 2

Я думаю, что это довольно сложная проблема, которая будет реализована на GPU с тягой.

Однажды я сделал аналогичную симуляцию, где мне пришлось интегрировать многие начальные условия одной и той же системы, но каждая интеграция прекращалась после разных шагов. Не только небольшие изменения, но и порядки величин, например, между 1000 и 10^6 ступенями. Я написал распараллеливание для этого, используя OpenMP, который работал на 48 ядрах. То, что я делал, было очень простым для OpenMP-распараллеливания: всякий раз, когда заканчивается одно начальное условие, начинается следующее. Это достаточно эффективно, если общее число траекторий намного больше, чем параллельные потоки. В принципе, вы можете реализовать его так же и на GPU. Как только одна траектория заканчивается, вы заменяете ее новым начальным условием. Конечно, вам нужно делать некоторые книги, особенно если ваша система зависит от времени. Но в целом это может сработать, я думаю.