Об эффективности параллельных алгоритмов при моделировании динамики лесных пожаров Текст научной статьи по специальности «Математика»

УДК 519.63

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

ОБ ЭФФЕКТИВНОСТИ ПАРАЛЛЕЛЬНЫХ АЛГОРИТМОВ ПРИ МОДЕЛИРОВАНИИ ДИНАМИКИ ЛЕСНЫХ ПОЖАРОВ

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

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

Известно негативное влияние природных пожаров на состояние окружающей среды. Так, одной из причин эффекта глобального потепления климата считают выбросы СО2 и других многочисленных газов, возникающих в результате пожаров [1-3]. Помимо влияния на климат Земли, лесные пожары наносят ущерб лесному фонду, а также способствуют возникновению таких катастроф, как большая задымленность больших городских территорий (например, в Москве в 2002 г.) или массовые городские пожары (например, в Лос-Аламосе (США) в 2002 г.).

Содержательный обзор исследований по проблеме моделирования процессов горения при лесных пожарах дан в работе [4]. Большой вклад в разрешение этой проблемы внесли Н. П. Курбатский, Э. Н. Валендик, М. А. Соф-ронов, А. М. Гришин, Г. Н. Коровин, Р. Ротхермел, M. E. Ллександер и другие ученые.

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

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

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

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

Постановка задачи. В качестве базовой модели процесса распространения лесного пожара взята модель, основанная на вычислении теплового баланса в лесных горючих материалах [3]. Массив горючего в общем случае представляет собой п параллельных однородных слоев горючего, расположенных один над другим. Произвольный і-й слой занимает по вертикали область 2І с координатами от ііН до ііК и содержит запас горючего материала ю(х, у, і, ґ) кг/ м3. Вертикальная координата середины і-го слоя обозначается і

іСр, а его толщина 5і. При этом

гіср

іН

ЧК

)

2

5,- - г1К -гш, I - 1,...,п.

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

Горючий материал в окрестностях точки С в некоторый момент времени может находиться в одном из трех состояний, описываемых функцией £ (х ^ z, I):

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

£ (х, у, X, ?)-• т. е. Ю (х, у, X, t )> 0, но горения не происходит,

1, если Ю (х, у, X, t) > 0 и горение происходит,

2, если Ю (х, у, X, t) — 0 и горение невозможно.

Области, соответствующие состояниям £ — 0, £ — 1, £ — 2, обозначаются соответственно О0 , О2 . Про-

екции областей , О, О2 на горизонтальную плоскость С будем обозначать соответственно £>0 , С , £>2, причем £>0 и С и С2 — С .

Уравнение нагрева для горючего в окрестности точки С — (х, у, 2), которое в момент времени t находится в состоянии £ (х, у, 2, t) — 0, имеет вид

дн] ( х, у, t) — дt

п

—X Я °г (хъ у^ г‘) (х - х^ у- у) ^1^1 +

,—1 С

+Фе (х, у, і) - к} (х, у) [Н} (х, у, і) - Н0} (х, у)], (1)

где Ну ( х, у, і ) = -

1

гуК

| Ну (х, у, г, і) ёг при на-

у -2Н г-

чальном условии Ну (х, у,0) = Н о у (х, у), (х, у) є Д

у - “0уЛл’У)’ \Л’У)^^0у,

/ = 1,..., п, здесь Иу (х, у, 2, t) - значение энтальпии в точке (х, у, г) в момент времени t; Ф , Фе - энергия, образующаяся при горении и поступающая от внешних источников соответственно, Вт/м3; (х - хl,у -У1) -

функция влияния пламени из точки (х1, У1) на точку (х,у) (функция Грина) из /-го слоя на у-й слой;

к (х, у, 2 ) = —, здесь а - коэффициент теплоотдачи,

Рс ,

Вт/(м2-град), 5 - удельная поверхность слоя, м-1, с - приведенная теплоемкость влажного материала, Дж/(кг • град).

Условие воспламенения горючего в у-м слое, т. е. перехода горючего в состояние (х, у, 2, t) = 1, имеет вид

/ * \ * / Ч У * у- /

Иу (х, у, tj ) > Иу (х, у), где tj = tj (х, у) - время воспламенения горючего в у-м слое в точке (х у )е Д0 у; и * - энтальпия начала газификации.

Уравнение расходования горючего будет следующим:

д_

ді

®у (х, у, і)

(2)

Ю0 у (х у )

Г-Гу при Юу (х, у, t) > 0.

10 при юу (х, у, t) = 0

с начальным условием юу (х, у,^) = Ю0у (х,у),

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

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

Зюу (х, у, t)

Ф у (х, У, і ) = -Ау (х, у)

ді

до;

Д

Д1і

Д

Д2і

моменты 1 = 0,1, 2,... с шагом Дt. Заменив в (1), (2) частную производную по времени разностными отношениями, а интеграл по области Дц суммой, получим следующее уравнение нагрева:

Ие (/,у,1 +1) = Ие (/,у,1) + Д1ДхДу х 2

хХ X Ф(т=п1 )х 4I (т -п -у') +

5=1 (т,п )еД]5

+Фе/ (^ у\ )Дt -Дк (^ у)х

(4)

х[И1 (^, Л <:)- И01 (^,})]

с начальными условиями И1 (/,у,0) = И01 (/, у) ,

(/, у ) е ДД01, 1 = Ъ2

Условие воспламенения горючего в I -м слое при

t

И (/,у, /)>И* (/,у) , узел (/,/) исключается из Дд0/ (/*) и присоединяется к

ДМ1 (t

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

Щ О, j, t + !) =

[ю/ (/,у,^-гД1 при Ю/ (/,у,t)>0, (5)

Iю/ (^ j,t) при щ (/, j,t)< 0

с начальным условием Щ ^ j, t ) = щ0/ (i, у), (/, у) е Ддц , / = 1,2, а уравнение тепловыделения - вид

ф/ ^, у\t) = -И/ О, у )х [щ/ (ь j,t)- Щ О, j,t -Д1 )]

Ді

(6)

(3)

(х у )е Д1 у, у = 1 где Иу - теплота сгорания горючего, Дж/кг.

Условие погасания (перехода в состояние

Б у ( х, у, t ) = 2)

Юу (х, у, 1г* ) = 0, (л; у)е Д2 у, у = 1,..., п.

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

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

^1 ^) и ^2 (t) .

Метод решения. Рассмотрим алгоритм построения горящей кромки, основанный на численном решении уравнений, описывающих распространение процесса горения. Введя в каждом из слоев прямоугольную сетку с шагами Дх и Ду по координатам х и у соответственно, заменим области Д0/ , Д1/ и Д2/ соответствующими сеточными областями Д

і = 1,2. Перейдем к дискретному времени, рассматривая систему в

(;,у)є Ддц, І =1,2 Условие погасания при і — і*

®І (i, у, і* ) = 0, ^ у ) є ДД1І, точка (і, у) исключается из Ддц (і*) и присоединяется к

ДД2І (і* ) •

Функции влияния Ь,іу (х - хь у - у1) зависят от характеристик леса, природн^іх и погодн^іх факторов и вычисляются по специальным подпрограммам [3]. Используем следующие исходные данные:

- область моделирования в виде двух множеств узлов

ДДІ у),; = 1,•••, п1,у = 1,•••, т1}, І = 1,2;

- начальный и конечный моменты времени іо и і у, временной шаг Ді;

- участки с одинаковыми характеристиками горючих материалов Цу, к = 1, •••, К, У Цу = Дді, І = 1,2;

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

Н0кІ (і), НкІ (і), ккІ (і), ЬкІ (і), ю0кІ, ФекІ (і),

юІ (і,у,і), І = 1,2, к = 1,^,К, і = і0,і0 +Ді, •••,іг;

- параметры функции, описывающей тепловое воздействие локального пламени Р0кі (і), а0кі (і), &кі (і),

Ьукі(і), фь (1) = ™ (1);

- скорости сгорания Гкі (/) для всех участков Ц в

каждый момент времени і = і0, і0 +Ді,^, іу,

к = 1,^, К;

- начальное состояние системы - области ДД0І (і0), Дд1І (і0 ) , ДД2І (і0 ) •

Используя информацию о лесном горючем и внешней среде, можно получить исходные данные для рассматриваемой модели. Наибольшую сложность представляет вычисление функции влияния пламени ^ (х - Х1, у - У1). При моделировании крупных многодневных пожаров оценку этой функции также можно получить путем анализа аэрокосмических снимков пожара, полученных в последовательные моменты времени (рис. 1). На основе этих снимков могут быть вычислены необходимые для расчетов значения функций ^ (х - хь у - у ) и скоростей. Поскольку конфигурация кромки реальных пожаров часто является достаточно сложной, то для упрощения расчетов целесообразно применять сглаживание границы контура на снимке пожара, например его аппроксимацию эллипсом [3]. Это позволяет оперативно получать необходимые для моделирования исходные данные.

Мясштяй

I____________________I______________________I

О 50 'ОС

' 16 ч 25 мин.

Рис. 1. Динамика контуров пожара

Распараллеливание вычислений. Дляраспараллели-вания процесса вычислений предлагается схема, вытекающая из физического содержания данной задачи. Расчет энтальпии для точки (г, ]) на (/ + 1)-м временном шаге происходит с использованием некоторого количества точек на им шаге. Численный расчет ведется итеративно: по имеющимся значениям /-го временного шага выстраивается (/ + 1)-й шаг и т. д. Таким образом, исходную задачу можно разбить на р подзадач для областей , Бд2г, Бдзг-, I = 1,2, к = 1,..., р, пересекающихся только по границе разбиений и независимых друг от друга на каждом расчетном шаге.

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

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

Сбалансированность вычислений и минимизация межпроцессорных обменов может быть обеспечена за счет выбора оптимального способа распределения данных и вычислений по процессорам. В рассматриваемом нами случае возможны два способа разбиения исходной области по вычислительным узлам: одномерное (на полосы) и двумерное (на квадраты).

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

Определим теперь потенциальное ускорение алгоритма. Будем оценивать время работы параллельной программы исходя из соотношения Sр = 7| / Тр . где Sр -ускорение; 7j - время вычислений на одном процессоре; Тр - время вычислений на р процессорах.

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

1D 3

^ comm = 2 - V • т при одномерном разбиении,

f ’2DComm = ~Т= ' ^ '1 при двумерном разбиении, у1Р '

где N - размерность задачи; р - количество вычислительных узлов; т - время пересылки одного числа. Алгоритм и его программная реализация являются масштабируемыми, если ускорение и производительность линейно зависят от количества используемых процессоров [6; 9]: Sp - O(p). На практике алгоритмы, для которых Sp - O(p/(Inp)), также считаются масштабируемыми.

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

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

Тестовая область W представляла квадрат W = [1,6л] х [1,6л]. Возвышение поверхности задавалось косинусоидой z(x) -1 + cos(x /150). Процесс распространения пожара инициировался из узла с координатами (0,5; 0,5). При запуске параллельных программ измерялось время их работы в секундах. На основе данных о времени работы параллельных программ определялись

другие их характеристики, такие как ускорение и эффективность распараллеливания (табл. 1). В вычислительных экспериментах было сделано по 100 шагов по времени.

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

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

1

. ■■ ■- ■

У Г.я

1 j 1 ' ' 1Ї ' 1 1 ■

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

—•—Ускорение при двумерном разбиении — ■ — ускорение при линейном разбиении Рис. 2. Зависимость ускорения от количества доступных процессоров

Также был проведен вычислительный эксперимент по выявлению зависимости ускорения вычислений от роста размерности задачи. Были рассмотрены случаи мелкой (10 000 узлов) и крупной (640 000 узлов) сетки. Тестовая область О представляла, как и ранее, квадрат О = [1,6л] х [1,6л], а возвышение поверхности задавалось

косинусоидой z( x) = i + cos( x /i5G). Характеристики горючего не изменялись. В вычислительных экспериментах было сделано также по iGG шагов по времени. Деком -позиция расчетной области - двумерная. Полученные значения ускорения и эффективности, показанные на кластере ИВМ СО РАН, приведены в табл. 2.

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

Значения ускорения вычислений с увеличением размерности задачи при использовании двумерной декомпозиции расчетной области возрастают логарифмически (рис. 3). Эту зависимость можно аппроксимировать

выражением y = G,3i5iln(x) - G,2G5 9 с достовернос-2

тью аппроксимации R = 0,882 1.

А

т"

ї I

м in II 'II II -in I nil ь in -hi : hi

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

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

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

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

Таблица 1

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

Характеристики параллельных программ Количество процессоров

i 2 4 8 i6

Двумерное разбиение расчетной области

Время выполнения основного цикла T, c 7GG,S6S 934 359,845 6 iSS,i42 GGi 81,975 424 53,299 i9

Ускорение S i,GGG GGG i,947 693 3,725 2i3 8,549 745 i3,i49 7i

Эффективность E = S / p i,GGG GGG G,973 S47 G,93i 3G3 i,G68 7i8 G,82i 857

Одномерное разбиение расчетной области

Время выполнения основного цикла T, c 7GG,S6S 934 359,8456 2G7,G77 7 iG5,259 7 87,268 G6

Ускорение S i,GGG GGG i,947 693 3,384 569 6,658 474 8,G3i 2i9

Эффективность E = S / p i,GGG GGG G,973 S47 G,S46 i42 G,832 3G9 G,5Gi 95i

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

В результате исследования были созданы следующие программы:

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

- программа для препроцессорной обработки данных - разрезания данных на отдельные файлы для многопроцессорных вычислений;

- программа для численного моделирования распространения кромки лесного пожара.

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

Библиографический список

1. Софронов, М. А. Огонь в лесу / М. А. Софронов,

А. Д. Вакуров. Новосибирск : Наука. Сиб. отд-ние, 1981.

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

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

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

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

6. Воеводин, В. В. Параллельные вычисления /

В. В. Воеводин, Вл. В. Воеводин. СПб. : БХВ-Петербург, 2002.

7. Message-Passing Interface Forum, MPI 2: Extensions to the Message-Passing Interface [Electronic resource]. Electronic data. Access mode: http://www.unix.mcs.anl.gov/ mpi/. Title from display.

8. RS/6000 SP: Practical MPI Programming [Electronic resource]. Electronic data. Access mode: http:// www.redbooks.ibm.com. Title from display.

9. Корнеев, В. Д. Параллельное программирование в MPI / В. Д. Корнеев. 2-е изд. Новосибирск : Изд-во Ин-та вычисл. математики и мат. геофизики Сиб. отд-ния Рос. акад. наук, 2002.

Таблица 2

Значения ускорения вычислений и эффективности распараллеливания в зависимости от размерности задачи

Характеристика параллельных программ Размерность задачи

iGGxiGG 2GGx2GG 4GGx4GG SGGxSGG

Время выполнения основного цикла на одном процессоре Г,с G,464 5i5 7,3GG 3SG i22,96i SGG i 532,G29 398

Время выполнения основного цикла на четырех процессорах Г, С G,iS5 772 2,i49 292 33,887 iiG 394,943 4GG

Ускорение £ 2,5GG 462 i 3,396 643 8 3,628 57i 5 3,879 iii

Эффективность Е = Б / р G,625 ii6 G,S49 i6i G,9G7 i43 G,969 778

M. S. Vdovenko, G. A. Dorrer

ON EFFICIENCY ANALYSIS OF THE PARALLEL ALGORITHMS FOR MODELING THE FOREST FIRE SPREADING

The numerical modeling offorest fire spread dynamics with parallel algorithms is considered. The various ways of decomposition of initial burning area are compared. The results ofthe preliminary theoretical analysis of the calculations acceleration and its comparing with results of computing experiment shows, that the greatest efficiency is achieved, at least, when using from 4 up 16 processors cluster to at two-dimensional splitting initial area.

Keywords: parallel algorithm, modeling of the forest fire, efficiency.