CC BY
18
УДК 519.688
К. П. Беляев, А. А. Кулешов2, И. Н. Смирнов3, К. А. С. Танажура4
ПАРАЛЛЕЛЬНОЕ УСВОЕНИЕ ДАННЫХ НАБЛЮДЕНИЙ В ГИДРОДИНАМИЧЕСКОЙ МОДЕЛИ ЦИРКУЛЯЦИИ ОКЕАНА*
Предлагается параллельная реализация метода ансамблевых фильтров Калмана для усвоения данных наблюдений в модели Мирового океана HYCOM. Усваиваются данные спутниковых наблюдений температуры поверхности океана и уровня океана как по отдельности, так и совместно. Проводятся численные эксперименты по коррекции модельных расчетов с помощью данных натурных наблюдений, и результаты этой коррекции сравниваются с модельными расчетами без усвоения. Подтверждена эффективность реализованного алгоритма распараллеливания.
Ключевые слова: усвоение данных наблюдений, модель циркуляции океана, параллельные вычисления, численные эксперименты.
1. Введение. В современных численных моделях динамики Мирового океана и изменений климата производится многопроцессорный расчет решения системы уравнений динамики океана и атмосферы и одновременное корректирование получаемого решения с помощью данных измерений. Разрабатываемые в ведущих мировых прогностических центрах модели Мирового океана высокого пространственного разрешения с усвоением непрерывно поступающих данных наблюдений позволяют осуществлять оперативный прогноз трехмерного состояния морской среды, аналогичный метеорологическим прогнозам погоды, а также являются важной составляющей в исследованиях изменения климата и системах мониторинга состояния окружающей среды.
Задача коррекции численного решения системы уравнений динамики океана и атмосферы данными измерений (задача усвоения данных) заключается в том, чтобы, с одной стороны, сохранить баланс физических параметров, априори задаваемый уравнениями модели, а с другой — приблизить решение этих уравнений к наблюдаемым значениям параметров. Решение таких задач невозможно без использования технологий распределенных вычислений в силу огромного объема получаемой и обрабатываемой информации. Количество ресурсов, необходимых для выполнения оперативных расчетов моделями высокого пространственного разрешения, исчисляется сегодня 102-103 вычислительных ядер для краткосрочных прогнозов и 104-105 — для средне- и долгосрочных. Особенно критичен вопрос времени усвоения данных в случае, когда численная модель океана и модуль усвоения данных функционируют в оперативном режиме для построения среднесрочных и краткосрочных прогнозов. Данные со спутников становятся доступными спустя 2-3 часа, высокое пространственное разрешение моделей позволяет вычислять и прогнозировать динамику вихревых структур, при этом усвоение спутниковых данных наблюдений способствует их своевременному выявлению. Как следствие, это позволяет предсказывать такие опасные природные явления, как штормы и тайфуны. Так в настоящее время работает система "Меркатор", разработанная в Метеорологической службе Франции [1].
Существует несколько алгоритмов усвоения данных, которые применяются в задачах прогноза погоды и в оперативной океанологии. Используемые подходы можно разделить на вариационные (3D-Var, 4D-Var) [2] и динамико-стохастические, преимущественно, ансамблевые фильтры Калмана (Ensemble Kaiman Filter, EnKF) [3-5]. В последнее время были предложены гибридные схемы усвоения, основанные на синтезе 3D-Var и EnKF [6-8]. Методы, основанные на EnKF, не требуют построения сопряженного оператора модели и решения обратной задачи, что для модели с большим числом параметров весьма затруднительно, и используют модель как "черный ящик". Такие методы хорошо распараллеливаются и вполне применимы для глобальных моделей, в то время как
1 Институт океанологии им. П. П. Ширшова РАН, вед. науч. сотр., д.ф.-м.н., e-mail: kbel55Qyahoo.com
2 федеральный исследовательский центр Институт прикладной математики им. М.В. Келдыша Российской академии наук, гл. науч. сотр., д.ф.-м.н., e-mail: andrew_kuleshovQmail.ru
3 Факультет ВМК МГУ, ст. преп., к.ф.-м.н., e-mail: ismirnovQcs.msu.ru
4 федеральный университет штата Байия, г. Сальвадор, Бразилия, проф., e-mail: clemente.tanajuraQgmail.com
* Работа выполнена при поддержке гранта Российского научного фонда, проект № 14-11-00434.
4Б-Уаг-метод в силу вычислительной сложности на сегодняшний день не используется ни в одной глобальной модели высокого пространственного разрешения.
Однако, в целом, вопрос оптимального усвоения данных, как в методическом, так и в чисто прикладном аспекте, еще не решен. Результаты применения того или иного метода зависят как от модели, так и от региона наблюдений и конкретных наблюдаемых параметров. Поэтому каждое новое корректно выполненное исследование в данном направлении имеет свою специфику и представляет несомненный интерес. Особенно важно разработать и практически использовать малозатратные, в смысле использованной памяти и времени счета, относительно простые алгоритмы, которые способны одновременно усваивать данные из различных источников, при этом учитывая корреляционные и функциональные связи между ними.
Цель данной работы заключается в параллельной реализации на многопроцессорных компьютерах алгоритма EnKF для одновременного усвоения в модели динамики океана данных наблюдений нескольких параметров, полученных из различных источников, и в анализе результатов проведенных численных экспериментов.
2. Краткое описание модели океана, алгоритма усвоения и источников данных наблюдений. В работе используется численная модель динамики океана HYCOM (Hybrid Circulation Ocean Model) [9, 10], разработанная в Университете Майами (США), которая применяется для исследования океанских термогидродинамических процессов в широком диапазоне пространственных и временных масштабов. Модель принадлежит классу 3D и основана на так называемом сигма-приближении, когда вся глубина от поверхности до дна разбивается на слои равной плотности заданной величины, поэтому координатами всех распределенных параметров являются координаты точки на плоскости и номер слоя плотности, толщина которого есть величина переменная и расчетная. На границе раздела воздух-вода применено нелинейное кинематическое условие свободной поверхности с явным описанием потоков массы воды и соли, импульса и тепла.
Уравнения решаются численным методом конечных разностей на разнесенных сетках. Обмен теплом, импульсом и влагой на границе атмосферы и океана задан так называемыми балк-формулами CORE [11, 12]. Поскольку основным содержанием настоящей работы является параллельная реализация метода усвоения данных наблюдений, а модель динамики океана используется только как инструмент, то более детальное описание модели и численных экспериментов по ее верификации приведено в работах [9, 10, 13].
Модель HYCOM сконфигурирована в области Атлантики от 79° ю.ш. до 55° с.ш., исключая Средиземное море, но включая Саргассово море и Мексиканский залив, с использованием горизонтальной широтно-долготной сетки с разрешением (1/4)°. По вертикали модель имеет так называемую изопикническую структуру, т.е. от поверхности до дна весь океан разбивается на заранее заданные уровни равной плотности (изопикны). В данной конфигурации используется 21 уровень. Разрешение модели по горизонтали составляет 480 х 760 точек. В модели вычисляется 109 независимых переменных, а именно: уровень океана; 3 баротропные переменные — горизонтальные компоненты вектора скорости, давление на поверхности океана и 5 бароклинных переменных для каждого из 21 выделенных уровней плотности; горизонтальные компоненты вектора скорости, толщина слоя, температура и соленость. Следовательно, число компонент вектора состояния модели составляет 480 х 760 х 109 = 39763200. В качестве начальных полей температуры и солености берутся среднегодовые значения из климатического атласа WOA09 [8]. На боковых границах задаются фиксированные значения для температуры и солености из [8], а также полагается постоянным поток массы: через пролив Дрейка — 110 Св (1 Св = 106м3/с); в районе Южной Африки и вдоль 20° в.д. — 10 Св и на боковой границе вдоль Антарктиды — 120 Св. На северной границе Атлантики поток полагается равным нулю.
Коррекция модельного расчета данными наблюдений (усвоение данных) производится по следующим формулам [5, 14]:
где Ха, Хь — векторы рассчитываемого в численной модели параметра соответственно после и до усвоения (analysis and background) размерности, равной количеству точек сетки п, которое
Ха = Xb + K(Yobs - НХЬ) К = Ш1' (II nil' + R)-1,
-1
(1) (2)
для поверхности океана с разрешением 0.25° имеет порядок 107; Y0bs — вектор наблюдений размерности, равной количеству точек наблюдения т, порядка 103; К — весовая матрица (Kalman gain matrix) размерности п х т; R — ковариационная т х m-матрица инструментальных ошибок наблюдений, которая имеет диагональный вид, так как предполагается, что инструментальные ошибки не коррелированны и Н — т х n-матрица проектирования значений модели в пространство наблюдений. Матрица В имеет специальное название ковариационной матрицы состояния модели. Ее строгое определение обычно не дается, вместо этого описывается метод ее расчета.
Основная идея метода EnKF заключается в том, что ковариационная матрица В не задается в явном виде или в виде функции, а получается из ансамбля различных расчетных векторов состояния модели (выборки). Пусть Х£п = (Х^.-.Х^) — матрица размерности п х еп, где еп — количество элементов ансамбля (обычно не более 100), столбцы которой равны значениям состояния модели минус среднее по ансамблю. Тогда матрица ковариации модели строится как В = (еп — 1)~1Х£п(Х%п)Т■ Для построения матрицы Н используется билинейная интерполяция. В настоящей работе используется вычислительно более дешевый метод EnOI, подробно описанный в [6, 13], в котором в качестве элементов ансамбля выступают состояния модели, полученные и архивированные в процессе счета при различных начальных и/или граничных условиях.
В качестве наблюдаемой информации использовались данные по уровню океана (УО) архива AVISO (Archiving , Validating and Interpolating Satellite Ocean data) и данные температуры поверхности океана (ТПО) архива OSTIA (Ocean Satellite Temperature Interpolating and Archiving data). Следует отметить, что в первом случае данные выбирались вдоль треков, лежащих в зоне Атлантики, а во втором случае данные задавались в сетке архива OSTIA. Несовпадение исходной информации по их локализации технически существенно усложняет процедуру их совместного усвоения.
3. Программная реализация алгоритма усвоения данных наблюдений и качественный анализ численных решений
3.1. Программная реализация алгоритма усвоения данных наблюдений. Усвоение данных наблюдений реализовано в форме отдельного программного модуля DAS (Data Assimilation Service). Расчетные значения параметров модели динамики океана поступают на вход программного модуля DAS, который использует другую процессорную декомпозицию расчетной области, что важно, поскольку размер трехмерных массивов параметров для модели океана с разрешением 0.25° составляет несколько гигабайт. Программный модуль реализован с использованием технологии MPI. При этом удалось добиться высокого параллелизма в силу независимости усвоения данных в каждой точке наблюдения по формулам (1), (2). Благодаря параллельной реализации метода усвоения корректировка вычисляемых параметров на 128 процессорных ядрах занимает около 20 с расчетного времени, вместо 18 мин на одном ядре, что было бы сравнимо со временем, затрачиваемым на расчет суточного модельного прогноза, а такое время усвоения неприемлемо.
3.2. Численные эксперименты и их результаты. В работе были проведены численные эксперименты с моделью HYCOM по следующей схеме: сначала модель "разгонялась" в течение 40 лет модельного времени, начиная с заданных климатических полей температуры и солености и нулевых начальных скоростей под воздействием климатических потоков тепла и ветра [12]. Затем поля вычисляемых в модели параметров за последние 10 лет архивировались на каждые календарные сутки и использовались в качестве значений ансамбля для вычисления матрицы В, как описано выше в п. 2. Далее поля вычисляемых параметров на каждом суточном шаге по времени корректировались с учетом данных наблюдений по формулам (1), (2) на расчетном интервале времени — 3 месяца 2011 г., начиная с 1 апреля. Данные по атмосфере, использованные как внешние условия на границе океан-атмосфера, задавались в соответствии с их календарным временем из архива CORE [12] и интерполировались на сетку. Наблюдения за температурой поверхности океана и уровнем океана также были взяты за те же даты из архивов OSTIA и AVISO соответственно. Отдельно проводилось три численных эксперимента: эксперимент 1, в котором в качестве наблюдаемой информации использовались только поля УО; эксперимент 2, в котором использовались только поля ТПО; эксперимент 3, в котором оба этих поля использовались совместно.
Ниже приводятся результаты за 1 день (2 апреля 2011 г.) до и после усвоения различных характеристик. Сравниваются аномалии уровня океана, т. е. поля УО минус среднее многолетнее
значение УО. Из сравнения полей на рис. 1, а, б хорошо видно, что после усвоения поле аномалий уровня становится менее хаотичным, амплитуды этих аномалий уменьшаются. Размер характерных вихрей становится меньше, при этом количество самих вихрей не возрастает. Все это говорит о том, что при усвоении и удалении среднего поля уменьшается дисперсия поля аномалий, в этом и состоит цель усвоения.
Рис. 1. Результаты расчетов аномалий УО: контрольный расчет (а); разность УО поело и до усвоения: для эксперимента 1(6); для эксперимента 2 (в); для эксперимента 3 (г)
При расчете поля аномалии уровня океана с усвоением данных температуры поверхности океана, влияющей на уровень океана (см. рис. 1, в), заметные аномалии остаются только в зоне Гольфстрима и в зоне Бразило-Мальвинского столкновения на юге Атлантики, т. е. там, где наблюдается наиболее сильная динамика океана и термическое расширение уровня не компенсируется средним потоком. Наконец, на рис. 1, г видно, что при совместном усвоении данных полей аномалии уровня океана и температуры поверхности океана остается аномалия только в зоне Гвинейского тече-
ния, где сильны именно прибрежные мелкомасштабные эффекты, не учитываемые моделью или методом усвоения.
3.3. Оценка ускорения и эффективности применяемых параллельных алгоритмов.
Были проведены расчеты ускорения и эффективности применяемых параллельных алгоритмов. В таблице представлено время вычислений £дгр в секундах в зависимости от количества используемых ядер Ир для описанного выше численного эксперимента 1 с усвоением данных наблюдений уровня океана.
Np 1 32 128 256 512 625
íjVp; С 1080 60 20 15 16 15
На рис. 2 представлен график ускорения параллельного алгоритма Б^р = ^л/^'Ир- На рис. 3 представлен график эффективности вычислений Ейдгр = Б^р/^р. На этих графиках видно, что при большом количестве (более 256) ядер ускорение становится величиной, близкой к постоянной, а эффективность падает при числе ядер более 32, что связано с увеличением числа обменов данными.
Рис. 2. Ускорение вычислений в зависимости от количества ядер
Рис. 3. Эффективность вычислений в зависимости от количества ядер
4. Выводы. Реализована параллельная схема усвоения данных в модели динамики океана HYCOM на основе стандартного фильтра Калмана. Схема позволяет оперативно усваивать большой объем данных и корректировать поля расчетных параметров в динамике на заданных шагах по времени. Показано, что предложенный алгоритм вычислительно эффективен и полученные в численных экспериментах результаты хорошо согласуются с данными наблюдений.
СПИСОК ЛИТЕРАТУРЫ
1. Valdivieso М.. Натев К.. Balmaseda М. An assossmont of air-sea heat fluxos from ocean and coupled геапа1увев // Climate DynamicB. 2015. N 10. Р. 1 26.
2. Агошков В.И.. ИпатоваВ.М.. Залесный В. Б. Задачи вариационной ассимиляции данных наблюдений для моделей общей циркуяции океана и методы их решения // Изв. РАН. Физика атмосферы и океана. 2010. 46. № G. С. 734 770.
3. Кауркин М.Н.. Ибраев Р. А.. Беляев К. П. Усвоение данных наблюдений в модели динамики океана высокого пространственного разрешения с применением методов параллельного программирования // Метеорология и гидрология. 2016. № 7. С. 47 57.
4. Кныш В. В.. Коротаев Г. К.. Мизюк А. И.. Саркисян А. С. Усвоение гидрологических полей Мирового океана // Изв. РАН. Физика атмосферы и океана. 2012. 48. № 1. С. 37 46.
5. Е ve lis en G. Data AsBimilation. the ЕпветЫе Kalman Filter. 2nd ed. Springer, 2009.
6. Беляев К.П., Танажура К. А. С., Тучкова Н.П. Сравнительный анализ экспериментов с усвоением данных дрифтеров АРГО // Океанология. 2012. 52. № 5. С. 643-653.
7. Tanajura С. A.S., Belyaev К. A. Sequential data assimilation method based on the properties of a diffusion-type process // Appl. Math. Model. 2009. 33. N 5. P. 2165-2174.
8. Kalnay E., Li H., Miyoshi Т., et al. 4D-Var or ensemble Kalman filter // Tellus. Series A: Dynamic Meteorology and Oceanography. 2007. 59. N 5. P. 758-773.
9. Bleck R., Boudra D.B. Initial testing of a numerical ocean circulation model using a hybrid quasi-isopycnal vertical coordinate // J. Phys. Oceanography. 1981. N 11. P. 755-770.
10. Bleck R. An oceanic general circulation model framed in hybrid isopycnic Cartesian coordinates // Ocean Modelling. 2002. 4. N 1. P. 55-88.
11. Griffies S.M., Biastoch A., Boning C., et al. Coordinated ocean-ice reference experiments (COREs) // Ocean Modelling. 2009. 26. N 12. P. 1-46.
12. Kalnay E. Atmospheric Modeling, Data Assimilation and Predictability. Cambridge: Cambridge University Press, 2003.
13. Xie J., Zhu J. Ensemble optimal interpolation schemes for assimilating Argo profiles into a hybrid coordinate ocean model // Ocean Modelling. 2010. 33. N 3. P. 283-298.
14. Belyaev K.P., Kuleshov A. A., Tanajura C.A.S. An application of a data assimilation method based on the diffusion stochastic process theory using altimetry data in Atlantic // Russ. J. Numer. Anal. Math. Modelling. 2016. 31. N 3. P. 137-148.
Поступила в редакцию 28.09.16
