36

Размещено на http://www.allbest.ru/

АВТОРЕФЕРАТ

диссертации на соискание ученой степени

доктора технических наук

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

Специальность: 05.13.18 - Математическое моделирование,

численные методы и комплексы программ

Пекунов Владимир Викторович

Иваново - 2009

Работа выполнена на кафедре "Высокопроизводительные вычислительные системы" Государственного образовательного учреждения высшего профессионального образования "Ивановский государственный энергетический университет имени В.И. Ленина"

Научный консультант:

доктор физико-математических наук, профессор Ясинский Федор Николаевич

Официальные оппоненты:

доктор технических наук, профессор Шувалов Сергей Ильич

доктор физико-математических наук Зибаров Алексей Владимирович

доктор физико-математических наук, профессор Григорьев Александр Иванович

Ведущая организация:

Институт прикладной математики

имени М.В. Келдыша РАН

Защита состоится "5" февраля 2010 г. в 11 часов в ауд. Б-237 на заседании диссертационного совета Д 212.064.03 при Ивановском государственном энергетическом университете.

С диссертацией можно ознакомиться в библиотеке ИГЭУ,

Автореферат разослан "___"_________ 2009 г.

Ученый секретарь диссертационного совета кандидат технических наук, доцент Шульпин А.

Общая характеристика работы

Актуальность темы

Существует множество экологических проблем, обусловленных загрязнением воздушной среды в окрестностях промышленных предприятий, а также в больших городах. Крупными источниками твердых и газообразных загрязнителей являются ТЭЦ и ТЭС, котельные установки, предприятия черной и цветной металлургии, транспорт. Поглощение загрязнителей (преимущественно оксидов азота и серы) каплями воды приводит к образованию кислотных дождей. Неблагоприятная экологическая ситуация часто является следствием принятия проектных решений, плохо обоснованных с точки зрения экологической безопасности. Проблема, обычно, заключается в использовании недостаточно точных оценок предполагаемых уровней загрязнений при анализе проектных вариантов. Ряд исследований подтверждает, что применяемые нормативные методики (ОНД-86 и другие) недостаточно учитывают ряд важных факторов, например, турбулентность и химическую кинетику, воздействие солнечного излучения на тепловые и фотохимические процессы, поглощение загрязнителей каплями воды.

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

Универсальные пакеты моделирования процессов аэрогидродинамики (FlowVision, FLUENT, GAS DYNAMICS TOOL, PHOENICS, Star-CD и многие другие) позволяют учесть все необходимые факторы и провести расчет требуемой сложности на многопроцессорных системах. Однако универсальность таких пакетов, фактически, заключается в попытках применения достаточно ограниченного набора моделей, алгоритмов и методов к множеству различных случаев. Очевидно, что программы, ориентированные на решение конкретного класса задач, потенциально способны решать данные задачи более эффективно.

Однако большинство известных автору отечественных и зарубежных специализированных программных продуктов (например, ADAM, CAL3QHC, Chensi, TASCflow, ISC_3, PANACHE, REMSAD, UAM_IV, ЭКОЛОГ, ПРИЗМА, VITECON, геоинформационные системы), предназначенных для расчета распространения загрязнений, ориентировано на однопроцессорные системы. В специализированных программных комплексах (например, ECOSIM и MAQSIP), адаптированных для работы на многопроцессорных системах, распараллелены лишь отдельные модули, что в ряде случаев не позволяет добиться высокой эффективности распараллеливания. Добавим, что в специализированных однопроцессорных разработках почти не рассматривается поглощение загрязнителей каплями воды (за исключением PANACHE), а в многопроцессорных - не учитывается перенос излучения и его влияние на тепловые и фотохимические процессы, что существенно снижает точность результатов.

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

Цель работы

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

Задачи исследования

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

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

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

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

Научная новизна

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

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

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

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

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

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

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

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

Практическая ценность

1. На базе предложенных модели и алгоритмов разработаны однопроцессорная и многопроцессорная версии программы моделирования загрязнений AirEcology-P. Многопроцессорный вариант оптимизирован для работы с МВС_1000, МВС-100, Power XPlorer, а также с кластерными и блочно-кластерными (с многоядерными блоками) системами на платформах Windows и Unix/Linux. Применение разработанных программ позволяет повысить качество решений при проектировании промышленных предприятий и объектов городской застройки за счет более точной и эффективной оценки предполагаемых уровней загрязнений.

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

3. На базе предложенных объектно-событийных моделей разработана система автоматизированного порождения программ PGEN++. Система используется для оперативной модификации и частичной верификации математической модели в программе AirEcology-P и для автоматизации разработки тестов в системе профильного тестирования ПРОФТЕСТ (компонент комплекса ГИПЕРТЕСТ).

4. Получены результаты параллельного численного моделирования процессов образования и распространения смога на улицах города, очистки атмосферы с помощью поливочных машин и распространения промышленных выбросов в окрестности предприятия. Даны рекомендации по применению результатов диссертации в САПР промышленных предприятий и САПР объектов городской застройки. Приведен пример итерационной процедуры поиска оптимального места размещения ТЭС с помощью разработанной подсистемы автоматизированной оценки предполагаемых уровней загрязнений.

Достоверность полученных результатов подтверждается итогами сравнения реальных данных с результатами, полученными в ходе численных экспериментов с применением предложенных модели, схем, алгоритмов, численных методов. Сравнение результатов моделирования переноса пассивной примеси, полученных с применением разработанной автором системы (AirEcology-P) и иных систем, показало, что по данному параметру AirEcology-P дает один из наиболее точных результатов.

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

Внедрение

Программа моделирования загрязнений AirEcology-P внедрена в учебный процесс кафедры "Безопасность жизнедеятельности" ИГЭУ по дисциплинам, предполагающим проведение экологических экспертиз;

в ЗАО "ВВП" (г. Нижний Новгород) для оценки состояния окружающей среды в районе предприятия в зависимости от режима его работы и климатических условий.

Имитаторы многопроцессорной системы Power XPlorer внедрены в учебный процесс кафедры "Высокопроизводительные вычислительные системы" ИГЭУ для обеспечения лабораторного практикума по дисциплине "Параллельные вычисления".

Система автоматизированного порождения программ в составе системы профильного тестирования ПРОФТЕСТ (компонент комплекса ГИПЕРТЕСТ) внедрена в ООО "ЛУКОЙЛ-Коми" (г. Усинск) для автоматизации разработки сценариев тестирования.

На защиту выносятся:

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

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

3. Быстрая и полная схемы выбора метода распараллеливания.

4. Методы и алгоритмы сокращения и оптимизации пересылок.

5. Алгоритм балансировки загрузки процессоров.

6. Объектно-событийные модели порождения программ.

Апробация работы

Материалы диссертации докладывались на следующих международных, всероссийских и межвузовских конференциях и семинарах:

II Международном симпозиуме "Математическое моделирование экологических процессов" (Иваново, 2000);

II Всероссийской научной конференции "Высокопроизводительные вычисления и их приложения" (Черноголовка, 2000);

Международной научно-технической конференции "Состояние и перспективы развития электротехнологий" (Иваново, 2001, 2003, 2005, 2007);

VIII Всероссийском съезде по теоретической и прикладной механике (Пермь, 2001);

IV Всероссийской научной internet-конференции "Компьютерное и математическое моделирование в естественных и технических науках" (Тамбов, 2002);

Всероссийской научно-технической конференции "Современные технологии и оборудование текстильной промышленности" (Москва, 2004);

Межвузовской научно-практической конференции "Экологические проблемы Ивановской области" (Иваново, 2005);

семинарах в Институте прикладной математики имени М.В. Келдыша, РАН (23.05.2002; 26.06.2003; 24.03.2005; 10.04.2008);

Межвузовском научно-техническом семинаре "Математическое моделирование и многопроцессорная вычислительная техника - 2006" (Иваново, 2006);

Региональной научно-технической конференции "Применение многопроцессорных суперкомпьютеров в исследованиях, наукоемких технологиях и учебной работе" (Иваново, 2007, 2008);

XXVII Российской школе по проблемам науки и технологий (Миасс, 2007);

семинаре в Нижегородском гос. университете (26 июня 2008);

семинаре в Институте вычислительной математики РАН (18 сентября 2008);

семинаре в Институте математического моделирования РАН (9 октября 2008).

Разработанные программа моделирования загрязнений AirEcology-P и система автоматизированного порождения программ PGEN++ зарегистрированы в Российском агентстве по патентам и товарным знакам (Роспатент), (свидетельства об официальной регистрации программ для ЭВМ №2006611068 и №2007614631).

Разработки, выполненные в рамках работы над диссертацией, были представлены на выставке авторских разработок "Зворыкинского проекта" в рамках международного форума "Дни Русских Инноваций" (Москва, 2009).

Работа выполнена при финансовой поддержке Минобразования и науки (грант РНП.2.2.1.1.7280) в совместной работе с Региональным научно-образовательным центром "Жидкие кристаллы" (г. Иваново).

По материалам диссертации опубликованы 47 печатных работ, в том числе монография, глава в коллективной монографии, 3 учебных пособия, 14 статей в ведущих рецензируемых научных журналах, перечень которых определен Высшей аттестационной комиссией, 5 статей в сборниках трудов вузов, 1 препринт, 17 тезисов докладов на конференциях и другие работы.

Структура работы

Диссертация состоит из введения, четырех глав, заключения, библиографического списка из 287 наименований и включает 274 страницы (из них 29 страниц библиографии), в том числе 69 рисунков,20 таблиц, а также приложение на 134 страницах.

Содержание работы

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

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

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

параллельное моделирование загрязнение воздушное

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

; j = 1, 2, 3;

; ; ;

;

,

где U - вектор скорости основной фазы, P - давление, _мол и _турб - молекулярная и турбулентная вязкости, T - температура, ₀, c, , b - теплопроводность, удельная теплоемкость, плотность и термический коэффициент расширения воздуха, - коэффициент, x_i, x_j, t - переменные Эйлера. Коэффициенты отражают влияние турбулентности на прочие процессы. Модель турбулентности K-E (RNG):

; ;

; ;

; ; ,

включая стандартные параметры: _K = 0,7194, = 0,7194, c₁ = 1,42, c₂₁ = 1,68, ₀ = 4,38, = 0,012, c = 0,0845. Модель K-E (RNG) показала хорошую точность и наибольшую адекватность передачи структуры потока в численных экспериментах с различными моделями (Абрамовича-Секундова, Спаларта-Аллмараса, четырьмя вариантами K-E: стандартным, Yap, Chen-Kim, RNG) при сравнении с опытными данными (см. приложение 1) для стандартной задачи об обтекании препятствия CEDVAL A-1.

;

;

; ,

где C - вектор концентраций веществ (газов и переносимых воздухом пылевых частиц), W_{j, -} скорость витания и коэффициент диффузии j_го вещества, N - число веществ, q - число химических реакций, R_ik и L_ik - стехиометрические коэффициенты при i-м веществе в правой и левой частях k-й реакции, A_k = A_k (T) - константа скорости k-й реакции. Для фотохимических реакций A_k 0 лишь в случае, если суммарная освещенность превышает заданный стартовый порог.

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

;

,

где - коэффициент диффузии пара, M_k - молярная масса воды. Уравнения переноса диффузного излучения в N_rad диапазонах (в последнем - теплового инфракрасного излучения) по модели K. N. Liou:

; ;

; ; ;

где - первый компонент разложения интенсивности излучения в j-м диапазоне (метод сферических гармоник), - интегральная освещенность (прямым солнечным излучением), и - начальная и конечная длины волн, B (T) - функция Планка, , , - коэффициенты ослабления, поглощения и рассеивания, ^j - альбедо однократного рассеивания, - фактор асимметрии. находится обратной трассировкой луча, что упрощает вычисления, особенно в областях сложной формы.

Из каждой ячейки трассируется несколько лучей (в направлении, обратном вектору падения солнечного излучения) до ячейки с известной величиной , для каждого луча рассчитывается ослабление по закону Бугера-Ламберта, результаты усредняются. Это повышает точность расчета для частично затененных ячеек. Расчет распараллелен (OpenMP) по лучам, в пределах блока области.

Уравнения для фазы тяжелых пылевых частиц:

; j = 1, 2, 3;

,

где U_p - вектор скорости пылевой фазы, _pмол - кинематическая вязкость пылевой фазы, g₁ = 0, g₂ = 0, g₃ = g, - сила сопротивления частицы потоку, _p - плотность пыли, D и - коэффициенты, введенные с целью придать единый вид параболическим уравнениям (в простейшем случае D = 1, = 0). Аналогичные коэффициенты в уравнениях для капельной фазы имеют тот же смысл по умолчанию.

Граничные условия. Для большинства переменных ставятся обычные условия первого/второго рода. Для энергии пульсаций K - условия второго рода, для скорости диссипации - первого рода: , где y - расстояние до стенки. Для диффузной интенсивности излучения () на открытых границах задается равенство потока освещенности извне известной величине . На закрытых границах для теплового инфракрасного излучения:

,

где _s - коэффициент черноты поверхности стенки, T_s - температура стенки, f_T (T_s) - интеграл функции Планка по длинам волн диапазона излучения. В иных диапазонах излучения:

,

где - альбедо стенки, ₀ - вектор падения лучей. Для температуры в верхнем слое стенки Ox_i решается двумерное уравнение теплопроводности с учетом кондуктивной и лучистой теплопередач:

,

где , c_s, _s, _s - температуропроводность, теплоемкость, плотность и коэффициент теплоотдачи материала стенки, - размер j-й ячейки стенки по Ox_i; F_s - поток лучистой энергии между стенкой и средой.

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

Предложена модель фазы капель из Z компонентов (то есть секций, к i-й секции относятся капли с диаметрами [d_i; d_i+1]):

;

;

; ; ;

,

межфазные () и межкомпонентные () переходы:

;

; ; ; ;

; ; ;

; ;

;

;

; ; ; ,

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

где , , - плотность, концентрация и вектор скорости i-го компонента, - концентрация j-го газа в каплях i-го компонента, , , и - коэффициенты, - сила сопротивления потоку с учетом деформации, - шаг интегрирования, L_i (D) и - потоки пара и j-го газа между i-м компонентом и средой, - плотность вещества i-го компонента, - концентрация пара на поверхности капли (с поправками Кельвина и Кехлера).

Новизна модели состоит в расчете параметров a_i и b_i кусочных распределений (линейных при полном заполнении [капли с диаметрами от d_i до d_i+1] и равномерных при частичном заполнении [с поиском начала x_i или конца y_i заполнения]) путем восстановления по значениям и . Сформулированы простые критерии для определения типа процесса в секции и ее заполненности. При полном заполнении решается система линейных уравнений относительно a_i и b_i:

а при частичном заполнении решается одно из двух уравнений:

; ,

соответственно, или ; причем , .

Подход имеет невысокую трудоемкость и особенно эффективен при активном переносе капель между ячейками со значительно варьирующимися и нестандартными распределениями капель. Такие случаи могут некорректно обрабатываться при использовании традиционных подходов (А.Е. Алояна, метода мод и других). Сравнение результатов моделирования популяции капель с применением разработанного подхода и лагранжевой модели показало, что при низких вычислительных трудозатратах (в 540 раз меньше по сравнению с лагранжевой моделью на 1000 капель) предложенный подход показал приемлемые величины относительной погрешности (до 218% по плотности капель и 1735% по загрязнителям), которая остается ограниченной. В практически важных случаях точность можно повысить, измельчив секции.

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

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

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

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

Представим подсистему уравнений химической кинетики в виде:

,

что соответствует подходу Рожкова. Перейдем к разностной форме:

; ,

h - шаг интегрирования по времени. Применив метод Адамса-Моултона и группируя выражения относительно , получим расчетную формулу компромиссной разностной схемы:

,

где - параметры метода Адамса-Моултона n-го порядка (n < 5). Экспериментально показано, что при высокой жесткости системы уравнений (например, при горении метана) совершенно устойчив расчет с периодическим изменением порядка метода: n = 1+ [s mod N]; N - максимальный порядок метода, N < 5.

При расчете с простым контролем погрешности (циклическим делением шага пополам, пока погрешность не достигнет заданного уровня) компромиссная схема с N = 4 при точности, сравнимой с точностью метода Гира, показала меньшую (приблизительно в 1,22 раза) трудоемкость при небольшом исходном шаге h и равную трудоемкость при большой исходной величине h.

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

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

а) если "вычислительная жесткость" [, где Q_j - трудоемкость интегрирования j_го уравнения, измеряемая в количестве арифметических операций] высока, то выбирается РПП, поскольку отношение времен выполнения итерации для методов РПФ и РПП

будет существенно больше единицы. Здесь N_ур - число уравнений [используется N_ур процессоров].

б) иначе анализируется структура орграфа взаимозависимостей переменных подсистемы динамических уравнений. Если есть подграф, являющийся пустым, или содержащий незначительное число дуг и не связанный с остальной частью графа, то для соответствующих переменных используется РПФ, для остальных - РПП.

Если такого подграфа нет, то оценивается отношение максимальных объемов пересылок для РПФ и РПП, приходящихся на долю одного процессора, в системах с индивидуальными каналами связи (коммутаторами):

; ,

где и - количества обменов данными на итерации при решении всех уравнений в блоке области (РПП) и при решении каждого уравнения во всей области (РПФ); P_i - степень i-й вершины орграфа взаимозависимостей; N_слоев - число слоев узлов, передаваемых за одну пересылку (РПП); N_z - число слоев узлов по оси, по которой разбивается область (РПП); M_i {0; 1} - флаг необходимости значения i-й переменной для предварительных вычислений. В системах с разделяемой коммуникационной средой (с общим коммутатором) оценивается отношение общих объемов пересылок для РПФ и РПП:

,

где N_П - число процессоров. При G₂ < 1 выбирается РПФ, иначе - РПП.

По предложенной методике был сделан предварительный выбор метода РПП.

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

Полная схема выбора оптимальной комбинации методов распараллеливания (РПФ, РПП, комбинированный) предназначена для систем с неоднородной средой передачи данных (например, кластерных систем из SMP-блоков, обменивающихся данными по коммуникационной сети). Минимизируются общие временные затраты на вычисления (то есть максимизируется полная полезная нагрузка процессоров) и пересылки данных. Соответствующая модель целочисленной нелинейной оптимизации с ограничениями содержит следующие переменные:

X_kj {0; 1}, выражающие факт отображения всего поля k-й переменной на j-й процессор (РПФ), , ;

F_k {0; 1}, определяющие РПП для поля k-й переменной;

G_j {0; 1}, определяющие участие j-го процессора в обработке переменных с применением РПП.

Целевая функция:

;

; ;

,

где

V - объем поля (в байтах) одной переменной, _r - скорость обработки (байт/с) поля r-й переменной, V⁺ - объем дополнительных данных, обрабатываемых процессором при использовании РПП, A_rf - матрица отношения зависимости f-й переменной от r-й переменной, Q_ij - матрица значений времени (с/байт) пересылки одного байта данных из i-го в j-й процессор (латентностью пренебрегаем, поскольку она существенно меньше времени передачи), - объем данных поля r-й переменной, пересылаемых между двумя соседними процессорами (в одном направлении) при использовании РПП.

Первый член целевой функции определяет общее время параллельного расчета полей переменных. Второй член выражает суммарное время передач данных. Функция определяет время передачи (в группе РПП-процессоров) в случае, когда для обработки r-й переменной используется РПП. Время передачи для прочих случаев (РПФ и РПФ, РПФ и РПП) вычисляется функцией .

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

; ;

; ;

.

Применение данной модели подтвердило обоснованность выбора, сделанного в пользу РПП для решаемой задачи при работе с 128-ядерной кластерной системой ИГЭУ (четырехъядерные блоки, связанные сетью Gigabit Ethernet) и системой МВС-1000/16 (Fast Ethernet). Комбинации РПФ/РПП эффективны при решении ряда задач специфической структуры (например, при наличии одного или нескольких обособленных или слабосвязанных с остальными уравнений).

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

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

Подстройка заключается в том, что для случайного набора из S узлов относительно q_i (t) градиентным методом решается задача: , где M - глубина предсказания, _jk<q_i (t) > - невязка в j-м узле для k-й точки, полученная при текущих значениях q_i (t). Для решаемой задачи подстройка снижает погрешность экстраполяции в 1,21,3 раза на наиболее критичных первых сотнях итераций, далее подстройка дает снижение не более чем в 1,05 раза.

Модель построения оптимальной асинхронной схемы обменов данными на стыках блоков. Совместим во времени вычисления и обмены данными на стыках. Расчет представим как серию процессов обработки групп переменных модели.

Совмещение вычислений и обменов ограничивается наличием зависимостей между группами.

Редуцируем орграф зависимостей переменных модели, объединяя узлы, образующие циклы. Получим сеть взаимозависимостей К групп переменных, по которой построим матрицу P отношений строгого следования между группами. Задача сводится к определению значений целочисленных переменных X_ij (i-я группа входит в блок групп, передаваемый j-м по счету и обрабатываемый (j+1) - м по счету; X_ij {0; 1}), при которых достигается минимум временных затрат:

;

; ,

где ; M - максимальное количество блоков групп; - табличная нелинейная функция времени обмена сообщениями (на j-й стадии) из N элементов; k - коэффициент замедления вычислений при асинхронном обмене; - доля объема данных блока поля, подлежащая передаче; _i - эмпирическая средняя скорость обработки поля одной переменной в i-й группе; V_i - объем данных в i-й группе (на одном процессоре). Введем отсутствие простоев при обработке блоков групп, вхождение группы в один блок и верный порядок обработки:

;

; .

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

Стратегии исключения и оптимизации пересылок апробированы при компоновке схем обменов для кластерных систем. Для системы МВС-1000/16 стратегия исключения обменов позволила повысить эффективность распараллеливания на 411% (при предсказании четырех новых точек по пяти предыдущим).

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

Новый алгоритм распределенной балансировки отличается минимальным числом пересылок. Алгоритм применим при медленном дрейфе "горячих" участков, когда , где h_min - минимальное расстояние между узлами расчетной сетки; - шаг интегрирования по времени; V_дрейф - максимальная скорость дрейфа "горячих" участков; k - коэффициент, k > 1. Для распределения нагрузки в начале итерации s + 1 с высокой степенью достоверности можно использовать информацию о времени обработки узлов на предыдущей итерации s. Считаем, что время прямо пропорционально числу проделанных итераций метода Гира. Алгоритм включает три этапа:

1. Каждый i_й процессор подсчитывает предполагаемое суммарное время L_i на обработку узлов, выделенных ему при блочном разбиении области, и передает значение L_i остальным процессорам. Подсчитывается предполагаемое среднее время выполнения итерации H.

2. Каждый i_й процессор определяет свой статус для данной итерации. Если больше допустимого дисбаланса, то i_й процессор будет участвовать в обменах данными. Если L_i > H, то он должен передать часть своих узлов другим процессорам. В противном случае процессор должен будет принять часть узлов от других процессоров. Число узлов и номера процессоров - партнеров по обмену рассчитываются по единой схеме всеми процессорами, в первую очередь распределяются наиболее длительно обрабатываемые узлы.

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

Эксперименты (см. приложение 2) показали, что эффективность распараллеливания кинетического блока на МВС_1000/16 при использовании предложенного алгоритма балансировки повышается с 91% до 98% для двух процессоров, с 91% до 100% для четырех, с 90% до 94% для восьми. Причем даже на восьми процессорах данный алгоритм показал лучшие результаты (на 12%) по сравнению с иными алгоритмами распределенной балансировки, что достигнуто за счет высокоэффективной организации обменов данными.

Обоснован выбор MPI в качестве основного программного средства распараллеливания по пространству. Также используются TCP Router и Emdedded Parix. При параллельной обработке фрагментов расчетной области многоядерными SMP-блоками используется комбинация MPI+OpenMP. Разработаны однопроцессорный и многопроцессорный варианты многоплатформенного (Windows, Linux, Parix) программного кода. Многопроцессорный вариант оптимизирован для работы с системами МВС_1000, Power XPlorer, Windows - и Linux-кластерами. В приложении 2 показаны высокая степень ускорения и хорошая общая эффективность распараллеливания на кластере ИГЭУ (95% на четырех, 75% на восьми, 81% на шестнадцати, 84% на 36 процессорах) и на МВС_1000/16 (88% на двух процессорах, 91% на четырех, 83% на девяти и 80% на одиннадцати). Комбинирование интерфейсов MPI/OpenMP создало предпосылки для более эффективного распараллеливания при использовании значительного количества процессоров в кластерной системе ИГЭУ с SMP-блоками, на которой и было получено ускорение в 30,4 раза на 36 процессорах, в 24,6 раза на 28 процессорах, в 15,6 раза на 20 процессорах, в 9 раз на 12 процессорах.

В приложении 3 описаны разработанные имитаторы Power XPlorer. Одномашинный имитатор реализует больший набор функций и обладает дополнительными отладочными возможностями по сравнению с аналогами. На базе PVM 3.4.2 разработан сетевой имитатор, позволяющий оперативно переносить параллельные программы (Embedded Parix) с Power XPlorer на Windows-кластеры. Имитаторы использованы при разработке программы моделирования распространения загрязнений.

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

Сделан обзор подходов [наиболее значимы: компилирующие шаблоны (Software Factories), прямое исполнение (Флора/FloraWare), прямое преобразование (Draco, TAMPR), решающие системы (ПРИЗ, IPGS), автоматный подход (SWITCH-технологии), специализированные системы (НОРМА, SciNapse, работы С.В. Востокина, В.Н. Нуждина) и, отчасти, ментальное программирование (IP)] к автоматизации программирования (алгоритмизации и кодирования). Принято решение о разработке нового, более гибкого подхода, позволяющего: а) легко интегрировать стратегии оптимального распараллеливания, б) производить сложное комбинирование и настройку алгоритмов, реализующих различные методики численного интегрирования.

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

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

Предлагается объектно-событийная модель (ОСМ), используемая как для описания задачи, так и для порождения программы. ОСМ представляется графом (P, E), где P - множество узлов (объектов различных классов), E - множество дуг (связей между объектами). На множестве P определим функцию принадлежности к одному из классов: class: , где C - множество классов.

Введем два основных типа связей. Основные связи определяют порядок срабатывания объектов и являются "каналами" передачи данных. Вспомогательные связи отражают структурные аспекты модели. Любая связь исходит из выходного контакта объекта и входит во входной контакт объекта. Граф модели может содержать циклы, но каждый цикл должен содержать вспомогательную связь. Подграф, построенный лишь на основных связях, является сетью.

Иерархия классов определена отображением parent: ,

где - пустая цепочка. Класс S - пятерка (N_с, I, O, F, M), где N_с - идентификатор класса; I = in (S), O = out (S) - множества входных и выходных контактов; F - поля; M - порождающие методы.

Контакты, поля и методы могут быть собственными, унаследованными и/или переопределенными. Если S = (, I^S, O^S, F^S, M^S), Q = (, I^Q, O^Q, F^Q, M^Q), причем Q = parent (S), то , , , . Допускается вызов методов-предков из M^Q в соответствующих переопределенных методах M^S.

Контакт класса S есть шестерка (N_t, T, L, Min, Max, D), где N_t - идентификатор контакта, T {входной, выходной} - тип контакта, L - множество допустимых выходных связей:

,

Min {0, 1} и Max {1, } - минимальная и максимальная степени контакта, D = D [N_k] - ассоциативный кортеж (буфер передачи данных), отражающий состояние контакта, N_k - имя k-й ячейки кортежа.

Множество полей F класса S представим ассоциативным кортежем F = F [N_j], хранящим текущее состояние объекта, принадлежащего классу S, причем N_j - имя j-й ячейки кортежа.

В множестве M класса S каждый метод является программным скриптом (на языке PHP), реализующим реакцию на некоторое событие, произошедшее в процессе интерпретации модели. Представим данное множество в виде вектора M = (M_S1, M_S2, … M_SN), где N - размер календаря событий. Реакция на событие может включать действия четырех типов: а) контекстную генерацию фрагмента кода на основании значений полей F и информации, поступившей в кортежи D входных контактов; б) планирование нового события (включение его в календарь событий); в) изменение значений в кортежах D выходных контактов; г) изменение значений в глобальном ассоциативном кортеже B - "почтовом ящике", предоставляющем дополнительные прямые средства обмена данными между объектами.

Для удобства представления модели в виде схемы существует специальный класс-контейнер S^* = (, I^*, O^*, , ), представляющий "обертку" для произвольного фрагмента модели (P^*, E^*). При этом

; ; .

При интерпретации модели контейнеры непосредственно заменяются фрагментами (P^*, E^*), включаемыми в общую модель (P, E).

Интерпретация модели осуществляется на двух взаимосвязанных уровнях: событийным и объектно-сетевым. На событийном уровне интерпретация управляется календарем событий C = (C₁,C₂, … C_N). В начале интерпретации автоматически планируются четыре обязательных события: размещение, инициализация, вызов, деинициализация. Произошедшее событие интерпретируется на объектно-сетевом уровне, после чего происходит следующее событие. Интерпретация модели заканчивается, если исчерпан календарь событий.

Для интерпретации i-го события (объектно-сетевой уровень) используется та же методика, что и для сетевого графика работ. В первую очередь срабатывают объекты без входных контактов (I = ), а также объекты, к входным контактам которых не подсоединены основные связи. При срабатывании j-го объекта активизируется его метод (M_j) _i, который генерирует фрагмент программы (символьную цепочку) (Z_j) _i, читает данные со входов, формирует данные на выходах и, возможно, планирует новые события. Все прочие объекты срабатывают тогда и только тогда, когда сработают все объекты, связанные по основным связям с их входными контактами.

Правила передачи данных по связям. Если ко входному контакту A = (, T^A, L^A, Min^A, Max^A, D^A) одного объекта подключены выходные контакты Bⁱ = (, Tⁱ, Lⁱ, Minⁱ, Maxⁱ, Dⁱ), , иных объектов, то в ячейку кортеж D^A [N_k] контакта A помещается копия результата слияния ячеек кортежей Dⁱ [N_k] контактов Bⁱ:

,

где Vⁱ - множество имен ячеек кортежей Dⁱ, знаком "?" обозначена операция слияния массивов.

Интерпретация события заканчивается, когда сработают все объекты модели. Последовательность срабатываний объектов может быть описана вектором Q, элемент q_k которого содержит номер объекта, сработавшего k-м по счету. Если несколько объектов сработали одновременно, то их номера расположены в векторе Q единой группой. По окончании интерпретации i-го события формируется цепочка

Y_i = ++…+,

где n - число объектов модели, знаком "+" обозначена операция слияния цепочек. Порожденная программа G определяется цепочкой

G = Y₁+Y₂+…+Y_N,

где N - общее количество произошедших событий.

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