CUDA: траектория Рунге-Кутты на каждом потоке графического процессора

Question

Резюме: Как избежать потери производительности, вызванной различной рабочей нагрузкой для разных потоков? (Ядро с циклом while на каждом потоке)

Проблема: Я хочу решить траектории частиц (описанные дифференциальным уравнением 2-го порядка), используя Рунге-Кутту для многих различных начальных условий. Траектории обычно имеют различную длину (каждая траектория заканчивается, когда частица попадает в какую-то цель). Кроме того, чтобы обеспечить числовую стабильность, размер ступени Рунге-Кутты устанавливается адаптивно. Это приводит к двум вложенным циклам while с неизвестным числом итераций (см. Последовательный пример ниже).

Я хочу реализовать процедуру Runge-Kutta для запуска на графическом процессоре с CUDA / C ++. Траектории не зависят друг от друга, поэтому в качестве первого подхода я просто буду распараллеливать различные начальные условия, чтобы каждый поток соответствовал уникальной траектории. Когда поток завершен с траекторией частиц, я хочу, чтобы он начинался с новой.

Если я правильно понимаю, однако, неизвестная длина каждого цикла while (траектория частиц) означает, что разные потоки будут получать разный объем работы , что может привести к серьезной потере производительности на GPU.

Вопрос: Можно ли простым способом преодолеть потери производительности, вызванные различной рабочей нагрузкой для разных потоков? Например, установить для каждого размера деформации только 1, чтобы каждый поток (деформация) мог работать независимо? r приведет ли это к другим потерям производительности (например, нет чтения объединенной памяти)?

Серийный псевдокод :

// Solve a particle trajectory for each inital condition
// N_trajectories: much larger than 1e6
for( int t_i = 0; t_i < N_trajectories; ++t_i )
{
    // Set start coordinates
    double x = x_init[p_i];
    double y = y_init[p_i];
    double vx = vx_init[p_i];
    double vy = vy_init[p_i];
    double stepsize = ...;
    double tolerance = ...;
    ...

    // Solve Runge-Kutta trajectory until convergence
    int converged = 0;
    while ( !converged )
    {
        // Do a Runge-Kutta step, if step-size too large then decrease it
        int goodStepSize = 0
        while( !goodStepSize )
        {
            // Update x, y, vx, vy
            double error = doRungeKutta(x, y, vx, vy, stepsize);

            if( error < tolerance )
                goodStepSize = 1;
            else
                stepsize *= 0.5;
        }

        if( (abs(x-x_final) < epsilon) && (abs(y-y_final) < epsilon) )
            converged = 1;
    }
}

Краткий тест моего кода показывает, что внутренний цикл while выполняется 2-4 раза в 99% всех случаев и> 10 раз в 1% всех случаев, прежде чем будет найден удовлетворительный размер шага Рунге-Кутты.

Параллельный псевдокод :

int tpb = 64;
int bpg = (N_trajectories + tpb-1) / tpb;
RungeKuttaKernel<<<bpg, tpb>>>( ... );

__global__ void RungeKuttaKernel( ... )
{
    int idx = ...;

    // Set start coordinates
    double x = x_init[idx];
    ...

    while ( !converged )
    {
        ...

        while( !goodStepSize )
        {
            double error = doRungeKutta( ... );
            ...
        }

        ...
    }
}

Answer 1

Я попытаюсь ответить на вопрос сам, пока кто-нибудь не найдет лучшего решения.

Подводные камни с прямым портированием серийного кода: Два цикла while приведут к значительному расхождению веток и потере производительности. Внешний цикл - это «полная» траектория, а внутренний цикл - это один шаг Рунге-Кутты с адаптивной коррекцией размера шага. Внутренний цикл: Если мы попытаемся решить Рунге-Кутту со слишком большим размером шага, тогда ошибка аппроксимации будет слишком большой, и нам нужно будет повторить шаг с меньшим размером шага, пока ошибка не станет меньше, чем наша терпимость. Это означает, что потокам, которым требуется очень мало итераций, чтобы найти подходящий размер шага, придется ждать потоков, которым требуется больше итераций. Внешний цикл: это отражает, сколько успешных шагов Рунге-Кутты нам нужно, прежде чем траектория будет завершена. Различные траектории будут достигать своей цели за разное количество шагов. Нам всегда придется ждать траекторию с наибольшим количеством итераций, прежде чем мы полностью закончим.

Предлагаемый параллельный подход: Мы замечаем, что каждая итерация состоит из одного шага Рунге-Кутты. Ветвление происходит из-за того, что нам нужно либо уменьшить размер шага для следующей итерации, либо обновить коэффициенты Рунге-Кутты (например, положение / скорость), если размер шага был в порядке. Поэтому я предлагаю заменить два цикла while одним циклом for. Первым шагом цикла for является решение Runge-Kutta, за которым следует оператор if, чтобы проверить, был ли размер шага достаточно мал или обновляется ли позиция (и проверяется ли полная сходимость). Все потоки теперь будут обрабатывать только один шаг Рунге-Кутты за раз, и мы обмениваем низкую занятость (все потоки ждут потока, которому нужно наибольшее количество попыток найти правильный размер шага) за стоимость расхождения ветвления одного, если -заявление. В моем случае решение Runge-Kutta стоит дороже по сравнению с оценками этого оператора if, поэтому мы сделали улучшение. Теперь проблема заключается в том, чтобы установить соответствующее ограничение цикла for и отметить потоки, которые требуют дополнительной работы. Этот предел устанавливает верхнюю границу для самого длительного времени, когда готовый поток должен ждать других. Псевдо-код:

int N_trajectories = 1e6;
int trajectoryStepsPerKernel = 50;
thrust::device_vector<int> isConverged(N_trajectories, 0); // Set all trajectories to unconverged
int tpb = 64;
int bpg = (N_trajectories + tpb-1) / tpb;

// Run until all trajectories are converged
while ( vectorSum(isConverged) != N_trajectories )
{
    RungeKuttaKernel<<<bpg, tpb>>>( trajectoryStepsPerKernel, isConverged, ... );
    cudaDeviceSynchronize();
}

__global__ void RungeKuttaKernel( ... )
{
    int idx = ...;

    // Set start coordinates
    int converged = 0;
    double x = x_init[idx];
    ...

    for ( int i = 0; i < trajectoryStepsPerKernel; ++i )
    {
        double error = doRungeKutta( x_new, y_new, ... );

        if( error > tolerance )
        {
            stepsize *= 0.5;
        } else {
            converged = checkConvergence( x, x_new, y, y_new, ... );
            x = x_new;
            y = y_new;
            ...
        }
    }

    // Update start positions in case we need to continue on trajectory 
    isConverged[idx] = converged;
    x_init[idx] = x;
    y_init[idx] = y;
    ...
}