CC BY
28Научный вестник НГТУ. - 2008. - № 4(33)
УДК 519.683.6
Об алгоритмах сортировки в методе частиц в ячейках*
В.А. ВШИВКОВ, Т.В. МАРКЕЛОВА, В.И. ШЕЛЕХОВ
Рассмотрена проблема сортировки частиц в методе частиц в ячейках. В задачах моделирования методом частиц в ячейках на каждом временном шаге необходимо проводить сортировку частиц по ячейкам прямоугольной трехмерной сетки. Сортировка не является полной, поскольку внутри одной ячейки частицы могут быть неупорядоченными. Предложен новый алгоритм неполной сортировки, основанный на перестановке каждой частицы в нужную ячейку. Получен эффективный алгоритм, дано его строгое и полное описание на языке предикатного программирования.
Ключевые слова: алгоритмы сортировки, методы частиц в ячейке, язык предикатного программирования, неполная сортировка, трехмерные модели, процесс коагуляции
ВВЕДЕНИЕ
Метод частиц в ячейках широко используется для решения многомерных задач астрофизики, физики плазмы [1, 2]. Суть метода частиц в ячейках заключается в следующем. Моделируемое вещество (газ, плазма) представляется набором дискретных частиц, которые под действием внешних сил (гравитация, электрическое и магнитное поля) перемещаются в пространстве. Необходимые для движения модельных частиц силы определяются в узлах пространственной сетки путем решения системы алгебраических уравнений, получающихся после аппроксимации дифференциальных уравнений, описывающих эти поля. Плотности и потоки вещества, используемые на этапе вычисления сил, определяют в узлах сетки по координатам и скоростям модельных частиц. Точность результатов, полученных с помощью метода частиц в ячейках, зависит от размеров ячеек сетки и количества модельных частиц. Чем меньше пространственные шаги, чем больше модельных частиц в ячейках, тем выше точность расчета. Поэтому естественным представляется использовать многопроцессорные ЭВМ, что не только повышает точность расчета, но и существенно уменьшает время счета.
Наиболее остро стоит вопрос увеличения точности и уменьшения времени счета при решении трехмерных нестационарных задач. В этом случае трехмерная область решения делится на непересекающиеся подобласти со своей частью дискретной пространственной сетки и соответствующей частью модельных частиц, находящихся в этой подобласти. Каждая такая подобласть (или несколько подобластей) размещается на своем процессоре. Минимальный размер подобласти - ячейка пространственной
* Получена 27 октября 2008 г.
Работа выполнена при поддержке программы президиума РАН «Происхождение и эволюция биосферы» № 18-2, программы «Происхождение и эволюция звезд и галактик», программы Рособразования «Развитие научного потенциала ВШ» (проект РНП.2.1.1.1969), гранты РФФИ 08-01-00615, 08-01-00622
сетки. При моделировании физического процесса модельные частицы за шаг по времени могут переходить из одних подобластей в другие, в том числе и в находящиеся на других процессорах. В связи с этим возникает проблема неполной сортировки частиц, т. е. сортировки по подобластям. Расположение частиц относительно друг друга внутри подобласти роли не играет.
Аналогичная проблема возникает при моделировании процессов парных взаимодействий частиц, когда в качестве парной частицы необходимо рассматривать только частицы, находящиеся в той же ячейке, что и первая частица. В качестве таких процессов можно принять коагуляцию, парные упругие столкновения, гравитационное взаимодействие с ближними частицами и т. п. Здесь также приходим к задаче неполной сортировки частиц по ячейкам.
Особенность метода частиц в ячейках в том, что частица за шаг по времени может переходить только в соседнюю ячейку. Это свойство необходимо использовать при создании алгоритма сортировки для уменьшения времени сортировки. Если в трехмерной области решения введена прямоугольная расчетная сетка, то каждая ячейка граничит с 26 окружающими ее ячейками, в которую могут перейти частицы. Данный случай рассмотрен в настоящей статье. Создание эффективного алгоритма сортировки затрудняет то, что ячейки невозможно пронумеровать таким образом, чтобы частицы за шаг по времени переходили в ячейки с соседними номерами, поэтому нельзя эффективно использовать простые методы сортировки.
Впервые задача создания экономичного алгоритма неполной сортировки возникла перед авторами при прямом моделировании процесса коагуляции. Коагуляция — процесс слияния двух сталкивающихся частиц за счет действия внутренних сил Ван-дер-Ваальса. Отличительная особенность прямого моделирования заключается в том, что процесс описывается на уровне взаимодействия частиц системы. То есть частицам вещества противопоставляются модельные частицы, и на каждом временном шаге необходимо находить наиболее вероятные для слияния пары, частицы которых лежат в одной и той же ячейке. Так как все характеристики частицы записываются в лагран-жевых координатах, массив частиц упорядочен только в начальный момент времени. Отметим характерные особенности этой задачи. Во-первых, внутри ячейки сортировать частицы нет необходимости, так как в итоговом массиве частицы не располагаются строго упорядоченно по координатам. Во-вторых, за один шаг по времени каждая частица перемещается не далее чем в соседнюю ячейку. Соответственно при использовании хорошо известных алгоритмов полной сортировки, например сортировки Хоара (так называемая быстрая сортировка [3]) время упорядочивания массива оказывается больше, чем время для расчета передвижения частиц. Поэтому была поставлена задача разработки нового алгоритма сортировки, учитывающего описанные выше особенности.
Кроме того, такой алгоритм может быть использован для распределения частиц по процессорам в многопроцессорных версиях программ моделирования движения частиц, а также для определения соседей для метода 8РИ.
Предложенный алгоритм реализует неполную сортировку методом перестановки частиц с использованием лишь информации о местоположении частиц на предыдущем временном шаге, когда частицы были отсортированы. Эта информация хранится в массиве а заголовков ячеек. Заголовок есть индекс первой частицы ячейки в массиве частиц. На предыдущем шаге частицы 1-й ячейки имели номера в диапазоне [а(0, а( + 1) — 1]. Для каждой частицы из этого диапазона определяется номер ячейки к, в которой находится частица после перемещения на данном шаге. Далее частица переставляется таким образом, что попадает в диапазон номеров для к-й ячейки. Такая
перестановка сопровождается подвижкой массива частиц между ячейками с номерами i и k. Оказывается, что для каждой промежуточной ячейки от i и k достаточно передвинуть только одну частицу и скорректировать заголовок ячейки.
Допустим, имеется трехмерная сетка размером ПхПхП. Пусть на очередном шаге частица перелетела в соседнюю ячейку по первому измерению. Тогда в процессе перестановки такой частицы на нужное место потребуется переставить ПхП частиц. При больших n алгоритм будет работать медленней, однако на порядок быстрее алгоритма Хоара.
В данной работе дается строгое и полное описание предложенного алгоритма сортировки. Описание алгоритма представлено в виде набора определений предикатов на языке P [4]. Используемый фрагмент языка P определен в разд. 1. Описание алгоритма на языке P проще и понятнее, чем на целевом языке программирования ФОРТРАН. Программа на языке ФОРТРАН, получаемая как результат оптимизирующей трансформации описания алгоритма на языке P, по эффективности лучше программы, запрограммированной вручную. Язык P использовался для детального описания ряда алгоритмов [5], среди них — нетривиальный алгоритм Маккрейта построения дерева суффиксов [6].
Формальная постановка задачи неполной сортировки частиц по ячейкам трехмерной прямоугольной сетки дается в разд. 2. Алгоритм сортировки методом перестановки частиц описывается в разд. 3. Эффективные программы на языке ФОРТРАН для предикатных программ алгоритма получены применением набора трансформаций и описаны в конце разд. 3. В разд. 4 дано сравнение алгоритмов сортировки при варьировании следующих параметров: времени очередного шага, размера сетки и числа частиц.
1. ЯЗЫК ОПИСАНИЯ АЛГОРИТМОВ СОРТИРОВКИ
В качестве языка публикаций алгоритмов используется фрагмент языка предикатного программирования P [4]. Предикатная программа есть замкнутый набор рекурсивных определений предикатов. Определение предиката имеет следующую структуру:
<имя предиката>(<описания аргументов>: <описания результатов>) = pre <предусловие> { <оператор> } post <постусловие>;
Параметрами предиката являются переменные — аргументы и результаты. Описание исходного или результирующего параметра предиката имеет следующий вид:
<тип><пробел><имя параметра>
Определение предиката описывает функцию, отображающую значения аргументов в значения результатов. Предусловие ограничивает область определения этой функции, а постусловие связывает значения аргументов и результатов. Предусловие и постусловие являются формулами исчисления предикатов высших порядков.
Операторами в языке P являются: оператор присваивания, вызов предиката, блок, параллельный оператор, условный оператор и др. Оператором присваивания является конструкция: <переменная> = <выражение>. Блок {<опера-тор1>; <Оператор2>; ...} определяет последовательное исполнение операторов. В параллельном операторе <Оператор1> || <Оператор2> результаты первого оператора не используются во втором, а второго - в первом; исполнение операторов может
проводиться параллельно. Отметим, что операторы могут быть легко преобразованы в эквивалентные логические формулы.
Среди операторов блока могут находиться описания локальных переменных: <тип> <пробел> <список имен переменных>. Возможно также описание переменной с инициализацией: <тип> <пробел> <имя переменной> = <выражение>.
Оператор может отсутствовать в определении предиката. В этом случае определение предиката играет роль спецификации предиката.
Ниже даны примеры описаний переменных. Текст текущей строки после пары символов "//" является комментарием.
nat i = 0; // описание натуральной переменной с инициализацией
seq int p; // p — последовательность целых чисел
Для обозначения типа, представленного изображением типа, некоторым именем используется описание типа. Например,
type Diap = 1..10; // тип диапазона целых чисел от 1 до 10
type SeqInt = seq int; // тип последовательности целых чисел
type Complex = struct real re, im end;// Complex—тип структуры (кортежа) с
полями re и im
type Ar= array 1..10, 1..5 of real; // двумерный массив вещественных с индексами от 1 до 10 и 1 до 5.
Для непустой последовательности s определены операции head(s) и tail(s), вычисляющие соответственно начальный элемент последовательности s и оставшуюся часть s. Пустая последовательность обозначается константой nil. Операция length(s) вычисляет число элементов последовательности s. В операции конкатенации a + b, где любой из аргументов a и b может быть последовательностью или элементом.
2. ЗАДАЧА СОРТИРОВКИ
Имеется конечное множество частиц внутри параллелепипеда в трехмерном пространстве. Каждая частица идентифицируется своими координатами x, y и z.
Описание трехмерной сетки. Параллелепипед покрыт трехмерной сеткой. Ячейка сетки имеет размеры Ax, Ay, Az в системе координат x, y и z. Левый ближний нижний угол параллелепипеда определяется координатами X0, Y0 и Z0:
real Ax, Ay, Az, X0, Y0, Z0;
Число ячеек, размещающихся в параллелепипеде по трем измерениям, nx, ny и nz соответственно:
nat nx, ny, nz;
Каждая ячейка сетки характеризуется некоторым индексом (i, j, k), где i = 1..nx, j = 1..ny, k = 1..nz. Левый ближний нижний угол ячейки имеет координаты xi, yj и Zk, где
xi = X0 + (i - 1) * Ax, yj = Y0 + (j - 1) * Ay, zk = Z0 + (k - 1) * Az.
Таким образом, мы определяем (/, j, k)-ячейку как множество
[xi .. xi + Ax) x [yj .. yj + Ay) x [zk .. zk + Az).
Упорядоченность ячеек естественным образом определяется лексикографической упорядоченностью их индексов. Эта упорядоченность представлена следующим отношением "<" для индексов:
(il, jl, kl) < (¡2, j2, k2) s il < ¡2 V il = ¡2 & (jl < j2 v ji = j2 & kl < k2).
Вычисление индекса произвольной точки параллелепипеда (x, y, z) реализуется функцией:
index(real x, y, z: nat i, j, k) = (1)
pre particleInside(x, y, z)
post i = entier((X0 - x)/Ax) + 1 & j = entier((Y0 - y)/Ay) + 1 & k = en-tier((Z0 - z)/Az) + 1;
Стандартная функция entier вычисляет целую часть вещественного аргумента. Предикат particleInside имеет следующее определение:
particleInside(x, y, z) s X0 < x < X0 + nx x Ax & Y0 < y < Y0 + ny x Ay &
Z0 < z < Z0 + nz x Az (2)
Вычисление по индексу (i, j, k) следующего индекса в соответствии с лексикографической упорядоченностью реализуется функцией
next(nat i, j, k: nat a, b, c) = pre inRange(i, j, k)
post if k = nz then c = l & if j = ny then a = i + l & b = l
else a = i & b = j + l end else a = i & b = j & c = k + l end;
inRange(i, j, k) s l < i < nx & l < j < ny & l < k < nz (3)
Отметим, что next(nx, ny, nz) = (nx + l, l, l). Ячейка с таким индексом находится за пределами параллелепипеда. Тем не менее данная особенность далее будет использована.
Для произвольных i и j, i = l..nx, j = l..ny определим секцию ячеек с индексом (i, j), называемую также (i,])-секцией. Это набор (i, j, ^-ячеек для всех k = l..nz. Секции с индексом (i, j) соответствует множество:
[xi .. xi + Ax) x [yj .. yj + Ay) x [Z0 .. Z0 + nz x Az). Введем понятие слоя с индексом i или -слоя. Это набор секций с индексами (i, j) для всех j = l.. ny. Слою с индексом i соответствует множество
[xi .. xi + Ax) x [Y0 .. Y0 + ny x Ay) x [Z0 .. Z0 + nz x Az). Последовательность частиц. Каждая частица представлена объектом типа PARTICLE с полями x, y и z, определяющими координаты частицы в пространстве. Тип PARTICLE имеет другие поля, которые не используются в алгоритме: type PARTICLE = struct real x, y, z, ... end;
Множество частиц внутри параллелепипеда представляется в виде последовательности следующего типа:
type SPART = seq PARTICLE;
Исходная последовательность частиц P, которую требуется отсортировать по ячейкам, определяется следующим описанием:
SPART P;
Общее число частиц число частиц N является исходным параметром задачи. ^ N
Частицы последовательности Р пронумерованы от 1 до N. Для частицы используется обозначение РЦ]; а j является номером частицы РЦ] в последовательности Р.
Частицы последовательности Р должны находиться внутри параллелепипеда, т. е. является истинным предикат:
т81с1е(Р) = vj=1..N. particleInside(P[j].x, РШ.у, Р[Л^) (4)
Сортированная последовательность частиц. Последовательность частиц Р является сортированной, если индексы частиц следуют в неубывающем порядке. Это определяется истинностью следующего предиката:
sorted(P) = vi, j = 1..№ (i < j ^ тСех(Рр].х, P[i].y, P[i].z) < index(P[j].x,
РШ.У, РШ^)). (4')
Рассмотрим свойства приведенного определения сортировки по ячейкам. Обозначим через ГО, j, к, Р) множество номеров частиц, находящихся в ячейке с индексом 0, j, к). Нетрудно доказать, что Г^, j, к, Р) определяет сплошной диапазон номеров: существуют такие И и ^ что Г^, j, к, Р) = [И..Ц. Для обозначения диапазонов будем использовать массивы Ир, j, к] и ф, j, к], где i = 1..ПХ, j = 1..пу, к = 1..nz. Таким образом, ГО, j, к, Р) = [Ир, j, k]..t[i, j, к]].
Если (и, jl, к!) < (i2, j2, к2) и ячейки с данными индексами не пусты, то Г(к, j!, к1, Р) < Г(i2, j2, к2, Р)*. Кроме того, если индекс (i2, j2, к2) непосредственно следует (в смысле лексикографической упорядоченности) за индексом j!, к1), т. е. О2, j2, к2) = next(i!, j!, к1), то
И[12, к2] = t[i!, j!, к1] + 1. (5)
Условимся, что соотношение (5) имеет место и тогда, когда 02, j2, к2)-ячейка окажется пустой; при этом
t[i2, j2, к2] = h[i2, j2, к2] - 1.
Кроме того, И[1, 1, 1] = 1.
Таким образом, сортированная последовательность Р разбита на диапазоны. Каждый диапазон определяет номера частиц в некоторой ячейке. Последовательность диапазонов соответствует лексикографической упорядоченности ячеек. С учетом соглашения (5) для пустых ячеек обеспечивается единственность определения массивов Ир, j, к] и t[i, j, к] для любой сортированной последовательности Р. В соответствии с соотношением (5) Ц^, j!, к1] может быть представлено через Ир2, j2, к2] следующего диапазона. Поэтому для сортированной последовательности достаточно рассматривать только массив Ир, j, к]. Значение h[i, j, к] будем называть заголовком 0, j, к)-ячейки.
Рассмотрим элемент РИ, 1 < S < N в отсортированной последовательности Р. Пусть (а, Ь, с) — индекс ячейки, в которой находится частица РИ. С учетом того что next(nx, ПУ, nz) = ^ + 1, 1, 1), определим + 1, 1, 1] = N + 1. Тогда истинно следующее соотношение:
И[а, Ь, с] < s < h[next(a, Ь, с)] (6)
* Отношение "<" определяется на диапазонах очевидным образом.
В итоге, связь отсортированной последовательности P с массивом заголовков ячеек h[i, j, k] определяется предикатом:
heads(h, P) = h[1, 1, 1] = 1 &
Vi=1..nx, j=1..ny, k=1..nz. if r(i, j, k, P) = 0 then h[next(i, j, k)] = h[i, j, k]
else h[next(i, j, k)] = max r(i, j, k, P) + 1 end;
Определим тип HEADCELLS для массива заголовков h[i, j, k].
type HEADCELLS = array 1..nx+1, 1..ny, 1..nz of nat;
Отметим еще одно свойство. Для частиц, находящихся в одной ячейке, порядок их следования в сортированной последовательности P может быть любым. Иначе говоря, внутри ячейки сортировка частиц не проводится.
Множество номеров частиц, находящихся в (i, j)-секции, представимо в виде некоторого диапазона, объединяющего диапазоны для (i, j, ^-ячеек. Заголовком секции назовем начало этого диапазона.
Свойства близости частиц. Задача сортировки частиц является частью задачи моделирования движения частиц. Сортировка реализуется после каждого акта перемещения последовательности частиц. Пусть Q — отсортированная последовательность частиц в начале акта перемещения, а P — в конце акта. Нашей задачей является сортировка последовательности P.
Считается, что в каждом акте перемещения частицы не могут разлететься слишком далеко: каждая частица может переместиться не далее соседней (смежной) ячейки*. Это условие формализуется как истинность следующего предиката:
near(P, Q) = Vl=1..N. nearIndex(index(P[l].x, P[l].y, P[l].z), index(Q[l].x, Q[l].y,
Q[l].z))
Близость ячеек с индексами (i, j, k) и (a, b, c) определяется предикатом:
nearIndex(i, j, k, a, b, c) = |a - i| < 1 & |b - j| < 1 &|c - k| < 1 (7)
В последовательности P рассмотрим произвольную частицу с номером s. Пусть (a, b, c) — индекс ячейки, в которой находится частица Q[s]. Пусть h[i, j, k] — массив заголовков для последовательности Q. Поскольку Q — отсортированная последовательность, имеет место соотношение (6) для Q. Нетрудно доказать, что из истинности near(P, Q) следует истинность следующего утверждения:
3a,b,c. nearIndex(index(P[s].x, P[s].y, P[s].z), a, b, c) & h[a, b, c] < s <
< h[next(a, b, c)] (8)
Будем говорить, что частица P[s] прописана в (a, b, ^-ячейке, если частица Q[s] находится в (a, b, c)-ячейке. Таким образом, в соответствии с утверждением (7) частица P[s] прописана в одной из смежных ячеек. Верно обратное утверждение. Частица P[s], прописанная в некоторой ячейке, реально находится в одной из смежных ячеек:
h[a, b, c] < s < h[next(a, b, c)] ^ nearIndex(index(P[s].x, P[s].y, P[s].z), a, b, c)
(9)
Утверждения (8) и (9) используются далее для построения эффективного алгоритма сортировки.
* Число смежных ячеек не более 26.
Внешняя спецификация программы сортировки. Программа реализует сортировку последовательности частиц P, являющейся результатом акта перемещения некоторой сортированной последовательности. Для нее имеется массив заголовков h. Результатами программы являются отсортированная последовательность P' и соответствующий ей массив заголовков h'. Обозначение h' определяет, что переменные h и h' должны быть склеены в реализации, т. е. в реализации h' заменяется на h. Программа представлена предикатом sort:
sort(HEADCELLS h*, SPART P*: ) = pre nearSorted(h, P)
post sorted(P') & perm(P, P') & heads(h', PO;
Мы используем запись h* и P* для обозначения того, что вместе с исходными параметрами предиката h и P имеются результирующие параметры h' и P', склеиваемые в реализации, соответственно с h и P; при этом из соображений компактности мы не указываем переменные h' и P' среди результатов предиката.
Предусловие nearSorted(h, P) определяет связь последовательности P с некоторой сортированной последовательностью до акта ее перемещения.
nearSorted(h, P) =
= i inside(P) & 3Q (inside(Q) & near(P, Q) & sorted(Q) & heads(h, Q))
Отметим, что предикат inside(Q) ограничивает число смежных ячеек, определяемых предикатом near(P, Q), для ячеек на границе параллелепипеда.
Последовательность P' должна быть перестановкой последовательности P, т. е. истинен предикат
perm(P, P') = 3 f: [1..N] ^ [1..N]. f — взаимно-однозначная & Vj = 1..N. P'[j] =
P[f(j)]
3. СОРТИРОВКА МЕТОДОМ ПЕРЕСТАНОВКИ ЧАСТИЦ
Задача сортировки частиц в ячейках может быть решена любым известным алгоритмом обычной сортировки. Для реализации поячеечной сортировки можно также адаптировать алгоритм обычной сортировки. Рассмотрим адаптацию алгоритма сортировки простыми вставками [3]. В этом алгоритме очередной элемент вставляется в нужное место отсортированной части последовательности. Аналогично, в предлагаемом нами алгоритме очередная частица вставляется в нужную ячейку с использованием старого массива заголовков. В алгоритме не используются отмеченные свойства близости (8) и (9).
Вместо трехмерного индекса удобнее использовать линеаризованный индекс, определяемый следующим образом:
linIndex(i, j, k: n) = pre inRange(i, j, k)
post n = (i - 1) * ny * nz + (j - 1) * nz + k;
Для линеаризованного индекса однозначно определяется соответствующий ему трехмерный индекс, т. е.
(vn=1...N+1) (3! i, j, k) inRange(i, j, k) & n = linIndex(i, j, k)
Таким образом, все условия, сформулированные для трехмерного индекса, могут быть использованы для линеаризованного индекса. Выполняется следующее соотношение:
linIndex(i, j, k) + 1 = linIndex(next(i, j, k)).
Введем тип LINHEADS, определяющий массив заголовков, индексируемый по линеаризованному индексу:
nat nn = nx * ny * nz;
type LINHEADS = array 1.. nn + 1 of nat;
Массив заголовков hl для линеаризованного индекса определим по массиву заголовков h для трехмерного индекса следующим образом:
hl[linIndex(i, j, k)] = h[i, j, k]
Преобразование заголовка от линеаризованного индекса к трехмерному определим предикатом
Lh3(hl: h) = h[i, j, k] = hl[linIndex(i, j, k)]
Программу сортировки методом перестановки частиц представим следующим определением:
sortL(LINHEADS* hl, SPART P*:) = pre inside(P) & 3Q. headsL(hl, Q)
{ sortD(hl, P, 1: hl', P') } post sorted(P') & perm(P, P') & headsL(hl', P');
headsL(hl', P') = heads(Lh3(hl'), P')
Предикат sortD представляет более общую задачу сортировки последовательности P начиная с ячейки, имеющей линеаризованный индекс n, в предположении, что все ячейки с меньшими номерами уже отсортированы и массив заголовков h соответствует отсортированному состоянию. Данное предположение фиксируется в предусловии предикатом sortedN(h, P, n):
sortD(LINHEADS h*, SPART P*, nat n: ) = pre inside(P) & 3Q. headsL(h, Q) & 1 < n < nn & sortedN(h, P, n) { sortSell(h, P, n: LINHEADS h1, SPART P1);
if n = nn then h' = h1 || P' = P1 else sortD(h1, P1, n + 1: h', P') end
}
post sorted(P') & perm(P, P') & headsL(h', P');
Предикат sortSell реализует сортировку частиц, прописанных в ячейке с линеаризованным индексом n, в предположении, что часть последовательности P (в ячейках с номерами от 1 до n-1) уже отсортирована:
sortSell(LINHEADS h*, SPART P*, nat n: ) = pre inside(P) & 3Q. headsL(h, q) & 1 < n < nn & sortedN(h, P, n) { sortDiap(h, P, n, h[n], h[n + 1]: h', P') }
post perm(P, P') & sortedN(h', P', n + 1);
Предикат sortDiap сортирует отрезок [u..s-1] последовательности P в ячейке с линеаризованным индексом n. Предикат sortedD(h, P, u) в предусловии определяет, что часть последовательности P для частиц с номерами от 1 до u — 1 уже отсортирована, причем массив заголовков h соответствует отсортированному состоянию:
sortDiap(LINHEADS h*, SPART P*, nat n, u, s: ) =
pre inside(P) & 3Q. headsL(h, Q) & 1 < n < nn & sortedD(h, P, u) &
h[n] < u < s = h[n+1] { if u = s then h' = h || P' = P
else sortParticle(h, P, u, s, n : hl, P1, nat u1, s1); sortDiap(h1, P1, n, u1, s1: h', P')
end
}
post perm(P, P') & sortedN(h', P', n + 1);
Предикат sortParticle сортирует частицу с номером u, прописанную в ячейке с линеаризованным индексом n. Частица перемещается в ячейку, индекс которой соответствует координатам частицы. Это реализуется цепочкой перемещений между двумя ячейками с корректировкой заголовков. Такие действия сохраняют сортированность начальной части (от 1 до u' — 1) последовательности P':
sortParticle(LINHEADS h*, SPART P*, nat u*,s*, n: ) =
pre inside(P) & 3Q. headsL(h, Q) & 1 < n < nn & sortedD(h, P, u) &
h[n] < u < s = h[n+1] { PARTICLE p = P[u];
nat m = linIndex(index(p.x, p.y, p.z));
if m = n then h' = h || P' = P || u' = u + 1 || s' = s
elsif m < n then moveBack(h, P, u, n, m, p: h', P') || u' = u + 1 || s' = s
else moveForward(h, P, u, n, m, p : h', P') || u' = u || s' = s - 1
end
}
post perm(P, P') & sortedD(h', P', u') & h'[n] < u'' < s' = h'[n+1] & (u' = u + 1 xor s' = s - 1);
Предикаты moveBack и moveForward перемещают частицу p в ячейку с номером m, где она и должна находиться в соответствии со своими координатами.
Предикат moveBack реализует перемещение частицы p с номером u, находящейся в ячейке с линеаризованным индексом n, в ячейку с номером m. Перемещение частицы p реализуется цепочкой элементарных перемещений. Каждое из них перемещает частицу p в предыдущую ячейку (рис.1): сначала частица p обменивается с начальной частицей ячейки q = P[h[n]], затем заголовок h[n] сдвигается на одну позицию вправо (h[n] := h[n]+1), в результате чего частица p оказывается в предыдущей ячейке. Если m = n — 1, то результат достигнут. Иначе реализуется элементарное перемещение частицы p из предыдущей ячейки и т. д., пока не будет достигнута ячейка с номером m.
h[n]_u_h[n+1]
1 q 1 1 Р 1
u h[n] ^ Г1 ^ h[n+1]
1 Р q I
Рис. 1. Перемещение частицы р в предыдущую ячейку; ячейка с номером П выделена серым цветом
В описанном алгоритме частица p последовательно обменивается с начальными частицами промежуточных ячеек. Очевидно, что можно избежать перемещения частицы p в промежуточные ячейки. Достаточно перемещать лишь начальные частицы во всех промежуточных ячейках, пока не будет достигнута ячейка с номером m. С учетом этого вместо P мы должны рассматривать сортировку последовательности:
P~ = (P[1]...P[u-1], p, P[u+1]...P[N]). (10)
Предикат sortedF(h, P, u, u0) в предусловии определяет, что часть последовательности P для частиц с номерами 1, ..., u — 1, u + 1, ..., u0 уже отсортирована, а массив заголовков h соответствует отсортированному состоянию. Значение u0 соответствует значению переменной u в предикате sortParticle: moveBack(LINHEADS h*, SPART P*, nat u, n, m, PARTICLE p: ) = pre inside(P) & 3Q. headsL(h, Q) & 1 < m < n < nn & sortedF(h, P, u, u0) & h[n] <
u < h[n+1] & m = linIndex(index(p.x, p.y, p.z)) { nat v = h[n];
P1[u] = P[v]; h1[n] = v + 1;
if m = n - 1 then P1[v] = p; { h' = h1 || P' = P1} else moveBack(h1, P1, v, n - 1, m, p : h', PO end
}
post perm(P~, P') & sortedD(h', P', u + 1) & P'[h'[m + 1] - 1] = p;
Значение последовательности P~ в постусловии определено формулой (10).
Предикат moveForward реализует перемещение частицы p с номером u, находящейся в ячейке с линеаризованным индексом n, в ячейку с номером m. Перемещение реализуется в виде цепочки элементарных перемещений. В каждом из них сначала частица p переносится в конец текущей ячейки, затем заголовок h[n + 1] смещается на позицию влево и частица p оказывается в следующей ячейке. В итоге частица p переносится в ячейку с номером m с сохранением содержимого проходимых ячеек: moveForward(LINHEADS h*, SPART P*, nat u, n, m, PARTICLE p: ) = pre inside(P) & 3Q. headsL(h, Q) & 1 < n < m < nn & h[n] < u < h[n+1] &
m = linIndex(index(p.x, p.y, p.z)) { nat v = h[n + 1] - 1;
P1[u] = P[v]; h1[n + 1] = v;
if m = n + 1 then P1[v] = p; { h' = h1 || P' = P1} else moveBack(h1, P1, v, n + 1, m, p : h', P') end
}
post perm(P~, P') & sortedD(h', P', u) & P'[h'[m]] = p;
Значение последовательности P~ в постусловии определено формулой (10).
Приведенные определения составляют полную предикатную программу сортировки методом перестановки частиц. Для получения эффективной программы к предикатной программе применен следующий набор трансформаций: 1) склеивание переменных, реализующее замену нескольких переменных одной; 2) замена хвостовой рекурсии циклом; 3) подстановка определения предиката на место его вызова; 4) кодирование последовательности частиц массивом.
4. РЕЗУЛЬТАТЫ СРАВНЕНИЯ АЛГОРИТМОВ
Описанные выше алгоритм был реализован на языке ФОРТРАН и сравнивался по времени счета с встроенной в пакет DDFPORT ФОРТРАНа процедурой быстрой сортировки qsort.
Анализ результатов запуска программ в разных режимах показал, что для частиц численностью более 10 000 qsort имеет гораздо меньшую скорость, чем разработанный алгоритм. В данных тестах замерялось время для выполнения процедуры сортировки на каждом шаге и вычислялось среднее время одной сортировки массива для нескольких десятков шагов. Соответственно путем варьирования шага по времени tau можно регулировать число частиц, вылетающих из своих ячеек. Чем меньше шаг, тем короче расстояние, пролетаемое частицей, а значит, и меньше число частиц, вылетающих из своих ячеек.
Первая серия расчетов проводилась на одной и той же сетке и с постоянным количеством частиц, менялось только tau (рис. 2).
Время счета, с
tau
сортировка методом перестановки частиц —■—быстрая сортировка
Рис. 2. Время счета на сетке 10x10x10 для 1 000 000 частиц
Для данного случая программа новой сортировки работает гораздо быстрее qsort.
Во второй серии экспериментов tau постоянное, число частиц постоянное, менялось количество ячеек, а соответственно и среднее число частиц в ячейке (рис. 3). Быстрая сортировка для этого случая, как и в предыдущий раз, работала медленнее.
В последней серии экспериментов варьировалось общее число частиц на двух видах сеток, для 1000 и 10 000 частиц обе сортировки показали сходное время, но с увеличением числа частиц время счета метода быстрой сортировки увеличивалось значительнее, чем у нового алгоритма (рис. 4 и 5).
Сравнение двух алгоритмов показало, что их использование оправдано для расчетов на подробных сетках с количеством частиц порядка одного миллиона и дает ускорение в три и более раза.
Время счета, с
1 10 100 10ОО
Число частиц в ячейке
сортировка методом перестановки частиц —■—быстрая сортировка
Рис. 3. Время счета для tau = 0.001 и 1 000 000 частиц
Время счета, с
Число частиц
сортировка методом перестановки частиц —■— быстрая сортировка
Рис. 4. Время счета для tau = 0.008 на сетке 10x10x10
Время счета, с
Число частиц
сортировка методом перестановки частиц —■—быстрая сортировка
Рис. 5. Время счета для tau = 0.008 на сетке 50x50x50
СПИСОК ЛИТЕРАТУРЫ
[1] Хокни Р., Иствуд Дж. Численное моделирование методом частиц. - М.: Мир, 1987.
[2] Березин Ю.А., Вшивков В.А. Метод частиц в динамике разреженной плазмы. - Новосибирск: Наука, 1980.
[3] Кнут Д. Искусство программирования для ЭВМ / пер. с англ. - М.: Мир, 1978. - Т.З.Сортировка и
поиск.
[4] Шелехов В.И. Язык предикатного программирования P. - Новосибирск, 2002. - (Препринт / ИСИ СО РАН; N 101).
[5] Методы предикатного программирования: сб. работ под ред. В.И. Шелехова. - Новосибирск: ИСИ СО РАН, 2004.
[6] Шелехов В.И. Разработка программы построения дерева суффиксов в технологии предикатного программирования. — Новосибирск, 2004. - (Препр. / ИСИ СО РАН; N 115).
[7] Hoare C. A. R. An axiomatic basis for computer programming // Communic. of the ACM. - 1969. -12(10). - P. 576-585.
[8] Kernighan B., Ritchie D. The C Programming Language. - N.Y.: Prentice Hall, 1988.
[9] Йенсен К., Вирт Н. Паскаль. Руководство для пользователя и описание языка. - М.: Финансы и статистика, 1982.
[10] Martinez C. Partial quicksort // Proc. of the 6th ACM-SIAM Work. Alg. Engng and Exp. (ALENEX) and the 1st ACM-SIAM Work. on Anal. Algorithm. and Combinat. (ANALCO). - P. 224-28.
Вшивков Виталий Андреевич, доктор физико-математических наук, профессор кафедры параллельных вычислительных технологий факультета прикладной математики и информатики Новосибирского государственного технического университета, заведующий лабораторией Института вычислительной математики и математической геофизики СО РАН. Основное направление научных исследований - математическое моделирование. Имеет более 200 публикаций, в том числе 8 учебных пособий и монографий.
Маркелова Тамара Валерьевна, аспирант кафедры математического моделирования механико-математического факультета Новосибирского государственного университета, младший научный сотрудник Института катализа СО РАН. Основное направление научных исследований - математическое моделирование.
Шелехов Владимир Иванович, кандидат технических наук, заведующий лабораторией Института систем информатики СО РАН. Основное направление научных исследований - методы трансляции и оптимизации программ, статический анализ программ, языки функционального программирования. Имеет более 60 публикаций.
V.A. Vshivkov, T.V. Markelova, V.I. Shelekhov
About Sorting Algorithms in Particle-in-Cell Method
The problem of particle sorting in solving 3D problems with PIC method is considered. Simulation with PIC method requires ordering of particle assignment to cells. It is important in parallelization of the method when each processor contains parts of the grid domain and corresponding particles. Particle order is required to search the closest particles. During simulation of a process particles can fly to neighbour cells. Thus it is necessary to perform particle sorting depending on the grid cells. The sorting algorithm must be fast enough. Sorting should not be complete because inside a cell particles may not be sorted. A new sorting algorithm is proposed. The algorithm is based on the procedure moving each particle to the necessary cell. The peculiarities of the PIC method enabled us to gain an efficient algorithm. The algorithm is thoroughly and strictly described in a predicate programming language.
Key words: sorting algorithms, .particle-in-cell methods, predicate programming, not full sorting, three-dimensional models, coagulation process
