Размещено на http://www.allbest.ru/

ФЕДЕРАЛЬНОЕ ГОСУДАРСТВЕННОЕ АВТОНОМНОЕ

ОБРАЗОВАТЕЛЬНОЕ УЧРЕЖДЕНИЕ

ВЫСШЕГО ПРОФЕССИОНАЛЬНОГО ОБРАЗОВАНИЯ

«НАЦИОНАЛЬНЫЙ ИССЛЕДОВАТЕЛЬСКИЙ УНИВЕРСИТЕТ

«ВЫСШАЯ ШКОЛА ЭКОНОМИКИ»

Московский институт электроники и математики им. А.Н. Тихонова

Выпускная квалификационная работа

по направлению 01.03.04 Прикладная математика

студента образовательной программы бакалавриата

«Прикладная математика»

Анализ перколяционных образцов случайной последовательной адсорбции k-меров на квадратной решетке

Слуцкий Михаил Георгиевич

Аннотация

Это исследование нацелено на изучение поведения плотностей перколяции и джамминга для образцов изотропной случайно последовательной адсорбции при помощи моделирования. Рассматриваются большие линейные -меры на квадратной решетке с периодическими граничными условиями. Для симуляции процесса разработан алгоритм, эффективный по использованию памяти и времени. Целью работы является подтверждение предыдущих результатов [1] и получение отношения плотностей перколяции и джамминга для -меров c длиной до

Abstract

This research is aimed at studying the behavior of percolation and jamming thresholds for isotropic random sequential adsorption samples by means of numerical simulations. Mostly large linear k-mers on a square lattice with periodic boundary conditions are considered. For simulating the process, an efficient in terms of memory and speed algorithm is presented. The goal is to verify previous results [1] and to obtain the ratio of percolation and jamming concentrations for lengths of k-mer up to .

Введение

перколяция джамминг адсорбция изотропный

В статистической физике теория перколяции описывает поведение кластеров частиц в какой-либо среде [7,8]. С помощью нее можно описать такие явления как протекание воды сквозь почву, протекание тока, распространение эпидемий и многие другие. В то время как перколяция на сравнительно небольших решетках довольно подробно изучена, с системами большого размера дело предстоит иначе. В силу несовершенства вычислительной техники и алгоритмов представлялось невозможным смоделировать систему для ее дальнейшего изучения. В данной работе будет описан алгоритм позволяющий моделировать системы с иглами

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

1. Основные понятия

1.1 Описание модели

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

Процесс будет происходить на дискретной квадратной решетке размера , где - длина стороны квадрата. Более того, рассматриваются периодические граничные условия, т.е. противоположные стороны квадраты “сшиты” друг с другом, превращая таким образом квадрат в тор. Каждая клетка решетки может быть либо занята, либо свободна. На торе будут случайным образом располагаться одинаково длинные иглы или -меры, которые на решетке будут являться прямоугольниками размера . Соответственно -меры на решетке могут иметь одну из 2 ориентаций: вертикальную или горизонтальную. Также в расчёт берется тот факт, что длина иглы не должна быть соизмерима с длиной стороны решетки, поэтому . Два -мера принадлежат одному кластеру, если существует такое множество граничащих друг с другом клеток, что в него входят клетки обоих-мера.

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

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

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

1.2 Обзор литературы

В недавних статьях [1,2], в которых рассматривалась аналогичная модель, были посчитаны плотности перколяции и джамминга для игл длиной , а также их отношение (красный график).

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

Помимо этого, в 2017 году вышла статья [4], где рассматривалась похожая система (отличие было в том, что рассматривались жесткие граничные условия и соответственно “spanning” перколяция), в которой было математически доказано, что для таких систем (вне зависимости от размеров иглы и решетки) на момент джамминга всегда существует перколяционный кластер.

Для теории перколяции представляет интерес проверить предположение Тарасевича и посмотреть, как ведет себя график при .

2. Описание алгоритма

2.1 Заполнение решетки

У каждой иглы есть самая верхняя левая клетка, и ее мы будем называть начальной клеткой -мера. Чтобы случайным образом расположить -мер на решетке, необходимо сгенерировать обе координаты начальной клетки, а также выбрать ориентацию -мера. После этого программа должна проверить, имеется ли возможность разместить k-мер, чтобы избежать наложений. если такая возможность есть, то -мер располагается на решетке, в противном случае выбираются новые координаты и ориентация. Но у этого метода есть существенные недостатки, которые приводят к необходимости его усовершенствования для больших решеток. Если решетка почти полностью заполнена, то программа будет совершать большое кол-во “холостых” итераций (т.е. неуспешных попыток размещения -мера). Например, если в решетке осталось место всего для одного к-мера (с учетом ориентации), то вероятность его расположить , т.е. невероятно маленькая. Более того, даже если решетка полностью заполнится, то программа об этом не узнает, и будет дальше пытаться размещать -меры, поэтому процесс размещения должен быть как-либо модифицирован.

В начале работы алгоритма, когда решетка практически пуста, достаточно рационально использовать случайное размещение, и “холостых” итераций будет незначительное количество. Поэтому имеет смысл сначала располагать иглы именно так, а в какой-то момент начать использовать другой алгоритм. Например можно создать список всех свободных клеток [10] и выбирать случайно уже из них, а после каждой итерации удалять из списка рассматриваемая клетку. В таком случае программа заканчивает свою работу, когда этот список становится пустым.

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

По сравнению с временем работы всей программы, составление списков правильных клеток занимает незначительное время, поэтому в программе используется именно этот алгоритм. Более того, заполнение горизонтального и вертикального списка можно запустить параллельно. Также экспериментальным путем было получено, что выгоднее всего составлять списки, когда решетка заполнена на 60%-65%.

2.2 Кластеризация -меров

Чтобы чем найти перколяционный кластер, необходимо знать, каким образом находить и маркировать кластеры из -меров. В [3,6] авторы описали union-find алгоритм для решения похожей проблемы. Несмотря на то, что их RSA сильно отличалась от рассматриваемой в этой статье, их идеи могут помочь с кластеризацией -меров. Основная идея заключается в том, чтобы в каждой клетке -мера хранить номер кластера, к которому он принадлежит. Union-find алгоритм состоит из двух частей.

Find: После добавления нового -мера проверяются соседние с ним клетки. Если все соседние клетки являются свободными, то этот -мер образует новый кластер.

Union: Если на стадии Find был найден хотя бы один кластер, то -мер присоединяется к нему. При нахождении следующего кластера эти кластеры сливаются в один.

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

2.3 Проверка перколяции

Для того, чтобы выделить из всего множества кластеров перколяционный, можно использовать алгоритм Machta [3,5]. Он подразумевает выделение из всего кластера одной клетки, которая будет называться “корневой”. Помимо того, для каждой клетки кластера нужно хранить два числа: расстояние до корневой клетки по соответствующей оси координат. Стоит отметить, что это расстояние может принимать отрицательные значения и при объединении двух разных кластеров, эти расстояния для каждой клетки должны изменяться вместе с номером кластера.

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

2.4 Экономия памяти

Основной проблемой при моделировании этого процесса является недостаток памяти. Алгоритм подразумевает, что в каждой клетке решетки нужно хранить как минимум 3 целых числа типа int (которые суммарно занимают 12 байт): номер кластера и расстояния до корневой клетки по каждой из осей. Так как количество клеток в решетке может достигать (для ), то это значит, что для работы программы может понадобиться порядка 12 терабайт оперативной памяти, что невероятно много, поэтому хранение информации требует серьезной оптимизации.

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

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

· 1 бит - клетка занята/свободна

· 1 бит - ориентация -мера, к которому принадлежит эта клетка

· 1 бит - является ли клетка начальной клеткой -мера

· 1 бит - дополнительная информация о -мере

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

-мера (двигаться по -меру в направлении влево или вверх в зависимости от его ориентации до тех пор, пока не встретится начальная клетка). Если -мер достаточно длинный (, то его первые 96 клеток можно разделить на 3 части по 32 клетки, как показано на Рисунке 3. Каждой части будет соответствовать одно из необходимых для работы программы чисел. Чтобы записать число, нужно разбить его на 32 бита, и записать их в соответствующие клетки, используя четвертый бит дополнительной информации. Соответственно, чтобы прочесть информацию из -мера, нужно проделать обратную операцию, т.е. прочесть 32 бита, а затем соединить их в одно число.

Такой подход к использованию памяти позволяет в каждой клетке решетки хранить всего 0.5 байта, что в 24 раза меньше, чем изначальный вариант, тем самым открывая возможности смоделировать решетки гораздо больших размеров.

2.5 Параллелизация с помощью OpenMP

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

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

Для того, чтобы разные нити не пытались расположить -мер в одном и том же месте, только одна нить должна выполнять эту функцию в каждый промежуток времени. Этого можно добиться используя критическую секцию в OpenMP. В ней программа должна выполнить следующие действия:

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

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

· расположить -мер

· пометить все окружающие кластеры как занятые в дополнительном массиве

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

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

Эту проблему можно решить немного изменив логику появления на решетке новых кластеров: при размещении -мера он сразу будет являться новым кластером (если возле него есть соседние -меры, то он впоследствии поменяет свой номер). Этот прием позволяет гарантировать отсутствие на решетке аномалий в виде нулевых кластеров.

Такой подход позволяет значительно ускорить вычисления.

Результаты

С помощью программы были смоделированы значения отношения плотностей перколяции и джамминга (Рисунок 4). Синим цветом обозначены значения, которые совпали с известными ранее значениями, полученными в [1,2]. Чтобы их получить была использована отдельная программа, с другим распределением памяти, так как вышеизложенный метод предполагает, что , и не подходит для работы с -мерами меньшей длины. Красным цветом обозначены новые, впервые полученные значения. Последняя точка на графике соответствует .

На Рисунке 5 плотности перколяции и джамминга представлены отдельно, красным и синим цветом соответственно.

Для каждой точки моделирование проводилось дважды, чтобы повысить точность вычислений. Сильно большее количество сильно повлияет на скорость работы, при том, что не даст сильного выигрыша в точности, потому что средняя погрешность незначительна. Для подтверждения этого факта решетка с была смоделирована 200 раз, при этом среднее квадратичное отклонение отношения плотностей оказалось равным 0.00511863.

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

Заключение

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

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

Список источников

[1] Tarasevich Yu, et al "Percolation of linear k-mers on a square lattice: From isotropic through partially ordered to completely aligned states." Physical Review E 86.6 (2012): 061116.

[2] Lebovka N., et al. "Random sequential adsorption of partially oriented linear k-mers on a square lattice." Physical Review E 84.6 (2011): 061603.

[3] Newman M., et al. "Fast Monte Carlo algorithm for site or bond percolation." Physical Review E 64.1 (2001): 016706.

[4] Kondrat G., et al. "Jammed systems of oriented needles always percolate on square lattices." Physical Review E 96.2 (2017): 022154.

[5] J. Machta et al., “Invaded cluster algorithm for Potts models”, Phys. Rev. E 54, 1332 (1996).

[6] Дональд Кнут, Искусство программирования, том 1. Основные алгоритмы

-- 3-е изд. -- М.: «Вильямс», 2006.

[7] D. Stauffer and A. Aharony, Introduction to Percolation Theory, 2nd edition,

Taylor and Francis, London (1992).

[8] Тарасевич Ю.Ю., Перколяция: теория, приложения, алгоритмы: Учебное пособие. М.: Едиториал УРСС, 2002.

[9] J. Talbot, G. Tarjus, P. R. Van Tassel, and P. Viot, From car parking to protein adsorption: an overview of sequential adsorption processes, Colloids Surf. A 165, 287 (2000).

[10] Тарасевич Ю.Ю., частное сообщение

Размещено на Allbest.ru