CC BY
22
УДК 519.63
МОДЕЛИРОВАНИЕ ЛЕСНЫХ ПОЖАРОВ С ИСПОЛЬЗОВАНИЕМ КЛАСТЕРНЫХ ВЫЧИСЛИТЕЛЬНЫХ СИСТЕМ
М.С. Вдовенко, Г.А. Доррер
ГОУ ВПО «Сибирский государственный технологический университет»
660049 Красноярск, пр. Мира, 82; e-mail: mail [email protected]
Рассмотрена возможность использования параллельных вычислительных систем и космических снимков для моделирования динамики лесных пожаров. Предложен способ оценки характеристик динамики пожара по космическим снимкам. Проанализированы различные способы декомпозиции сеточной области при численном моделировании динамики распространения лесного пожара. Приводятся оценка ускорения вычислений по данным эксперимента, выполненного на кластере ИВМ СО РАН.
Ключевые слова: моделирование лесных пожаров, параллельные алгоритмы, эффективность
The possibility of numerical modeling the forest fire dynamics with using the parallel computing systems and space photographs is considered. The method of estimation of the fire parameters on basis of space imagines is suggested. The various ways of decomposition of burning area are analyzed. The estimations of acceleration and its comparing with results of computing experiment on cluster of the Institute of Computing Modeling are implemented.
Keywords: modeling of the forest fires, parallel algorithms, efficiency
ВВЕДЕНИЕ
Известна роль лесных пожаров в формировании лесных биогеоценозов. Кроме того, они могут оказывать негативное влияние на состояние окружающей среды. Так, в числе причин эффекта глобального потепления климата считают выбросы СО2 и других многочисленных газов от пожаров (Софронов, 1981, Гришин, 2008, Доррер, 2008). Помимо глобального влияния на климат Земли, лесные пожары наносят ущерб лесному фонду, а также способствуют возникновению таких катастроф как большая задымленность больших городских территорий (Москва 2002 г.) или же массовых городских пожаров (США, Лос-Аламос, май 2002 г., Австралия, январь 2009г.).
Проблема математического моделирования процессов горения при лесных пожарах изучается уже в течение многих лет. Обзор результатов, полученных в этой области, приведен в работе (Weber, 1990). Большой вклад в решение данной проблемы внесли Н. П. Курбатский, Э. Н. Валендик, М. А. Софронов, А. М. Гришин, Г. Н. Коровин, R. Rothermel, M. E. Alexander и другие ученые.
Разработка математических моделей распространения пожара позволяет предсказать его поведение, что может помочь более эффективному проведению противопожарных мероприятий. Однако ключевой проблемой при этом является необходимость сбора большого количества информации об условиях горения и противопожарных мероприятиях. В последнее время в связи с созданием и вводом в эксплуатацию Информационной системы дистанционного мониторинга ИСДМ-Рослесхоз (Применение,2007), основанной на использовании спутниковой информации о пожарной обстановке в лесах, сложились благоприятные условия для разработки систем моделирования и прогнозирования лесных пожаров на всей территории России.
*
Работа поддержана РФФИ (проект 05-07-90201_в)
Следует отметить, что для решения задач моделирования крупных многодневных лесных пожаров требуются значительные вычислительные ресурсы, и использование кластерных вычислительных систем является одним из способов решения данной проблемы. При решении задач, связанных с исследованием природных процессов на больших территориях, общим подходом для получения равномерного распределения вычислительной нагрузки между процессорами является разделение вычислительной области на подобласти, количество которых совпадает с числом используемых процессоров, т. е. использование принципа геометрической декомпозиции (Воеводин, 2002).
Целью настоящей статьи является разработка параллельного алгоритма, основанного на разбиении моделируемой области на равновеликие подобласти и исследование его эффективности при использовании модели лесного пожара типа бегущей волны (Доррер, 2008).
АНАЛИЗ ЗАДАЧИ И ВЫЯВЛЕНИЕ ЕЕ ПОТЕНЦИАЛЬНОГО ПАРАЛЛЕЛИЗМА
Разработка параллельных программ практического уровня сложности представляет собою многоэтапный технологический процесс: постановка и анализ задачи, выбор модели программирования, декомпозиция задачи на параллельные процессы, анализ производительности и организации вычислительного эксперимента.
Постановка задачи. В качестве базовой модели процесса распространения лесного пожара нами взята модель, основанная на вычислении теплового баланса в лесных горючих материалах (Доррер, 2008). Массив горючего представляет собой в общем случае п параллельных однородных слоев горючего, расположенных один над другим. Произвольный I -й слой занимает по вертикали область Zi, с координатами от ZlH до ZiK и содержит запас
горючего материала соЩ^,у,г,і^ кг/ м3. Вертикальная координата середины і -го слоя обозначается ^ , а его толщина с),. При этом
2ісР = + 2ік 'І*2 , & і = гіК - гіН ,
ІСр ?
/ = 1,..., п •
Свойства горючего в пределах каждого слоя не зависят от 2 .
Горючий материал в окрестностях точки С в некоторый момент времени может находиться в одном из трех состояний, описываемых функцией: Б у, z,t
0 если в точке С в момент ? имеется ненулевой запас горючего (т. е. (О 4с, у, I > 0), но горения не происходит; если (О у, г, ^ > 0 и происходит горение;
если <х>4ь,у,2^ = 0 т. е. горение невозможно.
у, г, Г=
Области, соответствующие состояниям ^ = О , ^ = 1, Б = 2 обозначаются соответственно
о0 02 Проекции областей °о,йг ^2 на
горизонтальную плоскость Б будем обозначать соответственно £>о, А, /^2 причем В0 'и 'и В2 = В ■
Уравнение нагрева для горючего в окрестности точки С = 4с, у, г", которое в момент времени I находится в состоянии Б^,у,г,і^= О'.
----------= І II Ф/ V к-х1,у-у1 +
Ы '=Ц,- ;
(1)
где
н X, у, 0=
г — г тт
2
)к
\нл , у, г, t ^г ,
при начальном условии н ^ //о у
4і,у^^Во , / = 1,..., п ■ Здесь Н Щ(,у,г,і
зна-
чение энтальпии в точке ^,у,г , в момент времени I: Ф, Ф - соответственно, энергия, образующаяся при горении и поступающая от внешних источников, вт/м3; ^ ^ — Х^у — у1 - функция
влияния пламени из точки ^, уг на точку 4", у (функция Грина) из / -го слоя на / -й слой;
к С у 2 > где а - коэффициент теплоот-
’ ’ ^ рс'
дачи, вт/м2град; 8 - удельная поверхность слоя, м" с' - приведенная теплоемкость влажного материала, дж/кг град.
Условие воспламенения горючего в / -м слое, т.е. перехода горючего в состояние Sj у, г, ? 1 имеет вид:
Н] у, Г > Н* 4,У^, где ** = г* Щс, у - вре-
мя воспламенения горючего в у -м слое в точке Х,у~^ , II - энтальпия начала газифика-
ции.
Уравнение расходования горючего:
г ■ при со ■ у, ґ 3* 0 0 при со ■ у, ґ ^ 0
(2)
с начальным условием
<о1 ( у, г* }= со0]. <\ .у _• С/3е А/ •
Здесь О) ■ - активный запас горючего материала, в / -м слое, кг/м3; г. - относительная скорость сгорания / -го слоя, 1/с.
Уравнение тепловыделения в ] -м слое:
4,уУс>У
(3)
где - теплота сгорания горючего, дж/кг. Условие погасания (перехода в состояние
Б 4і,у,і^= 2):
к, У, t
= О •
X, у^
: /)
7 = 1
Модель теплопередачи из зоны горения к горючим материалам, задаваемая функцией влияния (функцией Грина) ^ь.^ — х1,у — у1 , подробно
рассмотрена в (Доррер, 2008).
Система уравнений (1) - (3) представляет собой модель процесса распространения горения по неоднородным слоям горючего типа бегущей волны. Особенностью рассмотренной модели является то, что часть входных и выходных переменных представляет собой множества ^ - -21 ^ - " ^2 ^
МЕТОД РЕШЕНИЯ
Алгоритм построения горящей кромки основан на численном решении уравнений (1) - (3), описывающих распространение процесса горения. В каждом из слоев вводится прямоугольная сетка с шагами по координатам х и у соответственно Ах- и Ау, области .О , ^ и ^ заменяются соответствующими сеточными областями £)д0., £>м., ВА21, I = 1,2. Перейдя к дискретному времени I = 0,1,2,... с
шагом At и заменив в (1) - (2) частные производные по времени разностными отношениями, а интеграл по области суммой, получим расчетные уравнения,
которые использовались при дальнейших исследованиях.
1
Функции влияния X — Л',, у — у, зависят
от характеристик леса, природных и погодных факторов и вычислялись по специальным подпрограммам. В качестве исходных данных использовались:
- область моделирования в виде двух множеств узлов
=1>■■■,”/,7 =!,■■■, >«/ // = 1,2;
- начальный и конечный моменты времени и tf9 временной шаг № ;
- участки с одинаковыми характеристиками горючих материалов С1к1, к = 1,..., К,
=Ач,- / = 1,2;
к
- теплофизические характеристики горючего для каждого из участков о в каждый момент
времени: нш < //;; <3 ки С К С
®он'Фей ^ С7,= 1,2. £ = 1,..., АГ,
/ = г0,г0 + ы,..., 1-г;
- параметры функции, описывающей тепловое воздействие локального пламени рои ^ .
«ои С 5и <3 < (рь <3 и’О
- скорости сгорания ги К. , для всех участков в каждый момент времени
г = г0,г0 + А?,..., «у. Ат = ■
- начальное состояние системы - области
-^до7 X) Аш <0 -®Д2/ ^0
Используя информацию о лесном горючем и внешней среде, можно получить исходные данные для рассмотренной модели. Наибольшую сложность представляет вычисление функции влияния пламени ^ у — уг". При моделировании крупных
многодневных пожаров оценку этой функции можно получить также путем анализа аэрокосмических снимков пожара, полученных в последовательные моменты времени. На основе этих снимков могут быть вычислены необходимые для расчетов значения функций ^ 7 у — у, и скоростей. Поскольку
конфигурация кромки реальных пожаров часто является достаточно сложной, то для упрощения расчетов целесообразно применять сглаживание границы контура на снимке пожара, например, аппроксимацию его эллипсом (Доррер, 2008). Пример такой оценки приведен на рисунке 1. Это позволяет оперативно получать необходимые для моделирования исходные данные.
Рисунок 1 - Аппроксимация контуров пожара эллипсом
Распараллеливание вычислений. Для распараллеливания процесса вычислений предлагается схема, вытекающая из физического содержания данной задачи. Расчет энтальпии для точки (/зу)на (7 + 1) -ом временном шаге происходит с использованием некоторого количества точек на ? -м шаге. Численный расчет ведется итеративно: по имеющимся значениям ? -го временного шага выстраивается (/ +1) -й, и т.д.
Таким образом, исходную задачу можно разбить на р подзадач для областей, /У^, £)* ,
I)дз , / = 1,2, к = 1,..., р •, пересекающихся только по границе разбиений, независимых друг от друга на каждом расчетном шаге. В случае технологии МР1 каждый процессор получает часть сетки, причем сетки соседних процессоров пересекаются по двум узлам. Это пересечение сеток позволяет частично продублировать вычисления на соседних процессорах, что сокращает межпроцессорные обмены.
Далее, в каждом процессоре происходят вычисления последовательно по шагам по времени на части сетки. На каждом шаге по времени соседние процессоры осуществляют обмен граничными значениями при помощи функций неблокирующих пересылок и приема данных МР1_ ^еМ и MPI_Irecv. По окончании расчетов каждый процесс обрабатывает свой массив вычисленных значений, записывая его в отдельный файл.
Для сбалансированности вычислений и минимизации межпроцессорных обменов ключевая роль отводится выбору способа распределения данных и вычислений по процессорам.
В рассматриваемом случае возможны два различных способа разбиения исходной области по вычислительным узлам - одномерное и двумерное разбиение. В обоих случаях исходная область включает взаимно перекрывающиеся подобласти, и пересчет значений на границах между данными областями предполагает согласно алгоритму суммирование энтальпий при обмене вычисленными значениями для граничных элементов.
При этом для перехода к следующей итерации необходимо согласование значений на границах расчетных подобластей. Пересылка данных осуществляется с использованием процедур библиотеки МР1 (Корнеев, 2002).
Определим теперь потенциальное ускорение алгоритма. Будем оценивать время работы параллельной программы исходя из соотношения Б = 1\ / 7 '/1 ■ где Б - ускорение, 7\ - время вычислений на одном процессоре, Т - время вычислений на р процессорах.
Для получения реалистических оценок помимо классической формулы Амдала (Воеводин, 2002) будем учитывать также время, затрачиваемое программой на обмены между процессами. Как следует из принятой нами схемы распределения данных, на каждом временном шаге требуется обмен границами.
Время пересылок для различных способов декомпозиции можно приблизительно выразить через количество пересылаемых данных (Воеводин, 2002):
нии,
V 2
= 2- N ■ т, при одномерном разбие-, при двумерном разбие-
нии, где N - размерность задачи, р - количество вычислительных узлов, г - время пересылки одного числа. Алгоритм и его программная реализация являются масштабируемыми, если ускорение и производительность зависят линейно от количества используемых процессоров (Воеводин, 2002, Корнеев, 2002): Б = ()(р) • На практике алгоритмы,
для которых Б = 0{р /{\пр)) также считаются масштабируемыми.
ПРОВЕДЕНИЕ ВЫЧИСЛИТЕЛЬНОГО ЭКСПЕРИМЕНТА И АНАЛИЗ РЕЗУЛЬТАТОВ
Параллельные программы тестировались на модельных задачах. Для модельной задачи был выбран один однородный слой горючего, расположенный на местности с уклоном, который изменяется по заданному закону. Расчеты производились по формулам (4) - (6) на кластерной системе МВС 1000 ИВМ СО РАН (г. Красноярск) на тестовой сетке 400x400 узлов при использовании от 1 до 16 процессоров. Тестовая область П представляла квадрат О = [1.6л] х [1.6л]. Возвышение поверхности задавалось выражением г(х) — 1 + соз(х /150). Процесс распространения пожара инициировался из узла с координатами (0.5, 0.5). При запуске параллельных программ измерялось время их работы в секундах.
Таблица 1 - Значение ускорения, показанного на кластере ИВМ СО РАН
Количество
процессоров
1
2
4
8
16
Время выполнения основного цикла, Т[с] Ускорение, Б Эфективность, Е = Б /р
Время выполнения основного цикла, Т[с] Ускорение, Б Эфективность, Е = Б / р__________
Разбиение расчетной области на квадраты (двумерное разбиение)
700.868934 359,8456 188,142001 81,975424
1 1,947693 3,725213 8,549745
1 0,973847 0,931303 1,068718
Разбиение расчетной области на полосы (одномерное разбиение)
700.868934 359,8456 207,0777
1,947693
0,973847
3,384569
0,846142
105,2597
6,658474
0,832309
53,29919
13,14971
0,821857
87,26806
8,031219
0,501951
На основе данных о времени работы программ (табл. 1) вычислялись другие характеристики параллельных программ, такие как ускорение и эффективность распараллеливания. В вычислительных экспериментах было сделано по 100 шагов по времени. Значения ускорения, показанного на кластере ИВМ СО РАН приведены в таблице 1.При применении двумерной декомпозиции и 8 процессоров значение эффективности больше единицы, что объясняется использованием в программе «динамических» массивов с подстраиваемыми под выделенное число процессоров размерами.
Таким образом, необходимы меньшие временные затраты на выборку обрабатываемых данных из оперативной памяти и передачу их через КЭШпамять. В случае использования 8 процессоров при данной размерности сетки весь массив помещается в КЭШе, что и определяет более быстрое выполнение вычислений за счет отсутствия необходимости обмена между оперативной памятью и КЭШем. На рисунке 2 приводятся полученные графики зависимости ускорения алгоритмов и эффективности распараллеливания в зависимости от количества используемых процессоров. Как видно из рисунка 2 результаты вычислительного эксперимента показали наличие хорошего ускорения при решении данной задачи.
- Ускорение при двумерном разбиении - -■ - - Ускорение при линейном разбиении
Рисунок 2 - Зависимость ускорения от количества доступных процессоров
Также был проведен вычислительный эксперимент, выявивший зависимость ускорения от роста размерности задачи. Были рассмотрены случаи мелкой (10000 узлов) и крупной (640000 узлов) сетки. Тестовая область О представляла, как и ранее, квадрат: □ =
11 _6тт | х 11,6тг|. а возвышение поверхности задавалось формулой г(х) = 1 + совСх /150). Характеристики горючего также не изменялись. В вычислительных экспериментах было сделано также по 100 временных шагов. Декомпозиция расчетной области - двумерная. Полученные значения ускорения показанного на кластере ИВМ СО РАН (г. Красноярск) приведены в таблице 2.
Таблица 2 - Значение ускорения, показанного на кластере ИВМ СО РАН
Размерность задачи 100x100 200x200 400x400 800x800
Время выполнения основного цикла на 1 0,464515 7,30038 122,9618 1532,029398
процессоре, Т, с Время выполнения основного цикла на 4 0,185772 2,149292 33,88711 394,9434
процессорах, Т, с Ускорение, Б 2,5004621 3,3966438 3,6285715 3,879111
Эффективность, Е = Б /Р 0,625116 0,849161 0,907143 0,969778
Выполненные расчеты показали, что с увеличением размерности задачи при использовании четырех процессоров наблюдается увеличение ускорения вычислений, по крайней мере, для задач размерностью не более 800^800. Из таблицы 2 видно, что при учете размера кэш-памяти можно добиться эффективности больше 1. Значения ускорения вычислений с увеличением размерности задачи при использовании двумерной декомпозиции расчетной области возрастают логарифмически. Эту зависимость можно аппроксимировать выражением у = 0,3151Ьп(х) - 0,2059, с достоверностью аппроксимации Н2 = 0,8821.
2 -------1------1------1------1------1-------1------
О 100000 200000 300000 400000 500000 600000 700000
Размерность задачи
Рисунок 3 - Зависимость ускорения от размерности задачи
На рисунке 3 показана зависимость ускорения вычислений от размерности задачи. При увеличении размерности задачи в 4 раза ускорение вычислений при использовании 4 процессоров возрастает. Средний уровень эффективности распараллеливания составляет около 83 %, что связано с неизбежными затратами времени на организацию межпроцессорных обменов и записи результатов в файл.
ЗАКЛЮЧЕНИЕ
Таким образом, при численном решении задачи моделирования динамики лесного пожара возможно применение как одномерной, так и двумерной декомпозиции. Однако результаты тестирования программ показали, что наиболее эффективной при использовании, по крайней мере, от 4 до 16 процессоров является двумерное разбиение исходной области. В результате исследования были созданы следующие программы:
программа, решающая обратную задачу - получения значений скоростей на основе аэрокосмических снимков контуров пожара в последовательные моменты времени; программа для препроцессорной обработки данных - «разрезания» данных на отдельные файлы для многопроцессорных вычислений; программа для численного моделирования распространения кромки лесного пожара.
Получены расчетные значения ускорений, позволяющие оценить масштабируемость алгоритма и его программной реализации. Эти результаты показывают, что алгоритм обладает значительным объемом потенциального параллелизма и хорошей, с точки зрения распараллеливания, структурой, что позволяет надеяться на получение ускорений, близких к линейным в зависимости от количества используемых процессоров для кластерных вычислительных систем. Значения ускорений, полученные в ходе вычислительного эксперимента, хорошо согласуются с теоретическими оценками.
БИБЛИОГРАФИЧЕСКИЙ СПИСОК
Воеводин, В. В. Параллельные вычисления. / В. В. Воеводин, Вл. В. Воеводин. - СПб.: БХВ-Петербург, 2002. -608 с.
Гришин, A. M. О математическом моделировании природных пожаров и катастроф. / А. М. Гришин // Вестник Томского государственного университета: Математика и механика. - Томск: Томский государственный университет, 2008 г. №2(3). С. 105-114.
Доррер, Г. А. Динамика лесных пожаров. / Г. А. Доррер. -Новосибирск: Изд-во СО РАН, 2008. - 404 с.
Корнеев, В. Д. Параллельное программирование в MPI. / В. Д. Корнеев. - 2-е изд. Новосибирск: Изд-во ИВМиМГ СО РАН, 2002. - 215 с.
Применение информационной системы дистанционного мониторинга «ИСДМ-Рослесхоз» для определения пожарной опасности в лесах Российской федерации: Учебное пособие - г. Пушкино (МО), ФГУ «Авиалесо-охрана», 2007 год. - 82 с.
Софронов, М. А. Огонь в лесу. / М. А. Софронов, А. Д.
Вакуров. - Новосибирск: Наука, 1981. 128 с. Message-Passing Interface Forum, MPI-2: Extensions to the Message-Passing Interface, 1997. [http://www.unix.mcs.anl.gov/mpi/], 13.03.2007 RS/6000 SP: Practical MpI Programming.
[www.redbooks.ibm.com], 11.08.2008 Weber, R. O. Modeling fire spread through fuel beds /R. O. Weber // Prog. Everg. Combust. Sci. 1990. V. 17. P. 65-82.
Поступила в редакцию 3 октября 2008 г. Принята к печати 8 июня 2009 г.
