Задача N тел или как взорвать галактику не выходя из кухни / forpes.ru

Главная
Задача N тел или как взорвать галактику не выходя из кухни

Задача N тел или как взорвать галактику не выходя из кухни +165

18.03.2019 19:25

AndrewSu 96 21800 Источник

Не так давно я прочёл фантастический роман «Задача трёх тел» Лю Цысиня. В нём у одних инопланетян была проблема — они не умели, с достаточной для них точностью, вычислять траекторию своей родной планеты. В отличии от нас, они жили в системе из трёх звёзд, и от их взаимного расположения сильно зависела «погода» на планете — от испепеляющей жары до леденящего мороза. И я решил проверить, можем ли мы решать подобные задачи.

Физика явления

Для понимания задачи необходимо разобраться с физикой явления. В рамках классической теории сила притяжения двух тел определяется законом Ньютона:

$\vec{F}(\vec{r}_1, \vec{r}_2)=-G m_1 m_2 \frac{\vec{r}_1 - \vec{r}_2}{|\vec{r}_1 - \vec{r}_2|^3},$

где

$\vec{r}_1, \vec{r}_2$ — положение тел в пространстве,

$inline$ — массы тел,

$inline$ — гравитационная постоянная.
В системе из

$inline$ тел на каждое из них будет действовать сила притяжения от остальных, что выражается уравнением:

$\vec{F}_n=-G \sum_{k \neq n} m_n m_k \frac{\vec{r}_n - \vec{r}_k}{|\vec{r}_n - \vec{r}_k|^3}.$

Воспользовавшись вторым законом Ньютона запишем ускорение для каждой частицы:

$\vec{a}_n = \vec{F}_n/m_n = -G \sum_{k \neq n} m_k \frac{\vec{r}_n - \vec{r}_k}{|\vec{r}_n - \vec{r}_k|^3}.$

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

$\frac{\partial^2 \vec{r}_n }{\partial t^2} = f_n =-G \sum_{k \neq n} m_k \frac{\vec{r}_n - \vec{r}_k}{|\vec{r}_n - \vec{r}_k|^3}.$

Здесь важно заметить, что сложность вычисления функции

$inline$ равна

$inline$ и сильно возрастает с увеличением количества взаимодействующих тел.

Математика

Первым и простейшим методом решения дифференциальных уравнений является метод Эйлера, который предназначен для решения уравнений вида:

$\frac{dy}{dx}=f(x,y).$

При переходе в дискретную область получим:

$y_{i}=y_{i-1}+h f(x_{i-1},y_{i-1}),\quad i=1,2,3,\dots ,m,$

где

$inline$ — шаг интегрирования, а

$inline$ — число шагов интегрирования. Таким образом, ести нам необходимо произвести вычисление положения тел на момент времени

$inline$ , то нам следует сделать

$inline$ шагов интегрирования. Тут видна первая проблема — если

$inline$ велико, то нам нужно сделать большое количество шагов интегрирования.

Для применения метода Эйлера к нашей задаче её следует свести к системе первого порядка. Для этого введём дополнительную переменную — скорость частицы:

$\frac{\partial \vec{v}_n }{\partial t} = f_n,\\ \frac{\partial \vec{r}_n }{\partial t} = \vec{v}_n.$

Второй проблемой в решении систем дифференциальных уравнений является точность решения и её контроль. Точность можно повышать двумя способами: уменьшением шага интегрирования и выбором метода с более высоким порядком точности. Оба способа ведут к увеличению вычислительной сложности, но разными путями. Например, можно использовать классический метод Рунге-Кутты четвертого порядка, он требует четырёх вычислений функции

$inline$ на каждом шаге, но имеет порядок точности

$inline$ (для сравнения, метод Эйлера имеет порядок точности

$inline$ и требует одного вычисления

$inline$ ). Контроль точности решения также можно осуществлять несколькими способами: сравнить с аналитическим решением, решить разними методами или с разным шагом и сравнить результаты, контролировать сторонние параметры и ограничения, которым должно соответствовать решение.
Также, у каждого из этих методов есть свои недостатки. Аналитические решения могут отсутствовать, или, даже в большинстве случаев, и вовсе отсутствуют. Например, для нашей задачи

$inline$ тел аналитическое решение есть только при

$inline$ , но даже этого достаточно для тестирования точности методов. Решение задачи двумя методами или с разным шагом увеличивает время вычислений, но этот подход возможно применять практически для любой задачи. Ограничения есть не у каждой задачи, но для нашей они есть: на каждом шаге интегрирования мы можем контролировать выполнение законов сохранения. Этот подход тоже увеличивает сложность вычисления, но здесь есть из чего выбирать, вычисление суммы импульсов или моментов импульса всех частиц имеет сложность порядка

$inline$ , в то время как вычисление полной энергии системы имеет сложность порядка

$inline$

Заметка про вычисление полной энергии

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

$inline$ слагаемых. Сложность в том, что все слагаемые имеют сильо разный порядок, и даже при вычислениях с двойной точностью не удаётся вычислить это значение с точностью, достаточной для сравнения на разных шагах. Для преодоления этой проблемы пришлось применить суммирование по алгоритму Кэхэна.

Рис 1. Пример эллиптической траектории.

Рассмотрим простой случай движения спутника по эллиптической орбите вокруг Земли. При приближении спутника к Земле его скорость увеличивается, а при удалении от Земли уменьшается, соответственно, напрашивается возможность уменьшать шаг интегрирования по времени на зелёном участке орбиты, и увеличивать на красном и синем без изменения точности решения. Попробуем сравнить более детально.

Таблица 1. Исследуемые методы решения дифференциальных уравнений

Обозначене	Порядок	Описание
adams	1-5	Метод Адамса-Башфорта
euler	1	Метод Эйлера
rk4	4	Классический метод Рунге-Кутты
rkck	5	Метод Каша-Карпа
rkdp	5	Метод Дормана-Принса
rkdverk	6	Метод Вернера 1) p. 181
rkf	7	Метод Фельберга 1) p. 180
rkgl	6	Неявный метод Гаусса-Лежандра
rklc	4	Неявный метод Лобатто
trapeze	2	Метод трапеций

Для выбора наилучшего метода для нашей задачи произведём сравнение нескольких известных методов. Для этого смоделируем столкновение двух систем тел

$inline$ и измерим относительное изменение полной энергии, импульса и его момента в конце моделирования (максимальное время моделирования

$T_{max} = 2.5$ ). При этом мы будем варьировать шаг и параметры методов интегрирования и замерять количество вызовов функции

$inline$ , соответственно, те методы, которые при меньшем числе вызовов приведут к меньшим потерям, будем считать более приемлимыми.


a)	b)
Рис 2. Относительное изменение энергии a), импульса b), в конце моделирования системы $inline$ тел различными методами в зависимости от количества вычислений функции $inline$ с двойной точностью (double)

Из графиков рисунке 2 видно, что наилучшее соотношение количества вычисления функции

$inline$ и относительного изменения энергии системы тел у методов Адамса пятого порядка и Дормана-Принса. Также видно, что для всех методов с ростом числа вычисления

$inline$ увеличивается относительное изменение импульса системы. Для относительного изменения энергии это также заметно, но только для нескольких методов, которые смогли достичь порога

$dE/E_0 < 10^{-12}$ . Этот эффект можно связать уже не с погрешностью методов, а с погрешностью вычислений, и дальнейшее увеличение точности возможно только совместно с увеличением точности вычислений с плавающей точкой.


a)	b)
Рис 3. Относительное изменение энергии a), импульса b), в конце моделирования системы $inline$ тел различными методами в зависимости от количества вычислений функции $inline$ с четверной точностью (__float128)

Из рисунков 3a и 3b видно, что применение вычислений с четверной точностью позволяет снизить относительные потери энергии вплоть до

$10^{-23}$ , но нужно понимать, что время вычислений по сравнению с двойной точностью увеличивается на два порядка. Как и в случае с вычислениями с двойной точностью наилучшим соотношением точности и количества вычислений функции

$inline$ обладают методы Адамса пятого порядка и Дормана-Принса.

Методы Дормана-Принса и Вернера относятся к классу вложенных методов и позволяют одновременно вычислить два решения с высоким и низким порядком точности (для метода Дормана-Принса порядки 5 и 4, а для метода Вернера порядки 6 и 5). Если эти два решения сильно отличаются, то мы можем разбить текущий шаг интегрирования на более мелкие. Что позволяет нам динамически изменять шаг интегрирования и уменьшать его только на тех участках, где это требуется.

Сравним методы Дормана-Принса, Вернера и Адамса пятого порядка более детально, на более длительном интервале моделирования нашей системы (

$T_{max} = 300$ ).

a) Относительное изменение энергии, импульса и его момента в процессе моделирования

b) Количество вычислений функции $inline$ на интервале моделирования $\Delta T=0.3$

Рис 4. Зависимости относительного изменения физических параметров системы и количества вычислений функции $inline$ от времени моделирования методом Дормана-Принса с переменным шагом $h=10^{-5}\dots 10^{-9}$

a) Относительное изменение энергии, импульса и его момента в процессе моделирования

b) Количество вычислений функции $inline$ на интервале моделирования $\Delta T=0.3$

Рис 5. Зависимости относительного изменения физических параметров системы и количества вычислений функции $inline$ от времени моделирования методом Вернера с переменным шагом $h=10^{-5}\dots 10^{-9}$

Рис 6. Относительное изменение энергии, импульса и его момента в процессе моделирования методом Адамса-Башфорта пятого порядка с шагом $h=10^{-6}$

Рис 7. Зависимости количества вычислений функции $inline$ для методов Адамса пятого порядка, Дормана-Принса и Вернера от времени моделирования

Видно, что на начальном этапе эволюции нашей системы (

$inline$ ) все три метода показывают схожие характеристики, но на более поздних этапах в системе происходят некие события, в результате которых резко подскакивают ошибки в основных параметрах системы (полной энергии, импульса и его момента). Но методы Дормана-Принса и Вернера справляется с этими изменениями существенно лучше за счёт возможности уменьшать шаг интегрирования на «сложных» участках, в результате чего возрастает число вычислений функции

$inline$ , что видно на рисунках 4b и 5b, но общее число вычислений

$inline$ у вложенных методов меньше, чем у метода Адамса-Башфорта, что видно на рисунке 7.

Интересно, что происходило с системой в эти моменты

Видео 1. Моделирование системы из 512 тел. Метод Дормана-Принса. Динамический шаг $h=10^{-5}\dots 10^{-9}$

Видео демонстрирует, что до момента времени

$inline$ движение относительно спокойное, а после происходит столкновение центров «галактик», которое приводит к резкому изменению траекторий и сильному увеличению скоростей некоторых частиц. При этом для сохранения точности решения необходимо уменьшать шаг интегрирования. Вложенные методы могут сделать это автоматически, на графиках видно, что на некоторых участках эволюции системы шаг интегрирования был уменьшнен почти на два порядка с

$10^{-5}$ до

$h=10^{-7}$ . При использовании метода Адамса и такого шага на всём интервале эволюции системы мы не дождались бы решения.

Итог

Для решения лучше использовать вложеные методы, которые позволяют динамически контролировать шаг интегрирования, и уменьшать его только на «сложных» участках траектории.

Не стоит гнаться за методами самого высокого порядка. Даже при использовании типа данных 'double' они не достигают своих потенциальных возможностей, использование же типов данных с большей точностью сильно увеличивает время решения задачи.

CPU реализация

Теперь, когда выбор метода решения уравнений определён, попробуем разобраться с вычислением силы взаимодействия для каждой частицы. Получим двойной цикл по всем частицам:

Код реализации 'simple'

for(size_t body1 = 0; body1 < count; ++body1)
{
    const nbvertex_t    v1(rx[ body1 ], ry[ body1 ], rz[ body1 ]);
    nbvertex_t          total_force;
    for(size_t body2 = 0; body2 != count; ++body2)
    {
        if(body1 == body2)
        {
            continue;
        }
        const nbvertex_t    v2(rx[ body2 ], ry[ body2 ], rz[ body2 ]);
        const nbvertex_t    force(m_data->force(v1, v2, mass[body1], mass[body2]));
        total_force += force;
    }
    frx[body1] = vx[body1];
    fry[body1] = vy[body1];
    frz[body1] = vz[body1];
    fvx[body1] = total_force.x / mass[body1];
    fvy[body1] = total_force.y / mass[body1];
    fvz[body1] = total_force.z / mass[body1];
}

Силы притяжения для каждого тела вычисляются независимо, и, для задействования всех ядер процессора, достаточно перед первым циклом написать директиву OpenMP:

Кусочек кода из реализации 'openmp'

#pragma omp parallel for
for(size_t body1 = 0; body1 < count; ++body1)

Т.к. каждое тело взаимодействует с каждым, то для уменьшения количества взаимодействий процессора с ОЗУ и улучшения использования кэша у нас есть возможность загружать в кэш часть данных и использовать их многократно:

Код реализации 'openmp+block'

#pragma omp parallel for
for(size_t n1 = 0; n1 < count; n1 += BLOCK_SIZE)
{
    nbcoord_t            x1[BLOCK_SIZE];
    nbcoord_t            y1[BLOCK_SIZE];
    nbcoord_t            z1[BLOCK_SIZE];
    nbcoord_t            m1[BLOCK_SIZE];
    nbvertex_t           total_force[BLOCK_SIZE];

    for(size_t b1 = 0; b1 != BLOCK_SIZE; ++b1)
    {
        size_t local_n1 = b1 + n1;

        x1[b1] = rx[local_n1];
        y1[b1] = ry[local_n1];
        z1[b1] = rz[local_n1];
        m1[b1] = mass[local_n1];
    }
    for(size_t n2 = 0; n2 < count; n2 += BLOCK_SIZE)
    {
        nbcoord_t            x2[BLOCK_SIZE];
        nbcoord_t            y2[BLOCK_SIZE];
        nbcoord_t            z2[BLOCK_SIZE];
        nbcoord_t            m2[BLOCK_SIZE];

        for(size_t b2 = 0; b2 != BLOCK_SIZE; ++b2)
        {
            size_t local_n2 = b2 + n2;

            x2[b2] = rx[local_n2];
            y2[b2] = ry[local_n2];
            z2[b2] = rz[local_n2];
            m2[b2] = mass[n2 + b2];
        }

        for(size_t b1 = 0; b1 != BLOCK_SIZE; ++b1)
        {
            const nbvertex_t    v1(x1[ b1 ], y1[ b1 ], z1[ b1 ]);
            for(size_t b2 = 0; b2 != BLOCK_SIZE; ++b2)
            {
                const nbvertex_t    v2(x2[ b2 ], y2[ b2 ], z2[ b2 ]);
                const nbvertex_t    force(m_data->force(v1, v2, m1[b1], m2[b2]));
                total_force[b1] += force;
            }
        }
    }

    for(size_t b1 = 0; b1 != BLOCK_SIZE; ++b1)
    {
        size_t local_n1 = b1 + n1;
        frx[local_n1] = vx[local_n1];
        fry[local_n1] = vy[local_n1];
        frz[local_n1] = vz[local_n1];
        fvx[local_n1] = total_force[b1].x / m1[b1];
        fvy[local_n1] = total_force[b1].y / m1[b1];
        fvz[local_n1] = total_force[b1].z / m1[b1];
    }
}

Дальнейшая оптимизация заключается в вынесении содержимого функции вычисления силы в основной цикл и исключении деления и умножения на массу тела m1[b1]. Кроме того, что мы немного сократили вычисления, компилятор на таком развёрнутом цикле сможет применить векторные инструкции процессора SSE и AVX.

Код реализации 'openmp+block+optimization'

#pragma omp parallel for
for(size_t n1 = 0; n1 < count; n1 += BLOCK_SIZE)
{
    nbcoord_t            x1[BLOCK_SIZE];
    nbcoord_t            y1[BLOCK_SIZE];
    nbcoord_t            z1[BLOCK_SIZE];
    nbcoord_t            total_force_x[BLOCK_SIZE];
    nbcoord_t            total_force_y[BLOCK_SIZE];
    nbcoord_t            total_force_z[BLOCK_SIZE];

    for(size_t b1 = 0; b1 != BLOCK_SIZE; ++b1)
    {
        size_t local_n1 = b1 + n1;

        x1[b1] = rx[local_n1];
        y1[b1] = ry[local_n1];
        z1[b1] = rz[local_n1];
        total_force_x[b1] = 0;
        total_force_y[b1] = 0;
        total_force_z[b1] = 0;
    }
    for(size_t n2 = 0; n2 < count; n2 += BLOCK_SIZE)
    {
        nbcoord_t            x2[BLOCK_SIZE];
        nbcoord_t            y2[BLOCK_SIZE];
        nbcoord_t            z2[BLOCK_SIZE];
        nbcoord_t            m2[BLOCK_SIZE];

        for(size_t b2 = 0; b2 != BLOCK_SIZE; ++b2)
        {
            size_t local_n2 = b2 + n2;

            x2[b2] = rx[local_n2];
            y2[b2] = ry[local_n2];
            z2[b2] = rz[local_n2];
            m2[b2] = mass[n2 + b2];
        }

        for(size_t b1 = 0; b1 != BLOCK_SIZE; ++b1)
        {
            for(size_t b2 = 0; b2 != BLOCK_SIZE; ++b2)
            {
                nbcoord_t        dx = x1[b1] - x2[b2];
                nbcoord_t        dy = y1[b1] - y2[b2];
                nbcoord_t        dz = z1[b1] - z2[b2];
                nbcoord_t        r2(dx * dx + dy * dy + dz * dz);
                if(r2 < NBODY_MIN_R)
                {
                    r2 = NBODY_MIN_R;
                }
                nbcoord_t        r = sqrt(r2);
                nbcoord_t        coeff = (m2[b2]) / (r * r2);

                dx *= coeff;
                dy *= coeff;
                dz *= coeff;

                total_force_x[b1] -= dx;
                total_force_y[b1] -= dy;
                total_force_z[b1] -= dz;
            }
        }
    }

    for(size_t b1 = 0; b1 != BLOCK_SIZE; ++b1)
    {
        size_t local_n1 = b1 + n1;
        frx[local_n1] = vx[local_n1];
        fry[local_n1] = vy[local_n1];
        frz[local_n1] = vz[local_n1];
        fvx[local_n1] = total_force_x[b1];
        fvy[local_n1] = total_force_y[b1];
        fvz[local_n1] = total_force_z[b1];
    }
}

Таблица 2. Зависимость времени вычисления (в секундах) функции $inline$ от количества взаимодействующих тел $inline$ для различных CPU реализаций

$inline$	2048	4096	8192	16384	32768
simple	0.0425	0.1651	0.6594	2.65	10.52
openmp	0.0078	0.0260	0.1079	0.417	1.655
openmp+block+optimization	0.0037	0.0128	0.0495	0.194	0.774

Параметры системы:

cистема: Debian 9, Intel Core i7-5820K (6 core)
компилятор: gcc 6.3.0

Хорошо видно, что версия с поддержкой OpenMP ускоряется в шесть раз, ровно по количеству ядер, а оптимизированная версия быстрее ещё немногим более чем в два раза. Так что, при оптимизации не стоит рассчитывать только на параллелизм. Интересно, что при вычислениях в один поток (simple) процессор работал на частоте 3.6 ГГц, в параллельной версии (openmp) сбросил частоту до 3.4 ГГц, а в параллельной и оптимизированной (openmp+block+optimization) сбросил ещё до 3.3 ГГц, но это не помешало ей работать в 13.6 раз быстрее. Также видно, что рост времени вычисления с увеличением размера задачи квадратичен, и дальнейшее увеличение

$inline$ делает задачу нерешаемой за разумное время.

GPU реализация

Но возникает желание производить вычисления ещё быстрее. Есть несколько доступных направлений для ускорения: вычисление на GPU, аппроксимация функции

$inline$ . Сначала для произведения вычислений на GPU я выбрал технологию OpenCL. Для более удобной работы была использована библиотека CLHPP. Главным достоинством OpenCL является то, что код можно запускать и на процессоре, и на GPU, что упрощает написание и отладку, а также расширяет список железа для запуска. В отладке помогает инструмент Oclgrind, который в рантайме показывает неверные вызовы OpenCL API и проблемы при обращении к памяти.

OpenCL

Для начала работы с OpenCL необходимо получить список доступных платформ. Самые распространённые платформы представляют компании AMD, Intel и NVidia.

Код

std::vector<cl::Platform>     platforms;
cl::Platform::get(&platforms);

Далее, после выбора платформы, нужно выбрать вычислительное устройство, которое эта платформа представляет:

Код

const cl::Platform&     platform(platforms[platform_n]);
std::vector<cl::Device> all_devices;
platform.getDevices(CL_DEVICE_TYPE_ALL, &all_devices);

И в заключение подготовительного этапа нужно создать контекст и очереди в рамках которых будет выделяться память и производиться вычисления. Например, контекст, объединяющий все вычислительные устройства выбранной платформы, создаётся следующим образом:

Код создания контекста и очередей

cl::Context                   context(all_devices);
std::vector<cl::CommandQueue> queues;
for(cl::Device device: all_devices)
    queues.push_back(cl::CommandQueue(context, device));

Для загрузки исходного кода на вычислительное устройство его необходимо скомпилировать, для этого предназначен класс cl::Program.

Код компиляции ядра

std::vector< std::string >   source_data;
cl::Program::Sources         sources;
for(int i = 0; i != files.size(); ++i)
{
    source_data.push_back(load_program(files[i]));//Загружаем из файла
    sources.push_back(std::make_pair(source_data.back().data(),
                                     source_data.back().size()));
}
cl::Program    prog(context, sources);
devices.push_back(all_devices);
prog.build(devices, options);

Для описания параметров функции (ядра), которая выполняется на вычислительном устройстве есть удолный шаблон cl::make_kernel.

Пример объявления ядра вычисления силы взаимодействия

typedef cl::make_kernel< cl_int, cl_int, //Block offset
                         cl::Buffer, //mass
                         cl::Buffer, //y
                         cl::Buffer, //f
                         cl_int, cl_int, //yoff,foff
                         cl_int, cl_int //points_count,stride
                         > ComputeBlock;

Далее все просто: объявляем переменную с типом ядра, передаём в неё скомпилированную программу и имя вычислительного ядра, можем запускать ядро почти как обычную функцию.

Код запуска ядра

ComputeBlock        fcompute(prog, "ComputeBlockLocal");
cl::NDRange         global_range(device_data_size);
cl::NDRange         local_range(block_size);
cl::EnqueueArgs     eargs(ctx.m_queue, global_range, local_range);
fcompute(eargs, ...все остальные аргументы); //Собственно, сам вызов ядра.

Само вычислительное ядро для OpenCL очень похоже на вариант 'openmp+block+optimization' для CPU, только в отличии от CPU версии управление первым циклом происходит при помощи OpenCL (диапазон цикла определяется переменной global_range из кода запуска ядра, а из ядра номер текущей итерации доступен с помощью функции get_global_id(0)). Сначала часть данных о телах загружается с локальную память, затем обрабатывается. Локальная память является общей для всех потоков в группе, поэтому загрузка происходит один раз для группы, а обрабатывается каждым потоком в группе и, т.к. локальная память существенно быстрее глобальной, вычисления происходят много быстрее.

Код ядра

__kernel void ComputeBlockLocal(int offset_n1, int offset_n2,
                                __global const nbcoord_t* mass,
                                __global const nbcoord_t* y,
                                __global nbcoord_t* f, int yoff,
                                int foff, int points_count, int stride)
{
    int        n1 = get_global_id(0) + offset_n1;
    __global const nbcoord_t*    rx = y + yoff;
    __global const nbcoord_t*    ry = rx + stride;
    __global const nbcoord_t*    rz = rx + 2 * stride;
    __global const nbcoord_t*    vx = rx + 3 * stride;
    __global const nbcoord_t*    vy = rx + 4 * stride;
    __global const nbcoord_t*    vz = rx + 5 * stride;

    __global nbcoord_t*    frx = f + foff;
    __global nbcoord_t*    fry = frx + stride;
    __global nbcoord_t*    frz = frx + 2 * stride;
    __global nbcoord_t*    fvx = frx + 3 * stride;
    __global nbcoord_t*    fvy = frx + 4 * stride;
    __global nbcoord_t*    fvz = frx + 5 * stride;

    nbcoord_t    x1 = rx[n1];
    nbcoord_t    y1 = ry[n1];
    nbcoord_t    z1 = rz[n1];

    nbcoord_t    res_x = 0.0;
    nbcoord_t    res_y = 0.0;
    nbcoord_t    res_z = 0.0;

    __local nbcoord_t    x2[NBODY_DATA_BLOCK_SIZE];
    __local nbcoord_t    y2[NBODY_DATA_BLOCK_SIZE];
    __local nbcoord_t    z2[NBODY_DATA_BLOCK_SIZE];
    __local nbcoord_t    m2[NBODY_DATA_BLOCK_SIZE];

    // NB! get_local_size(0) == NBODY_DATA_BLOCK_SIZE
    for(int b2 = 0; b2 < points_count; b2 += NBODY_DATA_BLOCK_SIZE)
    {
        int    n2 = b2 + offset_n2 + get_local_id(0);

        // Copy data block to local memory
        x2[ get_local_id(0) ] = rx[n2];
        y2[ get_local_id(0) ] = ry[n2];
        z2[ get_local_id(0) ] = rz[n2];
        m2[ get_local_id(0) ] = mass[n2];

        // Synchronize local work-items copy operations
        barrier(CLK_LOCAL_MEM_FENCE);

        nbcoord_t    local_res_x = 0.0;
        nbcoord_t    local_res_y = 0.0;
        nbcoord_t    local_res_z = 0.0;

        for(int local_n2 = 0; local_n2 != NBODY_DATA_BLOCK_SIZE; ++local_n2)
        {
            nbcoord_t    dx = x1 - x2[local_n2];
            nbcoord_t    dy = y1 - y2[local_n2];
            nbcoord_t    dz = z1 - z2[local_n2];
            nbcoord_t    r2 = (dx * dx + dy * dy + dz * dz);

            if(r2 < NBODY_MIN_R)
            {
                r2 = NBODY_MIN_R;
            }

            nbcoord_t    r = sqrt(r2);
            nbcoord_t    coeff = (m2[local_n2]) / (r * r2);

            dx *= coeff;
            dy *= coeff;
            dz *= coeff;

            local_res_x -= dx;
            local_res_y -= dy;
            local_res_z -= dz;
        }

        // Synchronize local work-items computations
        barrier(CLK_LOCAL_MEM_FENCE);

        res_x += local_res_x;
        res_y += local_res_y;
        res_z += local_res_z;
    }

    frx[n1] = vx[n1];
    fry[n1] = vy[n1];
    frz[n1] = vz[n1];
    fvx[n1] = res_x;
    fvy[n1] = res_y;
    fvz[n1] = res_z;
}

CUDA

Реализация для платформы NVidia CUDA немного проще, чем OpenCL, нам не нужно самим создавать контекст устройства и управлять очередью выполнения (по крайней мере, пока мы не захотим сделать multi-GPU реализацию). Как и в случае с OpenCL, нам необходимо выделить память на GPU, скопировать в неё наши данные, и потом можно запускать вычислительное ядро.

Детальнее о работе с CUDA можно прочитать здесь.

Код запуска CUDA ядра

dim3	grid(count / block_size);
dim3	block(block_size);
size_t	shared_size(4 * sizeof(nbcoord_t) * block_size);

kfcompute <<< grid, block, shared_size >>> (...параметры ядра...);

В отличие от OpenCL, в CUDA задаётся не полный диапазон итераций (в OpenCL реализации это global_range), а задаётся размер грида и размеры блока в гриде, соответственно немного меняется вычисление текущего номера тела, в остальном ядро очень похоже на OpenCL, за исключением других имён функций синхронизации и спецификатора для разделяемой памяти. Ещё полезной отличительной особенностью CUDA является то, что мы можем указывать необходимый размер разделяемой памяти при запуске ядра. Как и в OpenCL реализации, в начале каждого блока итераций мы копируем часть данных в разделяемую память и потом работаем с этой памятью из всех нитей блока.

Код CUDA ядра

__global__ void kfcompute(int offset_n2, const nbcoord_t* y, int yoff, nbcoord_t* f, int foff,
                          const nbcoord_t* mass, int points_count, int stride)
{
    int n1 = blockDim.x * blockIdx.x + threadIdx.x;

    const nbcoord_t*    rx = y + yoff;
    const nbcoord_t*    ry = rx + stride;
    const nbcoord_t*    rz = rx + 2 * stride;

    nbcoord_t    x1 = rx[n1];
    nbcoord_t    y1 = ry[n1];
    nbcoord_t    z1 = rz[n1];

    nbcoord_t    res_x = 0.0;
    nbcoord_t    res_y = 0.0;
    nbcoord_t    res_z = 0.0;

    extern __shared__ nbcoord_t shared_xyzm_buf[];

    nbcoord_t*    x2 = shared_xyzm_buf;
    nbcoord_t*    y2 = x2 + blockDim.x;
    nbcoord_t*    z2 = y2 + blockDim.x;
    nbcoord_t*    m2 = z2 + blockDim.x;

    for(int b2 = 0; b2 < points_count; b2 += blockDim.x)
    {
        int            n2 = b2 + offset_n2 + threadIdx.x;

        // Copy data block to local memory
        x2[ threadIdx.x ] = rx[n2];
        y2[ threadIdx.x ] = ry[n2];
        z2[ threadIdx.x ] = rz[n2];
        m2[ threadIdx.x ] = mass[n2];

        // Synchronize local work-items copy operations
        __syncthreads();

        nbcoord_t    local_res_x = 0.0;
        nbcoord_t    local_res_y = 0.0;
        nbcoord_t    local_res_z = 0.0;

        for(int n2 = 0; n2 != blockDim.x; ++n2)
        {
            nbcoord_t    dx = x1 - x2[n2];
            nbcoord_t    dy = y1 - y2[n2];
            nbcoord_t    dz = z1 - z2[n2];
            nbcoord_t    r2 = (dx * dx + dy * dy + dz * dz);

            if(r2 < NBODY_MIN_R)
            {
                r2 = NBODY_MIN_R;
            }

            nbcoord_t    r = sqrt(r2);
            nbcoord_t    coeff = (m2[n2]) / (r * r2);

            dx *= coeff;
            dy *= coeff;
            dz *= coeff;

            local_res_x -= dx;
            local_res_y -= dy;
            local_res_z -= dz;
        }

        // Synchronize local work-items computations
        __syncthreads();

        res_x += local_res_x;
        res_y += local_res_y;
        res_z += local_res_z;
    }

    n1 += foff;
    f[n1 + 3 * stride] = res_x;
    f[n1 + 4 * stride] = res_y;
    f[n1 + 5 * stride] = res_z;
}

Таблица 3. Зависимость времени вычисления (в секундах) функции $inline$ от количества взаимодействующих тел $inline$ для различных GPU реализаций и лучшей CPU реализации

$inline$	4096	8192	16384	32768	65536	131072
openmp+block+optimization	0.0128	0.0495	0.194	0.774	---	---
OpenCL+половинка NVidia K80	0.004	0.008	0.026	0.134	0.322	1.18
CUDA+половинка NVidia K80	0.004	0.008	0.0245	0.115	0.291	1.13

Где взять NVidia K80

OpenCL и CUDA реализации запускались на бесплатном сервисе Colab от Google, который предоставляет доступ к вычислителям NVidia K80. Подробнее о работе с этим сервисом можно прочитать на хабре.

В целом результат неутешительный, всего лишь в 5-6 раз быстрее, чем CPU реализация. Даже если мы будем вести расчёты на всей K80, то получим ускорение до 12 раз, но т.к. сложность задачи квадратична, то мы за разумное время сможем обрабатывать не 32768 взаимодействующих тел, а 131072, что только в 4 раза больше.

Аппроксимация функции $inline$

Если присмотреться к функции, которой задаётся сила притяжения двух тел, то видно, что она квадратично убывает с расстоянием. Поэтому мы можем точно вычислять силу взаимодействия между близкими телами, и приближённо между отдалёнными. Одним из известных подходов
является алгоритм treecode, предложенный Д. Барнсом и П. Хатом. В моделируемом пространстве строится октодерево, содержащее в своих листьях координаты и массы моделируемых тел. В родительских узлах содержится центр масс, общая масса дочерних узлов и радиус сферы, описанной вокруг тел дочерних узлов. Корень дерева содержит центр масс всех тел, их общую массу и радиус сферы, описанной вокруг них. При расчёте силы взаимодействия сначала считается расстояние до корня дерева, если отношение расстояния до узла к его радиусу больше некоторой константы

$\lambda_{crit}$ , то корень считается одним телом с координатами, равными координатам центра масс входящих в него тел, и массой, равной сумме масс дочерних узлов, если отношение меньше или равно

$\lambda_{crit}$ , то процедура рекурсивно повторяется для каждого дочернего узла. Если достигается лист дерева, то сила взаимодействия считается обычным способом. Таким образом, если одно тело сильно удалено от компактной группы других тел, эта группа представляется для него в виде одного тела, и сила взаимодействия рассчитывается не до каждого тела, а только до одного тела. Благодаря этому сложность алгоритма уменьшается с

$inline$ до

$O(N \log(N))$ ценой некоторой потери точности.

В своём подходе я решил использовать не октодерево, а kd-дерево т.к. оно проще в использовании и имеет более низкие накладные расходы на хранение по сравнению с октодеревом.

Вернёмся обратно к реализации на CPU. Узел kd-дерева можно представить в виде класса, содержащего указатели на левого и правого потомка и информацию о координатах и массе:

Узел kd-дерева

class node
{
    node*      m_left;  //!< Левый потомок
    node*      m_right; //!< Правый потомок
    nbvertex_t m_mass_center;  //!< Координаты центра масс узла
    nbcoord_t  m_mass; //!< Масса узла
    nbcoord_t  m_radius_sqr; //!< Квадрат радиуса описанной сферы, умноженный на lambda_crit
    nbvertex_t m_bmin; //!< Минимальные координаты описанного бокса
    nbvertex_t m_bmax; //!< Максимальные координаты описанного бокса
    size_t     m_body_n; //!< Номер тела, связанного с узлом
};

При таком способе хранения дерева у нас имеется два возможных варианта обхода дерева: либо использовать явную рекурсию, либо использовать стек самим. Я остановился на втором варианте.

Вычисление силы взаимодействия путем обхода дерева

nbvertex_t force_compute(const nbvertex_t& v1,
                         const nbcoord_t mass1)
{
    nbvertex_t  total_force;
    node*       stack_data[MAX_STACK_SIZE] = {};
    node**      stack = stack_data;
    node**      stack_head = stack;

    *stack++ = m_root;
    while(stack != stack_head)
    {
        node*            curr = *--stack;
        const nbcoord_t  distance_sqr((v1 - curr->m_mass_center).norm());

        if(distance_sqr > curr->m_radius_sqr)
        {//Узел достаточно далеко, вычисляем силу и пропускаем его детей
            total_force += force(v1, curr->m_mass_center, mass1, curr->m_mass);
        }
        else
        {// Узел слишком близко, запоминаем детей для последующей обработки
            if(curr->m_right != NULL)
            {
                *stack++ = curr->m_right;
            }
            if(curr->m_left != NULL)
            {
                *stack++ = curr->m_left;
            }
        }
    }
    return total_force;
}

Как и в cлучае «точной» CPU реализации, функция вычисления силы вызывается для каждого тела. Цикл по всем телам может быть легко распараллелен с помощью директив OpenMP.

Но соседние итерации цикла в этом случае будут обращаться к совершенно разным частям дерева, что не даёт эффективно использовать кэш процессора. Для преодоления этой проблемы обход всех тел нужно производить не по исходному порядку, а по порядку, по которому тела расположены в листьях kd-дерева, тогда соседние итерации будут происходить для тел, находящихся близко в пространстве, и обходить дерево по почти одинаковым путям.

Обход листьев дерева

template<class Visitor>
void traverse(Visitor visit)
{
    node*    stack_data[MAX_STACK_SIZE] = {};
    node**    stack = stack_data;
    node**    stack_head = stack;

    *stack++ = m_root;
    while(stack != stack_head)
    {
        node*                curr = *--stack;
        if(curr->m_radius_sqr > 0)
        {//Это не лист. Откладываем детей на стек.
            if(curr->m_left != NULL)
            {
                *stack++ = curr->m_left;
            }
            if(curr->m_right != NULL)
            {
                *stack++ = curr->m_right;
            }
        }
        else
        {// Это листовой узел. Вычисляем силу.
            visit(curr->m_body_n, curr->m_mass_center, curr->m_mass);
        }
    }
}

Эта реализация имеет другую проблему — нет универсального способа распараллелить такой обход дерева. Но мы можем полностью изменить способ хранения дерева в памяти — мы можем хранить все узлы в одном линейном масиве и полностью отказаться от хранения указателей на потомков, по аналогии с построением двоичной кучи. При начале нумерации узлов с

$inline$ , если узел находится в ячейке номером

$inline$ , то его левый потомок находится в ячейке

$inline$ , правый потомок в ячейке

$inline$ , а родитель в ячейке

$inline$ . Правый узел, соответствующий левому с номером

$inline$ , будет иметь номер

$inline$ . Сам массив будет иметь длину

$inline$ , а все листовые узлы будут расположены в последних

$inline$ ячейках, причём, близкие в пространстве узлы будут расположены в близких ячейках массива. Функция обхода дерева немного изменится, но основа остаётся прежней:

Вычисление силы путём обхода дерева в массиве

nbvertex_t force_compute(const nbvertex_t& v1, const nbcoord_t mass1)
{
    nbvertex_t total_force;

    size_t     stack_data[MAX_STACK_SIZE] = {};
    size_t*    stack = stack_data;
    size_t*    stack_head = stack;

    *stack++ = NBODY_HEAP_ROOT_INDEX;
    while(stack != stack_head)
    {
        size_t                curr = *--stack;
        const nbcoord_t       distance_sqr((v1 - m_mass_center[curr]).norm());

        if(distance_sqr > m_radius_sqr[curr])
        {
            total_force += force(v1, m_mass_center[curr], mass1, m_mass[curr]);
        }
        else
        {
            size_t    left(left_idx(curr));
            size_t    rght(rght_idx(curr));
            if(rght < m_body_n.size())
            {
                *stack++ = rght;
            }
            if(left < m_body_n.size())
            {
                *stack++ = left;
            }
        }
    }
    return total_force;
}

Но и это ещё не все возможности, которые открывает нам хранение узлов в массиве — мы можем отказаться от стека при обходе. Для этого в ветку кода, в которой мы переходим к детям узла, мы добавляем функцию вычисления следующего узла (

$next_{up}$ ), а в ветку, в которой вычисляем силу взаимодействия, мы добавляем вычисление следующего узла с пропуском текущего поддерева (

$skip_{down}$ ).

Для пропуска текущего поддерева нам нужно спускаться к корню (направление

$inline$ ), пока мы находимся в правом потомке, как только мы доходим до левого, то переходим в соответствующее ему правое поддерево (направление

$inline$ ), если мы попадаем в корень, то обход дерева завершается.

Код функции пропуска поддерева

index_t skip_down(index_t idx)
{
    // While index is 'right' -> go down
    while(is_right(idx))
    {
        index_t parent = parent_idx(idx);
        // We at root again. Stop traverse.
        if(parent == NBODY_HEAP_ROOT_INDEX)
        {
            return NBODY_HEAP_ROOT_INDEX;
        }
        idx = parent;
    }
    return left2right(idx);
}

Рис 8. Пропуск поддерева $skip_{down}$ .

Для перехода к следующему узлу надо, если возможно, переходить к левому потомку (направление

$inline$ ), а если потомка нет, то переходить к следующему узлу 'снизу' при помощи функции

$skip_{down}$ .

Код функции перехода к следующему узлу

index_t next_up(index_t idx, index_t tree_size)
{
    index_t left = left_idx(idx);
    if(left < tree_size)
    {
        return left;
    }
    return skip_down(idx);
}

Рис 9. Переходы к следующему узлу $next_{up}$ .

Может показаться, что рекурсию мы заменили на цикл

$inline$ в функции

$skip_{down}$ , и такая замена ничего не даёт, но давайте посмотрим, как определить, является ли узел с номером

$inline$ правым потомком. Это можно сделать, просто проверив его номер на нечётность (правый потомок находится в ячейке

$inline$ ), для этого достаточно вычислить

$i\&1$ . Т.е. мы делим число

$inline$ на два если младщий бит выставлен в единицу. Но это можно сделать и без цикла, во многих процессорах есть инструкция find first set, которая возвращает позицию первого установленного бита, воспользовавшись ей мы сворачиваем цикл в четыре инструкции процессора.

Код оптимизированной функции пропуска поддерева

index_t skip_down(index_t idx)
{
    idx = idx >> (__builtin_ffs(~idx) - 1);
    return left2right(idx);
}

После этого мы можем исключить стек из функции обхода дерева и заменить его на пару

$skip_{down} + next_{up}$ , после этого функция выглядит даже проще:

Вычисление силы путём обхода дерева в массиве без использования стека

nbvertex_t force_compute(const nbvertex_t& v1,
                         const nbcoord_t mass1) const
{
    nbvertex_t    total_force;
    size_t        curr = NBODY_HEAP_ROOT_INDEX;
    size_t        tree_size = m_mass_center.size();
    do
    {
        const nbcoord_t  distance_sqr((v1 - m_mass_center[curr]).norm());

        if(distance_sqr > m_radius_sqr[curr])
        {
            total_force += force(v1, m_mass_center[curr], mass1, m_mass[curr]);
            curr = skip_down(curr);
        }
        else
        {
            curr = next_up(curr, tree_size);
        }
    }
    while(curr != NBODY_HEAP_ROOT_INDEX);
    return total_force;
}

Всего у нас получилось шесть сочетаний обхода дерева и вычисления силы. Сравним эти подходы по времени вычисления и по качеству. Возъмём как меру качества относительное изменение полной энергии системы после 100 итераций. В качестве модельной системы возьмём две взаимодействующие «галактики», состоящих из

$inline$ тел каждая.

Таблица 4. Сочетания способа обхода дерева и вычисления силы

Тип обхода дерева / вычисления силы	Дерево со стеком	'Куча' со стеком	'Куча' без стека
Итерации по номеру тела	cycle+tree	cycle+heap	cycle+heapstackless
Обход листьев	nestedtree+tree	nestedtree+heap	nestedtree+heapstackless


a) Зависимость времени вычисления функции $inline$ от отношения критического расстояния до узла дерева к его радиусу ( $\lambda_{crit}$ ).	b) Зависимость относительного изменения энергии в системе от отношения критического расстояния до узла дерева к его радиусу ( $\lambda_{crit}$ ).
Рис 10. Результаты вычисления функции $inline$ различными способами на CPU (количество тел $inline$ )

Видно, что реализация 'nestedtree+tree' безнадёжно отстала по скорости, т.к. в ней отсутствует параллелизм. Лидируют же реализации с расположением узлов дерева в массиве и индексацией как в двоичной куче. Относительное изменение энергии пренебрежима мала для всех вариантов с

$\lambda_{crit} > 1$ . Также на рис. 10a видно, что при (

$\lambda_{crit} < 10$ ) все варианты вычисления функции

$inline$ существенно обгоняют по скорости самый оптимизированный вариант точного вычисления ('openmp+block+optimization'), при дальнейшем увеличении

$\lambda_{crit}$ реализации с деревом проигрывают точной версии.

Обход дерева на GPU

Обход дерева на GPU я пробовал реализовать как при помощи технологии OpenCL, так и CUDA. Вариант хранения узлов в виде дерева был сразу отброшен, и оставлены были только варианты с хранением дерева в массиве с индексацией как в двоичной куче. В целом, реализации вычисительного ядра не сильно отличаются от CPU версии.

OpenCL ядро для вычисления силы путём обхода дерева (обход по порядку нумерации тел)

__kernel void ComputeTreeBH(int offset_n1, int points_count, int tree_size,
                            __global const nbcoord_t* y,
                            __global nbcoord_t* f,
                            __global const nbcoord_t* tree_cmx,
                            __global const nbcoord_t* tree_cmy,
                            __global const nbcoord_t* tree_cmz,
                            __global const nbcoord_t* tree_mass,
                            __global const nbcoord_t* tree_crit_r2)
{
    int        n1 = get_global_id(0) + offset_n1;
    int        stride = points_count;
    __global const nbcoord_t*    rx = y;
    __global const nbcoord_t*    ry = rx + stride;
    __global const nbcoord_t*    rz = rx + 2 * stride;
    __global const nbcoord_t*    vx = rx + 3 * stride;
    __global const nbcoord_t*    vy = rx + 4 * stride;
    __global const nbcoord_t*    vz = rx + 5 * stride;

    __global nbcoord_t*    frx = f;
    __global nbcoord_t*    fry = frx + stride;
    __global nbcoord_t*    frz = frx + 2 * stride;
    __global nbcoord_t*    fvx = frx + 3 * stride;
    __global nbcoord_t*    fvy = frx + 4 * stride;
    __global nbcoord_t*    fvz = frx + 5 * stride;

    nbcoord_t    x1 = rx[n1];
    nbcoord_t    y1 = ry[n1];
    nbcoord_t    z1 = rz[n1];

    nbcoord_t    res_x = 0.0;
    nbcoord_t    res_y = 0.0;
    nbcoord_t    res_z = 0.0;

    int    stack_data[MAX_STACK_SIZE] = {};
    int    stack = 0;
    int    stack_head = stack;

    stack_data[stack++] = NBODY_HEAP_ROOT_INDEX;
    while(stack != stack_head)
    {
        int          curr = stack_data[--stack];
        nbcoord_t    dx = x1 - tree_cmx[curr];
        nbcoord_t    dy = y1 - tree_cmy[curr];
        nbcoord_t    dz = z1 - tree_cmz[curr];
        nbcoord_t    r2 = (dx * dx + dy * dy + dz * dz);

        if(r2 > tree_crit_r2[curr])
        {
            if(r2 < NBODY_MIN_R)
            {
                r2 = NBODY_MIN_R;
            }

            nbcoord_t    r = sqrt(r2);
            nbcoord_t    coeff = tree_mass[curr] / (r * r2);

            dx *= coeff;
            dy *= coeff;
            dz *= coeff;
            res_x -= dx;
            res_y -= dy;
            res_z -= dz;
        }
        else
        {
            int    left = left_idx(curr);
            int    rght = rght_idx(curr);
            if(left < tree_size)
            {
                stack_data[stack++] = left;
            }
            if(rght < tree_size)
            {
                stack_data[stack++] = rght;
            }
        }
    }

    frx[n1] = vx[n1];
    fry[n1] = vy[n1];
    frz[n1] = vz[n1];
    fvx[n1] = res_x;
    fvy[n1] = res_y;
    fvz[n1] = res_z;
}

В первом варианте обход дерева начинался по порядку нумерации тел в исходном массиве, поэтому соседние нити обходили совершенно разные участки дерева, что отрицательно сказывалось на производительности кэш памяти GPU. Поэтому во втором варианте был применён обход, начинающийся с вершины дерева, в этом случае соседние нити начинают обход дерева по одинаковому пути, т.к. соседние вершины дерева находятся поблизости и в пространстве. Также важно, что мы выбрали нумерацию в массиве узлов дерева не с нуля, а с единицы, в этом случае листья дерева хранятся во второй половине массива, и при количестве тел, равному степени двойки, у нас будет выровненный доступ к памяти по индексу tn1.

OpenCL ядро для вычисления силы путём обхода дерева (обход по порядку нумерации узлов дерева)

__kernel void ComputeHeapBH(int offset_n1, int points_count, int tree_size,
                            __global const nbcoord_t* y,
                            __global nbcoord_t* f,
                            __global const nbcoord_t* tree_cmx,
                            __global const nbcoord_t* tree_cmy,
                            __global const nbcoord_t* tree_cmz,
                            __global const nbcoord_t* tree_mass,
                            __global const nbcoord_t* tree_crit_r2,
                            __global const int* body_n)
{
    int        tree_offset = points_count - 1 + NBODY_HEAP_ROOT_INDEX;
    int        stride = points_count;
    int        tn1 = get_global_id(0) + offset_n1 + tree_offset;

    int          n1 = body_n[tn1];
    nbcoord_t    x1 = tree_cmx[tn1];
    nbcoord_t    y1 = tree_cmy[tn1];
    nbcoord_t    z1 = tree_cmz[tn1];

    nbcoord_t    res_x = 0.0;
    nbcoord_t    res_y = 0.0;
    nbcoord_t    res_z = 0.0;

    int stack_data[MAX_STACK_SIZE] = {};
    int    stack = 0;
    int    stack_head = stack;

    stack_data[stack++] = NBODY_HEAP_ROOT_INDEX;
    while(stack != stack_head)
    {
        int          curr = stack_data[--stack];
        nbcoord_t    dx = x1 - tree_cmx[curr];
        nbcoord_t    dy = y1 - tree_cmy[curr];
        nbcoord_t    dz = z1 - tree_cmz[curr];
        nbcoord_t    r2 = (dx * dx + dy * dy + dz * dz);

        if(r2 > tree_crit_r2[curr])
        {
            if(r2 < NBODY_MIN_R)
            {
                r2 = NBODY_MIN_R;
            }

            nbcoord_t    r = sqrt(r2);
            nbcoord_t    coeff = tree_mass[curr] / (r * r2);

            dx *= coeff;
            dy *= coeff;
            dz *= coeff;
            res_x -= dx;
            res_y -= dy;
            res_z -= dz;
        }
        else
        {
            int    left = left_idx(curr);
            int    rght = rght_idx(curr);
            if(left < tree_size)
            {
                stack_data[stack++] = left;
            }
            if(rght < tree_size)
            {
                stack_data[stack++] = rght;
            }
        }
    }

    __global const nbcoord_t*    vx = y + 3 * stride;
    __global const nbcoord_t*    vy = y + 4 * stride;
    __global const nbcoord_t*    vz = y + 5 * stride;

    __global nbcoord_t*    frx = f;
    __global nbcoord_t*    fry = frx + stride;
    __global nbcoord_t*    frz = frx + 2 * stride;
    __global nbcoord_t*    fvx = frx + 3 * stride;
    __global nbcoord_t*    fvy = frx + 4 * stride;
    __global nbcoord_t*    fvz = frx + 5 * stride;

    frx[n1] = vx[n1];
    fry[n1] = vy[n1];
    frz[n1] = vz[n1];
    fvx[n1] = res_x;
    fvy[n1] = res_y;
    fvz[n1] = res_z;
}

При обходе по порядку нумерации узлов дерева мы получили прирост производительости. Но этот вариант тоже можно улучшить. Глобальная память, в которой сейчас располагаются узлы дерева, оптимизирована для коллективного доступа, т.е. нити одной группы должны читать слова, расположенные в одном блоке памяти. В нашем случае, обход дерева начинается по одинаковым путям, и мы запрашиваем одинаковые данные всеми нитями группы, но чем дальше мы углубляемся в дерево, тем пути соседних нитей всё больше расходятся, и нам приходится запрашивать разные данные, что в несколько раз снижает производительность подсистемы памяти. Но узлы каждого поддерева, принадлежащие одному уровню, располагаются в относительно близких ячейках памяти. Т.е. при обходе оставшейся части дерева соседние нити вычислительного ядра обращаются не к одинаковым узлам дерева, но к близко расположенным в памяти. Для оптимизации такого доступа к памяти можно применить текстурную память. Но есть одна загвоздка. На текущий момент текстуры не поддерживают работу с данными двойной точности (мы же хотим точно вести расчёты). Но в CUDA есть функция __hiloint2double, которая собирает число двойной точности из двух целых чисел.

Код запроса чисел двойной точности из текстуры, хранящей целые числа

template<>
struct nb1Dfetch<double>
{
    typedef double4 vec4;
    static __device__ double fetch(cudaTextureObject_t tex, int i)
    {
        int2 p(tex1Dfetch<int2>(tex, i));
        return __hiloint2double(p.y, p.x);
    }
    static __device__ vec4 fetch4(cudaTextureObject_t tex, int i)
    {
        int        ii(2 * i);
        int4    p1(tex1Dfetch<int4>(tex, ii));
        int4    p2(tex1Dfetch<int4>(tex, ii + 1));
        vec4    d4 = {__hiloint2double(p1.y, p1.x),
                      __hiloint2double(p1.w, p1.z),
                      __hiloint2double(p2.y, p2.x),
                      __hiloint2double(p2.w, p2.z)
                  };
        return d4;
    }
};

При этом было сделано две реализации, в одной каждый элемент дерева (x, y, z, tree_crit_r2) запрашивался независимо, а во второй реализации эти запросы были объединены. Запрос массы узла происходит значительно реже, только при выполнении условия r2 > tree_crit_r2[curr], поэтому не имеет смысла объединять этот запрос с остальными. Ещё одной полезной функцией фреймворка CUDA является возможность регулирования соотношения размеров кэша L1 и размера разделяемой памяти (cudaFuncSetCacheConfig). В случае обхода дерева мы не используем разделяемую память, поэтому мы можем увеличить кэш L1 в ущерб ей.

CUDA ядро для вычисления силы путём обхода дерева (обход по порядку нумерации узлов дерева)

__global__ void kfcompute_heap_bh_tex(int offset_n1, int points_count, int tree_size,
                                      nbcoord_t* f,
                                      cudaTextureObject_t tree_xyzr,
                                      cudaTextureObject_t tree_mass,
                                      const int* body_n)
{
    nb1Dfetch<nbcoord_t>    tex;
    int        tree_offset = points_count - 1 + NBODY_HEAP_ROOT_INDEX;
    int        stride = points_count;
    int        tn1 = blockDim.x * blockIdx.x + threadIdx.x + offset_n1 + tree_offset;

    int          n1 = body_n[tn1];
    nbvec4_t     xyzr = tex.fetch4(tree_xyzr, tn1);
    nbcoord_t    x1 = xyzr.x;
    nbcoord_t    y1 = xyzr.y;
    nbcoord_t    z1 = xyzr.z;

    nbcoord_t    res_x = 0.0;
    nbcoord_t    res_y = 0.0;
    nbcoord_t    res_z = 0.0;

    int stack_data[MAX_STACK_SIZE] = {};
    int stack = 0;
    int stack_head = stack;

    stack_data[stack++] = NBODY_HEAP_ROOT_INDEX;
    while(stack != stack_head)
    {
        int            curr = stack_data[--stack];
        nbvec4_t    xyzr2 = tex.fetch4(tree_xyzr, curr);
        nbcoord_t    dx = x1 - xyzr2.x;
        nbcoord_t    dy = y1 - xyzr2.y;
        nbcoord_t    dz = z1 - xyzr2.z;
        nbcoord_t    r2 = (dx * dx + dy * dy + dz * dz);

        if(r2 > xyzr2.w)
        {
            if(r2 < NBODY_MIN_R)
            {
                r2 = NBODY_MIN_R;
            }

            nbcoord_t    r = sqrt(r2);
            nbcoord_t    coeff = tex.fetch(tree_mass, curr) / (r * r2);

            dx *= coeff;
            dy *= coeff;
            dz *= coeff;
            res_x -= dx;
            res_y -= dy;
            res_z -= dz;
        }
        else
        {
            int    left = nbody_heap_func<int>::left_idx(curr);
            int    rght = nbody_heap_func<int>::rght_idx(curr);
            if(left < tree_size)
            {
                stack_data[stack++] = left;
            }
            if(rght < tree_size)
            {
                stack_data[stack++] = rght;
            }
        }
    }

    f[n1 + 3 * stride] = res_x;
    f[n1 + 4 * stride] = res_y;
    f[n1 + 5 * stride] = res_z;
}

Анализ программы в профилировщике nvprof показал, что даже при использовании текстурной памяти для хранения дерева всё равно очень высокая нагрузка на глобальную память.

Действительно, в CUDA вся память ядра, которая адресуется по 'вычисляемым' адресам, хранится в глобальной памяти, и соответственно стек, который нужен для обхода дерева, расположен в глобальной памяти и 'отъедает' значительную часть пропускной способности чипов памяти, ведь стек есть у каждой выполняемой нити, а нитей очень много.

Но, к счастью, мы уже умеем обходить дерево без использования стека. Дополняя предыдущее вычислительное ядро функциями вычисления следующего узла дерева, получаем новое ядро, причём, более компактное.

CUDA ядро для вычисления силы путём обхода дерева без использования стека

__global__ void kfcompute_heap_bh_stackless(int offset_n1, int points_count, int tree_size,
                                            nbcoord_t* f,
                                            cudaTextureObject_t tree_xyzr,
                                            cudaTextureObject_t tree_mass,
                                            const int* body_n)
{
    nb1Dfetch<nbcoord_t>    tex;
    int        tree_offset = points_count - 1 + NBODY_HEAP_ROOT_INDEX;
    int        stride = points_count;
    int        tn1 = blockDim.x * blockIdx.x + threadIdx.x + offset_n1 + tree_offset;

    int          n1 = body_n[tn1];
    nbvec4_t     xyzr = tex.fetch4(tree_xyzr, tn1);
    nbcoord_t    x1 = xyzr.x;
    nbcoord_t    y1 = xyzr.y;
    nbcoord_t    z1 = xyzr.z;

    nbcoord_t    res_x = 0.0;
    nbcoord_t    res_y = 0.0;
    nbcoord_t    res_z = 0.0;

    int    curr = NBODY_HEAP_ROOT_INDEX;
    do
    {
        nbvec4_t     xyzr2 = tex.fetch4(tree_xyzr, curr);
        nbcoord_t    dx = x1 - xyzr2.x;
        nbcoord_t    dy = y1 - xyzr2.y;
        nbcoord_t    dz = z1 - xyzr2.z;
        nbcoord_t    r2 = (dx * dx + dy * dy + dz * dz);

        if(r2 > xyzr2.w)
        {
            if(r2 < NBODY_MIN_R)
            {
                r2 = NBODY_MIN_R;
            }

            nbcoord_t    r = sqrt(r2);
            nbcoord_t    coeff = tex.fetch(tree_mass, curr) / (r * r2);

            dx *= coeff;
            dy *= coeff;
            dz *= coeff;
            res_x -= dx;
            res_y -= dy;
            res_z -= dz;
            curr = nbody_heap_func<int>::skip_idx(curr);
        }
        else
        {
            curr = nbody_heap_func<int>::next_up(curr, tree_size);
        }
    }
    while(curr != NBODY_HEAP_ROOT_INDEX);

    f[n1 + 3 * stride] = res_x;
    f[n1 + 4 * stride] = res_y;
    f[n1 + 5 * stride] = res_z;
}

Производительность ядер, выполняемых на GPU, сильно зависит от размера блоков, на которые мы делим задачу. От этого размера зависит, сколько регистров, локальной памяти и прочих ресурсов будет доступно каждой вычислительной нити. Также нужно иметь в виду, что во время ожидания доступа к памяти в одной нити, другая нить может производить вычисления на шейдерном процессоре, таким образом, при достаточном количестве одновременно выполняемых нитей на одном процессоре время доступа к памяти будет скрыто за вычисленими. Поэтому перед сравнением производительности наших ядер нам надо для каждого из них вычислить оптимальный размер блока. Произведём сравнение на доступной нам половинке от NVidia K80.

Таблица 5. Зависимость времени вычисления (в секундах) функции $inline$ от размера блока для различных GPU реализаций при количестве тел $inline$ и $\lambda_{crit}=10$

Размер блока/ядро	8	16	32	64	128	256	512	1024
opencl+dense	5.77	2.84	1.46	1.13	1.15	1.14	1.14	1.13
cuda+dense	5.44	2.55	1.27	0.96	0.97	0.99	0.99	-
opencl+heap+cycle	5.88	5.65	5.74	5.96	5.37	5.38	5.35	5.38
opencl+heap+nested	4.54	3.68	3.98	5.25	5.46	5.41	5.48	5.31
cuda+heap+nested	3.62	2.81	2.68	4.26	4.84	4.75	4.8	4.67
cuda+heap+nested+tex	2.6	1.51	0.912	0.7	1.85	1.75	1.69	1.61
cuda+heap+nested+tex+stackless	2.3	1.29	0.773	0.5	0.51	0.52	0.52	0.52

Таблица 6. Зависимость времени вычисления (в секундах) функции $inline$ от размера блока для различных GPU реализаций при количестве тел $inline$ и $\lambda_{crit}=4$

Размер блока/ядро	8	16	32	64	128	256	512	1024
opencl+dense	366	179	89.9	69.3	70.3	69.1	68.9	68.0
cuda+dense	346	162	79.6	60.8	60.8	60.7	59.6	-
opencl+heap+cycle	16.2	18.2	20.1	21.2	21.2	21.3	21.2	21.1
opencl+heap+nested	10.5	7.63	6.38	8.23	9.95	9.89	9.65	9.66
cuda+heap+nested	8.67	6.44	5.39	5.93	8.65	8.61	8.41	8.27
cuda+heap+nested+tex	6.38	3.57	2.13	1.44	3.56	3.46	3.30	3.29
cuda+heap+nested+tex+stackless	5.78	3.19	1.83	1.21	1.11	1.10	1.11	1.13

Сложная ситуация, но, в отличие, от CPU версии обхода дерева, видно, что каждый шаг оптимизации приносит ощутимые плоды. Реализация 'opencl+heap+cycle' почти в 6 раз медленнее точного решения с полным вычислением функции

$inline$ . Реализация 'opencl+heap+nested', в которой обход дерева начинается с соседних узлов, быстрее предыдушей в 1.4 раза, т.к. лучше задействуется кэш память. В реализации 'cuda+heap+nested' увеличен кэш L1 в ущерб разделяемой памяти, что дало прирост скорости ещё в 1.4 раза, хотя не исключено, что в cuda реализации более оптимально скомпилированно вычислительное ядро (в версиях 'opencl+dense' и 'cuda+dense' ядра идентичные, а быстродействие у cuda версии выше в ~1.2 раза). Самый существенный прирост скорости вычислений (в 3.8 раза) достигается при расположении дерева в текстурной памяти и объединении запросов к элементам узла дерева. Реализация с обходом дерева без использования стека 'cuda+heap+nested+tex+stackless' быстрее ещё в 1.4 раза, т.к. в ней вся пропускная способность шины памяти используется только для доступа к данным об узлах дерева и не расходуется на стек. Таким образом, при

$\lambda_{crit}=10$ удалось достичь ускорения в два раза по сравнению с полным вычислением функции

$inline$ . Но

$\lambda_{crit}=10$ избыточно большое значение параметра, на графике зависимости относительного изменения энергии в системе от отношения критического расстояния до узла дерева к его радиусу для CPU реализации видно, что можно использовать более низкие значения

$\lambda_{crit}$ без видимой потери точности решения. Попробуем поварьировать

$\lambda_{crit}$ при оптимальных размерах блока, которые мы определили на предыдущем шаге.


a) $inline$	b) $inline$
Рис 11. Зависимость времени вычисления функции $inline$ от отношения критического расстояния до узла дерева к его радиусу ( $\lambda_{crit}$ ) для различных GPU реализаций обхода дерева

Видно, что при малых

$\lambda_{crit}$ все способы вычисления функции

$inline$ выходят на близкие значения, определяемые временем построения kd-дерева и подготовкой данных для GPU. Причём, время построения дерева вносит весомый вклад в общее время до

$\lambda_{crit}\leq 4$ , далее этим временем можно пренебречь. Интересно заметить, что при

$inline$ производительность опять улучшается при достижении

$\lambda_{crit}=1024$ , скорее всего это обусловлено тем, что все нити GPU обходят одни и те же вершины дерева одновременно, что улучшает использование кэша L1 и полностью исключает 'branch divergence'. Также видно, что наилучшее быстродействие у реализации без использования стека (cuda+heap+nested+tex+stackless), она обгоняет по скорости вариант со стеком примерно в

$inline$ раза. Прочие реализации ещё в несколько раз медленнее. Для закрепления результата проведём расчет времени на GPU с более новой архитектурой.

Результаты запуска на GeForce GTX 1080 Ti (одинарная точность)

Таблица 7. Зависимость времени вычисления (в секундах) функции $inline$ от размера блока для различных GPU реализаций при количестве тел $inline$ и $\lambda_{crit}=4$

Размер блока/ядро	8	16	32	64	128	256	512	1024
opencl+dense	47.8	23.4	11.6	11.59	12.8	12.8	12.8	12.8
cuda+dense	49.0	23.8	11.9	11.8	11.7	11.7	11.7	11.7
opencl+heap+cycle	7.48	8.26	7.73	7.36	7.33	7.27	7.25	7.26
opencl+heap+nested	1.33	1.20	1.41	2.46	2.53	2.49	2.44	2.47
cuda+heap+nested	1.23	1.10	1.31	2.28	2.29	2.28	2.23	2.25
cuda+heap+nested+tex	0.88	0.68	0.654	1.6	1.6	1.6	1.6	1.6
cuda+heap+nested+tex+stackless	0.71	0.47	0.45	0.43	0.43	0.42	0.41	0.395

Рис 12. Зависимость времени вычисления функции $inline$ от отношения критического расстояния до узла дерева к его радиусу ( $\lambda_{crit}$ ) для различных GPU реализаций обхода дерева

При использовании GeForce GTX 1080 Ti для вычисления, разница между реализациями обхода дерева со стеком и без стека достигает двух раз, при условии, что мы пренебрегаем временем построения дерева. Этот факт подталкивает нас к тому, чтобы в будующем перенестии на GPU и построение дерева. Из сравнения таблиц 5-7 видно, что нет единого оптимального размера блока для разного количества тел и тем более для разных архитектур GPU, причём, разница времени вычисления может достигать нескольких раз, даже если не брать в расчёт граничные значения. Таким образом, перед длительными расчетами имеет смысл определять оптимальный размер блока для каждой конфигурации.

Главное, чего мы достигли — это возможность моделировать попарное взаимодействие немногим более миллиона тел (

$2^{20}$ ) за разумное время на одном, не самом новом GPU. На более новых GPU (Tesla V100), видимо, можно будет обрабатывать уже около двух миллионов взаимодействующих тел за разумное время, т.к. их производительность примерно в четыре раза выше, чем у половинки Tesla K80. Хотя, это число ещё и несравнимо с количеством звёзд даже в карликовых галактиках, таких как Малое Магелланово Облако, но, на мой взгляд, является внушительным.

Заключение

Применение вложенных методов решения дифференциальных уравнений позволяет решать подобные задачи с высоким уровнем точности при сравнительно небольших временных затратах, а применение алгоритмов аппроксимации функции вычисления силы попарного притяжения позволяют обрабатывать колоссальные объёмы взаимодействующих тел. Таким образом, в отличие от инопланетян из «Задачи трёх тел», мы в состоянии решить задачу

$inline$ тел, и для небольшого числа тел и для целых звёздных скоплений, хотя и ценой некоторой потери точности.

Визуализация

Для тех, кто дочитал до конца приведу ещё несколько видеороликов с визуализацией процесса эволюции систем тел.

Моделирование столкновения двух галактик. Общее количество тел 60 тысяч.

Моделирование эволюции галактики, состоящей из миллиона звёзд. В центре тело с массой 99% от общей. Одиночные частицы практически неразличимы. Уже больше напоминает волны в капле жидкости. Раскрашено в соответствии с модулем скорости. Низкая скорость — синий цвет, средняя — зелёный, высокая — красный. Видно, что в центре скорость выше, и плавно убывает к краям, а самая низкая в экваториальной плоскости.

Более 'динамичное' моделирование эволюции галактики, состоящей из миллиона звёзд. В центре тело с массой 10% от общей. Центральное тело влияет на остальные существенно меньше. Сначала 'звёзды' разлетаются, потом собираются обратно в несколько скоплений, и в конце опять образуют одно большое скопление.

В процессе моделирования порядка 5% 'звёзд' покинуло начальную область безвозвратно.

На 10-й секунде очень напоминает существующую галактику Колесо телеги.

Код проекта можно найти на гитхабе.

Комментарии (96)

Wesha
18.03.2019 22:31
#19905490
+6
как взорвать галактику не выходя из кухни
Напомнило

akhalat
18.03.2019 22:59
#19905582
+3
В отличии от нас

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

А в целом: замечательная работа, но в реальности одной только теории ньютона недостаточно, приходится учитывать и эффекты теории относительности.

Да и в вычислительном плане тоже можно не просто любой метод решения дифф-уров брать с потолка, есть понятия «мягкой» и «жесткой» задач, чисел обусловленности, сходимости получившейся разностостной схемы и прочего подобного…
Но это для серьезных расчетов, а с тем чтобы просто поиграться — отлично.

Sdima1357
18.03.2019 23:55
#19905802
+2
Неплохо бы и привести сравнение с n-body из sdk Nvidia. Скорость, качество и тд. У них во всяком случае gpu-шная версия быстрее CPU на пару порядков. А так хорошая работа… Просто картина не полная.
1. AndrewSu Автор
  19.03.2019 21:47
  #19910960
  +1
  Посмотрел код примера из CUDA SDK. Вы этот имели ввиду?
  
  CPU код из SDK примерно соответствует моему на стадии 'openmp'. GPU код у нас организован похоже, только у них перед внутренним циклом есть директива 'unroll'.
  Таким образом, если сравнивать мои 'openmp' и СUDA реализации, то разница будет примерно в 20 раз. При сравнении с однопоточной версией CUDA реализация быстрее в 120 раз.
  
  В статье, сопропровождающей SDK, есть такие строки:
  
  The result is an algorithm that runs more than 50 times as fast as a highly tuned serial
  implementation (Elsen et al. 2006) or 250 times faster than our portable C implemen-
  tation.
  
  Непонятно что такое 'portable C implementation', но первая половина хорошо соотносится с моими результатами.
  
  На днях попробую собрать пример из SDK и сравнить более детально.
  1. Sdima1357
    19.03.2019 22:06
    #19911056
    Посмотрел код примера из CUDA SDK. Вы этот имели ввиду?
    
    Ага, именно ее.
    
    Непонятно что такое 'portable C implementation',
    
    Обычно под этим имеют в виду однопоточную и не заточенную под архитектуру конкретного процессора. Сейчас компиляторы значительно умнее тех что были в 2006.
1. AndrewSu Автор
  21.03.2019 00:23
  #19917832
  В последних версиях CUDA SDK появился модуль OptiX, предназначенный для трассировки лучей. В нём есть реализация обхода деревьев на GPU. Очень интересно сравнить с этой реализацией. Может кто-то уже сталкивался?

slovak
19.03.2019 00:09
#19905854
А как же черная материя и/или энергия? ;-)
1. agarus
  19.03.2019 11:04
  #19907318
  +1
  Тёмная она :) Или скорее мы. Пока.
1. Wizard_of_light
  19.03.2019 17:06
  #19909560
  +2
  Так это по сути эволюция темноматериальной системы и есть — учитывается только гравитационное взаимодействие. Чтобы в первом приближении прикинуть как это с тёмной материей выглядит, надо распределить равномерно по объему, 9/10 точек не отображать, а для отображаемых ввести небольшую потерю энергии при близком прохождении друг к другу.

bearad
19.03.2019 00:38
#19905922
Низкая скорость — синий цвет, средняя — зелёный, высокая — красный. Видно, что в центре скорость выше, и плавно убывает к краям, а самая низкая в экваториальной плоскости. В центре находится тело с массой 99.9% от общей

А как же темная материя и вот это все? Со схожими скоростями вещества по всему объему.
1. trapwalker
  20.03.2019 14:27
  #19915030
  +1
  Я так понимаю, что раз уж здесь не рассматриваются столкновения, а все массы — это математические точки, то иного взаимодействия, кроме гравитационного, нет. Это значит, что автор рассмотрел исключительно тёмную материю. Казалось бы, да? Но нет. Тёмная материя не собирается в кусочки именно потому, что не взаимодействует механически. Это значит, что в задаче появляется неимоверно огромное количество тел очень малой массы. Такое мы уж никак не смоделируем.

Daddy_Cool
19.03.2019 00:56
#19905966
Не так давно я прочёл фантастический роман «Задача трёх тел» Лю Цысиня

И я. ) Я тоже вдохновился и решил промоделировать с целью написать статью на Хабр когда-нибудь, но кажется слегка не успел. Вот моя простенькая симуляция:
youtu.be/dj1rm7rsqwI
Три тела равной массы, скорости заданы такими, чтобы суммарный импульс был равен нулю — чтобы звезды по возможности не улетели в другую галактику.
Хотя возможно я еще продолжу, меня интересуют более физические аспекты нежели программные.
Вообще я безумно впечатлён проделанной вами работой! Думается это вполне можно опубликовать в журнале типа «Вычислительные методы и программирование».
num-meth.srcc.msu.ru
Хотя если вы не занимаетесь непосредственно научной деятельностью — то смысла особого нет, думается на Хабре статью прочтет больше народа.

Одиночные частицы практически неразличимы. Уже больше напоминает волны в капле жидкости.

Я как-то видео работу — на конференции по неустойчивостям, так там так и делали — брали уравнение Навье-Стокса, добавляли гравитационную силу и получали спиральные рукава галактик, т.е. множество звезд моделировалось сплошной средой. (Хотя возможно это и распространенный метод в астрофизике, не знаю).
А какие у вас заданы абсолютные размеры/массы/времена, т.е. соответствует ли видео в каком-то виде реальности или не особо?
1. Wesha
  19.03.2019 01:26
  #19906038
  Три тела равной массы, скорости заданы такими, чтобы суммарный импульс был равен нулю
  А теперь и Вы напомнили...
  1. Daddy_Cool
    19.03.2019 10:10
    #19907038
    Каждому хочется почувствовать себя демиургом!
1. Fedorkov
  19.03.2019 14:40
  #19908518
  Хунта был великолепным таксидермистом. Штандартенфюрер, по словам Кристобаля Хозевича, — тоже. Но Кристобаль Хозевич успел раньше.
  Я сделал свою модель до того, как прочитал книгу, но мне это не помогло написать статью первым. :) github.com/dsdante/constel
  
  Ещё есть идея решить задачу методом самосогласованного поля.
  1. AndrewSu Автор
    20.03.2019 00:15
    #19911654
    +1
    Ещё есть идея решить задачу методом самосогласованного поля.
    Очень похоже на метод Ахмада-Коэна для решения задачи N-тел.
    Первое моделирование было проведено как раз с применением этого метода:
1. AndrewSu Автор
  19.03.2019 22:08
  #19911074
  +1
  Думается это вполне можно опубликовать в журнале типа «Вычислительные методы и программирование».
  
  Я однажды хотел опубликоваться в этом журнале, но был ответ, что публикация возможна только для тех, кто учится/работает в университете. А я уже закончил. Пришлось в другом месте публиковать. Сейчас, возможно, правила изменились.
  
  А какие у вас заданы абсолютные размеры/массы/времена, т.е. соответствует ли видео в каком-то виде реальности или не особо?
  
  Из физических констант задаю только G=1, а начальный диаметр 'галактики' задаю равным 100. Массы все одинаковые кроме центрального тела. К реальным величинам пересчитать возможно, но я этого пока не делал.

slava_k
19.03.2019 01:36
#19906058
Большое спасибо за интересную статью!

Alek_roebuck
19.03.2019 01:54
#19906094
+1
Если я правильно помню, метод Верле более или менее сохраняет энергию (по сравнению с методом Рунге-Кутта, несмотря на более высокую точность последнего).

Вообще, численное интегрирование задачи N тел — очень изученная тема, только в другом приложении, не астрономическом: в физике это принято называть молекулярной динамикой.
1. AN3333
  19.03.2019 08:56
  #19906694
  Совершенно верно. Есть куча методов сохраняющих энергию.
  И да, Верле самый знаменитый. Он ее сохраняет точно (видимо, с точностью до количества машинных знаков).
  
  И методы высоких порядков тоже имеются. Я пользуюсь VEFRL.
  У таких алгоритмов, не смотря на их дискретность, имеется эффективный гамильтониан. При увеличении шага падает точность, гамильтониан все дальше уходит о оригинала, а вот энергии деваться некуда.
  1. AndrewSu Автор
    19.03.2019 22:13
    #19911122
    +1
    Я пользуюсь VEFRL.
    
    Буду благодарен за ссылку на метод. Беглый поиск в гугле не даёт мне результата.
    
    AN3333
    20.03.2019 10:08
    #19912974
    https://www.sciencedirect.com/science/article/pii/S0010465502004514
    
    AndrewSu Автор
    20.03.2019 20:48
    #19917210
    Хабраблагодарю
1. AN3333
  19.03.2019 09:05
  #19906726
  Кстати, на счет молекулярной динамики. Оно вроде бы и изучено в хвост и в гриву, а хорошего решения до сих пор и не найдено.
  1. Daddy_Cool
    19.03.2019 10:12
    #19907056
    Проблема, что надо откуда-то брать потенциалы взаимодействий.
    Кстати из забавного — оказывается иногда хороший потенциал для своих нужд можно просто купить. Существует такой вариант э… платной науки.
    
    AN3333
    19.03.2019 10:27
    #19907148
    оказывается иногда хороший потенциал для своих нужд можно просто купить.
    
    Серьёзно… Не подскажете где?
    А то может и я начну потенциалами приторговывать.
1. Si010ver
  19.03.2019 10:07
  #19907004
  В астрофизике тоже. Только там он зовется SPH.
1. AndrewSu Автор
  19.03.2019 22:11
  #19911104
  На самом деле с сохранением энергии не всё так просто, даже у метода Верле. Но метод хороший, я его планировал реализовать, но пока руки не дошли.
  1. AN3333
    20.03.2019 10:25
    #19913072
    По вашей ссылке
    
    Symplectic integrators usually used in molecular dynamics, such as the Verlet integrator family, exhibit increases in energy over very long time scales, though the error remains roughly constant.
    
    Честно говоря, для меня это новость. Я не вижу повышения энергии, по крайнем мере за часы гоняния программы. Может быть тут проблема с ограниченной машинной точностью, а не в самом алгоритме, если этот эффект вообще есть?
    
    vlanko
    20.03.2019 15:11
    #19915370
    мой небольшой старый опыт с методами типа Рунге-Кутты показывает, что энергия только падает. Возможно, зависит от метода и железа.
    
    AN3333
    20.03.2019 15:30
    #19915466
    Забавно, а я думал что она только растет, для не симплектических интеграторов :)

TitovVN1974
19.03.2019 07:39
#19906454
{не проснувшись} для подобных задач используют методы частиц в ячейках.
Они снижают вычислительную сложность упуская детали локальной динамики парных взаимодействий.

Beshere
19.03.2019 08:46
#19906648
С одной, стороны подобную программу я написал еще в 10 классе на паскале, применяя куда как меньше математики и железа.

С другой — ваша «модель» не учитывает тёмную материю и энергию (или что там стоит за этими абстракциями) и потому — неверна.
1. Daddy_Cool
  19.03.2019 10:18
  #19907100
  Позволю себе вмешаться.
  В общем случае модель не может быть верной или не верной. Она может описывать результаты экспериментов с некой точностью в определенных условиях.
  Темная материя была введена/придумана как раз чтобы хоть как-то объяснить расхождения между имеющимися моделями и экспериментами.
  Кстати — а как качественно должны измениться результаты представленного моделирования если ввести тёмную материю? Я не специалист по астрофизике, но если вы ориентируетесь, то думаю на пальцах показать можете.
  1. AN3333
    19.03.2019 10:33
    #19907190
    Недавно читал, извините ссылку не помню, что при введении в модель темной материи как раз начинает получаться что-то похожее. Такая вот сетка с перемычками из спиральных галактик.
    За что купил, за то продаю :)
  1. rinaty
    19.03.2019 15:54
    #19909014
    Скорость звезд должна перестать убывать с радиусом
  1. Endeavour
    19.03.2019 18:00
    #19909884
    Никак, потому что это — симуляция системы гравитационно взаимодействующих тел. Темная материя — всего лишь вид материи, взаимодействующей исключительно гравитационно, поэтому данная симуляция полностью применима к ней. Для приближения к «реальности» достаточно выключить отображение 80% частиц.
    
    Большее приближение к реальности даст использование теории относительности вместо уравнения Ньютона. Темная энергия, по идее, автоматически при этом появится.

reykeycom
19.03.2019 10:07
#19907008
Большое спасибо за статью, интересно и поучительно. Интересно на сколько это далеко от реальности? Например из-за скорости распространения гравитации? Или из-за релятивистских эффектов. Интуиция подсказывает что картина должна быть кардинально иной, но интересно насколько?.. Не знаю как остальных, но но крайней мере меня бы очень порадовало релятивистское продолжение статьи.

HEKOT
19.03.2019 10:09
#19907034
А вот если добавить звездообразование, то можно спиральные рукава получить?
ошибочка
В центре тело с массой 99% от общей. Одиночные частицы практически неразличимы. Уже больше напоминает волны в капле жидкости. Раскрашено в соответствии с модулем скорости. Низкая скорость — синий цвет, средняя — зелёный, высокая — красный. Видно, что в центре скорость выше, и плавно убывает к краям, а самая низкая в экваториальной плоскости. В центре находится тело с массой 99.9% от общей.
1. vyo
  19.03.2019 13:01
  #19907932
  Для ошибочек ctrl+Enter ведь добавили. Не все ещё знают видать (или привычка многолетняя может делает своё дело).
  
  HEKOT
  20.03.2019 05:08
  #19912154
  Извиняюсь. На новостных сайтах пользуюсь, перлы малограмотных журналюг исправляю, а тут не проассоциировал :)
1. AndrewSu Автор
  19.03.2019 22:21
  #19911166
  +2
  Спиральные рукава и так есть. Просто они не так сильно 'визуально' выражены. В реальных галактиках они выглядят сильно ярче как раз из-за активного звёздообразования.
  
  Ошибочку поправил, спасибо.
  Но сообщать лучше с личку.
  
  HEKOT
  20.03.2019 03:58
  #19912068
  А вот откуда они взялись без звездообразования? Насколько я понимаю, без звездообразования их не то, что не видно, а из нет. Очевидно, что на периферии скорость вращения отличается от скорости врашения в центре, а значит, никаких протяжённых радиальных структур сохраняться не может, рукава быстренько закрутятся и потеряют свою красивую сприальную структуру. В данной модели в этом можно убедиться, воспользовавшись раскраской.

santa324
19.03.2019 10:12
#19907060
+1
Спасибо, интересная статья!
Скажите, Вы не пробовали оценивать не отклонение полной энергии и импульса системы по окончанию моделирования, а сравнить результаты у 2х моделирований, одно из которых имеет точность с огромным запасом? При этом интересно отклонение не средних значений системы и не общей энергии и импульса, а максимальное отклонение для отельного тела среди всех тел.
Вероятно, при вычислении с низким и высоким порядком точности именно это и сравнивается? Как я понимаю, там сравниваются координаты частицы полученные на следующем шаге для 2х разных точностей. Но не понятно какое отклонение считать малым…
Дело в том что если вся система «равномерная» но в ней есть компактная подсистема из всего нескольких тел но находящихся в очень жестком взаимодействии, то даже если эта подсистема полностью «разойдется» (так как для нее точность неприемлемо низкая) — то на общей энергии/импульсе системы это практически не скажется. Если же сравнивать координаты частиц при разных точностях, тогда не однозначно что является малым отклонением — незначительная погрешность на большом расстоянии от массивного тела может быть неприемлемо огромной в момент пролета вблизи него. Возможно можно сравнивать еще отклонения по импульсу и энергии отдельной частицы… Было бы очень интересно посмотреть на результаты таких экспериментов.
1. AndrewSu Автор
  19.03.2019 22:36
  #19911216
  Вероятно, при вычислении с низким и высоким порядком точности именно это и сравнивается? Как я понимаю, там сравниваются координаты частицы полученные на следующем шаге для 2х разных точностей.
  
  Да, именно так. У меня сравниваются не только координаты, но и скорости.
  
  Но не понятно какое отклонение считать малым…
  
  Этот параметр является настраиваемым, также как и другие.
  
  Дело в том что если вся система «равномерная» но в ней есть компактная подсистема из всего нескольких тел но находящихся в очень жестком взаимодействии, то даже если эта подсистема полностью «разойдется» (так как для нее точность неприемлемо низкая) — то на общей энергии/импульсе системы это практически не скажется.
  
  У мня реализованно так, что если часть системы разойдется, то шаг уменьшается для всех, и вычисления проводятся с новым шагом.
  Здесь, конечно, есть большой простор для оптимизации…
  При интегрировании системы законы сохранения никак не учитываются, это только дополнительный контролируемый параметр.
  1. santa324
    20.03.2019 10:02
    #19912948
    Этот параметр является настраиваемым, также как и другие.
    
    Вот как раз интересно как его настраивать, как понять что такое-то значение достаточно маленькое и не надо в настройках поставить еще меньше?
    
    AndrewSu Автор
    21.03.2019 00:39
    #19917872
    Нужно взять систему с известной эволюцией, опробовать различные параметры на ней, и определить границы применимости. Потом можно приступать к исследованию целевой системы. В нашем случае можно взять задачу Кеплера, попробовать разные траектории и сравнить с аналитическим решением.

paluke
19.03.2019 10:20
#19907112
Так дело не только в точности вычислений, но и в точности исходных данных. А при большом числе тел могут быть такие конфигурации, когда незначительное изменение исходных данных приводит к огромному расхождению по прошествии некоторого времени. ru.wikipedia.org/wiki/%D0%AD%D1%84%D1%84%D0%B5%D0%BA%D1%82_%D0%B1%D0%B0%D0%B1%D0%BE%D1%87%D0%BA%D0%B8
Возьмите например маленькое тело в точке Лагранжа L1 между двумя большими, и посмотрите на эволюцию системы, если начальное положение отклоняется от точки Лагранжа на 0,0001% в одну и в другую сторону.

rPman
19.03.2019 10:28
#19907160
Вопрос-предложение, где бы найти базу по звездам хотя бы нашей галактики, на сколько я понимаю она наиболее изученная в плане данных по координатам, скорости и массе звезд?

Чтобы подсунуть реальные данные в программу. Дело в том что начальные значения скорости в тестах в статье — нулевые, т.е. вся система испытывает невероятно сильные колебания и всплески именно из-за этого, но реалии таковы что звезды двигаются уже по относительно стабильным орбитам.

p.s. какое максимальное приближение звезд друг к другу наблюдается? интересно, на сколько высока вероятность, что пролетающая звезда сорвет одну из планет или изменит ее траекторию?
1. mobi
  19.03.2019 15:24
  #19908780
  Самые свежие и точные данные на текущий момент — это GAIA DR2: gea.esac.esa.int/archive (но всех звезд, естественно, не будет ни в одной базе до тех пор как человечество не начнет путешествовать по Галактике)
1. Wizard_of_light
  19.03.2019 15:50
  #19908980
  Эм, ну, реальные данные в принципе нарыть можно- воспользоваться, архивом телескопа Gaia, например. Но вот только есть нюанс — там данные о примерно 1,7 миллиардах звёзд (и это ещё предположительно порядка 1% от общего числа звёзд в нашей галактике). Правда, полные данные — координаты и скорости движения — есть только примерно для семи миллионов звёзд, но даже этого достаточно, чтобы при попытке прямого расчёта заставить любое доступное железо глубоко и надолго задуматься.
1. AndrewSu Автор
  19.03.2019 22:54
  #19911294
  Чтобы подсунуть реальные данные в программу. Дело в том что начальные значения скорости в тестах в статье — нулевые, т.е. вся система испытывает невероятно сильные колебания и всплески именно из-за этого, но реалии таковы что звезды двигаются уже по относительно стабильным орбитам.
  
  Если бы это было так, то всё звёзды просто упали бы на центральное тело. Скорости как раз выбраны так, что начальные орбиты 'почти' стабильные.
  
  В том моделировании это видно, все движутся по эллиптическим орбитам:

susel
19.03.2019 10:42
#19907230
Видео про симуляцию столкновения напомнило это: видео сравнения симуляции с фотографиями хаббла

huch
19.03.2019 10:50
#19907264
у меня сложилось впечатление, что вы за несколько минут моделируете процессы длительностью в сотни миллионов лет.
какие начальные условия вы задаете? массы тел, их скорости, расстояния между ними? сколько шагов интегрирования делаете?

Rhombus
19.03.2019 11:10
#19907356
Насколько я помню, в задаче трех (и более) тел решения являются неустойчивыми (aka динамический хаос), так что польза от этих вычислений сомнительная, как мне кажется.
1. Umpiro
  19.03.2019 11:49
  #19907562
  Насколько я помню, аналитического решения задача не имеет.
  1. Rhombus
    19.03.2019 11:54
    #19907584
    Я заметил, что это численное моделирование, ну и что? Это не отменяет того, что я написал.
    
    Немного оффтоп, но вспомнилась забавная цитата на тему точных решений:
    
    Разумной отправной точкой при рассмотрении проблемы многих тел мог бы служить следующий вопрос: сколько нужно тел, чтобы возникла проблема?
    
    Проф. Дж. Браун отмечал, что лица, интересующиеся точными решениями, могут найти ответ на этот вопрос, заглянув в историю.
    
    Для ньютоновской механики XVIII века задача трех тел была неразрешимой.
    
    С рождением общей теории относительности (где-то около 1910 г.) и квантовой электродинамики (1930 г.) стали неразрешимыми задачи двух и одного тела.
    
    В современной квантовой теории поля неразрешимой оказывается задача нуля тел (вакуума).
    
    Так что если мы интересуемся точными решениями, то ни одного тела — это уже слишком много.
    
    Из книги "Фейнмановские диаграммы в проблеме многих тел"
  1. Psychopompe
    19.03.2019 16:43
    #19909404
    Есть частные аналитические решения (и не одно).

amarao
19.03.2019 11:54
#19907586
Скажите, ваша точность достаточна для того, чтобы рискнуть жизнью и сказать «не сушитесь» иноплатенянам?
1. AndrewSu Автор
  19.03.2019 23:10
  #19911370
  Главное не «не сушитесь», а «засушиться» вовремя. Планету двигать они тоже не научились.
  Точности не хватит если система станет слишком жёсткой, а таковой она будет если планета будет пролетать слишком близко к одной из звёзд. В этом случае всё равно всё сгорит.
  Вообще удивительно, что трисолерианцы потеряли свою планету искусственным образом.
  Мне кажется, что при описанных в книге траекториях планету должно было вышвырнуть из системы или поглотить одно из солнц.
  1. amarao
    20.03.2019 10:56
    #19913284
    По книге сказано, что все остальные планеты так и погибли, это была последняя оставшаяся.

HEKOT
19.03.2019 12:16
#19907696
На 10-й секунде очень напоминает существующую галактику Колесо телеги.

Всё-таки, визуально, а не структурно. У Вас же в модели нет молодых ярких звёзд?

excoder
19.03.2019 12:51
#19907878
А как в романе объяснялась эта нестыковка — сверхвысокоразвитая цивилизация, способная захватывать другие планеты, неспособна численно устойчиво решать задачу трёх тел? Пусть и хаос, но — время Ляпунова и т.д., всё что можем вычислить — выглядит довольно просто даже для человечества. По мне это, честно говоря, ощущается как предпосылка сюжетного примитивизма.
1. TerraHominem
  19.03.2019 13:05
  #19907964
  На самом деле, описанное автором, не решает проблему для жителей Трисоляриса в системе звезд, которая в реальности является — Альфой Центаврой.
  В свое время Брунс и Пуанкаре доказали, что систему дифференциальных уравнений для трех тел невозможно свести к интегрируемой, динамические системы, судя по их открытию — не изоморфны. Получается, что эти интегрируемые системы допускают разложения на более мелкие, невзаимодействующие (из вики)
  1. excoder
    19.03.2019 20:25
    #19910540
    Почему не решает? В математике понятие интегрируемости ? решению системы в виде явного аналитического выражения. Для трёх тел это, действительно, в общем случае невозможно. Как и отсутствует символьное выражение для многих определённых интегралов. Отсюда и пошёл термин «интегрируемость», кстати. Не берутся некоторые интегралы, ну и интегралу всегда соответствует некоторый дифур. Ну и конечно, неингетрируемость системы, как и неинтегрируемость интегралов, в каком-либо явном символьном виде, не мешает решать их численно. Потому я и никак не избавлюсь от ощущения всего этого романа непростительной окологуманитарной ошибкой. Я уже и читать начал, и только пока укрепляюсь в этом впечатлении. Как, видимо, и комментаторы ниже.
  1. excoder
    20.03.2019 00:57
    #19911810
    Гипотетические жители Проксимы-Б тихо захватывали Землю, не обнаружив в земном интернете подобных вот знаний, которых у них не было у самих? core.ac.uk/download/pdf/25201586.pdf :)
    
    TerraHominem
    20.03.2019 09:11
    #19912662
    Интересненько, а как эти моделирования гипотетических устойчивых орбит помогли бы для жителей Проксимы? Ну хорошо, они пытались сдвинуть свою планету (OMG) на другую орбиту. Вы книгу полностью читали? На каком вы месте?
    До мотивов непосредственно экспансии Земли Вы уже дошли? Или вы рассуждаете на бегу?
    В любом случае, вы можете отложить книгу и вообще не прикасаться к научной-фантастике, либо берите из тех областей, где не будете слишком требовательны. Попробуйте фантастику с биологическими системами.
    
    TerraHominem
    20.03.2019 09:18
    #19912680
    И я еще раз повторю, на сколько мне известно — аналитического решения задачи трех тел на данный момент нет. Есть только частные решения при заранее заданных параметрах. В общем случае есть только приблизительное решение Вейерштрасса и все. Более того, в игре, в которой происходит попытка определить общее решение — идет описания истории трисолянцев (?), а не текущее научное развитие.

TerraHominem
19.03.2019 12:57
#19907912
Больше половины не понял, но за работу лайк поставил (Вообще такое делаю не часто, но всегда обращаю внимание на способ изложения материала и мыслей)

Весь цикл «В память о прошлом Земли» — просто великолепен. Очень интеллектуальная научная-фантастика. Прекрасно логично описана концепция «Темный лес» — объясняющая парадокс Ферми, а тема с субатомными частицами (софоны) — искажающие результаты исследований в фундаментальной физике и мгновенно передающие информацию от Земли (находясь во всех ее точках практически мгновенно) до Альфа Центавры — просто шикарна.
1. Polaris99
  19.03.2019 14:55
  #19908600
  Ничего общего книга с научной фантастикой не имеет. Да и персонажи картонные, если не сказать хуже.
  1. TerraHominem
    19.03.2019 15:41
    #19908904
    А Вы книгу то читали?
    
    Polaris99
    19.03.2019 16:24
    #19909252
    А Вы? Как показывает практика, в восторге от нее лишь люди с гуманитарным образованием, которым достаточно найти в тексте пару научных терминов для того, чтобы с полным на то основанием занести книгу в раздел строгой научной фантастики. Точно такая же история, кстати, и с Гиперионом Симмонса.
    
    TerraHominem
    19.03.2019 16:30
    #19909302
    Я как-то недавно наткнулся на разбор в Популярной Механике этой книги с т.зр. науки и знаете, судя по вашим мыслям — разбор в крупном журнале научпоп делал гуманитарий.
    Если следовать Вашим суждениям дальше, то получается для Вас научная фантастика — это раздел теоретической физики. Видимо мы с вами по разному понимаем жанр научная-фантастика.
    Приведите пример строго научной фантастики, на ваш взгляд.
    
    Polaris99
    19.03.2019 16:46
    #19909422
    Почитал я эту ересь в Популярной механике, даже комментировать лень. Попытка натянуть сову на глобус. Причем сова в данном случае — книга китайского автора.
    
    TerraHominem
    19.03.2019 16:56
    #19909488
    А пример предоставить, книжки то?
    
    TheYellingChives
    19.03.2019 17:52
    #19909836
    Да тролль же, очевидно. Знает что книгой восхищается большинство любителей НФ, вот и накидывает.

vlanko
19.03.2019 13:45
#19908188
Эволюция из миллиона частиц уже напоминает гипотезу создания Солнечной системы.

inetstar
19.03.2019 13:55
#19908240
Спасибо за статью! Как по мне — тянет на диссертацию.

Вы это сделали ради хобби? Или всё-таки для научной работы или за деньги?

Мне хотелось бы заниматься чем-то подобным за деньги и будет жаль, если вы всё это проделали ничего не получив, кроме статьи на Хабре.
1. AndrewSu Автор
  19.03.2019 23:28
  #19911454
  +2
  Вы это сделали ради хобби? Или всё-таки для научной работы или за деньги?
  
  Хобби. Для меня это просто интересно, даже не научная работа, не говоря уже про деньги.
  
  Мне хотелось бы заниматься чем-то подобным за деньги и будет жаль, если вы всё это проделали ничего не получив, кроме статьи на Хабре.
  Я получил разносторонний опыт. Но на основной работе проекты не менее интересные.
1. Bas1l
  20.03.2019 02:15
  #19911964
  +1
  Что-то такое, мне кажется, делают специальные команды разработчиков, занятые в суперкомпьютерных центрах. Приходят к ним физики, к примеру, и говорят: "распараллельте, пожалуйста, для CUDA наш код для моделирования гравитационного взаимодействи" (gpu сейчас стоят на многих суперкомпьютерах), и вот в идеале эти разработчики что-то такое должны сделать.

spam312sn
19.03.2019 13:59
#19908278
Восхитительно! Результаты моделирования завораживают

DimPal
19.03.2019 15:57
#19909036
А как изменится картина если предположить что гравитационные взаимодействия имеют скорость распространения? Свет и звук, например, имеют же скорость распространения…
1. SpaceEngineer
  19.03.2019 16:00
  #19909068
  Никак, скорость звёзд, даже самых быстрых — порядка 100 км/с, это в тысячи раз меньше скорости света.

SpaceEngineer
19.03.2019 15:59
#19909054
Моделирование эволюции галактики, состоящей из миллиона звёзд. В центре тело с массой 99% от общей.

Так это у вас не галактика, а зарождающаяся планетная система со звездой в центре и протопланетным диском. В реальных галактиках масса центральной чёрной дыры составляет 1-2%, а около половины — тёмная материя. Её надо моделировать теми же частицами (распределив изначально эллиптически), но не рендерить (или рендерить другим цветом).
Зарождение планет тоже интересно моделировать, но тогда нужно ввести возможность слипания частиц, и газодинамику (SPH например). Есть планы на это?

OnlySlon
19.03.2019 16:29
#19909298
Взрываем мою писанину в браузере.
http://tech-papa.com/hz.htm
Звезд поменьше чем у автора. Галактик побольше. :)

sim31r
19.03.2019 17:12
#19909594
В процессе моделирования порядка 5% 'звёзд' покинуло начальную область безвозвратно.

Интересно, а какие там ускорения возникают? Если звездная система с планетами, планеты останутся на орбитах вокруг своей звезды или они потеряют связь со звездой и звездная система разрушится? Что далее будет с потерянными планетами (какой-то процент будет потерян неизбежно), останутся свободными, или перейдут к другим звездам.
Думаю вычисления это не затруднит. Добавить к изначальной модели на 60 000 звезд, еще 1000 малых тел, которые ни на что влиять не будут, из-за малой гравитации, но сами будут подвержены влиянию гравитационных возмущений. Тела наподобие тел солнечной системы: Меркурия, Юпитера, Плутона, Цереры.

excoder
20.03.2019 01:01
#19911824
Благодарю за столь глубокую проработку темы. Если вы решите продолжить работать с темой, может быть, стоит попробовать вот этот подход: habr.com/ru/post/358352? Не исключено, что расширение горизонта прогнозирования на 1 порядок с помощью нейросетей позволит центаврянам отказаться от захвата землян в ближайшем будущем. :)

Bas1l
20.03.2019 03:03
#19912008
Это просто титанический труд! Спасибо за прекрасную статью!

У меня возникло несколько вопросов:
1. в моделировании гидродинамики есть несколько стандартных систем, с которыми сравнивают работу численных методов. Например конвекция Рэлея-Бенара и еще десяток других. У них необязательно есть полное аналитическое решение, но, к примеру, аналитически известна скорость диссипации энергии. Ну и всегда можно проверять на законы сохранения (притом известны типичные скорости расхождения с законами сохранения). Я не специалист в моделировании гравитации, но наверняка там есть что-то подобное. Вы не изучали этот вопрос и не делали сравнений? Первое, что приходит в голову, это однородное шарообразное облако точек. Точки, находящиеся на расстоянии R от центра испытывают притяжение как будто бы только от точек внутри сферы R (M = 4/3 pi R^3 \rho), гравитация от точек вне сферы компенсируется. Если точки двигаются с такой скоростью, чтобы точно компенсировать гравитацию, они должны бесконечно долго находиться на своих орбитах (может быть, точки можно расположить парами диаметрально противоположно, чтоб центр масс любой подсистемы радиуса R никогда не двигался). Скорость можно вывести из F = G m 3/3 pi R^3 \rho / R^2 = A R = a m = m v^2 / R. То есть точки должны вращаться с одинаковой по модулю угловой скоростью :-) (впрочем, не обязательно одинаковой по направлению). Качество симуляции может определяться тогда, к примеру, относительным отклонением орбиты от идеальной стационарной круговой орбиты за период (к примеру, < 10^-13 за период).
2. скажите, пожалуйста, откуда вы черпали вдохновение (кроме источников в конце, если были еще источники)? Особенно для представления Kd-дерева через "бинарную кучу" и stackless обхода геометрически близких узлов? спасибо!
3. если эти оптимизации (в таком виде/такой комбинации) не используются/не известны для моделирования задачи N тел (или используются, но не на GPU), то, мне кажется, можно запросто опубликовать отличную статью даже в каком-нибудь хорошем англоязычном журнале. Например, Journal of Computational Physics или Computer Physics Communications или даже более широкого профиля (или более специализированном). Аудитория у них меньше, но зато это люди, которым это реально может быть интересно. Англоязычные журналы могут быть более лояльны по поводу аффилиации с университетом.
1. AndrewSu Автор
  21.03.2019 00:06
  #19917778
  +2
  Качество симуляции я определяю только по основным законам сохранения, но дополнительно проверить на задачах с известным решением это здравая идея.
  
  Про обход деревьев на GPU я начал смотреть здесь. Kd-дерево — бинарное, хранение бинарного дерева по принципу кучи само напрашивается, даже не знаю где я это увидел. Основу идеи обхода дерева без стека я взял из формул для индексов элементов в куче. Но публикации на эту тему точно есть, но больше в области трассировки лучей. Про оптимизацию с использованием инструкции FFS я сам догадался, но что-то подсказывает, что такая инструкция появилась не зря, а возможно для подобных задач. До этого я про неё не слышал.
  
  Здесь нужно делать серьёзный обзор литературы, но без человека в теме я бы не стал. Я могу просто не знать слов для запроса в гугле, для нахождения искомых статей.
  1. Bas1l
    21.03.2019 01:43
    #19917970
    Спасибо большое!
    
    Я вот как раз тоже сам поискал и как раз похожие статьи про traversal нашел (но я просто k-d tree traversal погуглил): популярная со словом stackless (в ней как раз подходит картинка 2) и просто обзор.
    
    Ну и выше в комментах есь описание примера от nvidia, там просто много ссылок.
    
    Да, это, конечно, очень здорово, что вы такие идеи можете генерировать (про обход и FFS)!
    
    У меня еще несколько мелких вопросов:
    
    Когда вы пишете distance_sqr((v1 — m_mass_center[curr]).norm()), там не должен быть квадрат нормы? Обычно норма берет корень, а мы как раз хотим избежать этого и везде предполагаются квадраты
    
    В force_compute есть вычисление силы для относительно далеких узлов дерева. Но я все никак не нахожу вычисление силы для точек, которые близко, если условие distance_sqr > m_radius_sqr[curr] не выполняется
    
    Насчет обзора вы правы. Как идея, вы можете просто перевести эту статью на английский и выложить на хабре (и скинуть ссылку на реддит :-) ). Обидно будет, если такой труд будет доступен только русскоязычной части человечества..

santa324
20.03.2019 10:09
#19912978
А вы не думали вместо Barnes–Hut реализовать Fast multipole method?
По идее там достижим уровень в O(N) FMM_Tutorial
Хотя преимущества всего на миллионе частиц, может, и не будет.
1. AndrewSu Автор
  21.03.2019 00:12
  #19917794
  А вы не думали вместо Barnes–Hut реализовать Fast multipole method?
  Бегло просмотрел презентацию. Теперь думаю. Если нет подводных камней, то это заявка на миллиард частиц на одном GPU, хотя скорее в скорость доступа к памяти упрётся. Надо пробовать.
  1. AN3333
    21.03.2019 08:51
    #19918468
    Я правильно понимаю что Barnes–Hut работает для неоднородной среды, для конденсированной он не будет?
  1. AN3333
    21.03.2019 08:58
    #19918482
    Не знаю как на миллиард, а вот, если меньше миллиона то подводный камень FMM в том, что высоки начальные расходы, так что говорят, что до миллиона он бессмысленен несмотря на линейное поведение.

trapwalker
20.03.2019 15:49
#19915592
А разве суть проблемы в задаче трёх тел и одноименной книге не в том, что из-за неидеальной точности получившийся аттрактор в какой-то момент момент (в некотором бутылочном горлышке) начинает вести себя радикально иначе?
Попробуйте посчитать через какое время результаты разных способов расчета начинают расходиться радикально. Неужели все рассмотренные методы к концу симуляции получают одну и ту же картину?
Вот для случая с тремя звёздами и планетой эта симуляция нужна была, чтобы жители могли на каком-то достаточно продолжительном промежутке времени предсказать смену сезонов. Мне кажется основные проблемы там были не из-за постепенной «потери энергии», а из-за того. что модель переставала показывать правильно сезоны.

Vitter
21.03.2019 02:36
#19918076
исследования численных методов — это замечательно.
Только вот для симуляции галактик не хватает этак… 75% массы (тёмной материи), которые расположены примерно симметрично в гало (в сфере) этих галактик

Задача N тел или как взорвать галактику не выходя из кухни +165

Физика явления

Математика

Итог

CPU реализация

GPU реализация

OpenCL

CUDA

Аппроксимация функции

Обход дерева на GPU

Заключение

Визуализация

Комментарии (96)

AndrewSu Автор

AndrewSu Автор

AndrewSu Автор

AndrewSu Автор

AndrewSu Автор

AndrewSu Автор

AndrewSu Автор

AndrewSu Автор

AndrewSu Автор

AndrewSu Автор

AndrewSu Автор

AndrewSu Автор

AndrewSu Автор

AndrewSu Автор

AndrewSu Автор

Аппроксимация функции $inline$