О ВЫЧИСЛИТЕЛЬНОМ ТЕСТЕ ДЛЯ МОДЕЛИ АДИАБАТИЧЕСКОГО СЖАТИЯ ИДЕАЛЬНОГО БЕССТОЛКНОВИТЕЛЬНОГО ГАЗА Текст научной статьи по специальности «Физика»

УДК 531.3:519.6

О ВЫЧИСЛИТЕЛЬНОМ ТЕСТЕ ДЛЯ МОДЕЛИ АДИАБАТИЧЕСКОГО СЖАТИЯ ИДЕАЛЬНОГО БЕССТОЛКНОВИТЕЛЬНОГО ГАЗА

Д. А. Быковских1, В. А. Галкин2

1 Сургутский государственный университет, dmitriy.bykovskih@gmail. com 2Федеральный научный центр Научно-исследовательский институт системных исследований Российской академии наук, val-gal@yandex.ru

Работа посвящена исследованию тестовой задачи об адиабатическом сжатии идеального бесстолкновительного газа в одномерном пространстве. Аналитическое решение включает описание изменения во времени плотности распределения частиц, сгруппированных по скоростям. Представлено описание разработанного комплекса программ, позволяющих моделировать динамику идеального бесстолкновительного газа с подвижными границами в трехмерном пространстве методом прямого моделирования Монте-Карло. Проведены серии вычислительных экспериментов с различным числом частиц и скоростями движения границ, включая сравнение численных результатов с аналитическими решениями и оценку производительности разработанного комплекса программ.

Ключевые слова: идеальный бесстолкновительный газ, динамические системы, метод Монте-Карло.

ON COMPUTING TEST FOR ADIABATIC COMPRESSION MODEL OF IDEAL COLLISIONLESS GAS

D. A. Bykovskih1, V. A. Galkin2

1Surgut State University, dmitriy.bykovskih@gmail.com, 2System Research Institute, Russian Academy of Sciences, val-gal@yandex.ru

The article is devoted to computing test for adiabatic compression of an ideal collisionless gas in one-dimensional space. There is a description of particle density changes grouped by particle velocity over time for analytical solutions. The created software description is presented. The software allows simulating ideal collisionless gas dynamics with moving boundaries in the three-dimensional space by Monte Carlo method. Computing test series with different particle numbers and boundary velocities are run. The comparison between numerical results and analytical solutions is also made. The article contains software performance testing.

Keywords: ideal collisionless gas, dynamic system, Monte Carlo method.

Введение. Исследование динамических систем в изменяющейся во времени области (например, область, имеющая подвижные границы) является актуальной задачей в современных разделах естествознания. Решение такого сорта задач, как правило, связано с разработкой комплексов программ, позволяющих моделировать динамику поведения среды в изменяющейся во времени области. Промежуточным этапом может служить разработка тестовых задач, которые имеют аналитические решения и позволяют проверить корректность работы комплекса программ.

В работе рассматривается тестовая задача об адиабатическом сжатии идеального бесстолкновительного газа в одномерном пространстве. Такая задача относится к явлению параметрической неустойчивости, при котором нарастание энергии возмущения сопровождается непрерывным его сжатием с течением времени в пространстве. Это явление было подробно изучено и представлено А. И. Весницким в исследовании [1], связанной с колебанием струны и распространением волн в средах с подвижными границами, и А. Ф. Сидоровым в работе [2] о сжатии идеального газа поршнем.

1. Постановка тестовой задачи

Пусть в начальный момент времени ^ = 0 на отрезке [а, Ъ\ равномерно распределены частицы. Их плотность р0 = 1. Скорость каждой частицы |у0| = 1, а направление (т. е. знак

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

Пусть на концах отрезка [а, Ъ\ расположены границы. Одна из границ неподвижна и положение определяется уравнением Ъ = х. Другая граница движется со скоростью и е (0, у0 \, уменьшая первоначальную длину отрезка Ь = \Ъ — а|. Ее местоположение определяется уравнением а + и = х.

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

2. Аналитическое решение задачи

В начальный момент времени ^ = 0 существуют 2 группы частиц. Первая группа частиц имеет скорость V = V и функцию плотности распределения частиц /, а вторая - скорость у2 = —у0 и плотность /:

/ (х,0) = / (х,0) = ^ *(х, а, Ъ), (1)

, ч , ч , ч Г1, х е [а, Ъ\ где Их, а, Ъ) = в(х — а) в (Ъ — х) = \ г - ступенчатая функция.

[0, х £ [а, Ъ\

При взаимодействии частиц с границами изменяется их скорость, т. е. одновременно с исчезновением частицы из г -ой группы будет появляться частица из (г +1) -ой группы. Формула вычисления V скорости для г -ой группы частиц с учетом описанных ранее условий запишется в следующем виде:

=(" 'Г vo + 2u

i -1 2

Л

(2)

где г = 1, п.

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

i-i

L - u ^ t

t, =

,=0 (3)

+ (i - l)u

Траектория движения таких частиц с течением времени представлена на рис. 1. Функция (х, t) определяет принадлежность частиц i -ой группы в точке х в момент времени t:

¥г 0 = X(*, h, I-i)%(x,a + ut,b) + +X (t, t,.-i, tt) X (х, a + ut, a + uti-i + v, (t - t^)) (i mod 2) + (4)

+X (х, b + vi (t - tt-i), b) ((i + 1) mod 2)] +

+

X(л 1-1Д ) \Х (х, а + "I I + V, 0 - ^), 2>) (/' mod 2)

+

+% (х, а + ut, Ъ + V. (7 - ))((/ +1) mod 2)]. Плотность распределения для i -ой группы частиц с течением времени определяется

/; (х, 0=[ (t-t^zit, ti-i, U+fi-i

-(t-tjzit, tn I) ]((/-l)w + v0)/0(0^- (x, t) •

(5)

Рис. 1. Схема траекторий движения границ и первых и последних частиц в каждой группе с течением времени:

синие линии - траектории первых частиц, красные - траектории последних частиц

На рис. 2 представлен график изменения плотности распределения £(х, t) для частиц одной группы с течением времени на отрезке \а + п1, Ъ\. Изменение плотности распределения £(х, t) I -ой группы частиц с течением времени на отрезке \а + п1, Ъ\ представлено на рис. 2. На интервале [¿г_количество частиц увеличивается, а на [¿г, ^] их число уменьшается. На временном интервале , ^ их величина будет постоянной и равняться начальной плотности р0 /2.

Рис. 2. Схема изменения плотности распределения I -й группы частиц в зависимости от времени

С помощью плотности распределения частиц и их скоростей вычисляются остальные статистические оценки макроскопических параметров системы по следующим формулам [3, 4]:

n

p(x, t) = £ f (x, t),

(6)

i=1

V (х, г ) = (7)

р(х, г)

- V (х, г))2/(х, г)

е(х, г )=1 .м-„-, (8)

4 7 2 р(х, г)

2

Т (х, г ) = — е(х, г), (9)

р(х, г ) = - р(х, г) е (х, г). (10)

3. Описание комплекса программ

Траектория движения частицы в пространстве с учетом возможного взаимодействия с подвижными границами описывается следующими формулами:

к

х^ = х« +£Ч'+1)(г*+1 - гк), (11)

к=0

к = 0

"к = 1 «л Д^1)

vГ)ч , ;г+1) , , пр (12)

кОк-Л п(к), и(к)), к=1 к

где х - пространственные координаты частицы;

^ - скорость частицы после к -го взаимодействия с границей;

к к

2 г+1 - г* ) = 2 ^ = ^г - временной интервал;

"■к/ - 2'к

к=0 к=0

к - число столкновений частицы с границами; к - функция взаимодействия частицы с границей;

и (к) - скорость у -й границы в момент времени ^ на г -м шаге;

П(к) - нормаль к поверхности у -й границы в точке столкновения в момент времени гк на г -м шаге.

Пусть границы представляют собой плоскости или поверхности 1 -го порядка

(п, х - д(г)) = 0, (13)

)) = <

д(+1)= д()+ иАг, (14)

где п - нормаль к поверхности; д - координаты смещения границы.

Тогда момент времени, когда пересекается траектория движения частицы с движущейся плоскостью с учетом перехода в систему отсчета, где граница неподвижна, будет определяться как

г = (п х(г)- д)- шк-1) (15)

4 = (п, V')- и) . (15)

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

дится по другую сторону полупространства, ограниченного поверхностью, из-за вычислительной погрешности. Это приведет к повторному отражению и неверной траектории движения частицы (рис. 3). Решение заключается в добавлении поправочного коэффициента 5Мх, который будет сдвигать частицу по ее траектории назад на эту величину.

Рис. 3. Иллюстрация, поясняющая корректировку с помощью поправочного коэффициента 5 х :

жирная линия - граница; линия - траектория движения частицы с учетом взаимодействия с границей и поправочным коэффициентом 5Мх ; точечный пунктир - траектория движения частицы с учетом взаимодействия с границей без дополнительного поправочного коэффициента

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

,(<) -

(;+1)

П +1 =Р(П-Л П

](к ), и](к ).

) = ^)- 2п

I(к) V к

(г()- и

I (к), П](к)

(16)

На рис. 4 представлена схема расчета траектории движения частиц с учетом их возможного взаимодействия с подвижными границами. Такая схема позволяет ускорить расчеты за счет векторизации вычислений [7]. Хотя саму сортировку векторизовать невозможно, но за счет перестановки данных можно уменьшить число расчетов, переставляя в конец те частицы в группе, которые в дальнейшем на интервале Лt не взаимодействуют с границами.

Рис. 4. Схема расчета траектории движения частиц с учетом их возможного взаимодействия с подвижными границами

На рис. 5 представлена общая схема работы комплекса программ [8].

Рис. 5. Общая схема работы комплекса программ

Под установлением начального состояния модели подразумевается инициализация параметров частиц с помощью генератора псевдослучайных чисел, границ и элементарных областей. Каждая элементарная область накапливает суммы параметров частиц (общую массу, суммы компонентов скоростей и квадратов скоростей), принадлежащих определенной области в момент времени ^ . С помощью этих сумм вычисляются статистические оценки

макроскопических параметров по формулам, представленным в [9]. Сохранение промежуточных данных включает сохранение фазового пространства и сумм параметров частиц, которые используются на этапе обработки данных для визуализации процессов и вычисления статистических оценок макроскопических параметров соответственно. Также последний этап (обработка данных) может включать построение графиков изменения статистических оценок макроскопических параметров с течением времени и вычисление максимальных и относительных погрешностей, полученных в результате сравнения с точным решением.

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

4. Моделирование и анализ результатов

Были проведены серии вычислительных экспериментов при различном числе частиц и различных скоростях границ. На рис. 6 представлена визуализация процесса.

Рис. 6. Распределение частиц (N = 217) в различные моменты времени при и = 0,1:

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

На рис. 7 представлены графики изменений статистических оценок макроскопических величин при и = 0,1. Графики плотностей содержат плотности распределения г -й группы частиц в зависимости от времени. Из остальных графиков видно, что с увеличением числа частиц численное решение приближается к аналитическому. Следует отметить тот факт, что точное решение содержит осцилляции.

р Плотность Скорость вдоль оси X

р Адиабата т Температура

Рис. 7. Графики изменения статистических оценок макроскопических параметров в зависимости от времени при и = 0,1

В табл. 1 представлены результаты максимальных абсолютных и относительных погрешностей на временном интервале [0, 5], рассчитанные по формулам:

А¥ = тах| ¥ — ¥ I, (17)

| ап питр (17)

5F = х100%, (18)

| F | ( )

I an I

где Гап - аналитическое решение; ^пит - численное решение; ^ е [0; Т] - временной интервал.

Таблица 1

Максимальные абсолютные и относительные погрешности при u = 0,1 и различном числе частиц

N Av* AT Ap 5vx 5T 5p

102 4,23 х 10-1 1,68 х 10-2 2,13 х 10-1 5,10 х 105 1,85 х 101 1,85 х 101

103 1,18 х 10-1 3,96 х 10-3 4,77 х 10-2 2,44 х 104 5,18 5,18

104 3,16 х 10-2 7,01 х 10-4 9,06 х 10-3 3,29 х 104 7,82 х 10-1 7,82 х 10-1

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

Оценка производительности разработанного комплекса программ производилась на процессоре Intel Xeon E5-2690 V2 (240 GFlops DP для 10 ядер). Результаты показали 39 % от пиковой производительности на одно ядро (табл. 2) для основного расчетного блока (без вычисления статистических оценок макроскопических параметров) с учетом сохранения данных при следующих параметрах: шаг по времени At = 0.001, число итераций 5 000, число частиц N = 217. Ускорение оптимизированной версии по сравнению с неоптимизированной составило более чем 12 раз.

Таблица 2

Сравнение производительности оптимизированной и неоптимизированной

версии на Intel Xeon E5-2690 V2

t, с. Rtask, GFlops Rtask/Rpeak, %

Без оптимизаций 198,51 0,8 0,33

С оптимизациями 15,48 9,42 3,93

Заключение. Представлено описание тестовой задачи с подвижной границей для идеального бесстолкновительного газа. Был разработан комплекс программ, предназначенный для моделирования динамики идеального бесстолкновительного газа с подвижными границами. Проведены серии вычислительных экспериментов. Выполнено сравнение с точными решениями. На тестовых задачах об адиабатическом сжатии с применением декомпозиции и векторизации данных была достигнута производительность равная 39 % от пиковой для Intel Xeon E5-2690 V2.

Работа выполнена при поддержке гранта РФФИ №18-47-860004.

Литература

1. Весницкий А. И. Волны в системах с движущимися границами и нагрузками. М. : ФИЗМАТЛИТ, 2001. 320 с.

2. Сидоров А. Ф. Избранные труды: Математика. Механика. М. : ФИЗМАТЛИТ, 2001. 576 с.

3. Галкин В. А. Анализ математических моделей: системы законов сохранения, уравнения Больцмана и Смолуховского. М. : БИНОМ. Лаборатория знаний, 2011. 408 с.

4. Черчиньяни К. Математические методы в кинетической теории газов. М. : Мир, 1973. 246 с.

5. Цисарж В. В., Марусик Р. И. Математические методы компьютерной графики. Киев : ФАКТ, 2004. 464 с.

6. Гальперин Г. А., Чернов Н. И. Биллиарды и хаос. М. : Знание, 1991. 48 с.

7. Bik A. J. The Software Vectorization Handbook: Applying Intel Multimedia Extensions for Maximum Performance. Intel Press, 2004. 300 р.

8. Берд Г. Молекулярная газовая динамика / под ред. О. М. Белоцерковского, М. Н. Когана. М. : Мир, 1981. 319 с.

9. Быковских Д. А., Галкин В. А. О вычислительном тесте для одной модели бес-столкновительного идеального газа // Вестник кибернетики. 2017. № 3 (27). С. 119-127.