Научно-исследовательский вычислительный центр МГУ имени М.В. Ломоносова

Abstract

В работе изучается возможность эффективной программной реализации отдельных компонент модели общей циркуляции атмосферы на современных массивно-параллельных вычислительных архитектурах. Основное внимание уделяется методам решения эллиптических уравнений, возникающих при использовании полунеявных схем аппроксимации уравнений гидротермодинамики по времени. Результаты представлены для вычислительных архитектур на основе процессоров Intel Xeon Phi и графических процессоров Nvidia.

-

Эффективность реализации итерационных методов решения эллиптических уравнений в моделях общей циркуляции атмосферы на массивно-параллельных вычислительных системах* время счета до 80 раз за счет высокой оптимизации операций на сопроцессорах. Реализация динамической части гидростатической модели HILRAM (High Resolution Limited Area Model) на сопроцессорах для ускорения вычислений до 55 раз (относительно одного ядра центрального процессора) [ 19 ] потребовала разработки специализированных алгоритмов перекрытия вычислений с пересылкой данных между центральным процессором и графическим устройством. Высокая производительность в приведенных работах достигается за счет трудоемкой оптимизации – применение автоматических подходов к генерации программного кода, исполняемого на графических процессорах, снижает быстродействие алгоритмов в несколько раз [ 10 ].

Основной задачей на сегодняшний день применительно к климатическим исследованиями представляется разработка алгоритмов, обладающих высокой масштабируемостью на вычислительных кластерах, состоящих из сопроцессоров. В частности, несмотря на высокую производительность расчетов на одной видеокарте, программная реализация, предложенная в статье [ 19 ] масштабируется на 4 процессорах с ускорением в 2.7 раза при достаточно большой размерности вычислительной сетки (более чем 6×107 узлов).

В настоящей работе эффективность переноса вычислений на архитектуру сопроцессоров рассматривается на основе компонентов модели обшей циркуляции атмосферы Института вычислительной математики РАН [2]. Модель атмосферы основана на решении системы нелинейных уравнений гидротермодинамики на сфере с использованием  -координаты. Для численной аппроксимации рассматривается дискретизация уравнений по горизонтали на сдвинутой С-сетки [ 6 ], регулярной по долготе и широте. Разностная схема имеет второй порядок точности, а для интегрирования по времени используется полунеявный метод, основанный на выделении линейных членов, связанных с распространением гравитационных волн. Применение полунеявного подхода позволяет значительно ослабить ограничение на шаг по времени по сравнению с явной схемой, однако приводит к необходимости решения эллиптического уравнения на каждом шаге по времени. Численный метод нахождения решения данного уравнения является наиболее трудоемким компонентом расчета динамики атмосферы и, как правило, ограничивает масштабируемость модели на параллельных вычислительных системах [5]. 2. Численный метод решения эллиптического уравнения в модели общей циркуляции атмосферы Система уравнений для основных переменных (скорости, температуры и давления) в модели атмосферы можно свести к эллиптическому уравнению следующего вида: 1 2  1  2 ttd  t  Bh 

 2 ttd   RHS, где  tt (t)  (t  t)  2 (t)  (t  t) , t - шаг по времени, h - дискретный оператор Лапласа, d - горизонтальная дивергенция, B - матрица, полученная при объединении исходных уравнений, RHS - правая часть уравнения, учитывающая явные динамические тенденции. Трехмерное уравнение (1) с помощью метода разделения переменных и определения матрицы преобразования S преобразуется к набору двумерных уравнений схожих с уравнением Гельмгольца на каждом модельном уровне m  1, NL : 2hu  cu  g , c 

1 t 2 m

 0 , 1 где u  u( , )  2 ttd , g  g( , )  cS 1(RHS)S и m - собственные значения матрицы B . Предобусловленный метод сопряженных градиентов применяется для решения системы (2). В качестве предобусловливателя предлагается использовать геометрический многосеточный метод с так называемым V-циклом. Для сглаживания высокочастотной компоненты ошибки приближения на исходной и огрубленных сетках используется симметричный метод ГауссаЗейделя для красно-черного упорядочивания узлов сетки. Основной проблемой изложенного подхода для рассматриваемой задачи (2) является сильная анизотропия вычислительной сетки вблизи полюсов и усиление данного эффекта при переходе к огрубленной как по долготе, так и по широте сетке [ 18 ]. В настоящей реализации многосеточного метода используется следующий подход [ 8 ] к построению последовательности сеток – увеличение шага вдоль кругов широты в сочетании с условным огрублением вдоль меридианов вблизи экватора. Ширина огрубления вдоль меридианов при этом увеличивается при переходе на наиболее грубые сетки и оценивается по величине внедиагональных элементов матрицы. Процедура построения последовательности сеток в этом случае усложняется. Тем не менее поскольку сетка фиксирована то дополнительные вычислительные затраты приходятся только на этап инициализации расчетов. В данном случае оправдано применение точечного симметричного метода Гаусса-Зейделя при красно-черном упорядочивании узлов сетки. Данный подход представляется наиболее эффективным для реализации на параллельных архитектурах, в отличие от альтернативы «разрешения» анизотропии в многосеточном методе при дискретизации на сфере – применения линейных итерационных методов, предполагающих обращение трехдиагональных систем [ 18 ].

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

Дальнейшая модификация многосеточного метода состоит в разделении уровней по блокам в каждом из которых вычисления объединяются по вертикали. Ускорение сходимости итерационных методов на верхних уровнях приводит к тому, что сохранение фиксированной длины последовательности вложенных сеток для всех уровней в многосеточном методе является избыточным. Разделение вычислений на независимые блоки и последовательное уменьшение числа грубых сеток по вертикали позволяет использовать указанное свойство для уменьшения общего времени счета и числа обменов данными между параллельными процессами. Размерность блоков можно оценить по величине диагональных компонент, что обеспечивает близость числа итераций достаточных для сходимости. 3. Программная реализация итерационного метода на массивнопараллельных вычислительных системах Программная реализация на графических картах предполагает наличие двух уровней параллелизма в программе: распределение данных между отдельными графическими процессорами и внутренний параллелизм отдельных устройств. Библиотека MPI используется для организации обменов данными между устройствами, а для организации вычислений на графических процессорах применяется технология программирования CUDA [ 1 ].

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

Для оценки эффективности реализации итерационного метода на графических процессорах были проведены вычислительные эксперименты на системе Суперкомпьютерного центра МГУ «Ломоносов». На рисунке 1 приведено ускорение расчета на графической карте относительно четырех ядер центрального процессора при горизонтальном разрешении в 1 градус и 128 уровнях при изменении размерности вертикальных блоков. Рис. 1. Ускорение итерационного метода решения эллиптического уравнения на графическом процессоре относительно 4-х ядер центрального процессора. Рис. 2. Масштабируемость реализации итерационного метода на графических процессорах. Полученные оценки свидетельствуют о том, что для итерационного метода решения эллиптического уравнения, выполнение вычислений на архитектуре графических процессоров позволяет ускорить расчеты. За счет разделения переменных, и, как следствие, независимости систем уравнений на различных вертикальных уровнях удается перекрыть дополнительные «медленные» операции пересылки данных между центральным процессором и сопроцессором с нахождением решения дискретного двумерного эллиптического уравнения. Вместе с тем приведенная масштабируемость алгоритма представляется недостаточной и требует дальнейшей оптимизации обменов между сопроцессорами и их перекрытия с вычислениями.

В качестве архитектуры на основе процессоров Intel Xeon Phi, рассматривалась вычислительная система МВС-10П МП на базе модулей RSC Petastream. Программная реализация итерационного метода использует библиотеку MPI для организации обменов данными между устройствами и технологию OpenMP для распределения вычислений между ядрами (нитями) одной карты. К достоинствам вычислительной архитектуры можно отнести то, что реализация обменов данными между устройствами не требует «явного» копирования данных, а скорость передачи сообщений сопоставима с традиционной архитектурой [ 15 ]. Однако эффективное выполнение вычислений на Intel Xeon Phi требует тщательной оптимизации, и в частности возможности задействовать до 244 нитей одного устройства. На рисунке 3 показана масштабируемость многосеточного метода по числу ядер (4 нити на каждое ядро) одного процессора Intel Xeon Phi при различной размерности сетки. Для достижения производительности сопоставимой с 2-мя процессорами Intel Xeon E5-2690 размерность задачи на одной карте Intel Xeon Phi 7120D должна превышать 106 узлов. Рис. 3. Масштабируемость реализации итерационного метода на процессоре Intel Xeon Phi 7120D. 4. Заключение

В работе рассмотрена реализация итерационного метода решения эллиптического уравнения, возникающего за счет применения полунеявных схем аппроксимации уравнений термогидродинамики по времени в моделях общей циркуляции атмосферы, на массивнопараллельных вычислительных системах. Полученные оценки показывают, что для задач относительно большой пространственной размерности перенос вычислений на архитектуру ускорителей является предпочтительным и позволяет уменьшить время расчета. Решение климатических задач даже при высоком горизонтальном разрешении до 0.25 градуса, где время интегрирования системы по времени и масштабируемость модели являются основными параметрами, требует дальнейшей модификации численных методов и алгоритмов. Необходима разработка итерационных методов, в большей степени учитывающих особенности архитектуры массивно-параллельных систем и обеспечивающих оптимальную вычислительную сложность [ 12, 20 ]. Эффективное использование возможностей как графических карт, так и процессоров Intel Xeon Phi предполагает реализацию всех элементов модели атмосферы на данных устройствах. В этой связи более предпочтительными представляются перспективные системы, использующие многоядерную архитектуру на основе процессоров Intel Xeon Phi в качестве базового вычислительного элемента. Литература 3. Кулямин Д.В., Володин Е.М., Дымников В.П. Моделирование квазидвухлетних колебаний зонального ветра в экваториальной стратосфере // Известия РАН. Физика атмосферы и океана. 2009. № 1. С. 43-61. 4. Мортиков Е.В. Применение графических процессоров для численного моделирования течения вязкой несжимаемой жидкости в областях сложной конфигурации методом погруженной границы // Вычислительные методы и программирование. 2012. Т. 13, № 1. С. 177-191. The efficiency of the implementation of iterative methods for the solution of elliptic equations in atmospheric general circulation models on massively parallel systems Evgeny Mortikov Keywords: atmospheric general circulation model, massively parallel systems, Intel Xeon Phi, Nvidia GPU In this paper we study the possibility of the efficient implementation of atmospheric general circulation model and its particular components on modern massively parallel architectures. Main emphasis is placed on methods for solving elliptic equations, arising from applying semi-implicit in time approximations in thermo-hydrodynamic equations. Results are presented for computational architectures based on Intel Xeon Phi processors, as well as Nvidia graphic cards.

References

1. Боресков А.В. , Харламов А .А., Марковский Н .Д., Микушин Д .Н., Мортиков Е .В., Мыльцев А .А., Сахарных Н .А., Фролов В .А. Параллельные вычисления на GPU . Архитектура и программная модель CUDA: учебное пособие . М.: Изд-во Московского университета , 2012 . 336 C.

6. Arakawa A. , Lamb V.R. Computational design of the basic dynamical processes of the UCLA general circulation model // Methods Comput . Phys . 1977 . V. 17. P. 173 - 265 .

7. Bleichrodt F. , Bisseling R.H. , Dijkstra H.A. Accelerating a barotropic ocean model using a GPU // Ocean Modelling . 2012 . V. 41. P. 16 - 21 .

8. Buckeridge S. , Scheichl R . Parallel geometric multigrid for global weather prediction // Numerical Linear Algebra with Applications . 2010 . V. 17. P. 325 - 342 .

9. Demory M.-E. , Vidale P.L. , Roberts M.J. , Berrisford P. , Strachan J. , Schiemann R. , Mizielinski M.S. The role of horizontal resolution in simulating drivers of the global hydrological cycle // Clim . Dyn. 2014 . V. 42. P. 2201 - 2225 .

10. Govett M. Using GPUs for weather and climate models . Earth System Prediction Capability (ESPC) Workshop , New Orleans, USA. 2010 .

11. Kelly R. GPU computing for atmospheric modeling // Comp . Sci. Eng . 2010 . V. 12, N 4. P. 26 - 33 .

12. Li R. , Saad Y. GPU-accelerated preconditioned iterative linear solvers // Journal of Supercomputing . 2013 . V. 63, N 2. P. 443 - 466 .

13. Michalakes J. , Vachharajani M. GPU acceleration of numerical weather prediction // Parallel Processing Letters . 2008 . V. 18, N 4. P. 531 - 548 .

14. Roberts M.J. , Clayton A. , Demory M.-E . et al. J. Impact of resolution on the Tropical Pacific Circulation in a matrix of coupled modes / / J. Climate. 2009 . V. 22. P. 2541 - 2556 .

15. Semin A. , Druzhinin E. , Mironov V. , Shmelev A. , Moskovsky A. The performance characterization of the RSC PetaStream Module // Proc. of the 29th Supercomputing Int . Conf., ISC 2014 . 2014 . P. 420 - 429 .

16. Shimokawabe T. , Aoki T. , Ishida J . GPU acceleration of meso-scale atmospheric model ASUCA // Proc. of the 9th World Congress on Comp. Mech. and 4th Asian Pacific Congress on Comp. Mech. (WCCM/APCOM 2010 ). 2010 .

17. Thibault J.C. , Senocak I. CUDA implementation of a Navier-Stokes Solver on multi-GPU desktop platforms for incompressible flows . In: 47th AIAA Aerospace Sciences Meeting , 2009 .

18. Trottenberg U. , Oosterlee C.W. , Schüller A. Multigrid . London: Academic Press, 2001 . 631 P.

19. Vu V.T. , Cats G. , Wolters L. Graphics processing unit optimizations for the dynamics of the HIRLAM weather forecast model // Concurrency and Computation: Practice and Experience. 2013 . V. 25, N 10. P. 1376 - 1393 .

20. Yang X. , Mittal R . Acceleration of the Jacobi iterative method by factors exceeding 100 using scheduled relaxation // J. Comp. Phys. 2014 . V. 274. P. 695 - 708 .