CC BY
126
2009 Вычислительные методы в дискретной математике №3(5)
ВЫЧИСЛИТЕЛЬНЫЕ МЕТОДЫ В ДИСКРЕТНОЙ МАТЕМАТИКЕ
УДК 621.391.1:004.7
ДИСКРЕТНОЕ МОДЕЛИРОВАНИЕ ФИЗИКО-ХИМИЧЕСКИХ ПРОЦЕССОВ1
О. Л. Бандман
Институт вычислительной математики и математической геофизики СО РАН,
г. Новосибирск, Россия
E-mail: bandman@ssd.sscc.ru
В статье дано систематическое описание собственных результатов по исследованию дискретных моделей кинетики физико-химических процессов на микро- и наноуровнях. Модели являются расширениями классического клеточного автомата (КА) фон Неймана, отличаясь от него тем, что в них допускаются произвольные, в том числе вероятностные, функции переходов над подмножествами состояний клеток, а также асинхронные и смешанные режимы функционирования. Для математического описания моделируемых процессов используются формализмы «Алгоритма параллельных подстановок». Приводятся условия корректности и оценки эффективности параллельных реализаций для синхронных и асинхронных КА-моделей. Все представленные модели иллюстрируются результатами компьютерного моделирования.
Ключевые слова: дискретное моделирование, мелкозернистый параллелизм, асинхронный клеточный авт,ом,ат,, метод Монте-Карло, поверхностная химия, кинетика наносистем, параллельные вычисления.
1. Введение и неформальная постановка задачи
Исследование дискретных моделей пространственной динамики началось в 80-х годах прошлого столетия. Интерес к дискретному моделированию явился следствием нового взгляда на клеточный автомат (КА) фон Неймана [1], который стал рассматриваться как дискретное выражение пространственно-временной функции [2]. Во всех точках дискретного пространства (клетках) значения этой функции (состояния клеток) изменяются по одним и тем же правилам (правилам переходов), которые выражаются в виде булевых функций от состояний клеток из некоторой окрестности. Все клетки могут перейти в новые состояния одновременно или в любом порядке, изменив таким образом глобальное состояние КА. Если установить КА в начальное глобальное состояние, он начинает эволюционировать, переходя из одного состояния в другое, имитируя пространственно распределенный процесс, т. е. выполняя ту же задачу, что и дифференциальное уравнение в частных производных. Принципиальное различие между этими двумя моделями состоит в том, что дискретный характер правил переходов позволяет выразить присущую физико-химическим процессам нелинейность
хРабота поддержана Программой фундаментальных исследований Президиума РАН №2.15 (2009), Сибирским отделением РАН, Интеграционный проект 32 (2009).
и прерывность самым непосредственным образом, отображая события, составляющие моделируемый процесс: перемещения, химические взаимодействия и фазовые превращения некоторых абстрактных или реальных частиц, а в некоторых случаях молекул и атомов [3].
Такой «новый» взгляд на КА обусловлен, главным образом, двумя факторами. С одной стороны, появлением в нашем распоряжении супермощной вычислительной техники, которая позволяет моделировать КА с миллиардами клеток, что действительно требуется в реальных задачах, поскольку для моделирования поведения микро- и наночастиц нужны очень большие клеточные пространства. С другой стороны, выявившейся слабостью традиционных математических моделей, опирающихся на дифференциальные уравнения в частных производных, не способных описывать процессы и явления, которые на современном этапе развития науки и техники выходят на первый план, а именно: химические реакции, развитие живых организмов, поведение социума. Эти явления существенно нелинейны, диссипативны, синергетичны, обладают способностью к самоорганизации, самовоспроизводству и другими экзотическими свойствами. Поиск новых подходов к моделированию таких явлений пробудил интерес к дискретным методам моделирования вообще и особый интерес к КА как модели пространственной динамики.
Очевидно, что классический фоннеймановский КА с булевым алфавитом, булевыми правилами переходов и синхронным режимом функционирования не обладает достаточно широкими моделирующими способностями. Поэтому появилось множество модификаций КА, которые допускают целочисленный и символьный алфавиты, произвольные функции переходов, всевозможные структуры дискретных пространств и режимы функционировани, однако сохраняют два следующих свойства.
1) Функции переходов для всех точек пространства одинаковы и локальны, т. е. значение следующего состояния зависит только от состояний клеток из близкого окружения.
2) Алгоритм вычисления нового глобального состояния обладает внутренним параллелизмом: функции переходов клеток либо полностью независимы по переменным (синхронные КА, или СКА), либо частично независимы (асинхронные КА, или АКА), что делает их параллельную реализацию простой и эффективной.
Эти два свойства определяют широкий класс дискретных вычислительных моделей, объединенных понятием «мелкозернистый параллелизм» [5]. Дискретные модели физико-химических процессов, которым посвящена статья, составляют его подкласс, называемый иногда диффузионно-реакционными КА [5]. При этом диффузионная составляющая включает в себя имитацию пространственных перемещений частиц, а реакционная — имитирует химические преобразования (реакцию, диссоциацию, фазовые превращения). Значительную долю диффузионно-реакционных КА составляют асинхронные КА с вероятностными правилами переходов, которые иногда называют методом Монте-Карло.
Состояния клеток могут быть либо булевыми 0, 1, имитируя факт наличия или отсутствия частицы в клетке, либо символами C,C0,W,WC, если частицы характеризуются химическими свойствами, либо целыми числами 1, 2,..., указывая на число частиц в клетке. К диффузионно-реакционным КА не относятся КА-модели, называемые «решеточным газом» (Lattice-Gas Models) или КА-гидродинамикой [7], у которых состояниями служат векторы единичных скоростей. Этот класс КА-моделей наиболее развит. Он получил широкое применение с названием Lattice-Gas или Lattice-Boltzmann Models и составляет самостоятельное направление в моделировании. Исключив его
из нашего рассмотрения, мы можем называть интересующие нас КА-модели просто клеточными автоматами, или КА.
Практическое применение КА получили в научных исследованиях для имитации физико-химических процессов, таких, как химические реакции на катализаторах [8, 9], эпитаксиальный рост кристаллов [10], образование покрытий и эрозии металлов [11] и др., где они успешно заменяют методы молекулярной динамики, которые требуют слишком большого времени вычисления даже на самых мощных компьютерах. К сожалению, большинство известных применений выполнено без опоры на какую-либо теорию или методологию и носят самые разные названия: решеточная модель [12], эволюционный алгоритм [13], микроскопическая модель [14], метод Монте-Карло [9, 10] и др. Почти все эти модели имеют ограничения по времени из-за трудности параллельной реализации на кластерах и отсутствия каких-либо методов преобразований и оптимизации. А между тем понятие КА полностью покрывает их моделирующие возможности, что позволило сделать некоторые обобщения [6], применить к ним единую формальную модель [15] и на ее основе разработать условия вычислительной корректности [16] и методологию эффективного моделирования, в том числе для большеразмерных задач на суперкомпьютерах [17]. Результаты таких исследований разрознены и фрагментарны, что затрудняет их осмысливание и использование. Чтобы продвинуть изучение и, следовательно, применение КА-моделей, предпринимается попытка в рамках этой статьи объединить уже опубликованные и еще не известные результаты по теории и практике их применения.
Кроме введения и заключения, статья содержит еще три раздела. Во втором разделе даны определения понятий и формальное определение модели КА и условий его корректного функционирования. Третий раздел посвящен методологии моделирования сложных процессов, четвертый — многопроцессорной реализации КА-алгоритмов.
2. Формальное представление и свойства КА-модели
2.1. Основные понятия и определения
Для формального представления КА далее используется математическая модель, называемая «Алгоритм параллельных подстановок» (АПП) [15]. Эта модель, разработанная для описания архитектуры параллельных алгоритмов и структур на СБИС, оказалась удобной для представления КА-моделей по двум причинам: 1) в ней определены контекстные конфигурации, необходимые для задания вероятности и временных ограничений, и 2) в ней допускаются многоклеточные правила переходов (один оператор меняет состояния нескольких клеток).
Согласно принципам АПП, центральным понятием в КА является клетка, которая задается парой (а,т), где а Е А — состояние клетки из алфавита А, и т Е М — имя клетки из множества имен М. Значения состояний клеток обозначаются символами из первой половины латинского алфавита а,Ь,..., переменные состояний — буквами из его второй половины п,ь,и,х,у, г. Множество имен в общем случае счетное, но в конкретных примерах обычно оно конечно и состоит из векторов координат в Д-мерном дискретном пространстве. Поскольку на данном этапе исследованы АКА при Д ^ 2, то М = {(%,]) : I = 0,1,...,1; ] = 0, 1,...,3}. Символ т вместо (%,]) используется в общих формулах для краткости.
Конкретные значения А и М определяют класс клеточных массивов П(А, М) = = П С А х М, каждый представитель П(А, М) которого не может иметь двух клеток с одинаковыми именами. На множестве имен вводятся именующие функции <рк : М ^ М, к = 1, 2,...,д для некоторого д ^ 1. Их значения связывают с лю-
бой клеткой т Е М множество соседних клеток, совместно с которыми она образует локальную конфигурацию:
5 (т) = {(го,т), (^1,^1(т)),..., (ги, ^ (т)),..., (г,, (т))}. (1)
Пронумеровав и упорядочив соседей, локальную конфигурацию (1) можно задать двумя векторами:
и(5(т)) = (го,..., V,),
Т (5 (т)) = (т, ...,^и (т),...,<£, (т)). (2)
Первый вектор называется вектором состояний для конфигурации 5(т), второй — ее определяющим шаблоном.
Три локальные конфигурации
5(т) = |(^о, ^о(т)), (гь ^(т)),..., (г,, (т))},
5/(т) = {(г0, ^о(т)), К ^М^..^ (г, (m))}, (3)
5кк(т) = |(мо,^о(т)), (щь ^(т)),..., (м„,^„(т))},
где ^о(т) = ^о(т) = т, записанные в виде выражения
0(т) : 5(т) *5"(т) — 5к(т), (4)
составляют параллельную подстановку (далее просто подстановку). Эпитет «параллельная» означает, что она может быть применена ко всем т Е М параллельно.
Локальная конфигурация 5(т) в (4) называется базовой конфигурацией (или базой), 5кк(т) — контекстом, а 5к(т) — правой частью параллельной подстановки. База и правая часть подстановки имеют один и тот же определяющий шаблон. Шаблон контекста от них отличается. Состояния г'к в клетках правой части являются значениями функций перехода от значений состояний в клетках базы и контекста, т. е.
г'к = ¡и(го,г1,... ,г,,що,...,««). (5)
Локальная конфигурация 5(т) с переменными состояниями считается совместимой с клеточным массивом П(А, М), что обозначается как 5(т) С П, если
1) для всех (гк,^к(т)) Е 5(т), к = 0,1,...,д, и всех (щ,^(т)) Е 5кк(т), I =
= 0,..., п, переменные ги и щ определены в алфавите А;
2) одноименные клетки с константными состояниями (а, ^(т)) Е 5(т) и
(Ь, ^к(т)) Е П равны, т. е. а = Ь.
Параллельная подстановка 0(т) применима к клеточному массиву П(А, М), если в нем найдется непустое подмножество клеток Мк С М с именами тк Е Мк, для которых удовлетворяется условие
5(тк) 5кк(тк) С П. (6)
Применимая подстановка называется локальным оператором.
Применение локального оператора 0(т) к П в клетке тк Е Мк состоит в замене состояний в клетках (аи, (тк)) на значения функций ¡к, вычисленных по (5) (рис. 1).
Выполнение локального оператора происходит за время т, называемое тактом.
Применение 0(т) ко всем клеткам клеточного массива П(£) приводит к изменению
его глобального состояния, П(£) ——) П(£ + 1), и называется итерацией. Последовательность
Е(П) = П, П(1),..., П(*), П(* + 1),..., П(Т), (7)
Рис. 1. Графическое представление локального оператора
полученная в результате итеративного функционирования КА, называется эволюцией.
Эволюция КА является дискретным представлением пространственно-временной функции. Если итеративный процесс сходится, то эволюция КА имеет завершение, т. е. существует £ = Т, такое, что П(Т) = П(Т + 1) = П(Т + 2) = ... В иных случаях КА моделирует колебательный или хаотический процесс [18].
Левая часть локального оператора (4) разделена на базу и контекст для того, чтобы при применении его к любой конкретной клетке ти Е М отделить множества клеток, состояния которых она изменяет, от множества клеток 5кк(ти), состояния которых она не меняет, но использует в качестве переменных для вычисления функции переходов. При этом состояния клеток из контекста 5кк(ти) могут быть изменены на той же итерации применением этого же локального оператора к другой клетке тй, если
5кк(ти) П £(т5) = 0. (8)
Контексты подразделяются на контексты первого рода и контексты второго рода.
— Контексты первого рода содержат клетки вида (щи, (т)), в которых ■фи (т) опре-
делены в области эволюции КА и, чаще всего, находятся в непосредственной близости от базовых клеток.
— Контексты второго рода выполняют управляющие функции. Их локальные конфигурации содержат клетки, в которых состояния являются значениями предикатов применимости основного локального оператора, а имена не содержатся в пространстве моделирования и задаются отдельным подмножеством Мкк:
5кк(т) = {(а1 ,т1к),..., (аг , т")}, (9)
где {а1,... , аг} = Акк — контекстный алфавит, {тк/,... , т"} = Мкк — множество имен контекстных клеток и Мкк П М = 0.
Для вычисления состояний контекстных клеток необходимы дополнительные контекстные операторы вида
0£ : (хи ,тки) — (аи ,т4к), к = 1, 2,..., г. (10)
Состояния а и Е Акк являются предикатами, зависящими от внешних переменных хи, например значений генератора случайных чисел, или значений тактирующих импульсов, или указателей области значений имен клеток и др.
В диффузионно-реакционных КА обычно используютя следующие три типа управляющих контекстов, они все имеют вид (10) и различаются только функциями переходов а и (хи).
Тактирующий контекст применяется для организации п-тактного синхронного режима. В нем используется оператор с предикатом аТ, который зависит от номера шага или итерации:
Вероятностный контекст применяется в вероятностных КА, в его операторе вр (10) предикат ар зависит от случайного числа,
где rand — случайное число в интервале [0,1]; p — вероятность применимости подстановки.
Пространственный контекст применяется для выделения тех областей пространства моделирования, где соответствующий оператор применим. В его операторе 0s предикат зависит от имени клетки m Є M:
где М' С М.
Составной контекст применяется, когда сразу несколько условий необходимо выполнить, чтобы основной локальный оператор был применим. В частности, для случайного выбора координат клетки (г,^) Є М, к которой применяется локальный оператор при асинхронном режиме функционирования, нужен пространственновероятностный контекст с предикатом вида
Однако, поскольку для всех АКА контекст (14) одинаков, его можно убрать из выражений локальных операторов, оставив только индекс а при Н.
Управляющая роль контекстов второго рода распространяется также на распределение работы между несколькими основными локальными операторами в и (т),... ,
в, (т), участвующими в описании единого физико-химического процесса. Совместно с контекстными операторами они составляют систему локальных операторов
Пример 1. Вероятностный КА, моделирующий простую одномерную диффузию, имеет булев алфавит А = {0,1} и пространство моделирования, состоящее из 15 клеток М = {0,1... , г,... , 14}. Концентрация вещества моделируется количеством «единиц» в некоторой окрестности осреднения Аг(г). Пусть Аг(г) = {г — 1, г, г + 1}, |Аг(г)| = 3. Исходный клеточный массив П(0) = {(1, г) : г = 1,... , 4} и {(0, г) : г = 5,..., 15}. Процесс диффузии состоит в перемешивании вещества с целью получения равномерно распределенной концентрации. КА моделирует его путем обмена состояниями между соседними клетками. Поскольку у каждой клетки два соседа, то выбор левого (г — 1) или правого (г + 1) соседа выполняется операторами в1(г) или в2(г) с вероятностью р = 0,5. Кроме того, каждая итерация разбивается на три такта £о,£1,£2 таким образом, чтобы на к-м такте для каждого к операторы применялись к клеткам с именами г = k(mod 3). Для этого нужен тактирующий контекст, указывающий
1, если t mod n = 0, 0, если t mod n = 0.
(11)
если rand < p, если rand ^ p,
(12)
если m Є M', если m Є M',
(13)
если (i = randi ■ N) & (j = rand2 ■ N), в иных случаях.
(14)
0 = {0i(m),..., 9q(m), 0!',..., 0"}.
номер такта, и пространственный контекст для выделения подмножеств клеток, работающих на соответствующем такте. Таким образом, система локальных операторов
0 = {01 (т), 02(т), 0" ,0" ,0"} имеет вид
0і(і) : {(и, і), (г, і - 1)} * {(1, тт), (1,шр)(1,ш8)} - -»■ КМ),(и і -1)},
02(і) : {(и, і), (г, і + 1)} * {(1,тт), (0,шр)(1,ш8)} - і) і+1)},
0т'' : (¿,шт) — (ат (і),тт),
0Т : (гап^2,тр) — (ар(топ^2), тр),
0" : (і,ш8) — (а8(і),ш8),
где ат вычисляется по (11) с п = 3, ар(гап^) вычисляется по (12) с р = 0,5, и а8(к) — по (1.3) с М' = {і : і = fc(mod3)}.
Эволюция процесса (15) показана на рис. 2. Она не имеет завершения, а приходит к колебательному процессу вокруг глобального состояния при і > 26, в котором осредненное значение плотности (и(і)) — 1/3.
—► і
0 1 2 .... 14
36 | | ■ | ■ | ■ | ■ ■
Рис. 2. Эволюция КА одномерной диффузии (15)
Если контекст второго рода может быть выражен простым предикатом, то его можно заменить надписью над стрелкой в выражении локального оператора. Например, в (15) 01 (г) может иметь следующий вид:
01 (г) : {(и, г), (г, г - 1)} —> {(г, г), (и, г - 1)}. (16)
2.2. Режимы функционирования и условия корректности
Применение параллельных подстановок к клеткам клеточного массива может происходить в разных режимах. Режим определяет порядок применения локальных операторов к клеткам массива за время одной итерации. Основными режимами являются: синхронный, обозначаемый символом а, и асинхронный, обозначаемый символом а. Соответственно синхронный КА обозначается , а асинхронный — . Названные два
режима являются двумя крайними случаями. Известны также КА со смешанными и производными режимами, например блочно-синхронный КА [19].
При синхронном режиме изменение состояний происходят только после того, как все функции /(г0,... , ) для всех клеток вычислены, т. е. известны новые состояния
всех клеток массива. Если предположить, что все клетки реализованы аппаратур-но и выполняют вычисления функций перехода одновременно, реализуя потенциальный внутренний параллелизм КА, далее называемый клеточным параллелизмом, то время итерации равно одному такту. В обычном последовательном компьютере синхронная параллельная смена состояний клеток имитируется путем вычисления новых
состояний всех клеток последовательно при сохранении их в другом массиве, который служит глобальным состоянием для следующей итерации. Такая параллельность называется виртуальной, и вычисление нового глобального состояния занимает время Д* = |М| • т тактов.
При асинхронном режиме ни реальной, ни виртуальной одновременности смены состояний клеток не предполагается. На каждом такте происходит применение 0(т) к одной случайным образом выбранной клетке с именем т Е М, которая вычисляет функцию перехода от тех состояний соседей, которые на данный момент имеют место, и сразу производит замену текущего состояния на новое. Это значит, что функции переходов / (5) могут иметь в качестве аргументов как текущие, так и новые состояния, например (гк, (т)) € П(*), а (гг,^(т)) € П(* + 1). Одна итерация выполняется за
Д* = |М| • т тактов, т. е. за то же время, что и при виртуальном синхронном режиме. Внутренний параллелизм проявляется только в том, что смены состояний клеток происходят в произвольном порядке.
Для построения КА-моделей важным является тот факт, что два КА, отличающиеся только режимами, могут иметь совершенно разные эволюции.
Пример 2. Процесс агрегации частиц одной жидкости в смеси из двух жидкостей может быть представлен классическим КА, получившим название КА разделения фаз [20]. Пусть в начальном состоянии жидкость В (черная) размешана в жидкости С (белой). Концентрация В по объему равна рв(0). Как только размешивание прекращается, в смеси начинается агломерация частиц вещества В под действием меж-молекулярных сил связи. В двумерном приближении процесс выглядит как появление черных пятен, которые постепенно растут и меняют свою форму и размеры. Рост пятен постепенно замедляется, приближаясь к одному из двух устойчивых состояний: либо «глобальное черное», либо «глобальное белое». При этом существует некоторое критическое значение начальной концентрации £, такое, что если рв(0) > £, то процесс стремится к «глобальному черному», если рв (0) < £, то процесс стремится к «глобальному белому», и процесс застывает на каком-то образе при рв (0) = £. КА, моделирующий этот процесс, имеет алфавит А1 = {0,1}, множество имен М = {(г,]) : г,] = 0,... , N} и следующий локальный оператор
где (і,і) є М £''(і,і) = {(^(і,і)): ^к(і,і) = (і + + ^)}, є {-2,-1,12}
Исследование КА (17) показало, что его эволюции при синхронном и асинхронном режимах сильно различаются. На рис. 3 приведены четыре состояния из эволюции синхронного КА при рв(0) = 0,5. При і = 3000 глобальное состояние больше практически не меняется, т. е. критическая концентрация д равна 0,5. При асинхронном режиме д = 0,2245, и даже при вдвое меньшей начальной концентрации КА быстро эволюционирует к черному глобальному состоянию, показывая почти незаметные изменения на итерациях і > 200 (рис. 4).
Пример 2 служит иллюстрацией к тому факту, что КА-модель определена полностью, только если указан режим ее функционирования.
Индексом а обозначается синхронный режим, индексом а — асинхронный и индексом в — блочно-синхронный. Алгоритм функционирования КА Н полностью опреде-
(17)
и
1, если в < 24 или в = 25,
0, если в > 25 или в = 24,
22
в = ^2 ^2 гі+д,і+}і.
д=-2Н=-2
лен, если, кроме этого, задан исходный клеточный массив ^(0), т. е.
«к = (А,М, 0, ОД), к е|а,в,^}.
(18)
Рис. 3. Эволюция КА «разделения фаз» (18) при синхронном режиме функционирования
Рис. 4. Эволюция КА «разделения фаз» (18) при асинхронном режиме функционирования
2.3. Условия корректного функционирования
Как всякая математическая модель, КА должен удовлетворять условиям поведенческой и вычислительной корректности. Эти условия для КА просты и естественны, но легко могут быть нарушены при крупноблочном распараллеливании.
Условие детерминированности: локальный оператор 9(т) : £(т) ^ £'(т) в один и тот же момент т может изменять состояние не более чем одной клетки из каждого подмножества клеток £' (т) для всех т Є М.
Это условие гарантирует отсутствие коллизий, т. е. таких ситуаций, когда одной и той же клетке на одном и том же такте предписано сменить свое состояние разными функциями переходов /к (т) и / (^ (т)). Это значит, что одновременность (пусть виртуальная) применения 9 (т) допускается только к такому подмножеству клеток Мк С М, для которого выполняется условие
Утг,тЛ Є Мк (шд = т.,, ^ Т(5'(т„)) |~| Т(5'(тЛ)) = 0 .
(19)
Условие (19) не выполнено, например, в синхронном КА, если две клетки должны изменить состояния одновременно (рис. 5).
Очевидно, что (19) всегда выполняется для классических КА, так как их локальные операторы имеют одноклеточную базовую локальную конфигурацию, т. е. |£"(т)| = 1, и также в асинхронных КА, где одновременно может быть применен только один оператор и только к одной клетке.
Рис. 5. Графическое представление недетерминированного поведения. Одновременное применение в (г) и в (г+1) приводит к конфликту в клетках г, г + 1,г + 2
Условие равноправия клеток: на 1-й итерации для каждого £ локальный оператор 9(т) применяется ко всем клеткам (V, т) Є ^(£), причем к каждой клетке только один раз.
Это условие гарантирует, что все клетки участвуют в процессе эволюции в равной степени. Нарушение этого условия может быть вызвано несовершенством генератора случайных чисел, а также неправильной организацией взаимодействий между областями КА, если они расположены на разных процессорах. Очевидно, что синхронные КА всегда удовлетворяют условию равноправия, так как режим предусматривает детерминированный порядок применения локального оператора ко всем клеткам из ^(і).
3. Моделирование физико-химических процессов
3.1. Локальные операторы физико-химических процессов
Накопленный опыт по дискретному моделированию процессов на микро- и наноуровнях связан, в основном, с изучением химических реакций на каталитических поверхностях, причем в большинстве изученных моделей используется один и тот же набор элементарных действий. Пусть А = {а, Ь, с, 0}, где а, Ь, с обозначают наличие частицы определенного вещества, а 0 — свободное место, М = {(і, і) : і, і = 0,..., N}. Тогда типичные элементарные действия поверхностной химии можно отобразить в КА-моделях следующим образом.
1. Адсорбция моделирует попадание частицы на пустую клетку с заданной вероятностью ра, которая вычисляется исходя из физических данных (парциальное давление газа над поверхностью, скорость струи, направленной на поверхность, и др.):
9а : (0, (і,і)) -^ ^ (і,і)). (20)
Адсорбция иногда сопровождается диссоциацией, например, молекула Н2 или 02, попадая на поверхность катализатора, разлагается на два атома. Тогда ей требуются две соседние пустые клетки, которые могут быть смежными по оси і:
92а : {(0, (і,і)), (0, (і + 1, і))} {(Ь, (і,і)), (Ь, (і + 1,і))}, (21)
или по оси і:
92а : {(0, (і,і)), (0, (і,і + 1))} ^ {(Ь, (і,і)), (Ь, (і,і + 1))}. (22)
Десорбция, или сублимация, — действие, обратное адсорбции, обозначающее освобождение клетки. Поскольку простое исчезновение противоречит природе явлений, то обычно десорбция указывает на перемещение частицы из области моделирования во
внешнюю среду. Локальный оператор имеет вид (20), в котором левая и правая части поменялись местами, т. е.
0* :(а, (г,7)) (0, (^)). (23)
Химическая реакция — действие, в котором участвует несколько разных веществ. Здесь возможно несколько вариантов.
1) Если две частицы, вступающие в реакцию а + Ь ^ с, расположены в смежных клетках и продукт реакции остается в области моделирования, то он может занять любую из двух освобождающихся клеток равновероятно, что должно быть представлено четырьмя локальными операторами, например:
9іг
92г
93г
94г
{(а(і,і)),(Ь, (і +1,і))} {(а(і,і)), (Ь, (і +1,і))} {(Ь, (і,і^ (a, (і +1,і))} {(Ь, (і,і)), (а(і +1,і))}
Рг
1-Рг
Рг
1 Рг
{(0 (і,і))(с(і +1,і))}, {(с (і +1,і))}
{(0, (і,і))(c, (і +1,і))}, {(c, (і +1,і))},
(24)
где рг = 0,5.
Если реакция сопровождается сублимацией, то для ее описания достаточно одного оператора (24), у которого в правой части обе клетки пусты.
Диффузия. Процесс диффузии есть блуждание частиц в стремлении выровнять концентрацию вещества в пространстве. Наиболее естественным образом этот процесс отображается вероятностным АКА, впервые введенным в [19] и названным там «наивной диффузией». АКА наивной диффузии имеет булев алфавит, дополненный множеством контекстных символов Л" = {0,1,... , N}, и множество имен М = {(¿,,7) :
г,] = 0,1,... , N}, дополненное контекстным множеством М" = {ms,mi,mj}. Система локальных операторов содержит один основной локальный оператор
Ой: {(и (г,7')), К(г,7))} * С?',т-), (к,тк)} ^
^ {(ик, (г,7')), (и,^к(г,7'))}, один вероятностный контекстный для случайного выбора соседа <^>к(¿,7)
9'к : (х,тк) ^ (к,тк), к =
и два контекстных оператора для случайного выбора клетки
(25)
1, если 0 < гап^1 ^ 0,25,
2, если 0,25 < гап^1 ^ 0,5,
3, если 0,5 < гап^1 ^ 0,75,
4, если 0,75 < гап^1 ^ 1
(26)
9''
(27)
(ж,т^) ^ (г,т^), г = гап^ ■ N,
0;" : (х,т-) ^ (7,т-), 7 = гап^2 ■ N.
Если диффузия анизотропна, то интервалы вероятностей в (26) не одинаковы, а распределены в соответствии со значениями коэффициентов диффузии вдоль соответствующих направлений.
На рис. 6 показаны три глобальных состояния эволюции КА двумерной наивной КА-диффузии. Исходным состоянием служил заполненный нулями клеточный массив размерами 200 х 200, в центральной части которого имеется квадратная область высокой концентрации частиц:
Ш = {(1, (і,і)): і,і = 70,..., 130}.
Рис. 6. Наивная асинхронная КА-диффузия: глобальные состояния клеточного массива при і = 0, 250, 600
Из рис. 6 видно, что процесс диффузии постепенно превращает плотный квадрат в круг с уменьшающейся от центра концентрацией. Для получения модельного коэффициента диффузии производилось осреднение концентрации в каждой точке по окрестности с радиусом г = 8 клеток. Сопоставление осредненных клеточных массивов с решениями дифференциального уравнения Лапласа позволило вычислить модельный коэффициент диффузии Сто^, который должен быть равен безразмерному коэффици-
^ * т
енту в уравнении Лапласа С = , где d (м2/с) — физический коэффициент диффу-
Я2
зии, т (с) — шаг по времени, Я (м) — шаг по длине в дифференциальном уравнении [21]. Соотношение коэффициентов С и Стоа позволяет вычислить масштабы, связывающие дискретные модельные параметры с реальными физическими значениями.
3.2. Композиции локальных операторов
Построение КА-моделей сложных процессов требует использования в одном алгоритме нескольких взаимодействующих локальных операторов, которые составляют композиции разного типа [22]. При моделировании физико-химических процессов используются чаще всего суперпозиция локальных операторов и вероятностная последовательная композиция.
Оуперпозиция нескольких локальных операторов — это один локальный оператор вида
9(т) = 9|(9|_1 (... (91 (т)))),
(28)
который может применяться как в синхронном, так и в асинхронном режиме, причем компонентные операторы всегда применяются к клетке в заданном в (28) порядке и не могут быть разделены применением другого локального оператора к другой клетке. Вероятностная последовательная композиция множества операторов
@(т) = {91 (т), 92(т),... в\(т)}
(29)
предполагает применение компонентных операторов в случайном порядке к случайно выбранной клетке, причем в течение одной итерации каждый Ок(т) Е 0(т) применяется к каждой клетке т Е М по одному разу. Такой режим функционирования называется вероятностным асинхронным, или в [8, 9] методом Монте-Карло.
Пример 3. Упрощенная АКА-модель химической реакции окисления моноокиси углерода на катализаторе, называемая иногда ZBG-реакцией по инициалам авторов [23], подробно исследована в [8, 9]. Над поверхностью катализатора имеется смесь
двух газов: кислорода и моноокиси углерода с парциальными давлениями рь и ра соответственно. В АКА-модели поверхность катализатора, на котором происходит реакция, представлена клеточным массивом ^(Л, М), в котором М = {(¿,7) : ¿,7 = = 0,1,...^}, Л = {а, Ь, 0}, причем (а, (¿,7)), (Ь, (¿,7)) и (0, (¿,7)) обозначают клетки, в которых находятся молекулы СО, О, или клетка пуста соответственно. В алгоритме используется вероятностная последовательная композиция множества {О1 (г, 7), #2(г, 7), #3 (г, 7), #4(г, 7)} локальных операторов, моделирующих следующие элементарные действия.
1) Адсорбция СО из газа — О1 (¿,7): если клетка (¿,7) пуста, ее занимает молекула СО с вероятностью ра.
2) Адсорбция О2 из газа — #2 (г, 7): если клетка (г, 7) пуста и имеет пустого соседа, то с вероятностью рь в каждой их них оказывается молекула кислорода. Один из четырех соседей клетки (¿,7) выбирается с вероятностью рп = 0,25.
3) Реакция окисления СО (СО+О СО2) — #3(¿,7): если в клетке (¿,7) находится молекула СО и в соседней клетке находится молекула кислорода, то между ними происходит реакция. В результате получается молекула СО2, которая переходит в газ. Обе клетки опустошаются. Соседняя клетка выбирается с вероятностью рп = 0,25.
4) Та же самая реакция окисления СО (О+СО СО2) происходит, если в клетке (¿,7) находится молекула кислорода, а в соседней клетке имеется молекула СО —
#4 (¿7):
*1 (¿,7): {(0, (¿,7))} {(а, (¿,7))},
#2(г,7) : {(0, (г,7'))(0,^кз(г,7))} * (к2, т1) {(Ь (г,7'Ж (Ь (г,7))}, (30)
*з(г,7) : {(а, (г,7'))(Ь,^кэ(г,7))} * (кз, т2) —^ {(0, (г,7')), (0, (^кэ(г,7))},
*4(г,7) : {(Ь, (г,7'))(а,^к4(г,7))} * (к4,тз) —^ {(0, (г,7'Ж (0, (^к4(г,7))}.
В #2(¿,7),#3(¿,7),#4(¿,7) состояния к2, к3, к4 Е {1, 2, 3, 4} указывают на одного из четырех соседей клетки (г, 7), а именно: ^1 (¿,7) = (¿,7 + 1), ^2(¿,7) = (г + 1,7), ^3(¿,7) = = (г, 7 — 1), ^4(¿,7) = (г - 1,7). Значения к2, к3, к4 вычисляются соответствующими вероятностными контекстными операторами #2", #3", #47 вида (26).
На рис. 7 показаны три глобальных состояния, полученных при моделировании реакции на клеточном массиве размером 200 х 200 с периодическими граничными условиями.
Рис. 7. Три глобальных состояния эволюции АКА . Черные пикселы обозначают CO, серые — O, белые — пустые места
4. Многопроцессорные реализации КА-алгоритмов
Свойства мелкозернистости КА-алгоритмов предопределяют способ их многопроцессорной реализации с использованием декомпозиции области моделирования [17].
Пусть клеточный массив Q(A, M) имеет структуру прямоугольника, который разделен на n равных по количеству клеток доменов размерами |Dom| = |M|/n клеток. Эффективность распараллеливания выражается формулой
Q(n) = ШТ t , (31)
n * ldom + L * lex
где ln,ldom — времена, необходимые для выполнения одной итерации на всем КА и на одном домене соответственно, т. е. Iq = т • |M|; ldom = т • |Dom|; lex — время, затрачиваемое на один межпроцессорный обмен; L — число обменов за одну итерацию. Очевидно, что для получения высокой эффективности необходимо, чтобы
т • |Dom| >> L • tex. (32)
В синхронных КА обмен данными производится один раз за итерацию, т. е. L =1. Если принять 1ех/т ~ 103 (что соответствует реальности), то условие (32) легко выполняется при |Dom| ^ 104. Иными словами, высокая эффективность параллельных реализаций достигается без труда.
Для асинхронных КА межпроцессорный обмен должен выполняться всякий раз, когда изменяется состояние хотя бы одной клетки, находящейся на границе домена. Если это условие не выполняется, то не выполняются условия корректности (19), поскольку при стохастическом выборе клеток невозможно гарантировать, что новое состояние пограничной клетки (v,mg), mg G Domk, не понадобится на следующем же такте для вычисления нового состояния пограничной клетки соседнего домена (u, mi ), если mi G Domi, mi = ^JT1 (mg), ^(mg) G Ts. Отсюда следует, что для асинхронных КА за одну итерацию надо выполнить L = P пересылок, где P — периметр домена. Легко проверить, что при соотношении 1ех/т ~ 103 приемлемая эффективность (>70%)
достигается при |Dom| ^ 107, что делает распараллеливание нецелесообразным.
Поскольку препятствием для эффективного распараллеливания является асинхро-низм, то возникает соблазн ввести частичную синхронизацию режима функционирования, при этом
а) по возможности, сохранить эволюцию КА;
б) не нарушить условий корректности.
Это можно сделать, трансформировав асинхронный КА = (A, M, 0) в блочносинхронный КА = (A, M, 0) следующим образом [16].
1. На множестве имен M для каждого m G M определяется шаблон
B(m) = |^о(m), <^i(m),..., ^(m), <p,+i(m),..., ^(m)}, называемый блоком, который включает в себя определяющий шаблон Т(S(m)) базовой локальной конфигурации S(m):
Т(S(m)) с B(m). (33)
2. Строится разбиение M на независимые подмножества M1,..., Mr, для которых выполняются следующие соотношения:
|Mfc| = M, k = 1,... , r; (34)
r
У Mk = M, Vk,/ G {1,..., r} (k = / ^ Mk Q Mi = 0)
k=1
У В(ш) = М, Ушд, шг Є Ык (шд = шг ^ В(шд) ^ В(шг) = 0^ . (36)
тЄМ&
Иными словами, вводится г разбиений множества имен М блоками одного размера, причем каждый блок В (ш) включает в себя базовый шаблон Т (£ (ш)).
3. Смена глобального состояния П(і) ——> П(£ + 1) производится за г стадий. На каждой стадии производится применение всех локальных операторов ко всем клеткам из одного случайно выбранного независимого подмножества М^, і = 1, 2,... , г. Режим функционирования в пределах стадии — асинхронный, т. е. порядок выбора клеток в независимом подмножестве произвольный. Это не нарушает условия (19), так как функции переходов не пересекаются по переменным, что гарантируется условием (33). Следовательно, результат глобального перехода к следующей стадии — единственный и одинаков как для синхронного, так и для асинхронного режимов.
Тот факт, что при блочно-синхронном режиме на каждой стадии функции переходов разделены по переменным, позволяет производить обмен данными между доменами после завершения вычислений на каждой стадии. При этом ни одно условие корректности не нарушается, следовательно, допускается формирование пакета, в котором содержатся новые состояния клеток из обрабатываемого на этой стадии независимого подмножества. Объем V пакета равен с • Раот/г, где Раот — периметр домена, г = |В(ш)|. С учетом того, что стадии обрабатываются последовательно, условие эффективности распараллеливания (32) принимает вид
т • |Бош| > г • ¿еж. (37)
Из (37) следует, что чем больше размер блока г, тем больше время обменов. Следовательно, блок надо выбирать минимального размера, при котором удовлетворяются условия (19) и (33).
Пример 4. Процесс преобразования химической энергии в электрическую в водородном топливном элементе происходит следующим образом. Водород поступает на
пористый графитовый анод, активированный катализатором (платиной). Попадая на
платину, атомы водорода разлагаются на электроны и протоны:
Н2 ^ 2Н+ + 2е-. (38)
Протоны направляются к катоду, где, соединяясь с кислородом, образуют воду. Электроны поступают в электрическую цепь, создавая ток.
Диссоциация водорода на аноде определяет производительность топливного элемента, поэтому изучается особенно подробно. Моделирование его позволяет определить зависимость получаемого тока от парциального давления поступающего водорода и от структуры платиного покрытия в порах графитового анода.
Процесс моделируется АКА = (А, М, 0). В нём алфавит состояний есть А = (0, Н, е}, где 0 означает, что клетка свободна, Н — клетка содержит атом водорода, е — клетка содержит е электронов, е Є (0,1,...}. Множество имен клеток есть М = ((і,і) : і,і Є 0,1... , N} и (Е,ш8,шр} и соответствует клеткам на поверхности пор анода; Е — счетчик электронов; ш5,шр — контекстные клетки, предикаты в которых соответствуют вероятности ра попадания атомов водорода на платину, равной доле активированной поверхности, и вероятности выполнения диссоциации, равной парциальному давлению водорода рн. Множество 0 содержит пять основных локальных операторов (Ое, Он (к) : к =1, 2, 3, 4} и два контекстных вероятностных оператора вида (10):
°е : № (і,і)), (е,Е)}^((0, (м)), (е + 1,Е)}, ( )
Он (к): ((0, (і,і)), (0,^ (і,і))} {(Н, (і,і)), (Н,^ (і,і))}, к =1, 2, 3,4, ( )
где (i,j) G T(S(i,j)) — k-я клетка, смежная с (i,j); pH — парциальное давление водорода; ps — доля площади пор, активированной платиной.
Моделирование проводилось на клеточном массиве с размерами |М | = 10000 х 10000 клеток, что соответствует площади поверхности анода Sa = 10-8 мм2 при доле платинового покрытия ps = 0,5. Время выполнения 100 итераций составляло t = 64,4 с. Для того чтобы получить данные для площади диссоциации Sa ~ 1 мкм2, было выполнено моделирование на 36 процессорах с такими же размерами доменов. Для параллельной реализации АКА был преобразован в блочно-синхронный с блоками размера |B(m) | =3 х 3. Таким образом, итерация состояла из 9 асинхронных стадий с синхронной передачей данных (рис. 8).
ОД 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9 1
а б
Рис. 8. Моделирование процесса диссоциации водорода на аноде топливного водородного элемента: а — платиновое покрытие поверхности поры при ps = 0,5; б — зависимость количества полученных электронов от парциального давления поступающего водорода
Вычисления проводились на кластере MVS-1000 в ССКЦ СО РАН. Время вычисления 1000 итераций составило 1080 с, что соответствует эффективности распараллеливания Q = 0,6 и позволяет надеяться, что на 1000 процессорах за приемлемое время (несколько часов) можно моделировать получение тока с площадей около 1 см2.
Заключение
Представленные в статье дискретные методы компьютерного моделирования позволяют исследовать физико-химические процессы «изнутри», наблюдая химические и фазовые преобразования, а также перемещения частиц на микро- и наноуровне. Методы основаны на клеточно-автоматной модели природных явлений, которая является альтернативой дифференциальным уравнениям, и позволяют существенно расширить возможности математического и компьютерного моделирования. Поскольку в современных технологиях наблюдается стремление создавать материалы, полученные на основе превращений микро- и наночастиц, и разработка современной технологии невозможна без компьютерного моделирования, то можно надеяться, что КА-моделирование станет востребованным не только в науке, но и в промышленности.
ЛИТЕРАТУРА
1. Toffolli T. Cellular Automata as an Alternative to (rather than Approximation of) Differential Equations in Modeling Physics // Physica D.1984. V. 10. P. 117-127.
2. Wolfram S. Statistical mechanics of Cellular automata // Rev. Mod. Phys. 1993. V. 55.
P. 607-640.
3. Boon J. P., Dab D., Kapral R., Lawniczak A. Lattice-Gas Automata for Reactive Systems //
Phys. Rep. 1996. V. 273. P. 55-147.
4. Бандман О. Л. Мелкозернистый параллелизм в вычислительной математике // Программирование. 2001. №4. С. 1-18.
5. Бандман О. Л. Клеточно-автоматные модели пространственной динамики // Системная информатика: Сб. научн. тр. Новосибирск: Изд-во СО РАН. 2006. Вып. 10: Методы и модели современного программирования. C. 59-111.
6. Bandman O. Synchronous versus asynchronous cellular automata for simulating nano-systems kinetics // Bulletin of the Novosibirsk Computer Center. Series: Computer Science. Issue 27.
2006. P. 1-12.
7. Rothman B. H., Zaleski S. Lattice-Gas Cellular Automata. Simple Models of Complex Hydrodynamics. London: Cambridge Univ. Press, 1997. 320 p.
8. Makeev A. G. Coarse bifurcation analysis of kinetic Monte Carlo simulations: a
lattice-gas model with lateral interactions // J. Chem. Phys. 2002. V. 117. No. 18. P. 8229-8240.
9. Elokhin V. I., Latkin E. I., Matveev A. V., Gorodetskii V. V. Application of Statistical Lattice
Models to the Analysis of Oscillatory and Autowave Processes on the Reaction of Carbon
Monoxide Oxidation over Platinum and Palladium Surfaces // Kinet. Catal. 2003. V. 44. No. 5. P. 672-700.
10. Neizvestny I. G., Shwartz N. L., Yanovitskaya Z. Sh., and Zverev A. V. 3D-model of epitaxial growth on porous {111} and {100} Si surfaces // Comp. Phys. Commun. 2002. V. 147. P. 272-275.
11. Betz G., Husinsky W. Surface erosion and film growth studied by a combined molecular dynamics and kinetic Monte Carlo code // Izvestia of Russian Academy of Sciences. Phys. Ser. 2002. V.66. No. 4. P. 585-587.
12. Kovalev E. V., Resnyanskii E. D., Elokhin V. I., et al. Novel statistical lattice model for the supported nanoparticle. Features of the reaction performance influenced by the dynamically changed shape and surfaces morphology of the supported active particle // Phys.Chem.Chem.Phys. 2003. V. 5. P. 784-790.
13. Alba E., Troya J. M. Cellular Evolutionary Algorithms: evaluating the influence of ratio //
Lect. Not. Comp. Sci. 2000. V. 197. P. 29-38.
14. Malvanets A., Kapral R. Microscopic model for Fitz-Nagumo dynamics // Phys. Rev. E. 1997. V. 55. No. 5. P. 5657-5670.
15. AchasovaS., Bandman O., Markova V., Piskunov S.. Parallel Substitution Algorithm. Theory and Application. Singapore: World Scientific, 1994. 198 p.
16. Bandman O. Parallel Simulation of Asynchronous Cellular Automata Evolution // Lect. Not. Comp. Sci. V. 4173. 2006. P. 41-48.
17. Бандман О. Л. Параллельная реализация клеточно-автоматных алгоритмов моделирования пространственной динамики // Сибирский журн. вычислительной математики.
2007. №4. C. 45-361.
18. Wolfram S. A new kind of science Champaign, Ill., USA: Wolfram Media Inc, 2002. 2000 p.
19. Toffolli T., Margolus N. Cellular Automata Machine. USA: MIT Press. 1987, 284 p.
20. Vichniac G. Simulating Physics by Cellular Automata // Phys. D. 1984. V. 10. P. 86-112.
21. Bandman O. Comparative Study of Cellular automata Diffusion Models // Lect. Not. Comp. Sci. 1999. V. 1662. P. 395-399.
22. Бандман О. Л. Методы композиции клеточных автоматов для моделирования пространственной динамики // Вестник Томского госуниверситета. 2002. №9(1). С. 188-192.
23. ZiffR.M., Gulari E., Bershad Y. Kinetic phase transitions in irreversible surface-reaction model // Phys. Rev. Lett. 1986. V. 56. P. 553-2558.
