Как стать автором
Обновить

Комментарии 127

В отличии от нас

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

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

Да и в вычислительном плане тоже можно не просто любой метод решения дифф-уров брать с потолка, есть понятия «мягкой» и «жесткой» задач, чисел обусловленности, сходимости получившейся разностостной схемы и прочего подобного…
Но это для серьезных расчетов, а с тем чтобы просто поиграться — отлично.
Неплохо бы и привести сравнение с n-body из sdk Nvidia. Скорость, качество и тд. У них во всяком случае gpu-шная версия быстрее CPU на пару порядков. А так хорошая работа… Просто картина не полная.
Посмотрел код примера из 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 и сравнить более детально.
Посмотрел код примера из CUDA SDK. Вы этот имели ввиду?

Ага, именно ее.

Непонятно что такое 'portable C implementation',

Обычно под этим имеют в виду однопоточную и не заточенную под архитектуру конкретного процессора. Сейчас компиляторы значительно умнее тех что были в 2006.
В последних версиях CUDA SDK появился модуль OptiX, предназначенный для трассировки лучей. В нём есть реализация обхода деревьев на GPU. Очень интересно сравнить с этой реализацией. Может кто-то уже сталкивался?
Провёл сравнение, и не зря!
До рабочей станции не пошёл, проверял на ноутбуке.
Intel® Core(TM) i5-3230M CPU @ 2.60GHz
GeForce GT 740M
Ubuntu 16.04
gcc 5.4.0

Расчёты проводил с двойной точностью, как и в статье, для 65536 тел.
Результаты примера из CUDA-SDK
#./nbody -fp64 -numbodies=65536 -benchmark
4.8549 s
#./nbody -fp64 -numbodies=65536 -benchmark -cpu
57.7650 s

Как видим, слабый GPU, обгоняет слабый CPU в 12 раз. Стоит отметить, что по-умолчанию пример из SDK собирается без поддержки OpenMP, и многопоточность пришлось включать отдельно флагами компиляции.
Результаты моей реализации (у меня создаётся две 'галактики' по 32768 тел):
#./nbody-simulation --block_size=64 --engine=openmp --solver=euler --stars_count=32768 --max_time=3 --min_step=0.1 --max_step=0.1 --check_step=1 --verbose=1 --dump_step=0
25.012 s
#./nbody-simulation --block_size=64 --engine=block  --solver=euler --stars_count=32768 --max_time=3 --min_step=0.1 --max_step=0.1 --check_step=1 --verbose=1 --dump_step=0
10.1567 s
#./nbody-simulation --block_size=64 --engine=cuda   --solver=euler --stars_count=32768 --max_time=3 --min_step=0.1 --max_step=0.1 --check_step=1 --verbose=1 --dump_step=0
6.54931 s
#./nbody-simulation --block_size=64 --engine=cuda_bh_tex --solver=euler --stars_count=32768 --max_time=3 --min_step=0.1 --max_step=0.1 --check_step=1 --verbose=1 --dump_step=0 --tree_layout=heap_stackless --distance_to_node_radius_ratio=4
0.489003 s


Выглядит интересно.
Во первых, моя реализация взаимодействия каждый с каждым (--engine=cuda) в полтора раза медленнее реализации из CUDA-SDK. Я посмотрел внимательнее код примера и нашёл ещё одно отличие, вместо 1/sqrt(x) они применяют функцию rsqrt(x), которая немного экономит флопы ценой точности.
Отрывок из NVIDIA_CUDA_C_Programming_Guide_10.0.pdf из которого я делаю вывод, что rsqrt - аппроксимация, хотя в других местах про это явно не сказанно.
To preserve IEEE-754 semantics the compiler can optimize 1.0/sqrtf() into rsqrtf()
only when both reciprocal and square root are approximate, (i.e., with -prec-
div=false and -prec-sqrt=false ). It is therefore recommended to invoke rsqrtf()
directly where desired.

Я заменил у себя 1/sqrt(x) -> rsqrt(x) и время работы наших реализаций сравнялось.
#./nbody-simulation --block_size=64 --engine=cuda   --solver=euler --stars_count=32768 --max_time=3 --min_step=0.1 --max_step=0.1 --check_step=1 --verbose=1 --dump_step=0
4.7297 s


Во вторых, CPU реализация из CUDA-SDK, в примерно в два раза отстаёт от аналогичной моей (--engine=openmp) и почти в шесть раз от оптимизированной (--engine=block).
Особенно понравилость в CPU реализации из CUDA-SDK
T invDist = (T)1.0 / (T)sqrt((double)distSqr);

В коде для GPU используют rsqrt, а на CPU всегда считают с двойной точностью.


В третьих, реализация с вычисления силы взаимодействия, использующая kd-дерево (--engine=cuda_bh_tex), обходит любую из выше перешисленных 'плотных' реализаций на порядок.
Для раскрытия потенциала GPU необходимо использовать одинарную точность везде где только можно. В AMBER для этого применяется т.н. «hybrid single / double / fixed precision model», они статью по этой теме публиковали. Например, для 2080Ti выигрыш выходит 30 раз (float vs double).
Такой прирост неудивителен. Ведь с одинарной точностью 2080Ti почти 12 Тфлопс выдаёт, а с двойной только 0.36 Тфлопс. Игровая карта. Интересно статью посмотреть, были ли просадки в точности вычислений, и на каких операциях они экономили. Я у себя не пытался сравнивать по точности выходного результата. Может быть вы видели, можно ли под такие расчёты использовать новые RT ядра этих карточек?
Где это критично — точность не теряется благодаря использованию разной точности в разных местах, в т.ч. fixed precision arithmetics в некоторых местах. Про RX ядра не слышал, но если интересно — смотрите AMBER mailing list и публикации от авторов AMBER. Дело в том, что в разработке AMBER участвуют инженеры Nvidia, так что скорее всего именно они первыми эту тему попытаются раскрыть, если это вообще имеет смысл.
точность не теряется
Там не совсем это написано.
Согласно таблице 4, точность теряется. Авторы больше сравнивают модели вычисления силы с одинарной точностью и целочисленным 64 битным аккумулятором (SPFP), с аккумулятором двойной точности (SPDP) с эталонной моделью, которая всё вычисляет с двойной точностью (DPDP).
Обе модели SPFP и SPDP теряют в точности, по сравнению с DPDP, 3-5 десятичных знаков, но имеют прирост в производительности ~10 раз на карточках, рассматриваемых в статье. Модель SPFP на 60-80% производительнее чем SPDP, при схожих ошибках. Насколько эти ошибки критичны для результата моделирования особо не раскрывается.
Но всё равно интересно, спасибо!
Да, потому я написал про «где критично». Обоснования некритичности приводятся в двух других статьях.
we have already shown [9] that these errors are sufficiently small such that the SPDP GPU implementation generates trajectories with numerical stability that is equivalent to the CPU implementation. It is thus to be expected that numerically stable simulations are also possible with the SPFP precision model.
A thorough analysis showing that numerically stable MD simulations using the PME algorithm are possible with the SPDP and SPFP precision models will be presented elsewhere[19]
We thus recommend the GPU implementation of PMEMD to be used with the SPFP precision model to reduce thecomputational cost compared to the SPDP and DPDP precision models while maintaining numerical stability that is equivalent to the reference double precision CPU implementation

Generalised Born [9] model немного ближе к Вашей задачe, наверное.
Насколько эти ошибки критичны для результата моделирования особо не раскрывается.

Это тонкий вопрос.


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


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


Со слов одного из разработчиков пакета LAMMPS, для расчётов дальнодействующих сил методами, основанными на преобразовании Фурье — двойная точность может быть принципиальна, т.к. там идёт суммирование большого числа разнознаковых слагаемых. Ещё чувствительнее будет давление, т.к. там ещё одно суммирование разнознаковых слагаемых. Но даже в этом случае система будет вести себя "правдоподобно", т.е. энергия будет сохраняться, интегрирование будет устойчивое и т.п.

Все верно, как раз хотел то же самое написать. Добавлю, что под методом, основанными на преобразовании Фурье в MD подразумевается PME. У PME есть альтернатива — Reaction field method (RF). Наиболее быстрый известный мне GPU код, в котором реализован RF — OpenMM. В случае OpenMM RF вариант 1.5 раза быстрее чем PME, при том что PME там реализован на уровне AMBER, то есть очень оптимально. Там ещё есть бенчмарк, который показывает, что Titan X в 25 раз быстрее чем i7-7740X (4x4.3GHz) CPU. А 2080Ti будет ещё быстрее чем Titan X.
Правильно писать: чтобы получить существенное ускорение на игровой видеокарте, нужна задача, для которой не критична разница между двойной и одинарной точностью в вычислительно-интенсивной части.
Классическая молекулярная динамика с этим жить может, гидродинамика — не всегда, в квантовых расчётах — не ниже двойной точности.
Ещё одна особенность. В Makefile для примера nbody из CUDA-SDK отсутствовали флаги оптимизации. Я добавил -O3, стало сильно лучше.
#./nbody -fp64 -numbodies=65536 -benchmark -cpu
20.8 s
Вот они, прелести маркетинга. Заменить в бенчмарке видеокарты 1.0/sqrt(x) на rsqrt(x) и не включать флаги оптимизации на CPU — и двукратное ускорение превращается в порядок.
Хотя, если сравнивать с ноутбучной (и любой другой «игровой») видеокартой, честнее в одинарной точности смотреть — двойная же там урезана очень сильно.
А как же черная материя и/или энергия? ;-)
Тёмная она :) Или скорее мы. Пока.
Так это по сути эволюция темноматериальной системы и есть — учитывается только гравитационное взаимодействие. Чтобы в первом приближении прикинуть как это с тёмной материей выглядит, надо распределить равномерно по объему, 9/10 точек не отображать, а для отображаемых ввести небольшую потерю энергии при близком прохождении друг к другу.
Низкая скорость — синий цвет, средняя — зелёный, высокая — красный. Видно, что в центре скорость выше, и плавно убывает к краям, а самая низкая в экваториальной плоскости. В центре находится тело с массой 99.9% от общей


А как же темная материя и вот это все? Со схожими скоростями вещества по всему объему.
Я так понимаю, что раз уж здесь не рассматриваются столкновения, а все массы — это математические точки, то иного взаимодействия, кроме гравитационного, нет. Это значит, что автор рассмотрел исключительно тёмную материю. Казалось бы, да? Но нет. Тёмная материя не собирается в кусочки именно потому, что не взаимодействует механически. Это значит, что в задаче появляется неимоверно огромное количество тел очень малой массы. Такое мы уж никак не смоделируем.
Не так давно я прочёл фантастический роман «Задача трёх тел» Лю Цысиня

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

Я как-то видео работу — на конференции по неустойчивостям, так там так и делали — брали уравнение Навье-Стокса, добавляли гравитационную силу и получали спиральные рукава галактик, т.е. множество звезд моделировалось сплошной средой. (Хотя возможно это и распространенный метод в астрофизике, не знаю).
А какие у вас заданы абсолютные размеры/массы/времена, т.е. соответствует ли видео в каком-то виде реальности или не особо?
Три тела равной массы, скорости заданы такими, чтобы суммарный импульс был равен нулю

А теперь и Вы напомнили...

Каждому хочется почувствовать себя демиургом!
Хунта был великолепным таксидермистом. Штандартенфюрер, по словам Кристобаля Хозевича, — тоже. Но Кристобаль Хозевич успел раньше.
Я сделал свою модель до того, как прочитал книгу, но мне это не помогло написать статью первым. :) github.com/dsdante/constel

Ещё есть идея решить задачу методом самосогласованного поля.
Ещё есть идея решить задачу методом самосогласованного поля.
Очень похоже на метод Ахмада-Коэна для решения задачи N-тел.
Первое моделирование было проведено как раз с применением этого метода:
image
Думается это вполне можно опубликовать в журнале типа «Вычислительные методы и программирование».

Я однажды хотел опубликоваться в этом журнале, но был ответ, что публикация возможна только для тех, кто учится/работает в университете. А я уже закончил. Пришлось в другом месте публиковать. Сейчас, возможно, правила изменились.

А какие у вас заданы абсолютные размеры/массы/времена, т.е. соответствует ли видео в каком-то виде реальности или не особо?

Из физических констант задаю только G=1, а начальный диаметр 'галактики' задаю равным 100. Массы все одинаковые кроме центрального тела. К реальным величинам пересчитать возможно, но я этого пока не делал.
Большое спасибо за интересную статью!
Если я правильно помню, метод Верле более или менее сохраняет энергию (по сравнению с методом Рунге-Кутта, несмотря на более высокую точность последнего).

Вообще, численное интегрирование задачи N тел — очень изученная тема, только в другом приложении, не астрономическом: в физике это принято называть молекулярной динамикой.
Совершенно верно. Есть куча методов сохраняющих энергию.
И да, Верле самый знаменитый. Он ее сохраняет точно (видимо, с точностью до количества машинных знаков).

И методы высоких порядков тоже имеются. Я пользуюсь VEFRL.
У таких алгоритмов, не смотря на их дискретность, имеется эффективный гамильтониан. При увеличении шага падает точность, гамильтониан все дальше уходит о оригинала, а вот энергии деваться некуда.
Я пользуюсь VEFRL.

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

Серьёзно… Не подскажете где?
А то может и я начну потенциалами приторговывать.
в смысле, «хорошего»? Существующие решения достаточно хороши что бы весьма широко применяться.
В астрофизике тоже. Только там он зовется SPH.
На самом деле с сохранением энергии не всё так просто, даже у метода Верле. Но метод хороший, я его планировал реализовать, но пока руки не дошли.
По вашей ссылке
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.

Честно говоря, для меня это новость. Я не вижу повышения энергии, по крайнем мере за часы гоняния программы. Может быть тут проблема с ограниченной машинной точностью, а не в самом алгоритме, если этот эффект вообще есть?
мой небольшой старый опыт с методами типа Рунге-Кутты показывает, что энергия только падает. Возможно, зависит от метода и железа.
Забавно, а я думал что она только растет, для не симплектических интеграторов :)
Есть уже готовые универсальные реализации MD с портами на CUDA и OpenCL и mpi. Из наиболее мейнстримных — AMBER(Rutgers Uni, Nvidia и ко), GROMACS (Stockholm University и многие другие), OpenMM (Stanford). Эти движки используются в production на тысячах суперкомпьютеров и кластеров по всему миру, как в академии так и коммерческими компаниями. OpenMM и AMBER имеют примерно равно высокую производительность на GPU/CUDA, оба движка open source, но OpenMM ещё и полностью бесплатный. Для Вашей задачи, к сожалению, совсем без изменений эти движки использовать не получится. Масса в MD присутсвует, но это инерционная масса, гравитационное взаимодействие в MD обычно игнорируют. Можно задействовать «электростатическую» часть MD потенциала, у гравитации и кулона одинаковая форма потенциала, так что можно было бы гравитационное взаимодействие выразить через «заряд», но проблема в том, что у гравитационной массы нет знака и отталкивания, и MD движку это сложно «объяснить» — нужно лезть в код. Но зато там state of the art алгоритмы и супер оптимизированные их реализации, к примеру reaction field method. Правда там и много лишнего для данной задачи — bonded interactions, куча разных термостатов, SHAKE и десятки-сотни других функций. Вероятно, решения из астрономии типа Gravidy и Bonsai подойдут лучше, они и релятивистские поправки могут учитывать.
{не проснувшись} для подобных задач используют методы частиц в ячейках.
Они снижают вычислительную сложность упуская детали локальной динамики парных взаимодействий.
С одной, стороны подобную программу я написал еще в 10 классе на паскале, применяя куда как меньше математики и железа.

С другой — ваша «модель» не учитывает тёмную материю и энергию (или что там стоит за этими абстракциями) и потому — неверна.
Позволю себе вмешаться.
В общем случае модель не может быть верной или не верной. Она может описывать результаты экспериментов с некой точностью в определенных условиях.
Темная материя была введена/придумана как раз чтобы хоть как-то объяснить расхождения между имеющимися моделями и экспериментами.
Кстати — а как качественно должны измениться результаты представленного моделирования если ввести тёмную материю? Я не специалист по астрофизике, но если вы ориентируетесь, то думаю на пальцах показать можете.
Недавно читал, извините ссылку не помню, что при введении в модель темной материи как раз начинает получаться что-то похожее. Такая вот сетка с перемычками из спиральных галактик.
За что купил, за то продаю :)
Скорость звезд должна перестать убывать с радиусом
Никак, потому что это — симуляция системы гравитационно взаимодействующих тел. Темная материя — всего лишь вид материи, взаимодействующей исключительно гравитационно, поэтому данная симуляция полностью применима к ней. Для приближения к «реальности» достаточно выключить отображение 80% частиц.

Большее приближение к реальности даст использование теории относительности вместо уравнения Ньютона. Темная энергия, по идее, автоматически при этом появится.
Большое спасибо за статью, интересно и поучительно. Интересно на сколько это далеко от реальности? Например из-за скорости распространения гравитации? Или из-за релятивистских эффектов. Интуиция подсказывает что картина должна быть кардинально иной, но интересно насколько?.. Не знаю как остальных, но но крайней мере меня бы очень порадовало релятивистское продолжение статьи.
НЛО прилетело и опубликовало эту надпись здесь
Для ошибочек ctrl+Enter ведь добавили. Не все ещё знают видать (или привычка многолетняя может делает своё дело).
НЛО прилетело и опубликовало эту надпись здесь
Спиральные рукава и так есть. Просто они не так сильно 'визуально' выражены. В реальных галактиках они выглядят сильно ярче как раз из-за активного звёздообразования.
Ошибочку поправил, спасибо.
Но сообщать лучше с личку.
НЛО прилетело и опубликовало эту надпись здесь
Спасибо, интересная статья!
Скажите, Вы не пробовали оценивать не отклонение полной энергии и импульса системы по окончанию моделирования, а сравнить результаты у 2х моделирований, одно из которых имеет точность с огромным запасом? При этом интересно отклонение не средних значений системы и не общей энергии и импульса, а максимальное отклонение для отельного тела среди всех тел.
Вероятно, при вычислении с низким и высоким порядком точности именно это и сравнивается? Как я понимаю, там сравниваются координаты частицы полученные на следующем шаге для 2х разных точностей. Но не понятно какое отклонение считать малым…
Дело в том что если вся система «равномерная» но в ней есть компактная подсистема из всего нескольких тел но находящихся в очень жестком взаимодействии, то даже если эта подсистема полностью «разойдется» (так как для нее точность неприемлемо низкая) — то на общей энергии/импульсе системы это практически не скажется. Если же сравнивать координаты частиц при разных точностях, тогда не однозначно что является малым отклонением — незначительная погрешность на большом расстоянии от массивного тела может быть неприемлемо огромной в момент пролета вблизи него. Возможно можно сравнивать еще отклонения по импульсу и энергии отдельной частицы… Было бы очень интересно посмотреть на результаты таких экспериментов.
Вероятно, при вычислении с низким и высоким порядком точности именно это и сравнивается? Как я понимаю, там сравниваются координаты частицы полученные на следующем шаге для 2х разных точностей.

Да, именно так. У меня сравниваются не только координаты, но и скорости.
Но не понятно какое отклонение считать малым…

Этот параметр является настраиваемым, также как и другие.

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

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

Вот как раз интересно как его настраивать, как понять что такое-то значение достаточно маленькое и не надо в настройках поставить еще меньше?
Нужно взять систему с известной эволюцией, опробовать различные параметры на ней, и определить границы применимости. Потом можно приступать к исследованию целевой системы. В нашем случае можно взять задачу Кеплера, попробовать разные траектории и сравнить с аналитическим решением.
Так дело не только в точности вычислений, но и в точности исходных данных. А при большом числе тел могут быть такие конфигурации, когда незначительное изменение исходных данных приводит к огромному расхождению по прошествии некоторого времени. 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% в одну и в другую сторону.
Вопрос-предложение, где бы найти базу по звездам хотя бы нашей галактики, на сколько я понимаю она наиболее изученная в плане данных по координатам, скорости и массе звезд?

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

p.s. какое максимальное приближение звезд друг к другу наблюдается? интересно, на сколько высока вероятность, что пролетающая звезда сорвет одну из планет или изменит ее траекторию?
Самые свежие и точные данные на текущий момент — это GAIA DR2: gea.esac.esa.int/archive (но всех звезд, естественно, не будет ни в одной базе до тех пор как человечество не начнет путешествовать по Галактике)
Эм, ну, реальные данные в принципе нарыть можно- воспользоваться, архивом телескопа Gaia, например. Но вот только есть нюанс — там данные о примерно 1,7 миллиардах звёзд (и это ещё предположительно порядка 1% от общего числа звёзд в нашей галактике). Правда, полные данные — координаты и скорости движения — есть только примерно для семи миллионов звёзд, но даже этого достаточно, чтобы при попытке прямого расчёта заставить любое доступное железо глубоко и надолго задуматься.
Чтобы подсунуть реальные данные в программу. Дело в том что начальные значения скорости в тестах в статье — нулевые, т.е. вся система испытывает невероятно сильные колебания и всплески именно из-за этого, но реалии таковы что звезды двигаются уже по относительно стабильным орбитам.

Если бы это было так, то всё звёзды просто упали бы на центральное тело. Скорости как раз выбраны так, что начальные орбиты 'почти' стабильные.

В том моделировании это видно, все движутся по эллиптическим орбитам:
image
Видео про симуляцию столкновения напомнило это: видео сравнения симуляции с фотографиями хаббла
у меня сложилось впечатление, что вы за несколько минут моделируете процессы длительностью в сотни миллионов лет.
какие начальные условия вы задаете? массы тел, их скорости, расстояния между ними? сколько шагов интегрирования делаете?
Насколько я помню, в задаче трех (и более) тел решения являются неустойчивыми (aka динамический хаос), так что польза от этих вычислений сомнительная, как мне кажется.

Я заметил, что это численное моделирование, ну и что? Это не отменяет того, что я написал.


Немного оффтоп, но вспомнилась забавная цитата на тему точных решений:


Разумной отправной точкой при рассмотрении проблемы многих тел мог бы служить следующий вопрос: сколько нужно тел, чтобы возникла проблема?


Проф. Дж. Браун отмечал, что лица, интересующиеся точными решениями, могут найти ответ на этот вопрос, заглянув в историю.


Для ньютоновской механики XVIII века задача трех тел была неразрешимой.


С рождением общей теории относительности (где-то около 1910 г.) и квантовой электродинамики (1930 г.) стали неразрешимыми задачи двух и одного тела.


В современной квантовой теории поля неразрешимой оказывается задача нуля тел (вакуума).


Так что если мы интересуемся точными решениями, то ни одного тела — это уже слишком много.


Из книги "Фейнмановские диаграммы в проблеме многих тел"

Есть частные аналитические решения (и не одно).
Возможно, память Вас подводит. Есть достаточно устойчивые численные решения, активно применяются, работают: N-body simulations of gravitational dynamics. Кроме того есть ещё очень схожая задача молекулярной динамики, где тоже всё прекрасно решается численными методами.
Скажите, ваша точность достаточна для того, чтобы рискнуть жизнью и сказать «не сушитесь» иноплатенянам?
Главное не «не сушитесь», а «засушиться» вовремя. Планету двигать они тоже не научились.
Точности не хватит если система станет слишком жёсткой, а таковой она будет если планета будет пролетать слишком близко к одной из звёзд. В этом случае всё равно всё сгорит.
Вообще удивительно, что трисолерианцы потеряли свою планету искусственным образом.
Мне кажется, что при описанных в книге траекториях планету должно было вышвырнуть из системы или поглотить одно из солнц.
По книге сказано, что все остальные планеты так и погибли, это была последняя оставшаяся.
НЛО прилетело и опубликовало эту надпись здесь
А как в романе объяснялась эта нестыковка — сверхвысокоразвитая цивилизация, способная захватывать другие планеты, неспособна численно устойчиво решать задачу трёх тел? Пусть и хаос, но — время Ляпунова и т.д., всё что можем вычислить — выглядит довольно просто даже для человечества. По мне это, честно говоря, ощущается как предпосылка сюжетного примитивизма.
На самом деле, описанное автором, не решает проблему для жителей Трисоляриса в системе звезд, которая в реальности является — Альфой Центаврой.
В свое время Брунс и Пуанкаре доказали, что систему дифференциальных уравнений для трех тел невозможно свести к интегрируемой, динамические системы, судя по их открытию — не изоморфны. Получается, что эти интегрируемые системы допускают разложения на более мелкие, невзаимодействующие (из вики)
Почему не решает? В математике понятие интегрируемости ≡ решению системы в виде явного аналитического выражения. Для трёх тел это, действительно, в общем случае невозможно. Как и отсутствует символьное выражение для многих определённых интегралов. Отсюда и пошёл термин «интегрируемость», кстати. Не берутся некоторые интегралы, ну и интегралу всегда соответствует некоторый дифур. Ну и конечно, неингетрируемость системы, как и неинтегрируемость интегралов, в каком-либо явном символьном виде, не мешает решать их численно. Потому я и никак не избавлюсь от ощущения всего этого романа непростительной окологуманитарной ошибкой. Я уже и читать начал, и только пока укрепляюсь в этом впечатлении. Как, видимо, и комментаторы ниже.
Гипотетические жители Проксимы-Б тихо захватывали Землю, не обнаружив в земном интернете подобных вот знаний, которых у них не было у самих? core.ac.uk/download/pdf/25201586.pdf :)
Интересненько, а как эти моделирования гипотетических устойчивых орбит помогли бы для жителей Проксимы? Ну хорошо, они пытались сдвинуть свою планету (OMG) на другую орбиту. Вы книгу полностью читали? На каком вы месте?
До мотивов непосредственно экспансии Земли Вы уже дошли? Или вы рассуждаете на бегу?
В любом случае, вы можете отложить книгу и вообще не прикасаться к научной-фантастике, либо берите из тех областей, где не будете слишком требовательны. Попробуйте фантастику с биологическими системами.
И я еще раз повторю, на сколько мне известно — аналитического решения задачи трех тел на данный момент нет. Есть только частные решения при заранее заданных параметрах. В общем случае есть только приблизительное решение Вейерштрасса и все. Более того, в игре, в которой происходит попытка определить общее решение — идет описания истории трисолянцев (?), а не текущее научное развитие.
аналитического решения трех тел на данный момент нет
ну и что? Есть же численное, которое можно получить с любой требуемой точностью (в разумных пределах).
систему дифференциальных уравнений для трех тел невозможно свести к интегрируемой
Вы пропустили уточнение: «в виде конечных аналитических выражений», а численными методами всё решается, что и описывает автор.
Больше половины не понял, но за работу лайк поставил (Вообще такое делаю не часто, но всегда обращаю внимание на способ изложения материала и мыслей)

Весь цикл «В память о прошлом Земли» — просто великолепен. Очень интеллектуальная научная-фантастика. Прекрасно логично описана концепция «Темный лес» — объясняющая парадокс Ферми, а тема с субатомными частицами (софоны) — искажающие результаты исследований в фундаментальной физике и мгновенно передающие информацию от Земли (находясь во всех ее точках практически мгновенно) до Альфа Центавры — просто шикарна.
НЛО прилетело и опубликовало эту надпись здесь
А Вы книгу то читали?
НЛО прилетело и опубликовало эту надпись здесь
Я как-то недавно наткнулся на разбор в Популярной Механике этой книги с т.зр. науки и знаете, судя по вашим мыслям — разбор в крупном журнале научпоп делал гуманитарий.
Если следовать Вашим суждениям дальше, то получается для Вас научная фантастика — это раздел теоретической физики. Видимо мы с вами по разному понимаем жанр научная-фантастика.
Приведите пример строго научной фантастики, на ваш взгляд.
НЛО прилетело и опубликовало эту надпись здесь
А пример предоставить, книжки то?
Да тролль же, очевидно. Знает что книгой восхищается большинство любителей НФ, вот и накидывает.
НЛО прилетело и опубликовало эту надпись здесь
Эволюция из миллиона частиц уже напоминает гипотезу создания Солнечной системы.
Спасибо за статью! Как по мне — тянет на диссертацию.

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

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

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

Что-то такое, мне кажется, делают специальные команды разработчиков, занятые в суперкомпьютерных центрах. Приходят к ним физики, к примеру, и говорят: "распараллельте, пожалуйста, для CUDA наш код для моделирования гравитационного взаимодействи" (gpu сейчас стоят на многих суперкомпьютерах), и вот в идеале эти разработчики что-то такое должны сделать.

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

Восхитительно! Результаты моделирования завораживают

А как изменится картина если предположить что гравитационные взаимодействия имеют скорость распространения? Свет и звук, например, имеют же скорость распространения…
Никак, скорость звёзд, даже самых быстрых — порядка 100 км/с, это в тысячи раз меньше скорости света.
Моделирование эволюции галактики, состоящей из миллиона звёзд. В центре тело с массой 99% от общей.

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

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

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

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

Это просто титанический труд! Спасибо за прекрасную статью!


У меня возникло несколько вопросов:


  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. Качество симуляции я определяю только по основным законам сохранения, но дополнительно проверить на задачах с известным решением это здравая идея.
  2. Про обход деревьев на GPU я начал смотреть здесь. Kd-дерево — бинарное, хранение бинарного дерева по принципу кучи само напрашивается, даже не знаю где я это увидел. Основу идеи обхода дерева без стека я взял из формул для индексов элементов в куче. Но публикации на эту тему точно есть, но больше в области трассировки лучей. Про оптимизацию с использованием инструкции FFS я сам догадался, но что-то подсказывает, что такая инструкция появилась не зря, а возможно для подобных задач. До этого я про неё не слышал.
  3. Здесь нужно делать серьёзный обзор литературы, но без человека в теме я бы не стал. Я могу просто не знать слов для запроса в гугле, для нахождения искомых статей.

Спасибо большое!


Я вот как раз тоже сам поискал и как раз похожие статьи про 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] не выполняется

Насчет обзора вы правы. Как идея, вы можете просто перевести эту статью на английский и выложить на хабре (и скинуть ссылку на реддит :-) ). Обидно будет, если такой труд будет доступен только русскоязычной части человечества..

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

просто обзор
Спасибо, просмотрю.

там не должен быть квадрат нормы?
Там квадрат. Видимо, название функции не очень удачное.

Но я все никак не нахожу вычисление силы для точек, которые близко, если условие distance_sqr > m_radius_sqr[curr] не выполняется
Да, это совсем не очевидный момент. Но это условие выполняется только для близких точек. В оригинале оно выглядит так:
distance_to_node/node_radius < distance_to_node_radius_ratio
Тут я избавился от деления и квадратных корней, для этого все элементы m_radius_sqr домножены на квадрат distance_to_node_radius_ratio. У листьев дерева m_radius_sqr[curr] равен нулю, и это условие всегда для них истинно.
А вы не думали вместо Barnes–Hut реализовать Fast multipole method?
Бегло просмотрел презентацию. Теперь думаю. Если нет подводных камней, то это заявка на миллиард частиц на одном GPU, хотя скорее в скорость доступа к памяти упрётся. Надо пробовать.
Я правильно понимаю что Barnes–Hut работает для неоднородной среды, для конденсированной он не будет?
Авторы занимались астрофизикой, видимо, поэтому он больше там и распространён. Если я заменю гравитационный потенциал на Леннард-Джонса, то мой алгоритм не сломается. Но этот потенциал обычно обрывают на нескольких радиусах частицы, поэтому обход дерева надо немного иначе делать.
Не знаю как на миллиард, а вот, если меньше миллиона то подводный камень FMM в том, что высоки начальные расходы, так что говорят, что до миллиона он бессмысленен несмотря на линейное поведение.
А разве суть проблемы в задаче трёх тел и одноименной книге не в том, что из-за неидеальной точности получившийся аттрактор в какой-то момент момент (в некотором бутылочном горлышке) начинает вести себя радикально иначе?
Попробуйте посчитать через какое время результаты разных способов расчета начинают расходиться радикально. Неужели все рассмотренные методы к концу симуляции получают одну и ту же картину?
Вот для случая с тремя звёздами и планетой эта симуляция нужна была, чтобы жители могли на каком-то достаточно продолжительном промежутке времени предсказать смену сезонов. Мне кажется основные проблемы там были не из-за постепенной «потери энергии», а из-за того. что модель переставала показывать правильно сезоны.

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

Возможно, для преодоления проблем с не точностью нужно периодически повторно делать вычисления. И проблем не будет.

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

причем тут фильмы, если это в книге было?
проблема 3 тел именно в том и заключается, что долгосрочные прогнозы невозможны из-за очень быстро накапливающейся ошибки.
Повторные вычисления тут не помогут, это не даст возможности отдалить горизонт прогнозирования. Ошибка кроется даже не в точности вычислений, а в большей степени в точности измерений,и эту погрешность принципиально нельзя устранить.

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

Это если стандартные подходы применять. А нестандартными (вроде использования BigDecimal), думаю, никто не заморачивался, потому что Неуловимый Джо.

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

Думаю, что при правильном подборе измерительного оборудования вполне смогу предсказать результат подбрасывания монеты с вероятностью > 90%. Там ничего сложного, физическая модель элементарная — думаю, вполне достаточно замерить угловые и линейные скорости в момент броска по всем трём осям с точностью до пятого знака.

По классике же монету человек бросает, причём конкретным образом, снижающим устойчивость системы.

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

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

Насколько я понимаю проблема численного моделирования упирается в то, что сила (ускорение) изменяется пропорционально 1/r2, в то время как численное интегрирование линейно на конечном интервале. И каким бы маленьким мы его не выбрали, тела могут сблизиться на расстояние, где интервал численного вычисления начинает существенно влиять.
Аналитически проинтегрировать движение 3х точек тоже вроде не умеем.
Т.е. можно только более менее посчитать частные случаи типа вращения спутника на орбите Земли с учётом влияния Луны.
А прилёт ещё одной Луны на орбиту Земли посчитать уже не удастся.

Вы правы, но похоже солверы выше первого порядка решают проблему нулевого расстояния.
Для этих целей есть техника регуляризации уравнений движения (если интересно, можете поискать книгу «Gravitational N-Body Simulations: Tools and Algorithms», автор — Sverre Aarseth, там подробно изложены как возникающие проблемы, так и методы их решения).
Зарегистрируйтесь на Хабре, чтобы оставить комментарий

Публикации