Развитие климатических моделей на сегодняшний день связано как с уточнением используемых физических параметризаций, так и с повышением пространственного разрешения при фиксированных методах аппроксимации уравнений гидротермодинамики. Численные эксперименты [9,14] показывают, что увеличение горизонтального разрешения до 0.25 градуса в моделях атмосферы и 0.1 градуса для моделей Мирового океана позволяет существенно повысить точность моделируемых характеристик крупномасштабных циркуляций. При этом динамика морского льда, процессы перемешивания в океане, мелкомасштабная динамика атмосферы чувствительны к горизонтальному разрешению моделей Земной системы. Увеличение вертикального разрешения в моделях также значимо для климатических исследований. В частности, эксперименты [3] с моделью общей циркуляции атмосферы Института вычислительной математики РАН с 80 расчетными уровнями по вертикали позволили воспроизвести квазидвухлетние колебания зонального ветра в экваториальной стратосфере близкие к наблюдаемым. По этим причинам воспроизведение климата при высоком пространственно-временном разрешении представляет значительный интерес (например, можно отметить планируемую программу сравнения моделей высокого разрешения HighResMIP в рамках CMPI6 -Climate Model Intercomparison Project).
Численное моделирование климата Земли относится к задачам наиболее требовательным к вычислительным ресурсам. Возможность проведения экспериментов с высоким пространственным разрешением напрямую связана с эффективностью программной реализации климатической модели на современных параллельных вычислительных системах. Особый интерес представляет возможность ускорения расчетов за счет переноса вычислений на архитектуру перспективных сопроцессоров -графических карт и процессоров Intel Xeon Phi [11,13].
Важное направление исследований связано с изучением способов переноса расчета динамики крупномасштабных моделей на архитектуру ускорителей. В этой связи можно отметить и работы [4,17], в которых показана возможность ускорения до 30-50 раз и масштабируемость по сопроцессорам вычислительных методов решения гидродинамических задач -системы уравнений Навье-Стокса или осредненных уравнений, дополненных турбулентным замыканием. В статье [7] приведены схожие оценки ускорения для численного метода решения баротропного уравнения вихря. Перенос блока динамики негидростатической модели атмосферы на архитектуру графических процессоров в работе [16] позволил уменьшить время счета до 80 раз за счет высокой оптимизации операций на сопроцессорах. Реализация динамической части гидростатической модели HILRAM (High Resolution Limited Area Model) на сопроцессорах для ускорения вычислений до 55 раз (относительно одного ядра центрального процессора) [19] потребовала разработки специализированных алгоритмов перекрытия вычислений с пересылкой данных между центральным процессором и графическим устройством. Высокая производительность в приведенных работах достигается за счет трудоемкой оптимизации -применение автоматических подходов к генерации программного кода, исполняемого на графических процессорах, снижает быстродействие алгоритмов в несколько раз [10].
Основной задачей на сегодняшний день применительно к климатическим исследованиями представляется разработка алгоритмов, обладающих высокой масштабируемостью на вычислительных кластерах, состоящих из сопроцессоров. В частности, несмотря на высокую производительность расчетов на одной видеокарте, программная реализация, предложенная в статье [19] масштабируется на 4 процессорах с ускорением в 2.7 раза при достаточно большой размерности вычислительной сетки (более чем 6×10 7 узлов).
В настоящей работе эффективность переноса вычислений на архитектуру сопроцессоров рассматривается на основе компонентов модели обшей циркуляции атмосферы Института вычислительной математики РАН [2]. Модель атмосферы основана на решении системы нелинейных уравнений гидротермодинамики на сфере с использованием -координаты. Для численной аппроксимации рассматривается дискретизация уравнений по горизонтали на сдвинутой С-сетки [6], регулярной по долготе и широте. Разностная схема имеет второй порядок точности, а для интегрирования по времени используется полунеявный метод, основанный на выделении линейных членов, связанных с распространением гравитационных волн. Применение полунеявного подхода позволяет значительно ослабить ограничение на шаг по времени по сравнению с явной схемой, однако приводит к необходимости решения эллиптического уравнения на каждом шаге по времени. Численный метод нахождения решения данного уравнения является наиболее трудоемким компонентом расчета динамики атмосферы и, как правило, ограничивает масштабируемость модели на параллельных вычислительных системах [5].
Система уравнений для основных переменных (скорости, температуры и давления) в модели атмосферы можно свести к эллиптическому уравнению следующего вида:
Предобусловленный метод сопряженных градиентов применяется для решения системы (2). В качестве предобусловливателя предлагается использовать геометрический многосеточный метод с так называемым V-циклом. Для сглаживания высокочастотной компоненты ошибки приближения на исходной и огрубленных сетках используется симметричный метод Гаусса-Зейделя для красно-черного упорядочивания узлов сетки. Основной проблемой изложенного подхода для рассматриваемой задачи (2) является сильная анизотропия вычислительной сетки вблизи полюсов и усиление данного эффекта при переходе к огрубленной как по долготе, так и по широте сетке [18]. В настоящей реализации многосеточного метода используется следующий подход [8] к построению последовательности сеток -увеличение шага вдоль кругов широты в сочетании с условным огрублением вдоль меридианов вблизи экватора. Ширина огрубления вдоль меридианов при этом увеличивается при переходе на наиболее грубые сетки и оценивается по величине внедиагональных элементов матрицы. Процедура построения последовательности сеток в этом случае усложняется. Тем не менее поскольку сетка фиксирована то дополнительные вычислительные затраты приходятся только на этап инициализации расчетов. В данном случае оправдано применение точечного симметричного метода Гаусса-Зейделя при красно-черном упорядочивании узлов сетки. Данный подход представляется наиболее эффективным для реализации на параллельных архитектурах, в отличие от альтернативы «разрешения» анизотропии в многосеточном методе при дискретизации на сфере -применения линейных итерационных методов, предполагающих обращение трехдиагональных систем [18]. Для повышения масштабируемости численного метода системы (2) решаются совместно для всех уровней предобусловленным методом сопряженных градиентов. При этом по достижении сходимости на некотором уровне он исключается из расчета. В ходе работы алгоритма последовательно исключаются нижние уровни за счет монотонного уменьшения диагональных коэффициентов (и, как следствие, увеличения числа итераций). В итоге алгоритм сводится к выполнению операций для «двумерных» данных только на самом верхнем уровне. Данная методика позволяет объединить малые, как правило, неэффективные обмены данных между параллельными процессами для всех задействованных на текущем этапе уровней, а вместе с тем уменьшить накладные расходы, связанные с инициализацией большого числа коротких сообщений. Более того, алгоритм позволяет объединить и коллективные операции взаимодействия для расчета скалярных произведений. Предложенный подход позволяет значительно повысить масштабируемость на параллельных архитектурах и сохраняет арифметическую сложность алгоритма.
Дальнейшая модификация многосеточного метода состоит в разделении уровней по блокам в каждом из которых вычисления объединяются по вертикали. Ускорение сходимости итерационных методов на верхних уровнях приводит к тому, что сохранение фиксированной длины последовательности вложенных сеток для всех уровней в многосеточном методе является избыточным. Разделение вычислений на независимые блоки и последовательное уменьшение числа грубых сеток по вертикали позволяет использовать указанное свойство для уменьшения общего времени счета и числа обменов данными между параллельными процессами. Размерность блоков можно оценить по величине диагональных компонент, что обеспечивает близость числа итераций достаточных для сходимости.
Программная реализация на графических картах предполагает наличие двух уровней параллелизма в программе: распределение данных между отдельными графическими процессорами и внутренний параллелизм отдельных устройств. Библиотека MPI используется для организации обменов данными между устройствами, а для организации вычислений на графических процессорах применяется технология программирования CUDA [1].
На каждом шаге интегрирования по времени правая часть системы, расположенная в памяти центрального процессора, группируется и копируются в память графического устройства. Итерационный метод, реализованный на графической карте, используется для нахождения решения эллиптического уравнения. По достижении сходимости итерационного метода вектор решения, расположенный в памяти сопроцессора, копируется в основную память центрального процессора для продолжения расчета основной динамики. Блочная структура по вертикали итерационного метода позволяет выполнять дополнительные медленные обмены данными между центральным процессором и графической картой асинхронно с расчетами.
Для оценки эффективности реализации итерационного метода на графических процессорах были проведены вычислительные эксперименты на системе Суперкомпьютерного центра МГУ «Ломоносов». На рисунке 1 приведено ускорение расчета на графической карте относительно четырех ядер центрального процессора при горизонтальном разрешении в 1 градус и 128 уровнях при изменении размерности вертикальных блоков. Влияние перекрытия обменов данными между центральным и графическим процессором и вычислениями показано на рисунке пунктирной линией. Ускорение после начального этапа практически постоянно и незначительно увеличивается с числом уровней в блоке, что характерно для реализации многосеточного метода на графических процессорах в связи с последовательным уменьшением размерности задачи при огрублении сетки. На рисунке 2 показана масштабируемость реализации итерационного метода на графических процессорах при горизонтальном разрешении в 1 градус для фиксированного вертикального блока, содержащего 32 уровня.
В работе рассмотрена реализация итерационного метода решения эллиптического уравнения, возникающего за счет применения полунеявных схем аппроксимации уравнений термогидродинамики по времени в моделях общей циркуляции атмосферы, на массивнопараллельных вычислительных системах. Полученные оценки показывают, что для задач относительно большой пространственной размерности перенос вычислений на архитектуру ускорителей является предпочтительным и позволяет уменьшить время расчета. Решение климатических задач даже при высоком горизонтальном разрешении до 0.25 градуса, где время интегрирования системы по времени и масштабируемость модели являются основными параметрами, требует дальнейшей модификации численных методов и алгоритмов. Необходима разработка итерационных методов, в большей степени учитывающих особенности архитектуры массивно-параллельных систем и обеспечивающих оптимальную вычислительную сложность [12,20]. Эффективное использование возможностей как графических карт, так и процессоров Intel Xeon Phi предполагает реализацию всех элементов модели атмосферы на данных устройствах. В этой связи более предпочтительными представляются перспективные системы, использующие многоядерную архитектуру на основе процессоров Intel Xeon Phi в качестве базового вычислительного элемента.
The efficiency of the implementation of iterative methods for the solution of elliptic equations in atmospheric general circulation models on massively parallel systems
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.