CC BY
56
УДК 519.681.5 Е.А. Шельмина
Применение технологий параллельного программирования для решения обратных задач переноса примеси в атмосферном воздухе
Предложен алгоритм решения обратных задач переноса примеси в атмосферном воздухе с применением модифицированного метода Марчука. Для численного решения таких задач используются явные разностные схемы высокого порядка точности и технологии параллельного программирования. Для параллельной реализации численного решения обратной задачи переноса примеси используются функциональная, геометрическая декомпозиция и комбинация функциональной и геометрической декомпозиций.
Ключевые слова: обратные задачи, математическое моделирование, параллельная реализация.
К настоящему моменту разработано большое количество методов математического моделирования прямых и обратных задач математической физики, в том числе и задач переноса примеси. При использовании этих методов для постановки прямой задачи исследуемое физическое явление описывается дифференциальными уравнениями в частных производных с соответствующими начальными и граничными условиями, а решение поставленной задачи получают с использованием различных численных методов.
Решение обратных задач математической физики часто является элементом математического моделирования различных физических процессов. Это связано с тем, что далеко не все условия, описывающие процесс, заранее известны. На данный момент можно сказать, что методы решения обратных задач составляют важное направление исследований в области математической физики.
В настоящее время среди задач математической физики большое внимание уделяется задачам, которые описывают распространение примесей в атмосферном воздухе. Это связано с тем, что проблема охраны окружающей среды и ее восстановления является одной из важных задач науки, развитие которой стимулируется всевозрастающими темпами технического прогресса во всех странах мира. Поступление загрязняющих веществ в атмосферу избежать невозможно, но разумное использование природных ресурсов и постоянный контроль качества атмосферного воздуха позволят обеспечить безопасный уровень воздействия на атмосферу и избежать глобально негативных последствий. Большую помощь здесь может оказать применение методов математического моделирования, с помощью которых разрабатываются математические модели, позволяющие изучать атмосферные процессы. Основой этих моделей являются уравнения гидротермодинамики атмосферы, уравнения переноса примесей, уравнения, отвечающие за режим взаимодействия воздушных масс и примесей с неоднородной поверхностью Земли. Все эти уравнения записываются в дифференциальном виде и представляют собой законы сохранения энергии, массы и количества движения.
Следует отметить, что чаще всего проблема математического описания переноса примесей представляется двумя классами задач. Первый - это решение «прямых» задач, когда по известным параметрам источников примеси (месторасположение источника, время выброса) требуется найти значения концентраций выбрасываемой субстанции. Второй - решение «обратных» задач, когда известны значения концентрации примеси, измеренные в нескольких контрольных точках, требуется найти тип, координаты и мощность её источников.
Математическое описание «прямых» и «обратных» задач переноса примеси приводит к системам дифференциальных уравнений в частных производных. Но точное решение таких систем чаще всего невозможно получить в аналитическом виде, поэтому для решения таких задач используют численные методы. Но в силу сложных процессов распространения примеси такие задачи являются громоздкими и требовательными к вычислительным ресурсам. Перспективным способом решения этих проблем является использование современных высокопроизводительных многопроцессорных вычислительных систем, которые обеспечивают существенное ускорение получения результатов расчетов и повышение качества численного прогноза.
Целью данной работы является создание алгоритмов численного решения обратных задач переноса примеси с использованием явных разностных схем высокого порядка точности и высокопроизводительной вычислительной техники.
Физическая постановка прямой задачи переноса примеси. Распространение газообразной примеси в приземном слое атмосферы над промышленным центром и его окрестностями исследуется при следующих условиях. Рассматривается ограниченный в прямоугольнике участок территории, на котором в течение времени исследования происходит выброс газообразной примеси от различных источников. Мощность источников зависит от времени. Распределение по области исследования поступающей примеси зависит от метеорологических условий (поля ветра, осадков, интенсивности турбулентного обмена) и от характера взаимодействия примеси с подстилающей поверхностью. Предполагается, что эти необходимые для замыкания модели параметры известны. Кроме того, будем считать, что газообразная примесь инертна, т.е. не вступает в химические реакции с другими веществами и скорость ее переноса совпадает со скоростью движения окружающего воздуха, её молекулярный вес близок к молекулярному весу воздуха и примесь изотермична.
Математическая постановка прямой задачи переноса примеси. С учетом принятой физической постановки задачи адвективно-диффузионное уравнение, моделирующее перенос газообразной примеси в заданном потоке, представляется в следующем виде [1]:
” _ с
дх _
где С - концентрация примеси; и,У,№ - компоненты соленоидального вектора скорости атмо-
ди д¥ дМ
сферного воздуха, т.е.----1----1----= 0 ; Г,Кх - коэффициенты турбулентной диффузии; Q - ин-
дх ду дх
тенсивность поступления примеси от М рассматриваемых высотных источников,
М
Q =^ Qj (?)8(х - хо )8(у - уо )8(х - х0); Qj (?) - расход выброса из источника с координатами
] =1
дС дС дС дС д
—+и—+V—+М—+стС=— д? дх ду дх дх
гдС
дх
д
+
ду
гдС
ду
д
+
дх
Кх~
+Q, (1)
х
-0 ,.0 „0
],
у°,х] ; ст - интенсивность влажного осаждения (вымывания осадками) примеси.
Начальные и граничные условия для уравнения (1) имеют следующий вид:
? = 0: С (0, х, у, х) = С0( х, у, х);
0 дС 0 Т дС 0
х = 0: — = 0; х = Ьу.:— = 0;
дх х дх
0 дС 0 Т дС 0
у=0:^-=0; у=Ту :^~=0;
ду ^ ду
х = 0:Кх — = аС - Я; х = Тх :— = 0. (2)
дС С п Т дС — = аС - Я; х = Ь- : —
дх х дх
Я = Я(?,х,у) - интенсивность поступления газообразной примеси от наземных источников; а -скорость сухого осаждения газообразной примеси.
Физическая постановка обратной задачи переноса примеси. Требуется по известным метеорологическим параметрам атмосферы и результатам измерений концентрации примеси в N точках (Ы — количество постов наблюдения, на которых осуществляется измерение концентрации примеси), проводимых в течение некоторого временного отрезка Т, определить параметры (мощность, координаты и время срабатывания) источников примеси.
Математическая постановка обратной задачи переноса примеси. Чтобы записать математическую постановку обратной задачи переноса примеси, используем метод, который основан на решении уравнения, сопряженного с полуэмпирическим уравнением турбулентной диффузии (1), и двойственном представлении функционала от концентрации примеси [1]. Для этого адвективно-
*
диффузионное уравнение переноса примеси (1) умножается на сопряженные функции С£ (?,х,у,х)
* 2
{функции С£ еС [0, Ьх ]х[0,Ьу ]х[0, Ьх ] и определены в области исследования, £ = 1, N, N - количе -ство измерений концентрации С(?,х,у,х) в точках с координатами (х£,у£,х£) в момент времени ?£ }
и интегрируется по времени и пространству. После некоторых математических преобразований получаются сопряженные уравнения [1]:
дс* диС* дУС* дж* +аС*
ді дх ду ді С дх
начальные и граничные условия:
$ гдС* д 1 * С д І-Ч 1 д * [к С 1
дх ду ду д 2 де
= Рк; к = 1,...,#, (3)
С*(Т, х, у, і) = 0;
х = 0:иС* +ГС = 0; х = Ьх :иС* +ГС = 0; к дх х к дх
у = 0:УС* +ГС = 0; у = Ьу :УС* +ГС = 0;
і = 0:К
дС*
ду
= аС*; і = Ь
ду
.дС*
ді 1 ді
где Рк =8(х_хк)8(у_ус)8(і_ ік)8(і_ік); к = 1,...,#; 8(£) - дельта-функция (^єЖ ); N - количес-тво измерений концентрации С (і, х, у, і) в точках с координатами (хс , ук, ік) в момент времени ік. Здесь предполагаем, что ік >> 0 настолько, что справедливым является допущение
С*(0, х, у, і) = 0 (5)
С учетом (3)-(5) и осуществленных математических преобразований, в ходе получения сопряженных уравнений получаем
ТЬхЬуЬі ТЬхЬуЬі ТЬхЬу
= 0,
(4)
ПИ CPkdzdydxdі = ІШ С* Qdzdydxdі + I С*
0 0 0 0 0 0 0 0 0 0 0
і=0
Rdydxdі; к = 1,..., N.
где С* - решение к -й сопряженной задачи (3)-(4) с правой частью Рк . Вследствие принятого вида функции Рк последнее равенство перепишется как
ТЬхЬуЬ2 ТЬХ^у
\\\\ C*Qdzdydxdt + и I С* Rdydxdt = Ск ; к = 1,...,# (6)
0 0 0 0 0 0 0 z=0
и может рассматриваться в качестве двойственного представления функционала от концентрации примеси или некоторых дополнительных соотношений для оценки по измеренным значениям концентрации примеси Ск в точке с координатами (Хк,Ук,zk) в момент времени tk распределения и интенсивности источников выброса примеси в приземный слой воздуха. При этом нахождение решения системы уравнений (6) позволяет установить параметры источников загрязнения.
Таким образом, применение метода сопряженных уравнений позволяет построить алгоритм решения обратной задачи по нахождению параметров точечных источников по данным измерений концентрации примеси, полученных в ряде контрольных точек.
Численное решение. В данной работе для численной реализации обратной задачи переноса примеси (3)-(4), (6) использовались метод конечного объёма и явные разностные схемы высокого порядка точности: схема Ботта [2] и схема БКО [3] для аппроксимации адвективных членов. Тестирование разностных схем с целью определения их точности производилось на двумерном уравнении переноса примеси с известным аналитическим решением [4]. При этом было получено, что используемые разностные схемы высокого порядка точности дают результаты, достаточно близкие к аналитическому решению для случая постоянных значений скорости и коэффициента диффузии, что подтверждает теоретическое исследование их порядка точности.
При численном решении обратной задачи для проверки правильности её решения использовались результаты решения прямой задачи в качестве необходимых входных данных. Необходимые для расчетов параметры (скорость, направление приземного ветра, турбулентные характеристики атмосферного пограничного слоя) рассчитывались с использованием мезомасштабной модели ТГУ -ИОА СО РАН [5]. В качестве входных параметров использовались результаты измерения монооксида углерода на пяти постах наблюдения, расположенных в г. Томске. При моделировании рассматривался период времени, охватывающий двое суток (7-8 сентября 2000 г., Т = 172800 с), следова-
тельно, количество измерений {Ck,к = 1,...,N}, необходимых для решения обратной задачи, составляет N = 30. Область исследования имеет площадь размером 50х50 км, в ее центральной части находится городская территория. Верхняя граница расчетной области располагается на высоте 1000 м.
Параллельная реализация. При численном решении задач переноса примеси выделяют два основных способа параллельной реализации: распараллеливание по физическим процессам (функциональная декомпозиция) и геометрическая декомпозиция расчетной области.
Функциональная декомпозиция. Распараллеливание по физическим процессам метода численного решения задачи по определению параметров источников производилось с использованием принципа «master-slave». В этом случае управляющий master-процесс с частотой а передает каждому slave-процессу значения метеорологических параметров, необходимые для решения сопряженных задач. Подчиненные slave-процессы, получив данные, в свою очередь, ведут расчеты независимо друг от друга, и найденные на каждом шаге по времени приближенные решения своей сопряженной задачи (3)-(4) возвращают управляющему процессу, который ищет глобальный минимум функционала (6), что в итоге позволит определить параметры источников загрязнения.
В работе на основе теоретических оценок и в результате вычислительных экспериментов исследованы ускорение и эффективность данного параллельного алгоритма для задачи идентификации параметров мгновенного источника. Теоретические выкладки приводят к следующим результатам:
еФ1 ~ ^ ' 17Ф1 — 1* ~ " Г (П\
^Р ~ С , ЕР = ~ С , (7)
1 + Р 1+ 5%-ю
т т
Здесь р - число используемых процессов; т - число арифметических операций для вычисления одного значения сеточной функции; % - отношение времени межпроцессорной передачи одного
числа к времени арифметической операции; Ф1 - первый вариант функциональной декомпозиции.
Из (7) видно, что при таком способе организации параллельных вычислений не удается достичь идеального параллелизма, поскольку нет хорошей балансировки загрузки используемых в расчетах р процессоров. Это связано с тем, что master-процесс не участвует в проведении массовых вычислений. Поэтому с целью повышения эффективности алгоритма рассматривается другой подход, в котором master-процесс дополнительно к своей работе также решает одну из N сопряженных задач (3)—(4). В этом случае теоретическая оценка ускорения и эффективности может быть записана следующим образом:
Р-1 77Ф1 Sp 1-1/Р
8ф2 «_Р_, Еф2 — (8)
1+ 1+ 5%'ю т т
и разработанная параллельная программа будет иметь более высокие показатели по времени выполнения (Ф 2 - второй вариант функциональной декомпозиции).
Для подтверждения полученных оценок ускорения параллельного алгоритма был проведен ряд расчетов на кластере ТГУ СКИФ СуЬепа. Результаты представлены в табл. 1.
Таблица 1
Время счета (сек) параллельной программы для рассматриваемых вариантов функциональной
декомпозиции_______________________________________
1
Количество измерений,N Число используемых процессов/время счета для варианта Ф1 Число используемых процессов/время счета для варианта Ф2
5 6/312 5/323
10 11/316 10/337
20 21/322 20/354
Из таблицы видно, что оба рассматриваемых варианта имеют высокие показатели масштабируемости, что выражается в слабом изменении времени счета при добавлении процессорных элементов в связи с увеличением объема вычислений (N). Однако первый вариант функциональной декомпозиции алгоритма решения обратной задачи (3)—(4), (6) имеет небольшое преимущество во времени выполнения параллельной программы, что обусловлено появлением дополнительных возможностей совмещения работы активных процессов. В то же время соотношение эффективности
рассматриваемых параллельных программ еф =—Ф- <1 показывает лучшую масштабируемость
E ф 2
ЕР
второго варианта декомпозиции, что хорошо согласуется с теоретическими оценками (7)-(8), со-
î Ф\ Р
гласно которым (е )теор «------.
F p + 1
Геометрическая декомпозиция. При решении каждой сопряженной задачи (3)-(4) в данной работе также использовалась одномерная (вдоль оси OY) геометрическая декомпозиция сеточной области. В этом случае вычисления осуществляются следующим образом: N сопряженных задач (3)-(4) последовательно одна за другой решаются численно в подобласти вычислительной сетки, распределенной каждому процессу. Вычисления в течение промежутка времени [0,7’] проводятся
*
каждым процессом одновременно, но для корректного расчета сеточных значений функции Ck вблизи границ подобласти требуются некоторые сеточные значения с соседних по расположению подобластей. Это может быть обеспечено лишь путём межпроцессорной передачи данных с использованием, например, стандарта MPI.
Теоретическая оценка ускорения и эффективности данного алгоритма имеет следующий вид:
ЧГ і
S Г________Р E Г = Sp_________1 (9)
SP ~ Ґґ л V ЕР = ~ f с л v (9)
1+щ 5^+.4^ р 1+ц 5^+.4^
m ^ N Nx j m { N Nx j
где Nx - размер вычислительной сетки по оси Ox ; N - количество решаемых сопряженных задач.
Для подтверждения проведенных теоретических оценок были проведены вычислительные эксперименты, результаты которых представлены в табл. 2.
Таблица 2
Время счета параллельной программы для геометрической декомпозиции__________
Количество измерений, N Число используемых процессов, p Время, с
5 5 627
10 10 643
20 20 743
Сопоставляя результаты, приведенные в табл. 1 и 2, видно, что применение геометрической декомпозиции для рассматриваемой задачи идентификации параметров внезапного выброса в атмо -сферу по данным измерений существенно (более чем в два раза) уступает по быстродействию функциональной декомпозиции. Причиной этого являются коммуникационные затраты на межпроцессорную передачу данных, необходимую для корректного решения сопряженных задач (3) - (4) на всей сеточной области, разделенной между активными процессами.
Комбинированный способ распараллеливания. Рассмотренные выше подходы создания параллельных версий алгоритмов решения обратной задачи переноса примеси показали неплохие результаты по эффективности, однако количество используемых активных процессов в них ограничено: в случае функциональной декомпозиции числом проведенных измерений N, а при одномерной геометрической декомпозиции размером сетки Nx и выбранным сеточным шаблоном.
При использовании комбинированного подхода параллельной реализации предлагается для увеличения в расчетах количества используемых активных процессов совместить применение функциональной и геометрической декомпозиции. Все р активных процессов делятся на N групп (по числу измерений), каждая из которых решает численно одну сопряженную задачу (3)-(4), используя одномерную геометрическую декомпозицию. Для получения приближенного решения процессы отдельной группы обмениваются между собой рассчитанными значениями сеточной функции вблизи границ подобластей. Поиск минимума функционала выполняет один из р процессов, которому остальные передают результаты промежуточных расчетов.
Теоретические оценки такого способа параллельной реализации записываются в следующем виде:
комб
1+-Х
m
5ю +-
4 p
N ■ N.
^комб = _p p
комб
x У
1+Х
m
5ю +-
4 p
N ■ N.
N
= 12 3
(10)
x У
Сравнивая оценки параллельных реализаций (9) и (10), можно отметить, что эффективность комбинированного подхода не хуже, чем при геометрической декомпозиции, однако количество используемых активных процессов можно увеличить на порядок, что при имеющей место масштабируемости параллельной версии алгоритма должно привести к значительному сокращению времени счета.
Результаты вычислительных экспериментов для задачи идентификации характеристик источника примеси подтвердили предварительные теоретические оценки (табл. 3).
Таблица 3
Время счета параллельной программы для комбинации функциональной
Количество измерений, N Время счета последовательной программы, с Число используемых процессов / Время счета для комбинированного подхода, с
5 1626 100/33
10 3922 200/42
20 7799 400/59
Стоит отметить, что результаты численных расчетов при использовании комбинированного подхода показывают снижение времени счета с уменьшением вычислительной трудоемкости задачи, замедление ускорения параллельной программы с ростом числа используемых процессов и при сокращении количества сопряженных задач N. Комбинированный подход действительно на порядок позволяет увеличить число используемых процессов, при этом эффективность параллельной программы остается на уровне не менее 50%, что является неплохим для задач такого класса показателем. Однако, в целом функциональная декомпозиция имеет существенное преимущество, судя по (8)-(10).
Заключение. В данной работе приведены математические постановки прямой и обратной задач переноса примеси в атмосферном воздухе. Математические постановки обратных задач опираются на аппарат сопряженных уравнений и двойственное представление функционала от концентрации примеси. Для численного решения сопряженных задач используются явные разностные схемы высокого порядка точности. Для ускорения проведения вычислительных расчетов была разработана параллельная программа, реализующая комбинацию функциональной и геометрической декомпозиции.
Работа выполнена при поддержке гранта РФФИ № 12-01-31050 мол_а.
Литература
1. Марчук ГИ. Математическое моделирование в проблеме окружающей среды. - М.: Наука, 1982. - 315 с.
2. Bott A. Monotone Flux Limitation in the Area - preserving Flux - form Advection Algorithm // J. Monthly Weather Review. - 1992. - Vol. 120. - P. 2592-2602.
3. Chi-Wang Shu. Essentially non-oscillatory schemes for hyperbolic conservation laws // Preprint of Division of Applied Mathematics. -Brown University. - 1996. - 92 p. [Электронный реcурс]. - Режим доступа: http://www.ipam.ucla.edu/publications/pcatut2005/pcatut_5477_preprint.pdf, свободный (дата обращения: 24.02.2013).
4. Панасенко Е.А. Определение городских районов-загрязнителей атмосферного воздуха по данным наблюдений / Е.А. Панасенко, А.В. Старченко // Оптика атмосферы и океана. - 2009. - Т. 22, № 3. - С. 279-283.
5. Старченко А.В. Численное исследование локальных атмосферных процессов // Вычислительные технологии. - 2005. - №10, ч. 2. - С. 81-89 .
Шельмина (Панасенко) Елена Александровна
Канд. физ.-мат. наук, доцент каф. экономической математики информатики и статистики ТУСУРа
Тел.: 8 (383-2) 90-01-87
Эл. почта: [email protected]
Shelmina (Panasenko) E.A.
The numerical solution of inverse problems of pollutant transport using a finite difference schemes of high order and supercomputers
In the paper we suggest an algorithm for solving inverse problems of pollutant transport using a modified method of Marchuk. Independent conjugate equations are integrated for solving the inverse problem at each time step. Such conditions of numerical investigation allow to attract high-performance computing. There are the following ways for the parallel implementation of numerical solution of pollutant transport: geometrical decomposition, functional decomposition, the combination of functional and geometric decomposition.
Keywords: inverse problems, mathematical modeling, parallel implementation.
