Я пытаюсь интегрировать систему ODE с библиотекой odeint и параллельно работать над набором точек (это означает один и тот же ODE со многими различными начальными условиями). В частности, я использую алгоритм адаптивного размера шага runge_kutta_dopri5.
Есть некоторые моменты, в которых алгоритм терпит неудачу, уменьшая размер шага и чрезвычайно замедляя весь процесс интеграции.
Есть ли способ остановить процесс интеграции только для некоторых точек набора, если они не проходят определенный тест?
В моем конкретном случае, поскольку я интегрирую гравитационную задачу, я хотел бы остановить интегрирование, когда точка приближается к аттрактору, поэтому расстояние меньше определенного предела.
В последовательных вычислениях я думаю, что это можно сделать с помощью пользовательского цикла while с stepper.try_step
функция, как показано более или менее в идее позади этот вопрос.
Как это может быть выполнено в параллельных вычислениях с тягой?
Благодарю.
Как уже упоминалось, нет прямого решения этой проблемы.
Одним из возможных решений является предоставление
Эскиз для пользовательского средства проверки ошибок может быть (он скопирован из 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. Как только одна траектория заканчивается, вы заменяете ее новым начальным условием. Конечно, вы должны вести бухгалтерский учет, особенно если ваша система зависит от времени. Но в целом это может сработать, я думаю.