Параллельные алгоритмы моделирования процессов распространения лесных пожаров на основе математических моделей различных типов Текст научной статьи по специальности «Математика»

Вычислительные технологии

Том 18, № 1, 2013

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

М. С. Вдовенко, Г. А. Доррер, П. С. Шаталов Сибирский государственный технологический университет, Красноярск

e-mail: shatalovp@gmail.com

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

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

Введение

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

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

1. Математическая модель процесса

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

В зависимости от целей моделирования и уровня принимаемых решений можно выделить три уровня моделей: оперативные, тактические и стратегические [2]. На оперативном уровне производится достаточно подробное вычисление теплофизических параметров, как правило, одного крупного пожара на основе уравнений тепло- и массообме-на при горении растительных материалов в пологе леса [2]. Оперативные модели требуют численного решения систем дифференциальных уравнений в частных производных с соответствующими начальными и граничными условиями, вид которых зависит от конкретной задачи. Для получения дискретных аналогов используются численные методы (например, итерационно-интерполяционный, дифференциально-разностный, контрольного объёма и т. д.). В результате получаются сложные системы алгебраических уравнений, решаемые точными, а чаще всего приближёнными методами. Для ускорения этих расчётов целесообразно использовать технологии высокопроизводительных параллельных систем.

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

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

2. Физико-математические модели оперативного уровня

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

уровень распараллеливания уже каждого отдельно взятого подпроцесса, отвечающего какому-либо физическому процессу. Следует заметить, что весьма важным фактором является попытка отыскать в самой задаче так называемый естественный параллелизм, и результат этого поиска во многом зависит от того, какие стороны физического процесса (лесного пожара) выделяет исследователь. В качестве примера математической модели рассмотрим универсальную аэротермохимическую модель лесного пожара, разработанную А.М. Гришиным. Данная модель состоит из систем дифференциальных уравнений в частных производных:

f + f = * ^ = !• 2- 3, ^ = t - (D

J i=5

dvi dp dnj

= + PFi + - Qvi-

dt dxi dxj

3 3

1 X ^ / / / 4 2

2 / y v-i ~is/ ai 1 \vi vis)\vj vjs) aj i=l j=l

-PQS vi|v| + - Y^ \vi - vis)2 cos ai + Y1 \vi - vis)\vj - vjs) cos a

" = " + x " x" + 2"" x "ft, i,j = 1, 2, 3,

dT А Л dp д / дТ\ д ( дсЛ

^ ъ р^ = 2. К dt + 1 + дХ" (,Лэф" д^ + дХ" ^ °эфаСрадХ" -

¿=5 ¿=5 " у "' " \ а=1 " /

- (к + к(з)) [сЛя - В(Т)] + срз{Т8 - Т)(1 - ас)К1з + срТ{Т3 - Т)К2з + Ф2Я52+

+^53^53 + ^54^54 + ?2з (лв- - ^8+) . (3)

При выводе этой системы использовались следующие величины: £ — время, 1Г — радиус-вектор любой точки; Х" — декартовы координаты этой точки; к,к(з),кз — интегральные коэффициенты поглощения газа, дисперсных частиц и конденсированной фазы; В — функция Планка; Ип — интегральная плотность излучения; Сpi,Сpт, срз — теплоёмкости при постоянном давлении отдельных фаз, водяного пара и газообразных продуктов пиролиза; 1", |1| — пульсационные составляющие скорости потока и изгиб-ных колебаний элементов лесного горючего материала; pi — истинная плотность г-й фазы; = Я|+ - , Я|+, — массовые скорости конденсации паров и испарения свободной воды в газодисперсной фазе; 513,^2з,53з — теплоты пиролиза, испарения

связанной воды и горения коксика; q52,^53,^54 — теплоты горения СО, Н2, СН4; р —

(1) (2) ( \ давление газа; т" , т" — компоненты тензоров касательных (тангенциальных) напряжений для ламинарных и турбулентных течений; са — массовая концентрация; а — компонент в газодисперсной среде; N — количество компонентов в газодисперсной сре-де.Уравнение (1) представляет собой закон сохранения массы газодисперсной фазы; уравнение (2) — закон сохранения количества движения газодисперсной среды в проекциях на оси декартовой системы координат. В (2) входят члены, обусловленные силовым взаимодействием газодисперсного потока с костяком пористо-дисперсной среды. Уравнение (3) представляет собой закон сохранения энергии в газодисперсном потоке с учётом переноса энергии как конвекцией, так и излучением, а также выделения и поглощения тепловой энергии в результате различных физических и химических процессов. Данная часть уравнений даёт общее представление о сложности модели.

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

3. Метод геометрической декомпозиции расчётной области на тактическом уровне

В качестве базовой модели процесса распространения лесного пожара принята система уравнений теплового и материального баланса в массиве лесных горючих материалов, состоящем в общем случае из п параллельных плоских слоев [2, 3]. Горючий материал в каждом слое может находиться в одном из трёх состояний: 50 — исходное состояние материала, который способен загореться, 51 — материал находится в состоянии горения, 52 — горение материала невозможно. Проекции этих слоев на горизонтальную плоскость обозначаются соответственно D0j, П?, , где ] — номер слоя. Каждый слой по вертикали ограничен координатами ZjH, ZjK (^ = 1,..., п).

Усреднённое уравнение нагрева горючего материала в ]-м слое в окрестности точки С = (ж,у), которое в момент времени £ находится в состоянии Б(ж,у,£) = Б0, имеет вид

дН-(ж у £) ' Л Г Г

-' Фг(ж1У1*)^-(ж - Ж1У - У1)^МУ1 +

+ Фе(ж,У,£) - (ж,У)[Н? (ж,у,£) - Яс7 (х^у)].

где

Ъ К

Н?(ж^ о) = /

j j ъ Н

при начальном условии Н?(ж, у, 0) = Н0?(ж, у), (ж, у) € П0?, г = 1,..., п. Здесь Н?(ж, у, £) — значение энтальпии в ]-м слое в точке (ж, у) в момент времени Дж/м3; Ф^, Фе — соответственно энергия, образующаяся при горении и поступающая от внешних источников, Вт/м3; е? = (ж - ж1, у - у1) — функция влияния пламени из точки (ж1, у1) на точку (ж, у) (функция Грина) из г-го слоя на ]-й слой; к(ж,у^) = ав/рс', где а — коэффициент теплоотдачи, Вт/(м3- град); в — удельная поверхность слоя, м-1; с' — приведённая теплоёмкость влажного материала, Дж/(кг • град).

Условие воспламенения горючего материала в ]-м слое, т. е. перехода его в состояние б?(ж, = 51, имеет вид Н7(ж, у,£*) > Н*(ж,у), где = £*(ж,у) — время воспламе-

нения горючего в ^'-м слое в точке (ж, у) € П0?, Н* — энтальпия начала газификации.

Уравнение расхода горючего с начальным условием ш?(ж,у,£*) = ш0?(ж,у), (ж,у) € П?, имеет вид

д_

(ж,У,£) ш0? (ж,У)

при (ж,у,£) > 0, (5)

0 при ш?(ж,у, £) = 0.

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

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

дш' (х у ~Ь)

Фj(х,у,1) = — ^(х,у) е D1j■, 3 = 1,...,п, (6)

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

Условие погасания (перехода в состояние Sj(х, y,í) = $2) следующее:

шj(х,у,^*) = 0, (х,У) е ^, 3 = 1,...,п. (7)

Весь комплекс теплофизических параметров сосредоточен в функциях тепловыделения Фj(х,у,¿) и влияния (функция Грина) eij = (х — х1,у — у1), при этом функция Грина зависит от hf — высоты пламени, р0 — характерной длины влияния пламени, а/ (ш) — угла отклонения пламени от вертикали, — угла, определяющего направление ветра, и ш — скорости ветра [2]. Указанные параметры либо вычисляются с помощью специальных формул, либо оцениваются по данным наблюдений. При моделировании динамики крупных многодневных пожаров для этой цели используются данные космического мониторинга лесов.

На рис. 1 представлено сравнение расчётных результатов, полученных на основе модели (1)-(4), и экспериментальных данных, зафиксированных лабораторией лесной пирологии Института леса им. В.Н. Сукачёва СО РАН. Достаточно точное совпадение

Рис. 1. Сравнение расчётной конфигурации лесного пожара с экспериментальными данными

расчёта с данными дистанционной съемки конфигурации реального лесного пожара в последовательные моменты времени подтверждает достоверность модели процесса распространения лесного пожара [3].

4. Вычислительный алгоритм

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

• считывание исходных данных из соответствующих файлов;

• определение сфер влияния локального пламени, зависящих от направления и силы ветра, характеристик горючего материала, особенностей рельефа местности (склон, равнина, каньон и др.);

• пересылка информации (теплофизических характеристик горючего) о расчётных точках, расположенных на границе области вычисления данного процессора соседним процессором для обеспечения согласования вычислений на границах и балансировки вычислительной нагрузки;

• определение границ контура горения, включающее:

— параллельную обработку исходных данных (расчёт теплового баланса на данном временном шаге) по формуле (1);

— обмен данными между соседними процессорами (передача информации о величине энтальпии для всех граничных данных);

— вычисление значения энтальпий для каждой граничной точки согласно алгоритму суммирования, определение воспламенившихся и погасших точек по формулам (2)-(4);

• запись полученных результатов в выходной файл;

• освобождение занимаемой памяти.

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

5. Моделирование операций межпроцессорного обмена

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

Для моделирования параллельных процессов были использованы иерархические сети Петри — традиционный аппарат для исследования параллельных систем. Сеть имеет

Рис. 2. Модель обмена данными одного процесса с соседними; позиции: Po — разрешение обмена, Pg — завершение обмена, Pi — P4 — данные для обмена с соседними процессами, P5 — Pg — завершение обмена по направлениям, th, tt — начало и завершение операций обмена; обмены: т1 — up-dw, т2 — left-right, т3 — upleft-dwright, т4 — dwleft-upright

Рис. 3. Модель операций межпроцессорного обмена на примере операции т1: Р12—Р14 — данные, приготовленные для обмена, Рц — Р13 — данные, полученные при обмене, тц — Г12 — операции пересылки данных, К1, К2, К3, К4, К5 — маркеры готовности к обмену данными, — начало операции обмена, — конец операции обмена

два вида узлов: позиции и переходы. Первые моделируют наличие данных, вторые — операции. Полная модель отражает все этапы вычислительного процесса. Одна из операций — обмен данными с соседними процессорами — показана на рис. 2. Более детально схема операции обмена на примере операции т1 приведена на рис. 3. Остальные операции межпроцессорного обмена выполняются аналогично.

6. Исследование способов разбиения расчётной области

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

ет взаимно перекрывающиеся подобласти. Пересчёт значений энтальпий на границах Н2(х,у, ¿) между данными областями предполагает, согласно алгоритму, их суммирование при обмене вычисленными значениями для граничных элементов. При этом для перехода к следующей итерации необходимо согласование значений на границах расчётных подобластей.

Основными критериями качества параллельных реализаций алгоритма являются, во-первых, ускорение вычислений с ростом числа процессоров и, во-вторых, адекватность воспроизведения моделируемых явлений. Последнее подтверждается достаточно точным совпадением расчётных результатов с данными дистанционной съёмки конфигурации реального лесного пожара в последовательные моменты времени (см. рис. 1). Ускорение алгоритма определяется отношением времени счёта последовательного алгоритма ко времени счёта его параллельного варианта. Эффективность распараллеливания определяется отношением ускорения к количеству процессоров, на котором оно достигнуто. Стопроцентная эффективность соответствует идеальному ускорению, когда использование р процессоров в вычислениях дает ускорение в р раз. Для получения реальных оценок учитывалось время, затрачиваемое программой на обмены. Было проведено несколько серий вычислительных экспериментов на модельной задаче, в качестве которой рассматривался случай неоднородной среды (слой горючих материалов с возвышением). При запуске параллельных программ измерялось время работы основного цикла в секундах. На основе данных о времени работы программ вычислялись такие характеристики параллельных программ как ускорение и эффективность распараллеливания. Расчёты производились на кластерной системе ИВМ СО РАН (г. Красноярск) на тестовой сетке 400 х 400 при использовании до 16 процессоров.

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

Далее был проведён вычислительный эксперимент по установлению зависимости ускорения от роста размерности расчётной области (рис. 5) в случаях от крупной (10 000

16

О 2 4 Щ 8 10 12 14 16 18

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

Рис. 4. Зависимость ускорения алгоритмов от количества используемых процессоров: ■ — расчётное ускорение соответственно при двумерном и одномерном разбиении; прямые линии 1,2 — теоретическое ускорение соответственно при одномерном и двумерном разбиении

4 3.5

со

¡0 , а 3

о к

о р

и

2 -I-1111111

О 100000 200000 300000 400000 500000 600000 700000

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

Рис. 5. Зависимость ускорения от размерности расчётной области

узлов) до мелкой (640 000 узлов) сетки. Декомпозиция расчётной области двумерная. При увеличении размерности расчётной области в четыре раза ускорение вычислений возрастает, что указывает на наличие у параллельного алгоритма свойства масштабируемости [4].

7. Исследование влияния величины зоны перекрытия между расчётными подобластями

Независимо от выбранного способа декомпозиции расчётной области и количества расчётных узлов подобласти, полученные в результате разбиения, включают взаимные перекрытия (рис. 6). Величина перекрытия зависит от параметров функции влияния

<->

Перекрвггие

Рис. 6. Область влияния локального пламени

14,0000

12,0000

а ю.оооо

©

а

о И

8.0000

6.0000

щ 4.0000

2.0000

0.0000

/у_ г:

г

Г Л?^

•+■■ Перекрытие в 1 расчётную точку Перекрытие в 4 расчётные точки Перекрытие б 8 расчётных точек -в - Перекрытие в 13 расчётных точек

6 в 10

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

—Перекрытие в 2 расчётные точки —— Перекрытие в 5 расчётных точек Перекрытие в 10 расчётных точек - Перекрытие в 15 расчётных точек

12

14

16

1е

-■*- - Перекрытие в 3 расчётные точки Перекрытие в 6 расчётных точек —л— Перекрытие в 12 расчётных точек

Рис. 7. Зависимость ускорения вычислений от величины области перекрытия

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

Данные, полученные на основании серии расчётов, выполненных на разном количестве процессоров (от 1 до 16) при различной величине перекрытия, представлены на рис. 7. Результаты показывают, что время счёта колеблется от 36 с до 5 мин в зависимости от объёма пересылок. Рост последнего ведёт к увеличению времени счёта и относительного ускорения вычислений.

8. Примеры модельных расчётов

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

На рис. 8, а приведены графические изображения состояния области моделирования в последовательные моменты времени при горении по однородному слою при ветре вдоль направления оси V. На рис. 8, б представлено горение по неоднородному слою (два типа горючего) с препятствием в центре (которое определено на основе эллипса с использованием в качестве длины полуосей случайных чисел). Можно видеть, как пожар обходит несгораемое препятствие. Распространение пожара по двум слоям горючего (рис. 8, в) моделирует развитие верхового пожара, который представляет наибольшую опасность и моделирование которого наиболее сложно.

Рис. 8. Распространение процесса горения (примеры расчётов)

Заключение

При разработке параллельных алгоритмов моделирования динамики лесных пожаров на многопроцессорных вычислительных системах с распределённой памятью исследовались:

• способы геометрической декомпозиции расчётной области на подобласти [3];

• размерность расчётной сетки;

• величина перекрытия соседних подобластей.

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

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

Разработанный комплекс программ может быть использован ФГУ "Авиалесоохрана" и подразделениями МЧС РФ для прогнозирования массовых лесных пожаров.

Список литературы

[1] Ершов Д.В., Коровин Г.Н., Лупян Е.А. и др. Российская система спутникового мониторинга лесных пожаров // Современные проблемы дистанционного зондирования Земли из космоса: Физические основы, методы и технологии мониторинга окружающей среды, потенциально опасных объектов и явлений. Сб. науч. статей. М.: Полиграфсервис, 2004. С. 47-57.

[2] Доррер Г.А. Модель распространения фронта лесного пожара // Теплофизика лесных пожаров / Новосибирск: Ин-т теплофизики СО АН СССР, 1984. С. 86-98.

[3] Вдовенко М.С., Доррер Г.А. Об эффективности параллельных алгоритмов при моделировании динамики лесных пожаров // Вестник Сиб. гос. аэрокосм. ун-та. Математика, механика, информатика. 2009. Вып. 2(23). С. 152-156.

[4] Гергель В.П., Стронгин Р.Г. Основы параллельных вычислений для многопроцессорных вычислительных систем. http://www.software.unn.ac.ru/ccam/files/HTML Version/index.html

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

[6] MacDonald N., Minty E., Harding T., Brown S. Writing Message-Passing Parallel Programs with MPI. http://www.hpc.nw.ru/PS/mpi-course.zip

[7] McBryan O.A. An overview of message passing environments // Parallel Comput. 1994. Vol. 20. P. 417-441.

Поступила в 'редакцию 14 ноября 2012 г., с доработки — 3 декабря 2012 г.