Анализ параллельных алгоритмов при моделировании процесса распространения лесных пожаров

Выводы

В ходе данной работы проведено исследование процесса формирования страховой документации

000 «Надежда». Проанализированы различные программные средства разработки и предложены альтернативные варианты для создания программного комплекса решения автоматизации задачи.

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

Литература

1 Калянов, Г.Н. CASE: структурный системный анализ (автоматизация и применение) / Г.Н. Калянов. -

М.: ЛОРИ, 1996. - 180 с.

2 Розенберг, Д. Применение объектного моделирования с использованием UML / Д. Розенберг, К. Скотт. -

М.: ДМК, 2002. - 160 с.

3 Методология структурного анализа и проектирования / http://www.marathon.ru.

УДК 519.63 Г.А. Доррер, М.С. Вдовенко

АНАЛИЗ ПАРАЛЛЕЛЬНЫХ АЛГОРИТМОВ ПРИ МОДЕЛИРОВАНИИ ПРОЦЕССА РАСПРОСТРАНЕНИЯ ЛЕСНЫХ ПОЖАРОВ

Рассмотрена эффективность параллельных алгоритмов при моделировании динамики лесных пожаров. Проанализированы различные способы декомпозиции сеточной области. Приводится оценка ускорения вычислений по данным эксперимента, выполненного на кластере ИВМ СО РАН.

Ключевые слова: лесной пожар, эффективность параллельных алгоритмов, моделирование лесных пожаров.

G.A. Dorrer, M.S. Vdovenko ANALYSIS OF PARALLEL ALGORITHMS FOR THE FOREST FIRE PROPAGATION PROCESS MODELLING

The efficiency of the parallel algorithms in forest fire dynamics modelling is considered. Various ways of decomposition of the grid area are analyzed. The estimation of calculation acceleration on the basis of experiment data made on the Institute of Computing Modeling cluster is given.

Key words: forest fire, parallel algorithms efficiency, forest fires modelling.

Введение

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

Проблема математического моделирования процессов горения при лесных пожарах изучается уже в течение многих лет. Обзор результатов, полученных в этой области, приведен в работе [4]. Большой вклад в решение данной проблемы внесли Н.П. Курбатский, Э.Н. Валендик, М.А. Софронов, А.М. Гришин, Г.Н. Коровин, F.A. Albini, G.M. Byram, R.C. Rothermel, M.G. Cruz, M.E. Alexander и другие ученые.

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

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

1. Анализ задачи и выявление ее потенциального параллелизма

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

Постановка задачи. В качестве базовой модели процесса распространения лесного пожара нами взята модель, основанная на вычислении теплового баланса в лесных горючих материалах [4]. Массив горючего представляет собой в общем случае n параллельных однородных слоев горючего, расположенных один над другим. Произвольный /-й слой занимает по вертикали область Zi с координатами от zjH до zjK

и содержит запас горючего материала со у, z, t ^ кг/м3.

Вертикальная координата середины i -го слоя обозначается zicp, а его толщина 8 . При этом

z--" = %íh + 'к ^2, 8 — :iK — i =

юр *ЧН ' гК Л* ’ . ' ¡К ' Ш ’

Свойства горючего в пределах каждого слоя не зависят от г .

Горючий материал в окрестностях точки С в некоторый момент времени может находиться в одном из трех состояний, описываемых функцией Б 4, у, г, г^

[о, если в точке С в момент г имеется ненулевой запас горючего (т.е.

СО у, 2^ _> )), но горения не происходит;

= 1, если со [,у,г,г^> ) и происходит горение;

2, если со {, у, г, г = ), т.е. горение невозможно.

Области, соответствующие состояниям Л'= ), £ = , £ = I, обозначаются соответственно , Г2 ,0 . Проекции областей 0,0,0 на горизонтальную плоскость £> будем обозначать соответственно Д,, Вх, П2, причем £>0 и \ и \ - Э.

Уравнение нагрева для горючего в окрестности точки С = %,у,г ^ которое в момент времени г находится в состоянии £ 4с,у,г,г = ):

д 4,У,*.

д

* и л •

ъ

+ :^,у\^,у^- 10]-4,у2

= ЕЯФ \>У~ '13МУ1 +

-]К

}Н

|ЯУ 4;,у,г,^г , при начальном условии Н] 4[,у,0^= Ч0] 4(,,у ^

У _е оуг= ,•••,«.Здесь Н 4с, у,г, г - значение энтальпии в точке 4с, У,г ^ в момент времени Г;

Ф, Ф - соответственно, энергия, образующаяся при горении и поступающая от внешних источников, вт/м3; £ \,у- функция влияния пламени из точки 4сг,уг^ наточку 4,у ^(функция Грина) из

а8

/-го слоя на 7-й слой; к у, г - — , где а - коэффициент теплоотдачи, Вт/м2град; 8 - удельная по" Р

верхность слоя, м1; с' - приведенная теплоемкость влажного материала, Дж/кг град.

Условие воспламенения горючего в у-м слое, т.е. перехода горючего в состояние Sj 4ь,у,г^ ,

имеет вид Н. (¿,у^* ^ Ч* 4,У ^ где 7* = * 41, у ^ - время воспламенения горючего в у -м слое в точке (,^е >0 .; Я*-энтальпия начала газификации.

Уравнение расходования горючего:

д со 4;, у, t Л Г- г при со 4, у, I > О

Ш\_а_]4(,у^ [ 0 при со 4,У^^= О

с начальным условием со (¡,у,(* = о ¡4с,у ^ 4(,у ^

У

Здесь со - активный запас горючего материала, в у -м слое, кг/м3; г . - относительная скорость сгорания ] -го слоя, 1/с.

Уравнение тепловыделения в _/ -м слое:

, дсо 4,У^ . -

Ф, \,у,г = к - \,у ^^---------------г, %,у^ \г 7=

1 " > ^ д

где /г. - теплота сгорания горючего, Дж/кг.

Условие погасания (перехода в состояние Б]4(,у^^= 1)\ со ^у^г*= ), 4,У^ )2],

7= ,-,п.

Модель теплопередачи из зоны горения к горючим материалам, задаваемая функцией влияния (функцией Грина) £ 4(- \,у- подробно рассмотрена в [4].

Система уравнений (1)-(3) представляет собой модель процесса распространения горения по неод-

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

Метод решения. Алгоритм построения горящей кромки основан на численном решении уравнений (1)—(3), описывающих распространение процесса горения. В каждом из слоев вводится прямоугольная сетка с шагами по координатам х и у соответственно А и А , области Ц,г, £>1г и /)2г заменяются соответствующими сеточными областями ВА , Ба /= ,2. Перейдя к дискретному времени )Д,2,... с ша-

гом А и заменив в (1)-(2) частные производные по времени разностными отношениями, а интеграл по области £>1г суммой, получим расчетные уравнения, которые использовались при дальнейших исследованиях.

Функции влияния %_4(- \,у— зависят от характеристик леса, природных и погодных факторов и вычислялись по специальным подпрограммам [4].

1

г

В качестве исходных данных использовались:

область моделирования в виде двух множеств узлов ИА - = >1= >2;

начальный и конечный моменты времени /0 и временной шаг А ;

участки с одинаковыми характеристиками горючих материалов О ,к = = ЭА, / = ,2;

к

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

И()йС Н1Л¿к С ЛиС ;ф С  СлО/= ’2' * = 1= оА + Ь параметры функции, описывающей тепловое воздействие локального пламени р 1 аш С , <5 С ,

V- С 141 С;

скорости сгорания гк14,^ для всех участков О в каждый момент времени ? = 0,?0+ ^

начальное состояние системы - области БА , ^ , БА ^ ИА , ^

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

чить также путем анализа аэрокосмических снимков пожара, полученных в последовательные моменты времени [4]. На основе этих снимков могут быть вычислены необходимые для расчетов значения функций £ С- \,У~ ;1^ и скоростей. Поскольку конфигурация кромки реальных пожаров часто является достаточно сложной, то для упрощения расчетов целесообразно применять сглаживание границы контура на снимке пожара, например, аппроксимацию его эллипсом. Это позволяет оперативно получать необходимые для моделирования исходные данные.

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

(7 + ) -м временном шаге происходит с использованием некоторого количества точек на ?-м шаге. Числен-

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

Рис. 1. Последовательность вычислений при моделировании распространения пожара. Позиции:

P1 - данные для инициализации функций библиотеки MPI; P2 - данные о метеоусловиях и характеристиках горючего; P3 - отображение рассчитанной горящей кромки на карте местности.

Составные переходы: tl - операции по вводу данных; 12 - основной расчетный цикл

Таким образом, исходную задачу можно разбить на р подзадач для областей, D\ . , DA D\ it

l = ,2,k = пересекающихся только по границе разбиений, независимых друг от друга на каждом

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

Рис. 2. Модель основного расчетного цикла. Позиции: Р21 - данные о метеоусловиях, характеристиках горючего, горящих участках, начальном и конечном моментах времени горения, Р22 - данные о локальных функциях влияния, Р23 - данные об энтальпии горючего, Р24 - данные об энтальпии с учетом влияния соседних участков, Р25 - данные о состоянии горящих участков. Составные переходы: т - вычисление текущего временного шага (выполняется при условии t < final), т [-определение границы локальных массивов, т '• - вычисление значения энтальпии горючего, т  - обмен данными с соседними процессорами, т • - определение новых воспламенившихся и погасших точек, г > - вывод в файл данных о конфигурации пожара

На каждом шаге по времени соседние процессоры осуществляют обмен граничными значениями при помощи функций неблокирующих пересылок и приема данных MPIJsend и MPI_Irecv. Схема обмена с соседними процессорными узлами процессора показана на рисунке 3. По окончании расчетов каждый процесс обрабатывает свой массив вычисленных значений, записывая его в отдельный файл.

Рис. 3. Схема обмена процесса " prOC" с соседними процессами: операции обмена up - In,left- ight,upleft- Inright,upright - Inleft выполняются одновременно

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

ществляется с использованием процедур библиотеки MPI [7- 9]. Схема операций обмена данными (операция т [ на рисунке 2) показана на рисунке 4.

Рис. 4. Модель обмена данными с соседними процессами. Позиции: Р413 P42, Р43, Р44, Р45 - данные для обмена с соседними процессорами (согласно рисунку 2). Составные переходы: т Л - обмен up - Int, т 2 - обмен left - ight, т з обмен upleft — Inright, т 4 - обмен upright — Inleft

Более детально схема одной из операций обмена показана на рисунке 5.

Рис. 5. Модель операций межпроцессорного обмена на примере операции т 1 (рисунок 4), остальные операции выполняются аналогично: Р411, Р412 - данные, приготовленные для обмена; Р413, Р414 -данные полученные при обмене; т _и,т 12 - операции пересылки данных; ку,к2,к3,к4,къ - маркеры готовности к обмену данными; ^ - начало операции обмена; ^ - конец операции обмена.

Потенциальное ускорение алгоритма. Будем оценивать время работы параллельной программы исходя из соотношения £ = , где 5” -ускорение, Тх - время вычислений на одном процессоре, Тр

- время вычислений на р процессорах.

Для получения реалистических оценок помимо классической формулы Амдала [6] будем учитывать также время, затрачиваемое программой на обмены между процессами. Как следует из принятой нами схемы распределения данных, на каждом временном шаге требуется обмен границами. Время пересылок для различных способов декомпозиции можно приблизительно выразить через количество пересылаемых данных [6]:

V comm — 2 • N3 • т, при одномерном разбиении,

V2Dcomm= N3 • т, при двумерном разбиении,

VP

где iV3- размерность задачи, р- количество вычислительных узлов, т - время пересылки одного числа. Алгоритм и его программная реализация являются масштабируемыми, если ускорение и производительность зависят линейно от количества используемых процессоров [6, 9]: S = 1{р). На практике алгоритмы, для которых S = Э(р1(In р)) также считаются масштабируемыми.

2. Проведение вычислительного эксперимента и анализ результатов

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

Таблица 1

Значение ускорения, показанного на кластере ИВМ СО РАН

Количество процессоров 1 2 4 8 16

Разбиение расчетной области на квадраты

Время выполнения основного цикла T [с] 700,868934 359,8456 188,142001 81,975424 53,29919

Ускорение, S 1,000000 1,947693 3,725213 8,549745 13,14971

Эффективность, Е = ? / р 1,000000 0,973847 0,931303 1,068718 0,821857

Разбиение расчетной области на полосы

Время выполнения основного цикла T [с] 700,868934 359,8456 207,0777 105,2597 87,26806

Ускорение, S 1,000000 1,947693 3,384569 6,658474 8,031219

Эффективность, Е = ? / р 1,000000 0,973847 0,846142 0,832309 0,501951

Расчеты производились по формулам (4)-(6) на кластерной системе МВС 1000 ИВМ СО РАН (г. Красноярск) на тестовой сетке 400x400 узлов при использовании от 1 до 16 процессоров.

Тестовая область О представляла квадрат О = [1,6л] х [1,6л]. Возвышение поверхности задавалось выражением г(х) = + Ю8(.х7150). Процесс распространения пожара инициировался из узла с координатами (0.5, 0.5). При запуске параллельных программ измерялось время их работы в секундах. На основе данных о времени работы программ (табл. 1) вычислялись другие характеристики параллельных программ, такие как ускорение и эффективность распараллеливания. В вычислительных экспериментах было сделано по 100 шагов по времени.

Значения ускорения, показанного на кластере ИВМ СО РАН, приведены в таблице 1. При применении двумерной декомпозиции и 8 процессоров значение эффективности больше единицы, что объясняется использованием в программе «динамических» массивов с подстраиваемыми под выделенное число процессоров размерами. Таким образом, необходимы меньшие временные затраты на выборку обрабатываемых данных из оперативной памяти и передачу их через КЭШ-память. В случае использования 8 процессоров при данной размерности сетки весь массив помещается в КЭШе, что и определяет более быстрое выполнение вычислений за счет отсутствия необходимости обмена между оперативной памятью и КЭШем.

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

15

Количестко процессоров

—*— Ускорение при двумерном разбиении -    - Ускорение при линейном разбиении

Рис. 6. Зависимость ускорения от количества доступных процессоров

Как видно из рисунка 6, результаты вычислительного эксперимента показали наличие хорошего ускорения при решении данной задачи.

Также был проведен вычислительный эксперимент, выявивший зависимость ускорения от роста размерности задачи. Были рассмотрены случаи мелкой (100x100=10000 узлов) и крупной (800x800=640000 узлов) сетки. Тестовая область О представляла, как и ранее, квадрат: О = [1,6л] * [1,6л], а возвышение поверхности задавалось формулой о(х)= + ;оз(д:/150). Характеристики горючего также не изменялись. В вычислительных экспериментах было сделано также по 100 шагов по времени. Декомпозиция расчетной области - двумерная. Полученные значения ускорения, показанного на кластере ИВМ СО РАН (г. Красноярск), приведены в таблице 2.

Таблица 2

Значение ускорения, показанного на кластере ИВМ СО РАН

Размерность задачи 100x100 200x200 400x400 800x800

Время выполнения основного цикла на 1 процессоре Т, с 0,464515 7,300380 122,961800 1532,029398

Время выполнения основного цикла на 4 процессорах Т, с 0,185772 2,149292 33,887110 394,943400

Ускорение £ 2,5004621 3,3966438 3,6285715 3,879111

Эффективность Е = ? / р 0,625116 0,849161 0,907143 0,969778

Выполненные расчеты показали, что с увеличением размерности задачи при использовании четырех процессоров наблюдается увеличение ускорения вычислений, по крайней мере, для задач размерностью не более 800x800. Из таблицы 2 видно, что при учете размера кэш-памяти можно добиться эффективности больше 1. Значения ускорения вычислений с увеличением размерности задачи при использовании двумерной декомпозиции расчетной области возрастают логарифмически. Эту зависимость можно аппроксимировать выражением у = 0,3151Ьп(х)-0,2059, с достоверностью аппроксимации Я2 - ),8821.

4

2 т--------1---------1----------1---------1---------1---------1---------1

О 100000 200000 300000 400000 500000 600000 700000

Ра ¡мерность ¡адачи

Рис. 7. Зависимость ускорения от размерности задачи

На рисунке 7 показана зависимость ускорения вычислений от размерности задачи. При увеличении размерности задачи в 4 раза ускорение вычислений при использовании 4 процессоров возрастает. Средний уровень эффективности распараллеливания составляет около 83%, что связано с неизбежными затратами времени на организацию межпроцессорных обменов и записи результатов в файл.

Заключение

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

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

Получены расчетные значения ускорений, позволяющие оценить масштабируемость алгоритма и его программной реализации. Эти результаты показывают, что алгоритм обладает значительным объемом потенциального параллелизма и хорошей, с точки зрения распараллеливания, структурой, что позволяет надеяться на получение ускорений, близких к линейным в зависимости от количества используемых процессоров для кластерных вычислительных систем. Значения ускорений, полученные в ходе вычислительного эксперимента, хорошо согласуются с теоретическими оценками.

Литература

1. Применение информационной системы дистанционного мониторинга «ИСДМ-Рослесхоз» для определения пожарной опасности в лесах Российской Федерации: учеб. пособие / ФГУ Авиалесоохрана. -Пушкино (МО), 2007. - 82 с.

2. Софронов, М.А. Огонь в лесу / М.А. Софронов, А.Д. Вакуров. - Новосибирск: Наука, 1981. - 128 с.

3. Гришин, A.M. О математическом моделировании природных пожаров и катастроф / А. М. Гришин // Вестн. Томского гос. ун-та: Математика и механика. - Томск: Изд-во ТГУ, 2008. - №2(3). - С. 105-114.

4. Доррер, Г.А. Динамика лесных пожаров / Г.А. Доррер. - Новосибирск: Изд-во СО РАН, 2008. - 404 с.

5. Weber, R.O. Modeling fire spread through fuel beds / R.O. Weber// Prog. Everg. Combust. Sci. - 1990. -V. 17. - P. 65-82.

6. Воеводин, В.В. Параллельные вычисления / В.В. Воеводин, Вл.В. Воеводин. - СПб.: БХВ-Петербург, 2002. - 608 с.

7. Message-Passing Interface Forum, MPI-2: Extensions to the Message-Passing Interface, 1997. [http://www.unix.mcs.anl.gov/mpi/], 13.03.2007.

8. Rs/6000 SP: Practical MPI Programming. [www.redbooks.ibm.com], 11.08.2008.

9. Корнеев, В.Д. Параллельное программирование в MPI. - 2-е изд. / В.Д. Корнеев. - Новосибирск: Изд-во ИВМиМГ СО РАН, 2002. - 215 с.