автореферат диссертации по информатике, вычислительной технике и управлению, 05.13.11, диссертация на тему:Программный комплекс компьютерного исследования атмосферных процессов для многоядерных процессоров

САНКТ-ПЕТЕРБУРГСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ

На правах, рукописи

005010059

РАБА Никита Олегович

ПРОГРАММНЫЙ КОМПЛЕКС КОМПЬЮТЕРНОГО ИССЛЕДОВАНИЯ АТМОСФЕРНЫХ ПРОЦЕССОВ ДЛЯ МНОГОЯДЕРНЫХ ПРОЦЕССОРОВ

05.13.11 — Математическое и программное обеспечение вычислительных машин, комплексов и компьютерных сетей

АВТОРЕФЕРАТ

диссертации на соискание ученой степени кандидата физико-математических наук

1 С ОЕЗ I

Санкт-Петербург

2012

005010059

Работа выполнена на кафедре информатики математико-механического факультета Санкт-Петербургского государственного университета.

Научный руководитель: ’ кандидат физико-математических наук,

доцент АМПИЛОВА Наталья Борисовна.

Официальные оппоненты: доктор физико-математических наук,

профессор ФЛЕГОНТОВ Александр Владимирович (Российский государственный педагогический университет им. А.И. Герцена),

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

Ведущая организация: Санкт-Петербургский государственный

политехнический университет.

Защита диссертации состоится « 1 » мт,Ча 2012 года в 1^5# ■часов на заседании совета Д 212.232.51 по защите докторских и кандидатских диссертаций при Санкт-Петербургском государственном университете по адресу: 198504, Санкт-Петербург, Старый Петергоф, Университетский пр., 28, математико-механический факультет Санкт-Петербургского государственного университета, ауд. 405.

С диссертацией можно ознакомиться в Научной библиотеке им. М. Горького Санкт-Петербургского государственного университета по адресу: 199034, Санкт-Петербург, Университетская наб. 7/9.

Автореферат разослан «*7 » 2012 года.

Ученый секретарь дйссертацйоннс совета доктор физико-математичес наук, доцент ,

Кривулин Н.К.

ОБЩАЯ ХАРАКТЕРИСТИКА РАБОТЫ

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

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

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

Одним из способов ускорения вычислений является их распараллеливание. Современные графические процессоры (вРи — ОгарЫсБ

з

Processing Unit) имеют гораздо больше вычислительных ядер по сравнению с центральными процессорами (CPU — Central Processing Unit), поэтому еще один способ ускорения — перенос вычислений на GPU.

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

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

Основные задачи:

• разработать и реализовать программный комплекс для компьютерного исследования развития конвективного облака на основе расширенного представления;

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

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

• разработать модификации этих алгоритмов для графических процессоров.

Основные результаты:

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

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

• разработаны методы распараллеливания вычислений между ядрами процессора: распараллеливание по экспериментам и по пространству с использованием потоков; приведен способ балансировки нагрузки; для 4-х ядерного процессора распараллеливание позволило сократить время расчета в 2.7 раза;

• разработан многопоточный алгоритм расчета коагуляции для GPU с использованием технологии CUDA, представлен способ разбиения вычислений на несколько этапов, приведен и проанализирован способ распределения расчетов по потокам, описан метод предотвращения коллизий записи при обращении к памяти несколькими потоками одновременно, показаны способы оптимизации доступа к памяти; разработанный алгоритм позволил сократить время расчета в 15-20 раз.

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

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

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

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

Разработка алгоритма расчета коагуляции для GPU позволила сократить время вычислений в 15-20 раз и обеспечить оперативность расчетов. Этот алгоритм является универсальным и допускает применение в представлениях размерности 2 и 3.

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

Апробация работы. Результаты работы докладывались:

• на научной конференции институтов Росгидромета, посвящённой 50-летию отдела физики облаков Главной геофизической обсерватории им. Л.И. Воейкова «Теоретические и экспериментальные исследования конвективных облаков» (Санкт-Петербург, 18-20 ноября, 2008);

• на 5-ой Международной конференции «Dynamical Systems and Applications» (Констанца, Румыния, 15-18 июня, 2009);

• на International Conference on Computational Science, ICCS 2010 (Амстердам, Нидерланды, 31 мая - 2 июня, 2010);

• на International Conference on Computational Science and Its Applications, 1CCSA 2010 (Фукуока, Япония, 23-26 марта, 2010);

• на 9-ой Международной научно-практической конференции «Исследование, разработка и применение высоких технологий в промышленности» (Санкт-Петербург, 22-23 апреля, 2010);

• на International Conference on Computational Science and Its Applications, ICCSA 2011 (Сантандер, Испания, 20-23 июня, 2011);

• на 64-ой научной конференции «Герцеповские чтения» (Санкт-Петербург, 11-16 апреля, 2011).

Публикации. Основные результаты, полученные при выполнении диссертационной работы, опубликованы в 11 работах автора, перечисленных в прилагающемся списке. Работы [1-3] опубликованы в изданиях по перечню ВАК. В работах [3-10] автору принадлежит разработка и реализация представлений облака, алгоритмов распараллеливания на ядра процессора, тестирование, а соавторам — общие постановки задач.

Структура и объем работы. Диссертационная работа объемом 160 машинописных страниц состоит из введения, пяти глав, четырех приложений и списка литературы, включающего 64 наименования, и содержит 42 рисунка и 24 таблицы.

СОДЕРЖАНИЕ РАБОТЫ

Введение содержит обзор современного состояния предметной области и обоснование актуальности диссертационной работы. Сформулирована цель работы и представлены выносимые на защиту положения. Также рассмотрена структура и краткое содержание диссертации по главам.

В главе 1 приведено описание представлений конвективного облака, разработанных автором, основанных на полуторамерном представлении [12], подробно описанном в приложении Б.

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

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

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

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

В представлении 2 используется спектральное описание микрофизичсских процессов для жидкой фазы, что увеличивает точность компьютерного исследования. В каждом узле сетки хранится информация о распределениях по размерам (спектрах) капель и частиц аэрозоля. Весь диапазон масс капель делится на N интервалов. Каждый интервал характеризуется средней массой капли в нем (для г-го интервала масса капли т,). Спектр представляется в виде набора значений /, / от 1 до /V, где / — концентрация капель в г-ом интервале в единице объема. Чем больше размер капли, тем больше скорость падения. Маленькие капли, скорость падения которых почти равна пулю, образуют облако, крупные капли выпадают в виде осадков. Учитываются следующие микрофизические процессы: нуклеация (образование новых капель на частицах аэрозоля), конденсация, испарение, коагуляция (столкновение и слияние частиц друг с другом). Возникновение осадков связано главным образом с процессом коагуляции, т.к. благодаря коагуляции возникают капли таких размеров, при которых они способны падать со скоростью, большей, чем скорость восходящих потоков.

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

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

Сравнение разработанных представлений с традиционным и проверка их адекватности приведены в приложении В.

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

Уравнение для изменения концентрации / капель массы т (стохастическое уравнение коагуляции) имеет вид:

Э/ (т) / 3/ = \ | К(т-т', т)/(от - от')/(т')с!т' - /(от) К(т, от')/(т')с!т',

где К(т, тг) — ядро коагуляции (вероятность столкновения и слияния капель с массами от и от1).

Формула перераспределения капель, вызванного коагуляцией, для спектра в виде набора/ была предложена Коветцом и Олундом [13]:

«А/,л - Л", ад.

где К)к = К{тр тк),

{г] + гък -/;!,)/(л-3 - /;!,), если гД < г] + гък < г;3, в,ь = ■ (С* -'•/ -гк3)/('•,!, -/-/), если /-3 <г] +/•; <г^,

О, в остальных случаях,

г, — радиус капель массы т,. Для расчета всех новых значений /, / = 1, ..., А' такой алгоритм требует 0(тУ3) операций.

Автором был разработан следующий алгоритм. Вместо трехмерного массива В^, заполняется матрица весовых коэффициентов Л/к и матрица индексов 1ц\ А,к = \ - т,- тк) / (т,+] - от,) = В/к,, 1]к = г, если от, < т, + тк < т1+1

(если такого г не находится, то 1)к = 0). Заметим, что В/к1 = Л/к, если / = 1/к, В/к, =

1 -А]к, если / = /д + 1, и В/к, = 0 в остальных случаях.

За счет столкновения капель из интервалов у и к концентрация капель в интервалах у и к уменьшится, в интервалах /' и г + 1 увеличится (т, < т, + тк < /и,+ь т.е. / = 1)к\ капли, образовавшиеся при слиянии перераспределятся между интервалами / и / + 1 с весовыми коэффициентами А/к и 1 - Л/к): д//д1\,* = -КАД/к, д/кЩ,к = -К,к/,/к,

ЩЩ,к = К]к/,/кА]к, с[, {Щк = К)к/,/к (1 -А1к), д//дс\1к = 0, для х, не равных у, к, / и г I 1, где д//д1\1к — изменение концентрации капель в интервале у за счет

столкновения капель из интервалов у и к. Общая масса капель при столкновении не меняется.

Изменение концентрации капель в интервале .у после коагуляции:

Расчет изменения концентрации во всех интервалах после коагуляции (используется дополнительные массивы и Б' размерности /V):

1. Инициализировать £ и Л": $Г, := 0, 5", := О, для / от 1 до N.

2. Перебрать все пары (/', А:) интервалов сталкивающихся частиц, у < к, и для

каждой пары, для которой /,* Ф 0, выполнить: / := 1/к, Г := Л’, := 5^ + ^

5 * := + Я 5', := 5*, + 5 V1 := 5 V, + ^ (1 -ЛД -.

3. Вычислить д//д! := Л",-Л’ „ , для / от 1 до /V.

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

получающейся частицы зависит как от типов и масс сталкивающихся частиц,

гак и от окружающей температуры. Уравнение для изменения концентрации / частиц г-го типа массы т:

д/ (т) / 3/ = (Г) | (от - т\ т')/) (т- т)/к (т)с1т -

-/,("0Хи Г

где Л//-7’ — количество типов частиц, С,,к(Г) — коэффициент, который равен 1, если при слиянии частицы у-го типа с меньшей по массе частицей к-го типа при температуре Т образуется частица г-го типа, и равен 0 в противном случае, К,/т, т) — ядро коагуляции (вероятность столкновения и слияния частицы г-го типа с массой т и частицы у-го типа с массой т1).

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

X и X' размерности NrrxN^, — концентрация частиц типа п в интервале /; интервалы для частиц всех типов берутся одинаковыми):

1. Инициализировать 5' и 5'': 3 := О, Л’ ’ := 0, для п от 1 до Ыгг, / от 1 до N.

2. Перебрать все пары (/, к) интервалов сталкивающихся частиц, и для каждой пары, для которой 1]к Ф 0, перебрать типы сталкивающихся частиц т (тип частиц из интервалау), п (тип частиц из интервала к), т < п. Для каждой четверки (/л, п, у, к), для которой либо у < к, либо т Ф п, определить р — тип частицы, образующийся при слиянии (используется рассчитанная заранее матрица переходов), и выполнить: /' := 1)к, /<’ := К„т]к/т;/пк, Б „„ := Л’ „ч + 1\

3~„к:=5Г„к + Р, 5^,, := 5',,/ + Т7А$, Б¥рН\ := 5+/;/,1 + И( \ -А]к).

3. Вычислить д/„/д1 := Б - Б „„ для п от 1 до Nn■, / от 1 до N.

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

Доказано следующее

Утверждение. Сложность разработанного алгоритма расчета коагуляции для представления с жидкой микрофизикой 0(М2) (т.е. меньше, чем алгоритм Коветца-Олунда), для представления со смешанной микрофизикой 0(Д’2;.; • /V2).

В главе 3 приведены методы распараллеливания для разработанных представлений.

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

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

Запуск потоков, а также переключение между потоками, выполняемыми на одном ядре, требует некоторого времени, поэтому следует уменьшить их количество. Для А'-ядерпого процессора, ядра которого могут оптимизировать выполнение Т потоков (Т = 2 для современных процессоров с поддержкой технологии Hyper-Threading), оптимально запускать ТК потоков.

Для наилучшей загрузки ядер процессора необходимо примерно равное время выполнения потоков (балансировка нагрузки). Расчет параметров в различных узлах сетки требует разного времени (например, в тех узлах, в которых нет частиц и влажность меньше 100%, расчет микрофизических процессов проводить не нужно). Так как спектры частиц и влажность меняются незначительно от узла к узлу, то расчет параметров в соседних узлах сетки занимает примерно одинаковое время. Поэтому для балансировки нагрузки можно распределять узлы сетки по ТК потокам так, чтобы соседние узлы рассчитывались в разных потоках. В первом потоке рассчитываются все динамические и микрофизические процессы в узлах с номерами 1 + i-TK (i = 0,

1, 2,...), во втором — с номерами 2 + i-TK,..., в ТК— с номерами ТК + i-TK.

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

Современные графические процессоры (GPU) имеют гораздо больше вычислительных ядер (несколько сотен) по сравнению с CPU (2-8). Чтобы полностью реализовать преимущество GPU, нужны специальные «существенно параллельные» алгоритмы (их выполнение должно разбиваться па тысячи потоков, которые выполняют одни и те же действия над разными данными). Также ссть дополнительные ограничения: действия должны быть достаточно простыми, должно быть как можно меньше ветвлений, синхронизаций потоков и обращений к глобальной памяти.

В настоящее время существует несколько технологий для использования графических процессоров для вычислений общего назначения (GPGPU — General-Purpose computing on Graphics Processing Units): CUDA, ATI Stream, OpenCL, DirectCompute. Одной из самых популярных технологий является CUDA фирмы NVidia. Она позволяет наиболее полно использовать возможности GPU, а в ее основе лежит язык C/C++ с некоторыми расширениями и ограничениями, что упрощает написание программ.

Программа на С1ЮА задействует и СРи и вРИ На СРи выполняются последовательные вычисления, а на ОРи массивно-параллельные. Эти параллельные вычисления представляются в виде одновременно выполняющихся потоков.

Представлен следующий способ вычисления. Каждый поток обрабатывает столкновение частицы типа т из массовой категории ] с частицами разных типов (и > т) из массовой категории к, т.е. один поток обрабатывает максимум ЫГТ столкновений (или меньше, например, когда к </). Значения у, т и к зависят от индексов потока /ИгеасИсЬс и блока Ыоск/ск (в С1ЮА потоки объединяются в блоки, блоки в сетки, потоки и блоки имеют свои многомерные индексы, потоки внутри одного блока могут синхронизироваться и имеют доступ к быстрой разделяемой памяти).

Одновременно несколько потоков могут пытаться изменить значения в одной и той же ячейке памяти (например, в одной и той же ячейке массива 5'), при этом неизвестно, какое значение там будет в результате.

Для предотвращения таких коллизий записи был разработан следующий алгоритм. Для каждой пары интервалов (/, к) сталкивающихся частиц типов т и п, создаются 4 отдельно выделенные ячейки памяти: для уменьшения концентрации частиц типа т в интервале у и типа п в интервале к, и для увеличения концентрации частиц типар в интервалах / и / + 1 (/ = 1/к). Для этого используется «развернутый» массив МР (вместо массивов .V и 5 ). Для каждой четверки (т, п, у, к) хранятся 4 соответствующих индекса для этого массива: тфт,ф, МЬт,1к, Iпёгтщк, м<1)тп]к. Индексы рассчитываются таким образом, чтобы ячейки памяти, относящиеся к увеличению (или уменьшению) концентрации частиц одного типа находились рядом друг с другом (т.е. образовывали последовательные области в массиве МР). Более того, ячейки памяти, относящиеся к увеличению (или уменьшению) концентрации частиц одного типа в одном массовом интервале тоже должны находиться рядом друг с другом. Размеры таких областей (количество ячеек) должны быть кратны КВБ (значение этого параметра объяснено далее). Некоторые ячейки могут остаться незаполненными, для них может не быть соответствующего индекса тс1. Индексы концов этих последовательных областей: ро50п/ — индекс последнего элемента в массиве МР, относящегося к уменьшению концентрации частиц типа п в интервале г, роя— индекс последнего элемента в массиве МР,

относящегося к увеличению концентрации частиц типа п в интервале Приводится способ вычисления индексов ind и pos и размерности массива МР.

Расчет изменения концентрации происходит в 3 этапа. На первом каждый поток для каждого столкновения (соответствующего ему) рассчитывает F = Ктщк /щ fnk и заполняет четыре ячейки массива МР: MP[ind,)m„jk\ = -F, MP[indh,mik] = -F. MP[ind2„,„jk] = F A/h MP[indimnjk] = /* ’ (I - A,k). После расчета всех столкновений, массив МР оказывается полностью заполненным (он изначально заполняется 0; ячейки, для которых нет индекса ind, остаются 0). На втором и третьем этапах происходит суммирование отдельных областей массива, относящихся к одним и тем же интервалам (например, элементы массива МР с индексами от posliluA + 1 до pos0nl находятся в области, относящейся к уменьшению концентрации частиц типа п в интервале /). На втором этапе суммируются элементы блоков массива МР размерностью RBS. Это производится методом параллельной редукции, поэтому параметр RBS (reduction block size) изначально выбирается равным степени 2. Результат суммирования сохраняется в массиве МР в ячейках с индексами, кратными RBS. На третьем этапе рассчитывается окончательная сумма по областям (в каждой области суммируются только элементы, индексы которых кратны RBS, каждый поток производит суммирование в своей области). Пусть М — результат суммирования для области отpos0n,-\ + 1 доpos0„„ Р — для области от posx „, 1 + 1 до ро.Ч\„„ тогда dfjdt = М + Р.

Доказано следующее .

Утверждение. Количество обращений к памяти у приведенного алгоритма для GPU не более чем в 1.5 раза превышает количество обращений к памяти у алгоритма для CPU.

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

В главе 5 показаны результаты сравнения быстродействия различных алгоритмов вычисления коагуляции. В следующей таблице приведено время расчета коагуляции в зависимости от количества интервалов N (все остальные параметры одни и те же) на CPU (Intel Core 2 Duo E6400) и на GPU (NVidia GTX470), а также ускорение вычислений на GPU:

. ДГ : ' . 64 128 256 512

СРи (сек.) ; -3.34 . 13.72 55.48 225.02

вРи (сек.) . .. ■ 0.203 0.688 2.81 11.53

Ускорение (раз) 16.5 19.9 19.7 19.5

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

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

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

В приложении В проведено сравнение поведения облака в представлениях с одним и двумя цилиндрами. Приведенные результаты показывают присутствие фазы диссипации в жизненном цикле облака в представлении с двумя цилиндрами.

Приложение Г содержит описание интерфейса программного комплекса и входных данных.

СПИСОК ОПУБЛИКОВАННЫХ РАБОТ НО ТЕМЕ ДИССЕРТАЦИИ

Статьи в журналах, рекомендованных ВАК:

[1] Раба Н.О. Оптимизация алгоритмов расчета физических процессов в модели облака с детальным описанием микрофизики // Вест. СПбГУ. Сер. 10, Прикладная математика. Информатика. Процессы управления —2010,-Вып. 3С. 121-126.

[2] Раба Н.О. Разработка и реализация алгоритма расчета коагуляции в

модели облаков со смешанной фазой с использованием технологии CUDA // Вестн. СПбГУ. Сер. 10, Прикладная математика. Информатика. Процессы управления.-2011.-Вып. 4. - С. 94-104. •

[3] Раба И.О., Станкова Е.Н. Исследование влияния компенсирующего нисходящего потока на жизненный цикл конвективного облака с помощью численной полуторамерной модели с двумя цилиндрами // Труды Главной геофизической обсерватории им. А.И.Воейкова.-2009-Вып. 559-С. 192-209.

Другие публикации:

[4] Станкова Е.Н., Раба И.О., Ампилова Н.Б. Численное моделирование

микрофизических процессов в смешанных конвективных облаках. Сравнение результатов расчётов с данными натурных экспериментов // Научная конференция институтов Росгидромета «Теоретические и экспериментальные исследования конвективных облаков». - СПб., 2008. - С. 90-91. .

[5] Станкова Е.Н., Раба Н.О., Ампилова Н.Б. О влиянии алгоритма

определения высоты нижней границы облака . на прогнозирование конвективных явлений с помощью численных моделей // Научная конференция институтов Росгидромета «Теоретические и экспериментальные исследования конвективных облаков». - СПб., 2008. — С. 43—44.... ......

[6] Raba N.. Stankova Е„ Ampilova N. One-and-a-half-dimensional Model of

Cumulus Cloud with Two Cylinders. Research of Influence of Compensating Descending Flow on Development of Cloud // Ovidius University Annals. Series: Civil Engineering. - Constantza: Ovidius University Press, 2009. - Vol. 1, Special Issue 11.-P. 93-101. ‘

[7] Раба НО., Станкова E.H., Ампилова Н.Б. Распараллеливание вычислений при моделировании конвективных облаков // Высокие технологии, исследования, промышленность. Сборник трудов 9-ой международной научнопрактическая конференции «Исследование, разработка и применение высоких

А

l.V

технологий в промышленности». - СПб.: Изд-во Политехи, ун-та, 2010. - Т. 1. — С.213-216.

[8] Raba N, Stankova Е., Ampilova N. On Investigation of Parallelization Effectiveness with the Help of Multi-core Processors // Procedia Computer Science. -2010.-Vol. 1, Issue l.-P. 2757-2762.

[9] Raba N.. Stankova E. On the Possibilities of Multi-Core Processor Use for Real-Time Forecast of Dangerous Convective Phenomena II Lecture Notes in Computer Science. Computational Science and Its Applications - ICCSA 2010. -2010.-Vol. 6017.-P. 130-138.

[10] N. Raba, E. Stankova On the Problem of Numerical Modeling of Dangerous Convective Phenomena: Possibilities of Real-Time Forecast with the Help of Multi-core Processors // Lecture Notes in Computer Science. Computational Science and Its Applications - ICCSA 2011. - 2011. - Vol. 6786. - P. 633-642.

[11] Раба И.О. Использование технологии CUDA для расчета коагуляции в модели облаков // Некоторые актуальные проблемы современной математики и математического образования. Герценовские чтения - 2011. - СПб., 2011. -С. 198-205.

Цитированная литература

[12] Shiino J. A Numerical Study of Precipitation Development in Cumulus Clouds // Papers in Meteorology and Geophysics. - 1978. - Vol. 29, № 4. - P. 157— 194.

[13] Kovetz A., Olund B. The Effect of Coalescence and Condensation on Rain Formation in a Cloud of Finite Vertical Extent // J. Atm. Sci. - 1969. - Vol. 26. -P. 1060-1065.

Подписано к печати 19.01.2012г. Формат 60 х 90 1/16. Бумага офсетная. Гарнитура Times. Печать цифровая. Уел. печ. л. 1. Тираж 100 экз. Заказ________________________________.

Отпечатано в типографии «Адмирал»

199048, Санкт-Петербург, В.О., 6-я линия, д. 59 корпус 1, оф. 40 Тел.: (812) 716-5223

61 12-1/647

Санкт-Петербургский государственный университет

Раба Никита Олегович

программный комплекс компьютерного

исследования атмосферных процессов для многоядерных процессоров

05.13.11 — Математическое и программное обеспечение вычислительных машин, комплексов и компьютерных сетей

Диссертация на соискание ученой степени кандидата физико-математических наук

Научный руководитель кандидат физико-математических наук, доцент Ампилова Наталья Борисовна

Санкт-Петербург 2011

Оглавление

Введение..................................................................................................................4

Глава 1 Представления конвективных облаков..........................................25

1.1 Представление с 2-мя цилиндрами и параметризованной смешанной микрофизикой....................................................................................................25

1.2 Представление с 2-мя цилиндрами и спектральной жидкой микрофизикой....................................................................................................31

1.3 Представление с 2-мя цилиндрами и спектральной смешанной микрофизикой....................................................................................................41

Глава 2 Алгоритмы расчета коагуляции и их оптимизация.....................51

2.1 Алгоритм расчета коагуляции...................................................................51

2.2 Оптимизации алгоритма расчета коагуляции..........................................54

2.3 Алгоритм расчета коагуляции для представления со смешанной фазой...................................................................................................................62

Глава 3 Расчет с помощью многоядерных CPU...........................................71

3.1 Общие принципы распараллеливания......................................................71

3.2 Распараллеливание по экспериментам.....................................................73

3.3 Распараллеливание по пространству........................................................74

Глава 4 Расчет с помощью GPU.......................................................................77

4.1 Основные принципы расчетов на GPU, технология CUDA...................77

4.2 Алгоритм расчета коагуляции на GPU.....................................................82

4.2.1 Основные принципы расчета коагуляции на GPU...........................82

4.2.2 Способ предотвращения коллизий записи........................................85

4.2.3 Алгоритм вычисления индексов ind и pos.........................................90

4.2.4 Оптимизация суммирования...............................................................93

4.2.5 Сравнение алгоритмов расчета коагуляции на CPU и GPU..........101

Глава 5 Результаты тестирования.................................................................106

5.1 Алгоритмы расчета коагуляции...............................................................107

5.2 Распараллеливание по экспериментам...................................................111

5.3 Распараллеливание по пространству......................................................112

5.4 Расчет на GPU............................................................................................114

Заключение.........................................................................................................121

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

Приложение А Классификация представлений конвективных

облаков................................................................................................................132

Приложение Б Описание традиционного полуторамерного

представления....................................................................................................135

Приложение В Исследование влияния второго цилиндра на жизненный

цикл облака........................................................................................................142

Приложение Г Описание пользовательского интерфейса.......................151

Г. 1 Интерфейс программы для представления с параметрической

микрофизикой..................................................................................................151

Г.2 Интерфейс программы для представления со спектральной микрофизикой..................................................................................................157

Введение

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

Еще одной задачей является оценка эффекта активного воздействия на конвективное облако. К этой же задаче относится и определение оптимальных параметров воздействия для достижения необходимых результатов. В [26] приведены основные параметры воздействия: время воздействия относительно стадии жизни облака, место введения и дозировка реагента. Целями воздействия на облако могут быть как предотвращение, так и вызывание осадков, а также изменение характеристик облака (например, борьба с грозовыми облаками для обеспечения безопасности полетов).

Одним из наиболее эффективных инструментов изучения конвективных облаков и механизмов воздействия на них является компьютерное исследование, которое позволяет, не прибегая к дорогостоящим натурным экспериментам, провести анализ развития облака, проследить влияние различных способов воздействия и метеорологических условий на ход процессов облако- и осадкообразования. Такие исследования с помощью компьютерного моделирования описаны в [6] (изучение эволюции облака, выявление опасности обледенения или поражения самолета молнией, определение дальности видимости, турбулентности, вероятности града), [7] (исследование развития облака при взрыве или пожаре), [2, 10, 13, 16] (изучение влияния активных воздействий на предотвращение или увеличение осадков).

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

По размерности пространства различают одномерные, полуторамерные, двумерные и трехмерные представления. Для научных исследований, направленных, прежде всего, на детальное изучение динамических, микрофизических процессов в облаке и их взаимосвязи, необходимо использовать представления, наиболее полно отражающие природные процессы, т.е. двух- и трехмерные (см., например, [17, 36, 39, 49, 52]). Такие представления достаточно сложные и требуют больших вычислительных ресурсов.

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

Наилучшими, с этой точки зрения, являются полуторамерные представления. В основе современных полуторамерных представлений (например, [2]) лежит модель Шиино [50]. В ней облако образуется в области в форме цилиндра. Эта область разбивается на слои по высоте. Значение параметров, таких как температура, скорость потока, влажность и др., усредняются в каждом слое. В отличие от одномерных представлений, где скорость потока имеет только вертикальную составляющую, и взаимодействие происходит только в вертикальном направлении, в полуторамерных дополнительно учитывается радиальная составляющая скорости (в горизонтальном направлении; направлена внутрь цилиндра или во внешнюю среду) и взаимодействие с внешней средой через боковую поверхность цилиндра.

Полуторамерные представления позволяют достаточно быстро делать достоверные прогнозы на краткосрочный период и исследовать активные воздействия на метеорологические процессы. Описание исследования с применением таких представлений приведено в [7, 10, 13, 16], а одна из их модификаций прошла независимые испытания в аэропорту «Пулково» и рекомендована для использования в оперативной практике [6, 9].

Однако, при некоторых данных радиозондирования атмосферы (реальных и модельных) такое полуторамерное представление не может адекватно воспроизвести все стадии жизни облака, в частности, стадию диссипации (рассеяния). При этом облако проходит стадию развития и зрелости, затем фактически стабилизируется: ни динамические, ни микрофизические характеристики не претерпевают существенных изменений. Подобное поведение облака при использовании данного представления описано в [10, 13]. Это вызвано отсутствием нисходящего потока, который должен компенсировать развитие восходящего потока в облаке (роль нисходящих компенсирующих движений исследовалась в [27]).

Решение этой проблемы путем введения второго (внешнего) цилиндра было предложено Асаи и Касахарой [28]. Пространство в их модели разбивается на два концентрических цилиндра с разными диаметрами. В этом случае внутренний цилиндр соответствует области с восходящим потоком (в этом цилиндре и формируется облако), внешний — области с компенсирующим нисходящим движением. Значения параметров усредняются в каждом слое для каждого из цилиндров. Происходит взаимодействие в горизонтальном направлении между внутренним и внешним цилиндром. Однако, в отличие от модели Шиино, в полуторамерной модели Асаи и Касахары не учитывается взаимодействие с окружающей средой и не моделируются осадки. Поэтому актуальной является разработка представления, сочетающего преимущества этих двух моделей.

Многие современные представления облаков (не только полуторамерные, но и с большей размерностью) используют параметрическое описание микрофизических процессов (например, [2, 5, 7, 17, 50]). В таких представлениях форма спектров частиц не учитывается, фигурируют только интегральные характеристики от функции распределения облачных элементов по размерам. Предполагается, что существуют облачные капли (маленькие), которые не движутся относительно потока, дождевые капли (крупные) и градины, которые падают относительно потока (некоторые представления не учитывают градины, а иногда даже и дождевые капли, т.е. совсем не моделируют образование осадков, см. [28]). Переходы между паром, облачными каплями, дождевыми каплями и градинами (т.е. процессы конденсации, испарения, замерзания, таяния, нуклеации, коагуляции и др.) задаются параметрически и никак не учитывают реальные размеры частиц. Скорость падения дождевых капель и градин рассчитывается исходя из их концентраций. Такое описание микрофизических процессов является упрощенным.

Для увеличения точности прогноза желательно использовать не параметрическое, а наиболее полное спектральное описание микрофизических процессов. Представления со спектральной (детальной) микрофизикой (например, [15, 36-38, 40]) учитывают спектры частиц (используются кинетические уравнения для спектров, функция распределения присутствует в явном виде). Частицы разных размеров имеют разные устоявшиеся скорости падения относительно потока. Микрофизические процессы описываются более подробно и учитывают размеры, а иногда и форму частиц. Эти представления позволяют проследить за эволюцией спектров, т.к. форма спектра рассчитывается в течение моделирования. В [35] показано, что представления со спектральной микрофизикой способны наиболее подробно воспроизводить облачные процессы.

К настоящему времени было разработано несколько пресдставлений облака со спектральной микрофизикой. В представлении [40] уделено внимание только микрофизическим процессам, используется упрощенная динамическая часть, нет ледяных частиц. В полуторамерных представлениях [15, 43] не учитываются ледяные частицы и нисходящий поток. Двумерные представления [36-38] являются более совершенными, но слишком сложными с точки зрения вычислений и поэтому не подходят для оперативного прогноза. Также вычислительно сложным является и трехмерное представление [39], к тому же в нем отсутствуют ледяные частицы.

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

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

Различные методы и алгоритмы расчета коагуляции обсуждаются в [32, 35, 44]. В настоящий момент для решения уравнения коагуляции применяются: метод Коветца и Олунда [40] и его модификации [51], метод Берри и Реинхарда [30], метод моментов [53].

Метод Берри-Рейнхарда [30] применяется в представлениях [38, 39] и является одним из самых точных. В этом методе спектр дискретизируется (значения запоминаются только для узлов сетки по массам). Для решения уравнения коагуляции используется численное интегрирование. Для расчета значения функции распределения в точках, отличных от узлов сетки, использовуются 6-точечные формулы интерполяции Лагранжа. В [35, 44] описаны существенные недостатки этого метода: способ медленный, не сохраняет общую массу частиц (до 30% уменьшения массы в некоторых случаях). В [32] показано, что из-за этих недостатков нежелательно использовать этот метод в представлениях с детальной микрофизикой.

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

Метод Коветца-Олунда [40, 51] использует дискретизацию спектра (разбиение на интервалы и сохранение концентрации в каждом интервале). Частицы при слиянии пропорционально распределяются между двумя ближайшими смежными интервалами. Достоинства данного метода: достаточно быстрый (по сравнению с другими методами), сохраняет массу частиц. Есть у метода и недостаток — при небольшом количестве интервалов метод может приводить к расширению спектра. В [44] выявлено, что этот недостаток можно устранить, если увеличить число интервалов, вследствие этого заодно улучшится определение спектра. Метод Коветца-Олунда и его модификации используется во многих (в том числе и полуторамерных) представлениях облаков [14, 15, 44].

Как было показано выше, многие представления не учитывают наличие твердой фазы (т.е. ледяных частиц) в облаке. В реальном облаке вода присутствует не только в виде пара (газообразное состояние) и капель (жидкое состояние), но и в виде ледяных частиц, поэтому необходимо учитывать это. Такие частицы могут иметь различную форму, образуются при разных условиях. Классификация ледяных частиц и условия их образования приведены в [41]. При столкновении между собой ледяные частицы образуют более сложные составные частицы (снежинки, крупа, градины в зависимости от условий). Учесть все типы частиц очень сложно. В настоящий момент предложены представления [36-38] со смешанной (т.е. жидкой и твердой) фазой, в которых учитываются капли и 6 типов ледяных частиц.

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

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

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

Большинство современных ПК имеют многоядерные (2-8 ядер) процессоры (CPU — Central Processing Unit), допускающие параллельное выполнение расчетов, поэтому одним из способов ускорения вычислений является их распаралл�