close

Вход

Забыли?

вход по аккаунту

?

Оптимизация алгоритма расчета коагуляции в модели облака со спектральной микрофизикой.

код для вставкиСкачать
Сер. 10. 2010. Вып. 3
ВЕСТНИК САНКТ-ПЕТЕРБУРГСКОГО УНИВЕРСИТЕТА
УДК 004.021+004.051
Н. О. Раба
ОПТИМИЗАЦИЯ АЛГОРИТМА РАСЧЕТА КОАГУЛЯЦИИ
В МОДЕЛИ ОБЛАКА СО СПЕКТРАЛЬНОЙ МИКРОФИЗИКОЙ
1. Введение. Развитие конвективных облаков в определенных условиях влечет
за собой такие опасные природные явления как гроза, град, шквал, смерч. Прогноз
возникновения и эволюции данного типа облаков – одна из важнейших задач, стоящих
перед исследователями. К наиболее эффективным инструментам изучения конвективных облаков относится численное моделирование. Существует множество моделей развития таких облаков, отличающихся как размерностью, так и степенью детализации
описания микрофизических процессов.
Для целей оперативного прогноза, например в метеорологических центрах аэропортов, где для расчетов могут использоваться не суперкомпьютеры, а обычные ПК,
необходимы простые модели (с меньшей размерностью). Они не требуют больших вычислительных ресурсов, однако способны воспроизводить те характеристики облака,
которые могут стать предикторами для существующих методов прогноза опасных конвективных явлений, практически в режиме реального времени. Наилучшими, с этой
точки зрения, являются полуторамерные модели [1, 2]. Для увеличения точности прогноза в таких моделях желательно использовать не параметрическое, а наиболее полное
(спектральное) описание микрофизических процессов [3, 4].
При применении таких моделей основная часть компьютерных ресурсов тратится
на вычисление микрофизических процессов (в зависимости от параметров в 40–600 раз
больше по сравнению с расчетом динамических процессов). В настоящей статье показаны методы оптимизации алгоритмов расчета микрофизических процессов, позволяющие ускорить вычисления. Также представлены модификации этих алгоритмов,
реализующие параллельное выполнение расчетов.
2. Краткое описание физической части модели. Используется нестационарная полуторамерная модель конвективного облака с детальным описанием микрофизических процессов для жидкой фазы. Область конвективных течений моделируется двумя концентрическими цилиндрами. Внутренний цилиндр соответствует области
с восходящими потоками (облачной области), внешний – окружающей области с нисходящими потоками (безоблачная область). Все характеристики облака (скорость, температура, давление, плотность, содержание пара, облачных капель, дождевых капель,
града) усредняются по сечению цилиндров. Вводится сетка по высоте с шагом Δh.
Ранее эта модель была подробно рассмотрена [5], но с параметрическим описанием
Раба Никита Олегович – аспирант кафедры информатики математико-механического факультета Санкт-Петербургского государственного университета. Научный руководитель: кандидат физикоматематических наук, доц. Н. Б. Ампилова. Количество опубликованных работ: 8. Научное направление: компьютерное исследование динамических систем. E-mail: no13@inbox.ru.
c Н. О. Раба, 2010
121
микрофизических процессов. Поскольку здесь применяется более детальное описание
таких процессов, остановимся на нем подробнее.
В каждом узле сетки хранится информация о спектре капель и аэрозолей. Учитываются следующие микрофизические процессы (в скобках кратко описывается их действие и указывается влияние на различные параметры):
нуклеация (образование новых капель на частицах аэрозоли; меняются спектр капель, спектр аэрозолей, содержание пара, температура);
конденсация/испарение (рост/уменьшение размеров капель из-за конденсации/испарения пара на/с их поверхности, общее количество капель при этом не меняется; изменяются спектр капель, содержание пара, температура);
коагуляция (столкновение и слияние капель различных размеров друг с другом,
меняется спектр капель, общая масса капель остается постоянной; благодаря данному
процессу происходит формирование капель значительных размеров).
На расчет коагуляции тратится наибольшая часть времени вычисления, потому
в статье описан процесс оптимизации алгоритма вычисления именно такого процесса.
2.1. Коагуляция. Приведем уравнение для изменения концентрации f капель массы m:
! ∞
!
1 m
∂f (m)
=
K(m − m , m )f (m − m )f (m )dm − f (m)
K(m, m )f (m )dm . (1)
∂t
2 0
0
Первое слагаемое в левой части (1) соответствует столкновению и слиянию капель,
суммарная масса которых равна m (слияние капель приводит к увеличению концентрации капель массы m), второе слагаемое – столкновению и слиянию капель массы m
с другими каплями (их слияние вызывает уменьшение концентрации капель массы m);
K(m, m ) – ядро коагуляции (вероятность столкновения и слияния капель с массами m
и m ), которое рассчитывается по формуле
K(m, m ) = π(r(m) + r(m ))E(r(m), r(m ))|V (m) − V (m )|,
(2)
где E(r, r ) – эффективность захвата; V (m) – скорость (установившаяся) падения капель массы m.
3. Алгоритмы.
3.1. Расчет коагуляции. Алгоритм перераспределения капель, вызванного коагуляцией, был разработан Коветцом и Олундом [3]. Он заключается в следующем.
Спектр капель представляется в виде зависимости концентрации капель от радиуса.
Спектр делится на N промежутков. В каждом из них запоминается концентрация капель (с радиусами в данном промежутке).
Уравнение (1) приобретает следующий вид:
i−1
N
i
∂fi =
Bjki Kjk fj fk − fi
Kij fj ,
∂t
j=1
j=1
(3)
k=j+1
здесь Kjk = K(mj , mk ) = K(m(rj ), m(rk )) – ядро коагуляции, fi = f (ri ) – концентрация капель с радиусом, находящимся в промежутке с центром ri ;
⎧ 3
3
r
−rj3 −rk
3
⎪
, если ri3 < rj3 + rk3 < ri+1
,
⎪ i+1
3
ri+1
−ri3
⎨
3
3
3
rj +rk −ri−1
Bjki =
3
, если ri−1
< rj3 + rk3 < ri3 ,
3
3
⎪
⎪
⎩ ri −ri−1
0,
в остальных случаях.
122
Для расчета коагуляции в данной модели применялся модифицированный метод
Коветца–Олунда [4]. Вместо разбиения по категориям радиусов (ri ) было выбрано разбиение по массам капель (mi ). Рассчитывается масса капель в каждой категории масс
(Mi = fi mi ), а не концентрация (fi ). Предполагается, что последний промежуток соответствует каплям с бесконечно большой массой и концентрация капель в нем равна 0,
т. е. fN = 0. При домножении на mi уравнение (3) преобразуется:
i−1 N
−1
i
∂Mi
= mi
Bjki Kjk fj fk − Mi
Kij fj .
∂t
j=1
j=1
(4)
k=j+1
Можно выразить Bjki через категории масс таким образом:
⎧ m −m −m
i+1
j
k
⎪
, если mi < mj + mk < mi+1 ,
⎨
mi+1 −mi
mj +mk −mi−1
Bjki =
, если mi−1 < mj + mk < mi ,
mi −mi−1
⎪
⎩
0,
в остальных случаях.
Представим правую часть уравнения (4) в виде Si+ − Si− :
Si− = Mi
N
−1
Kij fj =
j=1
N
−1
Kij Mi Mj /mj =
j=1
N
−1
−
Sij
,
j=1
−
= Kij Mi Mj /mj ,
где Sij
Si+ = mi
i−1 i
Bjki Kjk fj fk =
j=1 k=j+1
i−1 i
(mj + mk )Bjki Kjk fj fk .
j=1 k=j+1
Учтем, что коэффициенты Bjki симметричны относительно j и k, Kjk – тоже симметричны относительно j и k, Kjj = 0 (см. уравнение (2)). Тогда
Si+ =
i−1 i
mj Bjki Kjk fj fk +
j=1 k=j+1
=
i−1
i−1
mj
j=1
i
mj
i
i
j=1
Bjki Kjk fj fk +
k=j+1
Bjki Kjk fj fk +
mj
i
k=1
i
mk
k=2
i
mj
j=2
k=j+1
=
mk Bjki Kjk fj fk =
k=2 j=1
j=1
=
i k−1
Bjki Kjk fj fk =
j−1
k−1
Bjki Kjk fj fk + 0 =
j=1
Bjki Kjk fj fk +
k=1
i i
j=1 k=1
i
mj Bjjk Kjj fj fj =
j=1
Bjki Kjk Mj Mk /mk =
i i
−
Bjki Sjk
.
j=1 k=1
3.2. Оптимизация расчета коагуляции.
Базовый алгоритм. Поскольку массы mi фиксированные, то матрица Kij (ядра
коагуляции) заполняется в начале и не пересчитывается на каждом шаге (для каждого
узла сетки по высоте своя матрица, так как скорость падения капель и эффективность
захвата зависят от высоты). Коэффициенты Bijk также вычисляются только в начале
(не зависят от высоты). Считаем новые значения Mi по алгоритму:
123
для всех i от 1 до N − 1
для всех j от 1 до N − 1
−
Sij
= Kij Mi Mj /mj
для всех i от 1 до N − 1
Si− = 0
для всех j от 1 до N − 1
−
Si− = Si− + Sij
для всех i от 1 до N − 1
Si+ = 0
для всех j от 1 до i
для всех k от 1 до i
−
Si+ = Si+ + Bjki Sjk
(∗)
для всех i от 1 до N − 1
Mi = Mi + (Si+ − Si− )dt.
−
Очевидно, что для определения массива Si+ (при вычисленном Sij
) по этому алгорит3
му требуется O(N ) операций ((N − 1)N (2N − 1)/6 раз выполнится присваивание (∗)
−
−
в самом внутреннем цикле), для Sij
и Si− (при вычисленном Sij
) – по O(N 2 ) операций.
Оптимизация 1. Как и в базовом алгоритме, матрица Kij рассчитывается в начале (для каждого узла сетки). Вместо 3-мерного массива Bijk заполняются матрица
весовых коэффициентов Aij и матрица индексов Iij :
Aij = Bijk , Iij = k, если mk < mi + mj < mk+1 .
З а м е ч а н и е. Если mk < mi + mi < mk+1 , то
mi + mj − mk
mk+1 − mi − mj
=1−
= 1 − Bijk .
mk+1 − mk
mk+1 − mk
Обратное преобразование:
⎧
если k = Iij ,
⎨ Aij ,
1 − Aij , если k = Iij + 1 (см. замечание),
Bijk =
⎩
0,
в остальных случаях.
Bijk+1 =
Теперь расчет новых значений Mi выглядит так:
обнулить массивы Si− и Si+
для всех i от 1 до N − 1
для всех j от 1 до N − 1
−
= Kij Mi Mj /mj
Sij
k = Iij
если 0 < k < N, то
(5)
−
Si− = Si− + Sij
−
Sk+ = Sk+ + Aij Sij
+
+
−
Sk+1 = Sk+1 + (1 − Aij )Sij
для всех i от 1 до N − 1
Mi = Mi + (Si+ − Si− )dt.
Такой алгоритм потребует выполнения только O(N 2 ) операций.
Оптимизация 2. Производится расчет столкновений лишь в области с ненулевой концентрацией капель. На каждом шаге (и для каждого узла сетки) определяются
124
индексы imin , imax , такие что imin < imax , Mimin = 0, Mimax = 0, Mj = 0, для j < imin
и для j > imax . Это потребует O(n) операций. Далее в алгоритме (5) индексы i и j
меняются не в диапазоне от 1 до N − 1, а в диапазоне от imin до imax .
Естественно, такая оптимизация дает различное ускорение вычислений в зависимости от спектрального распределения капель. В худшем случае (когда imin = 1,
imax = N − 1) этот алгоритм потребует такого же количества операций, как и алгоритм (5).
3.3. Распараллеливание вычислений. Так как большинство современных ПК
имеют многоядерные процессоры, допускающие параллельное выполнение расчетов,
то еще одним способом ускорения вычислений является их распараллеливание. Это
можно сделать путем расчета параметров в разных узлах сетки на различных ядрах
процессора (распараллеливание по пространству).
Распараллелить можно не только микрофизическую часть, но и динамическую (так
как параметры на текущем шаге времени зависят только от параметров на предыдущем
шаге и не зависят от параметров в соседних узлах на текущем шаге). Хотя в каждом
узле последовательно определяются динамические, а потом и микрофизические процессы, расчеты в нескольких узлах можно выполнять параллельно.
Реализовать параллельное вычисление можно с помощью потоков. Создаются потоки, которым передаются ссылки на данные, полученные на предыдущем шаге. Каждый
поток выполняет вычисления в определенных узлах сетки. Переход к следующему шагу
происходит, когда все потоки выполнятся.
3.4. Оптимизация распараллеливания. Запуск потока требует некоторого времени, потому следует уменьшить количество потоков. Это можно сделать, рассчитывая в отдельном потоке параметры не для одного узла сетки, а сразу для нескольких.
Для N -ядерного процессора оптимально запускать N потоков (либо 2N , если ядра могут оптимизировать выполнение двух потоков, например, при поддержке технологии
Hyper-Threading).
Поскольку на каждом шаге происходит ожидание завершения всех потоков,
то для оптимальной загрузки процессора необходимо примерно равное время их выполнения (т. е. все потоки завершатся почти одновременно). Определение параметров
в различных узлах сетки требует разного времени (например, в тех узлах, в которых
нет капель и влажность меньше 100%, расчет микрофизических процессов проводить
не нужно). Так как спектры капель и влажность меняются плавно от узла к узлу,
то расчет параметров в соседних узлах сетки занимает примерно одинаковое время.
Поэтому для достижения равного времени выполнения потоков можно распределять
узлы сетки по N потокам следующим образом. В первом вычисляются все динамические и микрофизические процессы в узлах с номерами, сравнимыми с 0 по модулю N
(т. е. n ≡ 0(mod N )); во втором – с 1, . . . ; в N -м – с N − 1. То есть при таком распределении соседние узлы рассчитываются в разных потоках, что позволяет равномернее
распределить нагрузку на все ядра.
Надо заметить, что некоторые фрагменты программы все же выполняются в однопоточном режиме (например, создание и запуск потоков, учет граничных условий
и т. п.).
3.5. Тестирование скорости вычислений. Проводились измерения скорости
вычислений следующих экспериментов: Э-50 – 51 категория масс по спектру капель,
33 категории масс по спектру аэрозолей, шаг по времени – 1 с, шаг по высоте –
100 м, радиус внутреннего столба – 3 км, радиус внешнего столба – 10 км, коэффициент бокового перемешивания – 0.1, коэффициент вертикальной турбулентности – 100,
125
моделировалась эволюция облака в течение 1 ч 30 мин; Э-70 – 71 категория масс по спектру капель, Э-100 – 101 категория, Э-150 – 151 категория, Э-250 – 251 категория, остальные параметры совпадают. Для вычисления использовался ПК в следующей конфигурации: процессор Core 2 Quad Q8200 (4 ядра, 2.33 ГГц), оперативная память 2.5 Гб.
Тестировались полуторамерные модели конвективного облака с такими алгоритмами вычисления коагуляции: базовый алгоритм (модифицированный Коветца и Олунда); оптимизированный алгоритм 1; оптимизированный алгоритм 2. Дополнительно тестировалась модель, в которой были распараллелены вычисления динамических и микрофизических процессов (для расчета коагуляции использовался оптимизированный
алгоритм 2). Результаты тестирования приведены в таблице. Отметим, ускорение вычисления оптимизированного алгоритма при многопоточном выполнении по сравнению
с однопоточным от 2 до 3.1 раз, а не в 4 раза (так как часть ресурсов тратится на создание и запуск потоков, разные потоки выполняются за различное время, и часть
вычислений проводится в однопоточном режиме).
Результаты тестирования (время расчета в секундах)
Алгоритм
Базовый
Оптимизированный 1
Оптимизированный 2
Параллельный
Э-50
402.57
155.52
86.05
42.51
Э-70
918.62
254.53
119.82
50.58
Э-100
2569.82
472.79
182.03
65.75
Э-150
7940.01
1038.17
281.48
93.84
Э-250
33673.81
2604.76
537.18
171.96
4. Заключение. Была проведена оптимизация расчетов микрофизических процессов (при сохранении точности) и достигнуто существенное ускорение вычислений. Произведено тестирование работы алгоритмов (базового и оптимизированного) при одинаковых параметрах. При достаточно большом количестве массовых категорий (150)
оптимизированный алгоритм затратил почти в 30 раз меньше времени, чем базовый.
Также проверялась эффективность параллельных вычислений. При расчете
на 4-ядерном процессоре скорость многопоточного алгоритма по сравнению с однопоточным оказалась в 3 раза быстрее (а если сравнивать с базовым алгоритмом,
то в 80 раз). Это позволяет делать прогноз развития облака более оперативно.
Такая оптимизация расчетов может быть применена и в моделях с большей размерностью, и в моделях со спектральной микрофизикой, включающей твердую фазу.
Литература
1. Asai T., Kasahara A. A Theoretical Study of the Compensating Downward Motions Associated with
Cumulus Clouds // J. of the Atmospheric Sciences. 1967. Vol. 24. P. 487–497.
2. Shiino J. A Numerical Study of Precipitation Development in Cumulus Clouds // Papers in
Meteorology and Geophysics. 1978. Vol. 29, N 4. P. 157–194.
3. Kovetz A., Olund B. The Effect of Coalescence and Condensation on Rain Formation in a Cloud of
Finite Vertical Extent // J. of the Atmospheric Sciences. 1969. Vol. 26. P. 1060–1065.
4. Stankova E. N., Zatevakhin M. A. The Modified Kovetz and Olund Method For the Numerical
Solution of Stochastic Coalescence Equation // Proc. 12th Intern. Conference on Clouds and Precipitation.
Zurich, 1996. P. 921–923.
5. Раба Н. О., Станкова Е. Н. Исследование влияния компенсирующего нисходящего потока на жизненный цикл облака с помощью численной полуторамерной модели с двумя цилиндрами // Труды Гл. геофиз. обсерватории им. А. И. Воейкова. СПб., 2009. Вып. 559. С. 192–209.
Статья рекомендована к печати член-корр. РАН, проф. Г. А. Леоновым.
Статья принята к печати 1 апреля 2010 г.
Документ
Категория
Без категории
Просмотров
10
Размер файла
228 Кб
Теги
коагуляция, алгоритм, оптимизация, спектральная, облака, расчет, микрофизикой, модель
1/--страниц
Пожаловаться на содержимое документа