CC BY
109
2017 Дискретные модели реальных процессов №35
ДИСКРЕТНЫЕ МОДЕЛИ РЕАЛЬНЫХ ПРОЦЕССОВ
УДК 681.3.06: 681.323
КЛЕТОЧНО-АВТОМАТНЫЕ МОДЕЛИ ЕСТЕСТВЕННЫХ ПРОЦЕССОВ И ИХ РЕАЛИЗАЦИЯ НА СОВРЕМЕННЫХ КОМПЬЮТЕРАХ
О. Л. Бандман
Институт вычислительной математики и математической геофизики СО РАН,
г. Новосибирск, Россия
Представлены результаты анализа моделирующих способностей и вычислительных свойств методов клеточно-автоматного (КА) моделирования нелинейных пространственно распределенных процессов. Работа преследует две цели: 1) показать соответствие свойств КА-моделей современным тенденциям развития параллельных многопроцессорных архитектур (дискретность представления данных, локальность взаимодействий) и 2) раскрыть возможности КА-методов для компьютерного моделирования естественных существенно нелинейных, диссипативных процессов, не поддающихся традиционным методам математического моделирования. Обобщается опыт КА-моделирования, полученный в ИВМиМГ СО РАН. Работа содержит формальное представление КА-моделей, методы их построения, а также результаты реализации ряда тестовых и реальных задач на суперкомпьютерах.
Ключевые слова: математическое моделирование, дискретная математика, параллельные вычисления, клеточно-автоматные модели пространственной динамики, клеточно-автоматная гидродинамика, реакционно-диффузионные процессы.
DOI 10.17223/20710410/35/9
CELLULAR-AUTOMATA MODELS OF NATURAL PROCESSES, IMPLEMENTATION ON SUPERCOMPUTERS
O. L. Bandman
Institute of Computational Mathematics and Mathematical Geophysics SB RAS, Novosibirsk,
Russia
E-mail: bandman@ssd.sscc.ru
Conventional mathematical models based on differential calculus are sometimes not capable to simulate nonlinear dissipative processes on micro- or nano-level of resolution. This fact stimulates the development of new approaches to spatial dynamics simulation. Among them, cellular automata (CA) modeling is one of the promising methodologies, due to CA large simulation capability and compatibility with modern trends in supercomputer architecture. Although CA simulation is intensively studied and used in different fields, a few attention is paid to studying the parallel implementation peculiarities of large scale CA-models on supercomputers. Just this aspect of
CA-simulation is the subject of the paper aiming to analyse CA-simulation methods adaptiveness to supercomputing implementation based on validity conditions requirements for different modes of CA operation. For this purpose, the concept of the operation mode well known for simple CA (having only one transition rule) is expanded for composed CA (containing many transition rules). The new concept determines a CA-transition rules execution order, which in turn determines the behavioral properties of CA-model and their influence on the simulation performance. The obtained results are illustrated by some examples, which show CA methods at work by simulating essentially nonlinear and dissipative processes: superposition of asynchronous CAs for simulation of water permeating through porous medium and parallel composition of two CAs simulating pattern formation on a heated plate. Basically, the paper generalizes CA computer simulation theoretical results and experience obtained by researchers from the Institute of Computational Mathematics and Mathematical Geophysics of Siberian Branch of Russian Academy of Sciences.
Keywords: mathematical modeling, parallel computing, cellular automata models, Lattice-Gas hydrodynamics, reaction-diffusion processes.
Введение
Клеточный автомат, введённый фон Нейманом в середине прошлого века [1], в настоящее время рассматривается как математическая модель пространственно распределённого естественного процесса, в которой все компоненты (пространство, время и состояния) дискретны. Теория клеточных автоматов в её современном виде [2] составляет часть теории динамических систем. Созданная Пуанкаре в позапрошлом веке теория динамических систем опиралась на математические описания явлений в виде дифференциальных уравнений. Появление компьютеров и численных методов распространило её положения на конечно-разностные представления [3]. Наконец, развитие КА-моделирования привело к появлению раздела, изучающего полностью дискретные динамические системы [4]. Однако понятие КА возникло не в результате постепенной эволюции теории динамических систем к дискретному виду, а в результате развития идеи фон Неймана о возможности построения вычислительной модели устройства, воспроизводящего себе подобного. Это создавало впечатление о КА как об абстрактном представлении некоторого биологического сообщества, состоящего из множества «клеток», образующих регулярную пространственную структуру. Каждая клетка является элементарным вычислителем, который может находиться в одном из двух состояний — 0 и 1 (белое или черное) и изменять это состояние в зависимости от состояний клеток ближайшего своего окружения, называемого соседством. Алгоритм вычисления следующего состояния у всех клеток одинаковый. Все клетки выполняют переход в новое состояние одновременно, изменяя таким образом глобальную черно-белую картину распределения состояний по пространству.
Интерпретация и осмысливание концепции КА происходили вместе с развитием вычислительной науки и технологии и в эпоху современных суперкопьютеров приняли форму нового направления в математическом и компьютерном моделировании. Это вызвано тем, что стало ощущаться несоответствие между идеологией, на которой основана работа компьютера, и идеологией, на которой построена вычислительная математика. Противоречия между ними сводятся к следующим двум положениям.
1. Компьютер работает с дискретными данными в дискретном времени, тогда как вычислительные модели основаны на непрерывной математике. Это противоречие имеет два последствия. Во-первых, при переходе к вычислительным алгоритмам диффе-
ренциальные уравнения аппроксимируются разностными схемами, что приводит к потере точности. Во-вторых, непрерывность обрабатываемых данных требует работы с вещественными числами, т. е. использования арифметики с плавающей точкой, которая отбрасывает малые разряды в процессе счёта. Это приводит к неточным результатам и, что ещё более важно, вызывает неустойчивость вычислительного процесса. Последнее является сильным ограничением для выбора шага по времени при использовании явных схем дискретизации.
2. Тенденции развития современных суперкомпьютеров направлены в сторону увеличения количества параллельно работающих процессоров, тогда как вычислительная математика остаётся на фундаменте последовательных алгоритмов. Возникает проблема распараллеливания, которая решается не всегда эффективно.
Поиск новых моделей, свободных от этих противоречий, заставил пересмотреть концепцию КА. Новый взгляд на КА как на модель пространственной динамики, изложенный в фундаментальных работах [5, 6], привёл к осознанию того факта, что, несмотря на простоту каждой клетки, их кооперативная работа может моделировать сложные и разнообразные процессы, которые иногда невозможно (или, по крайней мере, неизвестно как) описать другим способом. И, что немаловажно, имея одного и того же «родителя», КА и компьютер близки по «духу и букве», что улучшает их совместимость. Настоящим потрясением основ в моделировании пространственной динамики стало появление клеточно-автоматной гидродинамики, названной «Lattice-Gas Cellular Automata» (LGA — модель вязкой жидкости), для которой была доказана эквивалентность уравнению Навье — Стокса [7]. Вслед за LGA появились КА-модели процессов диффузии [6], разделения фаз [8], реакционно-диффузионных процессов [9, 10], знаменитой реакции Белоусова — Жаботинского [11], движения солитонов [12] и др. Теория и методология КА-моделирования интенсивно развиваются: проводятся ежегодно 3-4 международных конференции, выпускаются два специальных международных журнала («Complex Systems», США, и «Cellular Automata», Великобритания). В результате оформилось новое направление компьютерного моделирования, в котором понятие «клеточного автомата» фон Неймана играет роль идеологической основы. К сожалению, в России только отдельные разрозненные группы занимаются КА-мо-делированием.
Цель работы — показать моделирующие способности КА-методов и их адаптируемость к современной параллельной вычислительной технике, иллюстрируя их результатами, полученными в течение ряда лет группой исследователей Института вычислительной математики и математической геофизики СО РАН. В п. 1 даны формальные определения основных понятий и классификация КА-моделей; п. 2 посвящён КА-моде-лям гидродинамики и сопоставлению их с непрерывными аналогами. В п. 4 представлены композиции КА, моделирующие сложные процессы. В заключении формулируются выводы из изложенного и проблемы дальнейшего развития.
1. Клеточно-автоматные модели пространственной динамики 1.1. Формальное представление клеточно-автоматной
модели
Формально, КА-модель определяется четвёркой понятий [13, 14]
N = (A,X, 0,р),
где X — множество клеток, представленных своими именами; A — алфавит состояний клеток в X; 0 — множество локальных операторов клеток, называемое также глобаль-
ным оператором; р — режим функционирования модели. Состояниями клеток могут быть числа, булевы векторы или символы, именами клеток — натуральные числа или векторы координат точек дискретного пространства. Клетка х Е X в состоянии а Е А записывается как пара (а,х). Таким образом, А х X есть множество всех клеток во всевозможных их состояниях. Любое его подмножество П = {(а^х^,... , (а8,х8)}, где {х1,... , х8} = X, называется клеточным массивом КА-модели, а набор состояний клеток в нём (а1,... , а8) —его и КА-модели (глобальным) состоянием.
Для определения локального оператора вводятся именующие функции ^ : X ^ X, г = 1, 2,... , д, такие, что для каждого х Е X все клетки х, ^1(х),... , (х) различные. Их подмножества
N(х) = {х,^1(х),...,^га(х)}, Е(х) = {^га+1(х),...,^,(х)}, 0 ^ п ^ д,
называются соответственно соседством и окрестностью клетки х, а их подмножества вместе с некоторыми их состояниями
5(х) = {(ио,х), («1, ^1(х)),..., (м„, р«(х))}, С(х) = {(м„+1, ^га+1(х)),..., («, Рд(х))}
— соответственно локальной конфигурацией и контекстом для этой клетки. Множества всевозможных локальных конфигураций и контекстов для клетки х обозначаются Б (х) и С(х) соответственно. Таким образом,
Б(х) = {5(х) = {(«о,х), («1, (х)),..., (м„, р«(х))} : «о,«1,... ,«п Е А}, С(х) = {С(х) = {(«п+1, ^«+1 (х)),... , («, Рд(х))} : м„+1,... Е А}.
Локальные операторы в в ставятся в соответствие клеткам так, что каждой клетке соответствуют один или более таких операторов. На множестве в(х) операторов, сопоставленных клетке х, может быть задано распределение вероятности, в соответствии с которым может определяться порядок их выбора для применения в процессе функционирования КА-модели в том или ином режиме. Произвольный локальный оператор в в в(х) обозначается в(х) и называется оператором клетки х. Это есть отображение в(х) : Б(х) х С(х) ^ Б(х), определяемое следующим образом: результат Б'(х) = в(х)[Б(х), С(х)] его действия над произвольной парой [5(х), С(х)] Е Б(х) хС(х) получается преобразованием локальной конфигурации Б (х) клетки х к её новой локальной конфигурации
Б'(х) = в(х)[Б(х), С(х)] = {(«0,х), («1, Р1(х)),..., (и«,^„(х))}
путём замены состояний м0,м1,...,мга клеток соседства N(х) новыми состояниями «0, «1,... ,«« соответственно, вычисленными как
«к = /к(«0, «1,..., ), к = 0,1,..., п,
где /к : Ад+1 ^ А и векторная функция / = (/0/1... /) называется функцией переходов.
Далее в случае пустого контекста, т.е. когда д = п, вместо в(х)[Б(х),С(х)] пишем в(х)[Б(х)]. Кроме того, локальную конфигурацию Б(х) и контекст С(х) под знаком оператора в(х) будем относить не только к клетке х, но и к оператору в(х), называя их его локальной конфигурацией и контекстом соответственно.
Нетрудно заметить, что данное определение локального оператора в(х) фактически совпадает с определением конечного автомата с входным алфавитом Ад, множеством
состояний А, выходным алфавитом Ап и функциями переходов — : А9 х А ^ А и выходов : А9 х А ^ Ап, удовлетворяющими соотношениям — (п^ ... пппп+1... п9, п0) = «0 и ^(п1п2 ... пппп+1... п9, п0) = п1 п2 ... пП, где п0 и «0 являются соответственно текущим и следующим состояниями клетки х; п1; п2,... , пп и п1, п'2,... , «п — соответственно текущие и новые состояния соседних для х клеток; пп+1, пп+2,... , п9 — текущие состояния клеток из контекста клетки х. Именно поэтому рассматриваемые здесь модели реальных процессов носят название клеточно-автоматных.
При заданном начальном клеточном массиве П(0) функционирование КА-моде-ли К происходит в дискретном времени £ = 0,1,... , каждый промежуток которого называется итерацией. Результатом £-й итерации является клеточный массив П(£ + 1) = в(П(£)). Он получается применением локальных операторов в(х) Е в ко всем клеткам (п,х) Е П(£) в порядке, определяемом заданным режимом, и тем самым совершается глобальный переход П(£) ^ П(£ +1) от текущего клеточного массива к следующему. Применение в(х) к клетке (п,х) Е П(£) означает замену состояний клеток из Б(х) в П(£) на состояния одноименных клеток из Б'(х) = в(х)[Б(х), С(х)].
Здесь уместно обратить внимание на роли локальной конфигурации Б(х) и контекста С(х) как частей клеточного массива П(£) в определении локального оператора в(х). Вместе они образуют аргументы, по значениям которых вычисляется значение Б'(х) на выходе в(х). С другой стороны, в результате выполнения оператора значение контекста С(х) не изменяется, а на место конфигурации Б(х) записывется новая — Б'(х).
Говорят, что глобальные переходы в КА-модели совершаются корректно [15], а вместе с ними корректна и КА-модель, если в процессе вычислений П(£ +1) по П(£) для всех £ ^ 0 не происходит ни одной коллизии, т. е. ни одной попытки изменить состояние одной и той же клетки более одного раза в один и тот же момент времени Формально это условие выражается так: для любых вк(х),в,(у) Е в, х,у Е X и Б(х) Е Б(х), Б (у) Е Б (у)
х = у ^ (Ж (в* (х)[Б (х), С (х)]) п N (в, (у)[Б (у), С (у)]) = 0), (1)
где N(Б'(х)) = {х, ^(х),..., <£п(х)}.
Результат итеративного применения глобального оператора в к массивам П(£) описывает эволюцию моделируемой системы в форме последовательности клеточных массивов П(0), П(1)... , П(£),... Если этот итеративный процесс сходится, т.е. существует такое, что П(£) = П(£ + 1) = ..., то система приходит к устойчивому состоянию. В иных случаях КА моделирует колебательный или хаотический процесс [5].
Реальные процессы, как природные, так и рукотворные, обычно состоят из множества различных элементарных действий. Соответственно в их КА-моделях используется много разных локальных операторов [14]. В их обозначениях буквой в последняя в рамках одной модели помечается далее разными числовыми индексами, как это сделано в (1). Для указания типа КА-модели, к которой относится данный оператор, или режима, в котором он применяется в ней, эта буква помечается символами, указывающими на эти тип модели или режим. Аналогичным образом метятся и буквы, обозначающие разные глобальные операторы и клеточные массивы в рассматриваемых моделях.
Различают три уровня сложности КА-моделей: 1) простые КА-модели, имитирующие одно элементарное действие (диффузия, конвекция, фазовые переходы, химические реакции и т.д.), используют один локальный оператор, |в| = 1;
2) сложные КА-модели, в них глобальный оператор работает с множеством локальных операторов, |в| > 1;
3) композиции глобальных операторов вг, такие, как последовательная (суперпозиция глобальных операторов) и параллельная (система глобальных операторов). В них все в г работают последовательно или параллельно, обновляя предписанные им подмассивы Пг С П и используя в качестве контекстов наборы клеток из всего массива П.
Во всех случаях итерацией считается исполнение глобального оператора в(П(£)), т.е. процедура, состоящая из IX| • |в| шагов применения всех вг(х) Е в(х) ко всем («, х) Е П(£). Порядок исполнения этих шагов определяется режимом функционирования р, который, в случае сложной модели, имеет две составляющие: пространственную рх, определяющую порядок выбора клеток, и операторную , определяющую порядок выбора локальных операторов. В случае композиций глобальных операторов каждый оператор работает в своем режиме, а взаимодействие между ними определяется типом их композиции: последовательная или параллельная.
1.2. Режимы функционирования Режим работы простой КА-модели характеризуется только пространственной составляющей рх, обозначаемой индексами из множества греческих символов {а, а, к, в}, которыми индексируются операторы, исполняемые в соответствии со следующими алгоритмами.
• В синхронном режиме (режим а) применение вст(х) ко всем клеткам (и,х) Е П(£) происходит в любом порядке. Новые значения (и',х) записываются во вспомогательный массив П', который после заполнения играет роль П(£ + 1). Для синхронных КА-моделей условие (1) обеспечивается применением к каждой клетке х операторов ва(х), не изменяющих состояний других клеток, что в классических синхронных КА выполняется по определению. При синтезе КА-модели этот факт является ограничением, зато при реализации синхронной КА-модели на многих процессорах параллельно обмен состояниями граничных клеток производится один раз за итерацию, что обычно гарантирует высокую эффективность параллельной реализации.
• В асинхронных режимах операторы ва(х) или вк(х) применяются к клеткам («, х) Е П(£), выбираемым последовательно в случайном (режим а) или заданном (режим к) порядке, причём текущие состояния клеток («•, ( (х)) Е 5(х) заменяются на новые («, ((х)) Е Б'(х) немедленно. При реализации на одном процессоре условие (1) выполняется автоматически, независимо от количества клеток в Б(х), изменяющих свои состояния одновременно. Однако при параллельной реализации на кластере асинхронный режим требует выполнения межпроцессорных обменов после обновления состояния каждой клетки, имеющей соседа на другом процессоре. Это делает процесс распараллеливания практически бессмысленным. Решение проблемы состоит в преобразовании асинхронного режима в блочно-синхронный (режим в), что, как показано на ряде примеров [16, 17], не искажает результатов моделирования, зато обеспечивает приемлемую эффективность параллельной реализации.
• Блочно-синхронный (режим в) [16] является промежуточным между асинхронным и синхронным. На множестве X строится разбиение П^) = {X1,X2,... ^т}, такое, что
ш > |N(Б'(х))|, (2)
VXfc Е П(Х)Ух,у Е X* [(х = у) ^ (Ж(Б'(х)) П N(Б'(у)) = 0)]. (2)
Итерация разбивается на ш стадий, исполняемых последовательно в любом (возможно, случайном) порядке, причем на к-й стадии операторы в в (х) применяются синхронно ко всем клеткам с именами в X*, к = 1, 2,... , ш. Режим позволяет преобразовывать локальные конфигурации любого размера п > 1, отсутствие коллизий при этом гарантируется условиями (2). При реализации вв на параллельных системах требуется ш межпроцессорных обменов за итерацию, что вполне приемлемо по эффективности.
Режим функционирования сложной модели характеризуется не только пространственной, но и операторной составляющей и обозначается парой символов (рх, ре), причём ре Е {$, 7} и указывает на один из двух разных способов назначения оператора в(х) Е в(х) для применения к выбранному х Е X: детерминированный (режим $) и вероятностный (режим 7). Наиболее известны следующие режимы сложных КА-моделей:
• Стохастический режим (р = (а, 7)). Он является максимально рандомизированным представлением поведения частиц на микроуровне. Глобальный оператор ва,7 выполняется следующим образом:
1) клетка (п,х) Е выбирается случайно (режим а);
2) локальный оператор в(х) выбирается среди операторов в в(х) по распределению частоты имитируемых ими действий и применяется к (п, х) немедленно.
Итерация состоит из IX| • |в| повторений п. 1 и 2. Стохастический режим используется в моделировании сложных химических реакций [17, 18].
• Блочно-синхронный детерминированный (р = ($, в)). Разбиение X на подмножества X!,... ,Xm выполняется согласно (2), с той лишь разницей, что
п
ш ^ 10 N (в,- (х)[Б (х), С (х)])|.
3 = 1
Итерация разбивается на ш стадий, которые выбираются в детерминированном порядке, на к-й стадии все локальные операторы применяются в режиме (а, $).
• Блочно-синхронный стохастический режим (р = (а,в)). Отличается от режима ($, в) тем, что стадии выбираются в случайном порядке. Применяется при параллельных реализациях КА-моделей, изначально заданных в стохастическом режиме (а, 7).
Композиция глобальных операторов объединяет в единую модель несколько глобальных операторов, каждый из которых может быть простым или сложным и работать в своём заданном заранее режиме. Наиболее известны два типа композиций:
• Последовательная композиция глобальных операторов ^(в1,... , в^}, работающих каждый в своем режиме, причём исходным клеточным массивом для в3 служит П3-1(£) = в3-1(П(£)). Итерация состоит из вычислений Ь промежуточных глобальных переходов:
+1) = вь(вь-1(... (в1 (ВД)))).
При реализации на суперкомпьютере последовательную композицию целесообразно размещать на кластерах, декомпозируя клеточный массив на подмассивы и организуя обмен состояниями их граничных клеток.
• В параллельной композиции [19] Ф(в(1),..., в(ь)) множество клеток разделено на Ь равномощных подмножеств (слоев): X = X(1) и ... и X(ь), XО) = {хО) : г = = 1, 2,..., N}, ] = 1, 2,..., Ь, IX| = N • Ь. Соответственно разделено и множество глобальных операторов: П = {П(1) и ... и П(ь)}. Все подмножества ПО) Е П обновляются одновременно соответствующими глобальными операторами вО") , при этом контексты локальных операторов в(г,Е в могут быть в любых других подмножествах из П.
При реализации КА-модели на суперкомпьютере глобальные операторы целесообразно размещать на ядрах компьютера с общей памятью, где каждый КА ведёт вычисления на своем ядре, пользуясь состояниями клеток всех КА из общей памяти [20].
1.3. Классификация КА-моделей Основными характеристиками всех КА-моделей, унаследованными от классического КА, являются дискретность пространства, времени и состояний, локальность межклеточных взаимодействий и итеративность глобальных переходов. Однако такие параметры КА-моделей, как алфавит состояний клеток, режим функционирования, тип и размер соседства, выбираются в соответствии со свойствами моделируемого явления и могут быть разного типа. На рис. 1 схематически показаны классы совре-
Рис. 1. Классы современных КА-моделей
менных КА-моделей [14]. Каждый исходящий из вершины путь к нижнему уровню соответствует определённому классу приложений. Классический КА [8] соответствует пути КА ^ булев ^ синхронный ^ детерминированный; КА-модели гетерогенных каталитических реакций [17, 18] лежат на пути КА ^ символьный ^ асинхронный ^ стохастический; целочисленная LGA-модель газового потока [21] соответствует пути КА ^ целочисленный ^ асинхронный ^ последовательная композиция.
С точки зрения методологии моделирования пространственной динамики КА-мо-дели целесообразно подразделить на две группы: КА-модели, имеющие аналоги в традиционной (непрерывной) математике, и КА-модели, для которых таких аналогов нет. Представителями первой группы являются КА-модели газовых и жидкостных потоков (Lattice-Gas Cellular Automata [5, 7, 22]). Представителями второй — многочисленные КА-модели химических, биологических, социальных явлений, фазовых переходов.
2. КА-модели, имеющие аналоги в непрерывной математике
2.1. КА-модели диффузии
Самой востребованной простой КА-моделью, имеющей аналог в виде дифференциального уравнения, является КА-модель диффузии. Её популярность объясняется тем, что она присутствует как компонент в процессах типа «реакция — диффузия». Она существует в нескольких видах: асинхронный [6] и синхронный [23] варианты, причём каждая модель имеет булеву и целочисленную версии. Асинхронный двумерный КА, называемый «наивной диффузией», —это простейший КА, в котором А = {0,1}; X = = {(г,3) : г = 0,1..., I; 3 =0,1,..., 3}; локальные операторы
в1 (г,3ЖСп (г,3)), Кк(г,3))}] = , (г>3)), (п,№3))}, к = 0,1, 2, 3, (3)
где в каждом применении оператора в,(г,3) соседняя клетка (г,3) выбирается с вероятностью рДг, 3) в интервале к/4 < рДг,3) < (к + 1)/4, к = 0,1, 2, 3. КА моделирует процесс изотропной диффузии с постоянным безразмерным коэффициентом диффузии Б = 0,5 [24].
2.2. КА-модели потока вязкой жидкости Известен ряд моделей газовых и жидкостных потоков, базовой моделью для которых служит КА-модель К = (А, X, в, (а, 5)}, получившая название РИР-модели (по первым буквам фамилий авторов [7]). Модель является представлением кинетического процесса движения и столкновения абстрактных частиц в гексагональном дискретном пространстве. Расстояния между центрами соседних гексагонов (клеток) равны единице (рис.2). Частицы могут перемещаться между центрами клеток с единичной скоростью. При столкновениях частицы меняют направления своего движения так, чтобы законы сохранения массы и импульса не нарушались (рис. 3). Множество имён клеток X = {х : х = (г,3), г = 0,... 1,3 = 0,... 3}. Соседство клетки х содержит, кроме неё самой, еще шесть её ближайших соседей N(х) = {х, ^1(х),..., ^6(х)} (рис. 2). Алфавит состояний есть А = : к = 1,..., 26}, где = (п1,... ,п6), п Е {0,1}, г = 1,..., 6. Если п = 1, то клетка (и>*, х) Е П содержит частицу, единичный вектор скорости е^ которой направлен к центру клетки ^¿(х) Е N (х).
Рис. 2. Клетка (ад,ж), ад = (010101), и её шесть соседей
Рис. 3. Две функции столкновения: трёхчастичная детерминированная (вверху) и двухчастичная вероятностная (внизу)
Глобальный оператор перехода ва,г (П) работает в режиме детерминированной суперпозиции синхронных глобальных операторов — сдвига вМ (П) и столкновения
©C(П). Оператор сдвига применяется ко всем (wk, x) G Q(t) синхронно. Он перемещает каждую частицу на одну клетку в направлении вектора её скорости. Это может быть выражено посредством локального оператора как
дм(x)[{(w, x)}, {(w'', ^(x)),..., (W', ^a(x))}] = (K x)}, (4)
где компоненты вектора нового состояния w' = (м1,... ,uj,... ,u'6) равны компонен-
II t (i) (i) (i)4
там векторов состояний контекстных для x клеток wi = (v1 ,... , vj ,... , v6 ) со сдвигом, сохраняющим направление движения соответствующей частицы, т. е. uj =
= (^(j + 2) mod 6 + 1) = Vj .
Глобальный оператор ©C (П) состоит из одного бесконтекстного локального оператора дс (x), заменяющего состояние одной клетки x на значение вероятностной функции перехода, которая меняет направления векторов движения частиц, имитируя их столкновение:
дс(x)[{(ui.. .U6,x)}] = {(u1.. .u6,x)}, (u1... u6) = f (ui.. .Мб). (5)
Функция перехода f удовлетворяет законам сохранения массы и импульса:
6 6 6 6
£ м = Е Е uiei = Е uiei. (6)
i=1 i=1 i=1 i=1
В трёхмерном случае в [22] предложено использовать дискретное пространство, состоящее из ромбододэкаэдров, имеющих 12 граней и, следовательно, соседство из 12 клеток (RD-модель). Сложность вычисления (табличного представления) булевых функций перехода в RD-модели вполне приемлема. Величину отклонения от полной изотропии аналитически не удалось получить. Однако вычислительные эксперименты по моделированию потока в трубе показали, что параболы Пуазейля в разных продольных разрезах трубы симметричны и неразличимы в пределах принятой в экспериментах точности. Отсюда сделан вывод, что модель применима там, где малое нарушение изотропии не играет роли, например при моделировании потока в пористой среде [25].
Преобразование результатов моделирования к привычным вещественным значениям скорости производится путём осреднения полученных значений по пространству и времени. Эта процедура сглаживает полученные зависимости, уменьшая или уничтожая так называемый «автоматный шум» (изломанный характер кривых, обусловленный дискретностью значений переменных). Для параллельных реализаций создан программный комплекс [26], работающий на кластере Межведомственного суперкомпьютерного центра РАН (г. Москва). На тестовых задачах эффективность параллельных реализаций при количестве процессоров n < 500 составила > 90 %, уменьшаясь до 70 % при увеличении числа процессоров до 1600.
2.3. Целочисленная КА-модель газового потока
Целочисленная LGA-модель потока (сначала названная автором мультичастич-ной) [21] предложена для смягчения негативных проявлений булевых представлений переменных: автоматного шума, неспособности моделировать турбулентные потоки, невозможности моделировать движение твёрдого тела в потоке и, следовательно, сжатия газа. Целочисленная модель отличается от булевой алфавитом состояний, который содержит не булевы, а целочисленные векторы: в двумерном случае Wk = (u1,... ,м6),
где п Е Ъ. Глобальный оператор, как и в булевой модели, является суперпозицией глобальных операторов сдвига (4) и столкновения (5). При этом для каждого состояния клетки правой части (5) функция переходов равна одному из всех возможных выбранных равновероятно значений т' = (п1.. .п'6), удовлетворяющих (6). Целочисленность алфавита позволяет моделировать движение твёрдых объектов в газе, имитируя сжатие газа. Чтобы реализовать это свойство, разработаны специальные локальные операторы взаимодействия частиц с перемещающимся твёрдым телом. Моделирование движения поршня в трубе и вдвигающейся задвижки, а также формирования порошковой струи под действием взрыва показало качественное сходство с моделируемыми процессами [21]. Тестовые вычислительные эксперименты по моделированию потоков в трубе показали, что число Рейнольдса И,е имеет значения 800 ^ 1000 и за препятствием хорошо видны отрывающиеся вихри.
2.4. Сравнение с непрерывными моделями
Модели КА-гидродинамики имеют аналоги в непрерывной математике. Поэтому имеет смысл провести сопоставление их свойств и возможностей с другими методами моделирования жидкостных потоков.
Трудности трёхмерной реализации и, главное, непривычность к булевому выражению физических величин привело к появлению более приближённой к традиционной математике модели, так называемой ЬаШсе-Сав-Во^шапп (ЬВМ) модели [27]. В ней дискретное состояние клеток заменено функцией распределения плотности, а функции столкновения — непрерывным оператором. По сути, модель получилась равной дискре-тизированной форме уравнения Больцмана. Идея построения вычислительного метода для моделирования гидро- и газодинамики в обход дифференциальных уравнений в частных производных реализована также в кинетически согласованных разностных схемах (КСРС) [28], где из кинетического описания среды на микроуровне выполняется непосредственный переход к разностным уравнениям. В отличие от ЬОА-моделей, которые находятся на стадии исследования, эти две модели хорошо развиты. Их распространение подтверждает целесообразность создания вычислительных алгоритмов моделирования пространственной динамики непосредственно из кинетических описаний. Это сокращает путь от представления моделируемого явления до его имитации на компьютере.
Действительно, для традиционных методов гидродинамики полный путь моделирования состоит из следующих этапов: 1) вывод дифференциального уравнения из законов сохранения; 2) дискретизация времени и пространства; 3) построение вычислительного метода; 4) программирование (чаще параллельное) и 5) компьютерная реализация (рис.4). На каждом этапе происходит потеря точности. При этом на пер-
Рис. 4. Схема построения моделей пространственной динамики
вом этапе выполняется переход от дискретного микропредставления к непрерывному макропредставлению, а на втором — обратная операция дискретизации пространства
и времени. На третьем этапе надо сделать выбор между эффективно распараллеливаемыми явными вычислительными методами с ограничением Куранта и неявными методами с неэффективной суперкомпьютерной реализацией. Наконец, на последнем этапе компьютер аппроксимирует вещественные числа булевыми векторами заданной разрядности. Для моделей КСРС и ЬБМ путь начинается с кинетического описания процесса, от которого, в каждом случае по-своему, совершается непосредственный переход к конечно-разностному представлению и явному методу программирования с более слабыми условиями Куранта. Наконец, ЬОА-модель, исходя из кинетического описания явления, выраженного булевыми функциями, совершает переход сразу к программе в двоичных кодах, что позволяет избежать преобразования чисел с плавающей точкой к двоичному коду. Это последнее исключает погрешности вычислений, вносимые округлением до заданной разрядности, и снимает проблему обеспечения вычислительной точности и вычислительной устойчивости. Так, в [29] показано, что при увеличении размеров области вычисления значений простых функций (вт(ж)) растёт размер ошибки округления. Поэтому при вычислениях на мощных суперкомпьютерах ЬОА-метод может оказаться предпочтительным.
3. Клеточно-автоматные модели, не имеющие непрерывных аналогов
3.1. Типичные компоненты сложных процессов КА-модель процесса пространственного разделения смесей двух веществ на компоненты («разделения фаз») имеет булев алфавит А = {0,1}, множество клеток X = {(г,3) : г,3 = 0,...,М} —в виде декартовой решётки, контекст клетки (г,3) представлен квадратом с центром в этой клетке и содержит д клеток
^ 3 Ж К (г, 3))} {(иъ 3)),..., (ия, <^я& 3))}] = (г, 3))} (7)
9 9
1, если Е ик > (д +1)/2 или Е ик < (д — 1)/2, „, _ , к=0 к=0 ^0 = \ я Я
0, если Е ик < (д — 1)/2 или Е ик < (д + 1)/2.
к=0 к=0
Режим работы КА-модели — синхронный. Из рис. 5, где показано четыре глобальных состояния из эволюции при исходном П(0), равном случайному равномерному распределению единиц в клеточном массиве с плотностью 0,5, видно, что процесс диссипа-тивный (энтропия уменьшается).
Рис. 5. Четыре глобальных состояния из эволюции КА «разделения фаз»
Ещё одна типичная простая КА-модель имитирует процесс формирования устойчивых образов [30]. В ней тоже используется булев алфавит и плоская декартова решётка.
Модель работает как в синхронном, так и в асинхронном режимах, хотя при этом эволюционирует совершенно по-разному. Локальный оператор имеет вид (7) с функцией перехода, заданной выражением
1, если пи> 0,
^о = < и=0
0, если пктк ^ 0, к=0
где ти —вещественные числа, интерпретируемые как «веса» состояний контекстных клеток (если ти < 0, то клетка — ингибитор, если ти > 0, то активатор) (рис. 6).
Рис. 6. Четыре глобальных состояния из эволюции асинхронной КА-мо-дели процесса формирования устойчивого образа
3.2. КА-модель гетерогенной химической реакции
на катализаторе Репрезентативным примером КА-моделирования сложного процесса может служить гетерогенная реакция окисления моноокиси углерода CO ^ CO2 на поверхности катализатора [30, 31]. В связи с реконструкцией кристаллической решетки платины под воздействием адсорбированных молекул CO структура поверхности может быть гексагональной (hex) или кубической (1 х 1). Элементарные события (абсорбция, реакция, десорбция, перестройка структуры поверхности, диффузия, реакция и др.) описываются девятью химическими уравнениями. КА-модель Н = (A,X, в, (а, 7)) этого процесса имеет следующие характеристики: алфавит состояний
A = {*1хЪ *hex, COLXs\ COhdX, °LXs\ OhdX}
содержит символы химических элементов, адсорбированных на поверхности двух типов; X = {(i,j) : = 0,1,...,M}, e(i,j) = {^1(i,j),..., 0g(i,j)}; локальные операторы в в моделируют все элементарные действия процесса реакции. Глобальный оператор в(П) работает в стохастическом режиме. Четыре глобальные состояния из эволюции процесса моделирования показаны на рис. 7. Моделирование выявило области параметров реакции, при которых процесс приходит к одному из возможных устойчивых состояний: стационарное однородное покрытие или автоколебания.
Реализация КА-модели на клеточном массиве размерами 8000 х 8000 выполнялась в блочно-синхронной версии на суперкомпьютере МВС-100К МСЦ РАН. Наибольшую эффективность распараллеливания (0,73 в сильном смысле) удалось получить при параллельном режиме на 128 ядрах.
tVM Фф'Щ&.Л f* ~ Ч 1 V ^ PL.fi; -dHwi « ■ ЧИЬьЬУ-Ч'ЧЖ^ G&tv;е- -Ч*'-41*^' V — О"
■ -Oads H"COads Q-*1x1 □-"hex
Рис. 7. Изменение состояний поверхности в ходе реакции COads ^ Oads на Pt100. Слева направо: образование островков COads, которые затем покрывают всю поверхность, на которой появляются островки Oads, распространяющиеся по всей поверхности
4. Композиции глобальных операторов клеточных автоматов 4.1. Последовательная композиция Ярким примером последовательной композиции (суперпозиции) глобальных операторов может служить КА-модель = (А, X, ^(в)) увлажнения почвы [32]. Алфавит состояний А = {м, 5, д, К] содержит символы, интерпретируемые как вода (м), твердая порода (в), газ (д) и разбухшее гидрофильное включение (К). Почва задана трёхмерным клеточным массивом П = {(и, (г, к)) : и € А, (г, к) € X, г,] = 0,1,..., 700, к = = 0,1,... , 1480] (рис. 8). В начальном состоянии на поверхности почвы — слой воды. Моделируется постепенное проникновение этой воды в пористую структуру почвы (рис. 9) с одновременным испарением и с учётом влияния гидрофильных включений, которые запирают поры при разбухании и открывают их при высыхании.
700
Рис. 8. Структура клеточного массива Рис. 9. Разрез фрагмента клеточного КА-модели увлажнения почвы массива при Ь = 10000 итераций
Глобальный оператор построен как суперпозиция пяти простых глобальных операторов р(в) = в£(в|(в§(в*(в£(П(0)))))), где — вс — оператор конвекции, содержащий локальный оператор
0с(г,^,к)[{(^, (г,^, к)), (д, (г,^,к + 1))]] = {(д, (г,^, к)), (м, (г,^,к + 1))], (8)
который моделирует движение частицы воды вниз под действием гравитации. Режим выполнения — асинхронный выбор клеток с вероятностью рс, упорядоченный вдоль вертикальной оси к вниз;
©к -
©E
оператор испарения, отличающийся от (8) взаимозаменой состояний т и д, значением вероятности ре и направлением упорядочивания асинхронного выбора клеток (вдоль оси к вверх); — в^ — оператор разбухания гидрофильной стенки, содержит локальный оператор
ds(i, j, k)[{(h, (i, j, k)), (w, M(i, j, k)) : l = 1, 2, 3, 4}] = {(h, (i,j,k)), (h,Mi(i,j,k)): l = 1, 2, 3, 4},
(9)
где Mz(i, j, k) = (i+ai, j+bi, k) и ai, bi для каждого l = 1, 2, 3,4 определяются как первые два числа в тройке, выбранной с равной вероятностью из четырёх последних троек в соседстве N(i,j,k) = {(i, j, k), (о, 1,k), (1,0,k), (0,-1,k), (-1,0,k)};
— ©Q — оператор высыхания моделирует выделение воды из разбухшей стенки при высыхании, он отличается от (9) взаимозаменой состояний w и h;
— ©D —оператор диффузии моделирует выравнивание поверхности воды в порах, кавернах и на поверхности, для полного выравнивания поверхности выполняется 20 раз на каждой итерации.
Поскольку КА-модель ориентирована на параллельную реализацию, операторы диффузии и гидрофильного влияния применяются в блочно-синхронном режиме с девятью стадиями.
Важным свойством предложенной КА-модели является тот факт, что параллельный алгоритм вычислений не требует выполнения обменов в трёх измерениях. При декомпозиции клеточного массива на вертикальные параллелепипеды плоскостями (k, i) и (k,j) обмены данными между подмассивами происходят только в горизонтальных плоскостях. Всего необходимо 9 • 22 = 198 обменов на каждой итерации. Клеточный массив разрезался на 7 • 7 = 49 подмассивов, размещаемых на 49 ядрах кластера НКС-30 ССКЦ СО РАН (Intel Xeon E5540, 2,53 ГГц). Для выполнения 50 000 итераций потребовалось 10 часов при эффективности распараллеливания примерно 70%.
5. Параллельная композиция глобальных операторов
Известно несколько КА-моделей с глобальным оператором в виде параллельной композиции Ф(©1(П),..., ©L(H)) [19, 33]. Большинство известных КА-моделей [30, 34] используют два глобальных оператора. Один из них моделирует исследуемое явление, а второй является контекстным, так как моделирует влияющие на него изменяющиеся внешние условия.
Типичным примером служит КА-модель процесса формирования узора на поверхности остывающей металлической пластины H = (A, X, Ф(©)}, где © = = {©а, ©а} и соответственно X = Xu U Xv с клетками (i, j)u в Xu и (i, j)v в Xv и П = Qu U Qv. Оба оператора простые, функционируют асинхронно, имеют булев алфавит состояний клеток. Оператор ©а имитирует распространение тепла по пластине от горячей области, где u = 1, к холодному квадрату в центре (u = 0) и к тёплым краям пластины (рис. 10) путём применения локальных операторов диффузии в1 (i, j) (3). Эволюция этой модели имитирует изменение распределения температуры пластины.
Простой глобальный оператор ©а выполняет основную функцию формирования устойчивого образа (узора) путём применения к Qv (t) локального оператора тотали-стического типа [8] $tot(i, j), у которого контекст C((i,j)v) состоит из двух частей:
Рис. 10. Начальное состояние Qv (0)
C'((i,j )v) = {(Vi+a,3+b, (i + a,j + b)v) : a,b = -r, -r + 1,..., 0,1,...,r + 1} \ {(i,j)}, C ((i,j )u) = {(Mi+a.j+b, (i + a, j + b)u) : a, b = -r, -r + 1,..., 0,1,...,r + 1}, r < |X | -радиус соседства 0tot(i, j)[{(vo, (i, j)v)},C'((i, j)v) U C((i, j)„)] = {(v0, (i, j)v)}, где v0 вычисляется функцией переходов как
r+1
1, если ^ Vi+a,j+b • w(i + a, j + b) ^ 0,
/ _ J a,b=-r
V0 = ^ r+1
0, если ^ Vi+a,j+b • w(i + a,j + b) < 0.
a,b=-r
Здесь w(i + a, j + b) = -0,2ui+j • (-1)(a+b).
В начальном массиве (0) значения u =1 равномерно распределены с плотностью 0,5.
На рис. 11 показаны три глобальных состояния Qv (t). Моделирование выполнялось на двух ядрах процессора Intel Core i-4770K. Параллельная работа глобальных операторов обеспечивалась с помощью библиотеки OpenMP. Устойчивое состояние на клеточном массиве размером 2 х 400 достигалось за 2,5 мин.
Рис. 11. Три глобальных состояния О(¿), полученных при моделировании процесса формировния узора на остывающей металлической пластине
Заключение
Клеточно-автоматные модели пространственной динамики — это одно из направлений в математическом и компьютерном моделировании, которое согласуется с современными тенденциями параллелизации архитектуры вычислительных систем. Основные свойства КА-моделей — локальность взаимодействий вычисляющих операторов и отказ от вещественных переменных — избавляют от острых проблем вычислительной математики: проблемы распараллеливания, проблемы потери точности и, наконец, проблемы сложности программирования. Более того, вычислительные свойства КА-моделей — способность моделировать нелинейные явления (самоорганизация, автоволны, образование фрактальных структур) — согласуются с насущной потребностью развивать методы моделирования химических и биологических процессов. Автор сочтёт, что его цель достигнута, если эта работа убедит кого-то из специалистов в необходимости исследовать и развивать КА-моделирование.
ЛИТЕРАТУРА
1. Фон Нейман Дж. Теория самовоспроизводящихся автоматов. М.: Мир, 1971. 384 с.
2. Wolfram S. Computation theory of cellular automata // Commun. Math. Phys. 1984. V. 96. P. 15-57.
3. Sandefur J. T. Discrete Dynamical Systems: Theory and Applications. Oxford: University Press, 1990. 445 p.
4. Donienunzio A., FormentiE., and Kurka P. Cellular automata dynamical systems // Handbook of Natural Computing. V. 1. / eds. G. Rozenberg, T. Back, J.N. Kok. Berlin: Springer, 2012. P. 24-76.
5. Wolfram S. Statistical mechanics of cellular automata // Rev. Mod. Phys. 1983. V. 35. P. 601-643.
6. Toffolli T. and Margolus N. Cellular Automata Machines: A New Environment for Modeling. USA: MIT-Press, 1987. 260 p.
7. Frish U., Hasslacher B., and Pomeau Y. Lattice-gas automata for Navier — Stokes equation // Phys. Rev. Let. 1986. V. 56. P. 1505-1508.
8. Wolfram S. A New Kind of Science. Champaign, Ill. USA: Wolfram Media Inc., 2002. 1197 p.
9. Weimar J. Cellular automata for reaction-diffusion systems // Parallel Computing. 1997. V. 23. No. 11. P. 1699-1715.
10. Ванаг В. К. Диссипативные структуры в реакционно-диффузионных системах. Ижевск: Ин-т компьют. исслед., 2008. 300 c.
11. Madore B. and Freedman W. Computer simulation of the Belousov — Zhabotinski reaction // Science. 1983. V.222. P. 615-618.
12. Park J. K., Steiglitz K., and Thurston W. P. Soliton-like behavior in automata // Physica D. 1986. V. 19. P. 423-432.
13. Achasova S., Bandman O, Markova V., and Piskunov S. Parallel Substitution Algorithm. Theory and Application. Singapore: World Scientific, 1994. 198 p.
14. Бандман О. Л. Клеточно-автоматные модели пространственной динамики // Системная информатика. 2006. Вып. 10. C. 59-111.
15. Ачасова С. М., Бандман О. Л. Корректность параллельных вычислительных процессов. Новосибирск: Наука, 1990. 256 с.
16. Bandman O. Parallel simulation of asynchronous cellular automata evolution // ACRI-2006. LNCS. 2006. V. 4173. P. 41-47.
17. Sharifulina A. and Elokhin V. Simulation of heterogeneous catalytic reaction by asynchronous cellular automata on multicomputer // PaCT-2011. LNCS. 2011. V.6873. P. 210-216.
18. Шарифулина А. Е. Параллельная реализация каталитической реакции (СО + О2 ^ СО2) с помощью асинхронного клеточного автомата // Труды ПаВТ-2012. Челябинск: Изд-во ЮрГУ, 2012. С. 325-335.
19. Bandman O. Cellular automata composition techniques for spatial dynamics simulation // Simulating Complex Systems by Cellular Automata / eds. A. G. Hoekstra, J. Kroc, P.M.A. Sloot. Berlin: Springer, 2010. P. 81-116.
20. Bandman O. Using multi core computers for implementing cellular automata systems // PaCT-2011. LNCS. 2011. V.6873. P. 45-58.
21. Медведев Ю.Г. Многочастичная клеточно-автоматная модель потока жидкости FHP-MP // Вестник Томского государственного университета. Управление, вычислительная техника и информатика. 2009. №1(6). С. 33-40.
22. Медведев Ю. Г. Моделирование трехмерных потоков клеточными автоматами // Вестник Томского государственного университета. 2002. Приложение. №1 (II). С. 236-240.
23. Медведев Ю. Г. Многочастичные клеточно-автоматные модели диффузионного процесса // Новые инф. технологии в исслед. сложных структур. Тез. докл. Томск: Изд-во НТЛ, 2010. С. 21-22.
24. Bandman O. Cellular automata diffusion models for multicomputer implementation // Bull. Nov. Comp. Center. Ser. Comp. Sci. 2014. V. 36. P. 21-31.
25. Медведев Ю. Г. Применение клеточно-автоматной модели потока вязкой жидкости в исследовании трёхмерных пористых сред // Автометрия. 2006. Т. 42. №3. С. 21-31.
26. Медведев Ю. Г. Программный комплекс клеточно-автоматного моделирования газопорошковых потоков // Труды ПаВТ-2012. Челябинск: Изд-во ЮрГУ, 2012. С. 732.
27. Wolf-Gladrow D. Lattice-Gas Cellular Automata and Lattice Boltzmann Models. Berlin: Springer, 2000. 310 p.
28. Четверушкин Б. Н. Кинетически-согласованные схемы в газовой динамике: новая модель вязкого газа, алгоритмы, параллельная реализация, приложения. М.: Изд-во МГУ, 1999. 232 с.
29. Numerical Algorithms and Round-off Errors. http://www.sml.ee.upatras.gr/ uploadedfiles/roundofferror.pdf.
30. Kireeva A. Two-layer CA model for simulating catalytic reaction at dynamically varying temperature // ACRI-2014. LNCS. 2014. V.8751. P. 166-175.
31. Бандман О. Л., Киреева А. Е. Стохастическое клеточно-автоматное моделирование колебаний и автоволн в реакционно-диффузионных системах // Сиб. ЖВМ. 2015. Т. 18. №3. С. 255-274.
32. Bandman O. Cellular automata model of fluid permeation through porous material // PaCT-2013. LNCS. 2013. V. 7979. P. 295-316.
33. Bandman O. Parallel composition of asynchronous cellular automata simulating reaction diffusion processes // ACRI-2010. LNCS. 2010. V.6350. P. 395-398.
34. Витвицкий А. Клеточные автоматы с динамической структурой для моделирования роста биологических тканей // Сиб. ЖВМ. 2014. Т. 17. №4. С. 315-325.
REFERENCES
1. Von Neumann J. Teoriya samovosproizvodyashchikhsya avtomatov [The Theory of Self-Reproducing Automata]. Moscow, Mir Publ., 1971. 384p. (in Russian)
2. Wolfram S. Computation theory of cellular automata. Commun. Math. Phys., 1984, vol.96, pp. 15-57.
3. Sandefur J. T. Discrete Dynamical Systems: Theory and Applications. Oxford: University Press, 1990. 445 p.
4. Donienunzio A., Formenti E., and Kurka P. Cellular automata dynamical systems. Handbook of Natural Computing. V. 1. / eds. G. Rozenberg, T. Back, J.N. Kok. Berlin: Springer, 2012, pp. 24-76.
5. Wolfram S. Statistical mechanics of cellular automata. Rev. Mod. Phys., 1983, vol.35, pp.601-643.
6. Toffolli T. and Margolus N. Cellular Automata Machines: A New Environment for Modeling. USA: MIT-Press, 1987. 260 p.
7. Frish U., Hasslacher B., and Pomeau Y. Lattice-gas automata for Navier — Stokes equation. Phys. Rev. Let., 1986, vol.56, pp. 1505-1508.
8. Wolfram S. A New Kind of Science. Champaign, Ill. USA: Wolfram Media Inc., 2002. 1197 p.
9. Weimar J. Cellular automata for reaction-diffusion systems. Parallel Computing, 1997, vol. 23, no. 11, pp.1699-1715.
10. Vanag V. K. Dissipativnye struktury v reaktsionno-diffuzionnykh sistemakh [Dissipative Structures in Reaction-Diffusion Systems]. Izhevsk, Institute of Computer Science Publ., 2008. 300 p. (in Russian)
11. Madore B. and Freedman W. Computer simulation of the Belousov — Zhabotinski reaction. Science, 1983, vol.222, pp. 615-618.
12. ParkJ.K., Steiglitz K., and Thurston W.P. Soliton-like behavior in automata. Physica D., 1986, vol. 19, pp. 423-432.
13. Achasova S., Bandman O, Markova V., and Piskunov S. Parallel Substitution Algorithm. Theory and Application. Singapore: World Scientific, 1994. 198 p.
14. Bandman O. L. Kletochno-avtomatnye modeli prostranstvennoy dinamiki [Cellular automata models of spatial dynamics]. System Informatics, 2006, no. 10, pp. 59-111. (in Russian)
15. Achasova S. M., Bandman O. L. Korrektnost' parallel'nykh vychislitel'nykh protsessov [The Correctness of Parallel Computing Processes]. Novosibirsk, Nauka Publ., 1990. 256 p. (in Russian)
16. Bandman O. Parallel simulation of asynchronous cellular automata evolution. ACRI-2006, LNCS, 2006, vol.4173, pp. 41-47.
17. Sharifulina A. and Elokhin V. Simulation of heterogeneous catalytic reaction by asynchronous cellular automata on multicomputer. PaCT-2011, LNCS, 2011, vol.6873, pp.210-216.
18. Sharifulina A. E. Parallel'naya realizatsiya kataliticheskoy reaktsii (SO + O2 ^ SO2) s pomoshch'yu asinkhronnogo kletochnogo avtomata [Parallel implementation of the catalytic reaction (SO + O2 ^ SO2) using an asynchronous cellular automaton]. Proc. PaVT-2012, Chelyabinsk, SUSU Publ., 2012, pp. 325-335. (in Russian)
19. Bandman O. Cellular automata composition techniques for spatial dynamics simulation. Simulating Complex Systems by Cellular Automata / eds. A. G. Hoekstra, J. Kroc, P.M. A. Sloot. Berlin, Springer, 2010, pp. 81-116.
20. Bandman O. Using multi core computers for implementing cellular automata systems. PaCT-2011, LNCS, 2011, vol.6873, pp. 45-58.
21. Medvedev Yu. G. Mnogochastichnaya kletochno-avtomatnaya model' potoka zhidkosti FHP-MP [An extension of the cellular-automaton FHP-I flow model to the FHP-MP multi-particle model]. Vestnik Tomskogo gosudarstvennogo universiteta. Upravlenie, vychislitel'naya tekhnika i informatika, 2009, no. 1(6), pp. 33-40. (in Russian)
22. Medvedev Yu. G. Modelirovanie trekhmernykh potokov kletochnymi avtomatami [The 3D fluid flow simulation]. Vestnik Tomskogo gosudarstvennogo universiteta. Prilozhenie, 2002, no. 1(II), pp. 236-240. (in Russian)
23. Medvedev Yu. G. Mnogochastichnye kletochno-avtomatnye modeli diffuzionnogo protsessa [Many-cellular automata models of diffusion process]. Novye inf. tekhnologii v issled. slozhnykh struktur. Tez. dokl. Tomsk, NTL Publ., 2010, pp. 21-22. (in Russian)
24. Bandman O. Cellular automata diffusion models for multicomputer implementation. Bull. Nov. Comp. Center. Ser. Comp. Sci., 2014, vol.36, pp.21-31.
25. Medvedev Yu. G. Primenenie kletochno-avtomatnoy modeli potoka vyazkoy zhidkosti v issledovanii trekhmernykh poristykh sred [Applying the cellular automaton model of a viscous fluid flow in investigating 3D porous media. Avtometriya, 2006, vol.42, no.3, pp. 21-31. (in Russian)
26. Medvedev Yu. G. Programmnyy kompleks kletochno-avtomatnogo modelirovaniya gazoporoshkovykh potokov [The software package for cellular automata modeling of gas-powder flow]. Proc. PaVT-2012. Chelyabinsk, SUSU Publ., 2012, p. 732. (in Russian)
27. Wolf-Gladrow D. Lattice-Gas Cellular Automata and Lattice Boltzmann Models. Berlin, Springer, 2000. 310 p.
28. Chetverushkin B. N. Kineticheski-soglasovannye skhemy v gazovoy dinamike: novaya model' vyazkogo gaza, algoritmy, parallel'naya realizatsiya, prilozheniya [Kinetically consistent schemes in gas dynamics: the new model of a viscous gas, algorithms, parallel implementation, application]. Moscow, MSU Publ., 1999. 232 p. (in Russian)
29. Numerical Algorithms and Round-off Errors. http://www.sml.ee.upatras.gr/ uploadedfiles/roundofferror.pdf.
30. Kireeva A. Two-layer CA model for simulating catalytic reaction at dynamically varying temperature. ACRI-2014, LNCS, 2014, vol.8751, pp. 166-175.
31. Bandman O. L. and Kireeva A. E. Stochastic cellular automata simulation of oscillations and autowaves in reaction-diffusion systems. Numerical Analysis and Applications, 2015, vol.8, no. 3, pp. 208-222.
32. Bandman O. Cellular automata model of fluid permeation through porous material. PaCT-2013, LNCS, 2013, vol.7979, pp. 295-316.
33. Bandman O. Parallel composition of asynchronous cellular automata simulating reaction diffusion processes. ACRI-2010, LNCS, 2010, vol. 6350, pp. 395-398.
34. Vitvitsky A. A. Cellular automata with dynamic structure to simulate the growth of biological tissues. Numerical Analysis and Applications, 2014, vol. 7, no. 4, pp. 263-273.
