Прекратить интеграцию с odeint используется с тяги
Я пытаюсь интегрировать систему ODE с библиотекой odeint и параллельно работать над набором точек (это означает один и тот же ODE со многими различными начальными условиями). В частности, я использую алгоритм адаптивного размера шага runge_kutta_dopri5.
Есть некоторые моменты, в которых алгоритм терпит неудачу, уменьшая размер шага и чрезвычайно замедляя весь процесс интеграции.
Есть ли способ остановить процесс интеграции только для некоторых точек набора, если они не проходят определенный тест?
В моем конкретном случае, поскольку я интегрирую гравитационную задачу, я хотел бы остановить интегрирование, когда точка приближается к аттрактору, поэтому расстояние меньше определенного предела.
В последовательных вычислениях я думаю, что это можно сделать с помощью пользовательского цикла while с stepper.try_step функция, как показано более или менее в идее позади этот вопрос.
Как это может быть выполнено в параллельных вычислениях с тягой?
Благодарю.
Решение
Как уже упоминалось, нет прямого решения этой проблемы.
Одним из возможных решений является предоставление
- Вектор целых, сообщающий, что ODE уже остановил. Таким образом, если i-й элемент этого вектора равен нулю, это означает, что i-й ODE все еще работает.
- Пользовательская проверка ошибок в runge_kutta_dopri5, которая устанавливает ошибку в 0, если текущая система уже остановлена. Таким образом, вы избегаете того, чтобы эта ошибка влияла на механизм управления размером шага, по крайней мере, она не влияла отрицательно на управление размером шага. Ниже приведен эскиз, как это можно реализовать
- Функция проверки, если интеграция остановлена. Это может быть в принципе помещено в наблюдателя.
Эскиз для пользовательского средства проверки ошибок может быть (он скопирован из 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;
};
Вы можете использовать такой шашки в контролируемом Rutge 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. Если хотите, вы можете проверить ветку cuda_controlled_stepper в репозитории odeint и / или отправить мне сообщение для получения дальнейших инструкций.
Другие решения
Я думаю, что это довольно сложная проблема для реализации на GPU с упором.
Однажды я провел похожую симуляцию, где мне пришлось интегрировать множество начальных условий одной и той же системы, но каждая интеграция остановилась после разного количества шагов. Не просто небольшие вариации, а порядки величин, например, от 1000 до 10 ^ 6 шагов. Я написал параллелизацию для этого, используя OpenMP, который работал на 48 ядрах. То, что я сделал, было очень просто для распараллеливания OpenMP: при выполнении одного начального условия запускается следующее. Это достаточно эффективно, если общее число траекторий намного больше, чем у параллельных нитей. В принципе, вы можете реализовать это и на GPU. Как только одна траектория заканчивается, вы заменяете ее новым начальным условием. Конечно, вы должны вести бухгалтерский учет, особенно если ваша система зависит от времени. Но в целом это может сработать, я думаю.
detector