Использование технологии CUDA для решения задачи переноса теплового и солнечного излучения в атмосфере Земли Текст научной статьи по специальности «Математика»

DOI: 10.25702/KSC.2307-5252.2018.9.5.183-191 УДК 551.593.5.531.521.3

К. Г. Орлов, И. В. Мингалев, Е. А. Федотова

ИСПОЛЬЗОВАНИЕ ТЕХНОЛОГИИ CUDA ДЛЯ РЕШЕНИЯ ЗАДАЧИ ПЕРЕНОСА ТЕПЛОВОГО И СОЛНЕЧНОГО ИЗЛУЧЕНИЯ В АТМОСФЕРЕ ЗЕМЛИ

Аннотация

Быстрое развитие вычислительной техники позволяет проводить расчеты не только на многоядерных серверах или кластерах, но и графических ускорителях с использованием технологии США. В данной работе приведено описание реализованных авторами на видеокартах алгоритмов вычисления коэффициента молекулярного поглощения атмосферного газа на сетке с высоким разрешением по частоте, построения параметризации молекулярного поглощения и решения в земной атмосфере уравнений переноса собственного и солнечного излучений. Также приведена оценка быстродействия алгоритмов, реализованных с использованием графических ускорителей (технология CUDA) и центрального процессора (технология OpenMP).

Ключевые слова:

собственное и солнечное излучение атмосферы, параметризация молекулярного поглощения, технология CUDA.

K. G. Orlov, I. V. Mingalev, E. A. Fedotova

THE CUDA TECHNOLOGY AS APPLIED TO THERMAL AND SOLAR RADIATION TRANSFER IN THE EARTH'S ATMOSPHERE

Abstract

Rapid development of computers allows us to make calculations not only on multicore servers or clusters, but also on graphics accelerators by the CUDA technology. The paper presents the algorithms used by the authors on video cards to calculate the molecular absorption coefficient of the atmospheric gas on a grid of high frequency resolution, to construct molecular absorption parameterization, and to solve transport equations of the thermal and solar radiation in the Earth's atmosphere. The speed of algorithms has been evaluated and the accuracy of the calculations results has been analyzed.

Keywords:

thermal and solar atmosphere radiation, molecular absorption parameterization, CUDA technology.

Введение

Расчет поля собственного ИК-излучения атмосферы необходим для решения различных физических задач, в частности, для интерпретации данных дистанционного зондирования атмосферы, для расчета нагрева атмосферы собственным ИК-излучением при моделировании общей циркуляции атмосферы Земли. Для достижения точности 1% и лучше при расчетах интенсивности излучения разрешение по частоте должно составлять примерно 0,001 см-1. Расчеты с таким высоким разрешением по частоте называются эталонными расчетами (Line-by-Line). Они требуют очень больших вычислительных затрат и

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

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

Отметим, что на высотах 0-70 км нужно учитывать изменение газового состава атмосферы с высотой. Ниже высоты 15 км вклад водяного пара в коэффициент молекулярного поглощения является существенным, а вклад озона мал. Выше высоты 20 км снижается роль водяного пара и возрастает вклад озона. Поэтому спектры поглощения на малых и больших высотах не коррелируют. Авторами предложен новый алгоритм построения параметризации молекулярного поглощения, который отличается от метода k-корреляции, учитывает изменение газового состава атмосферы с высотой, не требует проводить подгоночные расчеты для каждого модельного канала и относительно прост в программной реализации по сравнению с другими алгоритмами. Созданная авторами программа, реализующая этот алгоритм, позволяет менять число модельных каналов параметризации в широких пределах, а предложенная параметризация молекулярного поглощения в атмосфере Земли обладает хорошей точностью в диапазоне высот 0-76 км, как при отсутствии облачных слоев, так и при их наличии.

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

1. Вычисление коэффициента молекулярного поглощения атмосферного газа и других оптических параметров для узких спектральных каналов.

2. Группировка узких спектральных каналов в широкие модельные каналы и определение оптических параметров для этих модельных каналов.

3. Численное решение уравнений переноса собственного и солнечного излучения в каждом модельном канале.

Ниже приведено описание алгоритмов решения этих подзадач. Эти алгоритмы были разработаны авторами данной работы и реализованы на графических ускорителях фирмы NVIDIA с использованием технологии CUDA.

Алгоритм вычисления коэффициента молекулярного поглощения

В данной работе коэффициент поглощения для узких спектральных каналов с шириной 0.001 см"1 в атмосфере Земли вычисляется как сумма коэффициентов поглощения семи газов, среди которых метан, углекислый газ, водяной пар, азот, озон, кислород и окись азота по формуле:

где а - сорт молекул, па - концентрация молекул сорта а , Оа (V) - сечение

поглощения молекулы сорта а . Коэффициенты молекулярного поглощения этих семи газов рассчитываются с использованием спектроскопической базы данных НИКЛЫ 2012 [1] по стандартной теории. Согласно этой теории на каждой частоте вклады в коэффициент поглощения различных линий поглощения суммируются при условии обрезания крыльев линий на расстоянии 25 см"1 от центра линии. Для водяного пара и углекислого газа к этой сумме прибавляется дополнительное слагаемое, называемое сечением континуального поглощения. Зависимость сечения континуального поглощения от частоты излучения, а также давления и температуры воздуха задают с помощью эмпирических моделей, построенных в результате сравнения экспериментальных данных и теоретических расчетов. Мы использовали эмпирическую модель МТ СКО [2]. Главная трудность вычисления коэффициента молекулярного поглощения состоит в том, что в заданном частотном интервале число линий поглощения для некоторых газов может превышать 100 000, причем вклад каждой линии поглощения на каждой высоте (с учетом обрезания крыльев линий) необходимо вычислять в 50 000 узлах сетки по частоте. Результатом этих расчетов является двумерный массив коэффициентов молекулярного поглощения, одно измерение которого - по частоте, а другое - по высоте. Размерность по частоте может составлять более 4 000 000 узлов сетки, в зависимости от рассчитываемого частотного диапазона, а размерность по высоте - около 400 узлов сетки. Проведение таких расчетов без использования технологий распараллеливания требует больших временных затрат.

Авторами представленной работы было реализовано несколько алгоритмов расчета коэффициента молекулярного поглощения с использованием параллельных вычислений на графических ускорителях. Среди этих способов наиболее эффективным по скорости вычислений оказался тот, в котором каждая вычислительная нить (СИБЛ - ядро видеокарты) рассчитывает вклад в коэффициент поглощения одной линии поглощения заданного газа в интервале частот шириной 50 см"1 на заданной высоте. Таким образом, вклады всех линий в коэффициент молекулярного поглощения для одного газа на заданной высоте рассчитываются параллельно. Этот способ предполагает создание внешнего цикла по высоте, на каждой итерации которого вызывается функция вычисления сечений поглощения молекулы сорта а, реализованная на видеокарте (СИБЛ-функция). Время, которое тратится на вызов СИБЛ-функции мало по сравнению со временем работы этой функции, поэтому многократные вызовы СИБЛ-функции незначительно влияют на скорость работы программы.

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

(1)

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

Каждая вычислительная нить для своей линии поглощения сначала вычисляет интенсивность этой линии по формуле:

Sa (T) = S* (Tref )

Qa(Tref )exp(-C2Ean /T)(1 -exp(-C2Ecd /T)) Qa (T) exp(-C2Ean /Tf )(1 - exp(-C2Ecd /Tref ))

(2)

где Тг^ = 296К - нормальная температура, <2а (Т) - произведение вращательной и колебательной статистических сумм молекулы сорта а , Еот - энергия (в см-1) нижнего уровня перехода этой молекулы, Еш. - энергия (в см-1) перехода, соответствующая линии поглощения с номером I этой молекулы, С = ке / = 1.438769 см ■ К - вторая радиационная постоянная, в которой к -постоянная Планка, е - скорость света, кв - постоянная Больцмана. Далее вычисляются полуширина линии контура Лоренца и параметр доплеровской ширины линии по формулам:

(

aL =

7c

self

fp Л

a

\ P r l V ref j

г p - p ^ гy'

+ 7°

/ a

P- P

P r

v ref J J

ref

T

V 1

aD =

V

ai С

2RT

Uc

(3)

где Уа^ - коэффициент самоуширения линии с номером i молекул сорта а , УаГ - коэффициент уширения воздухом этой линии, /За - коэффициент

температурной зависимости этой линии,

P,

- парциальное давление молекул

сорта а, р = 1 атм, Р - полное давление атмосферного газа, Я - универсальная

газовая постоянная, ца - молярная масса молекул сорта а , Уа - частота центра линии поглощения молекул сорта а .

Рис. 1. Схематичное изображение контуров Фойгта (красные кривые). Черными стрелками изображены центры линий поглощения.

После этого вычисляются значения контура Фойгта (рис. 1) в узлах сетки по частоте. Для расчета контура Фойгта используется функция, описанная в работе М. Кунтца [3] и модернизированная авторами данной работы для проведения параллельных расчетов на видеокарте. Контур Фойгта рассчитывается в диапазоне от

Vo-25 см"1 до Vo+25 см"1, где Vo - центр линии поглощения. Каждая нить выполняет цикл от V0 — 25 см-1 до v0 + 25 см-1 с шагом 0.001 см-1, рассчитывая значения контура

Фойгта в 50 000 узлах сетки по частоте.

Точность расчета коэффициента поглощения на видеокарте оказалась сравнимой с точностью расчета на центральном процессоре (с использованием технологии распараллеливания Open MP). При этом скорость расчетов на GPU в десятки или сотни раз (в зависимости от типа видеокарты) превышает скорость расчета на CPU.

Алгоритм построения параметризации молекулярного поглощения и других оптических параметров

Для построения параметризации весь участок спектра разбивается на интервалы шириной от 50 см-1 до 4000 см-1 (в зависимости от рассчитываемого частного интервала), которые далее будем называть интервалами осреднения (рис. 2). В каждом интервале осреднения узкие спектральные каналы по различным алгоритмам объединяются в широкие модельные каналы, которые еще называют носителями резонансов [4]. В данной работе представлен новый алгоритм объединения узких каналов в широкие модельные каналы, позволяющий учитывать изменение газового состава атмосферы с высотой в диапазоне высот 0 - 76 км. Считается, что ниже высоты 15 км вклад водяного пара в коэффициент молекулярного поглощения является существенным, а вклад озона - мал, а выше высоты 20 км снижается роль водяного пара и возрастает вклад озона.

Основная идея нового алгоритма состоит в том, чтобы построение модельных каналов проводить в два этапа. На первом этапе выбирается высота первой сортировки в диапазоне 5 - 17 км, чтобы учесть линии поглощения водяного пара. Все узкие

каналы из интервала осреднения разбиваются на N групп так, чтобы коэффициенты молекулярного поглощения узких каналов внутри каждой группы были достаточно близки между собой на этой высоте, а также на высотах 0-20 км. На втором этапе выбирается высота второй сортировки в диапазоне 40-55 км, чтобы учесть линии поглощения озона. Каждая полученная после первой сортировки группа узких каналов

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

узких каналов внутри каждой подгруппы были достаточно близки между собой на этой высоте и на высотах 0-76 км. Узкие каналы, вошедшие в одну подгруппу, объединяются в один модельный канал. В итоге получается NN модельных каналов

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

Расчет значений оптических параметров, среди которых коэффициенты полинома Лежандра, альбедо однократного рассеяния, коэффициент экстинкции и функция источника первичного солнечного излучения, в модельных каналах требует большого объема вычислений. Для ускорения расчетов авторами данной работы был разработан и реализован алгоритм вычисления значений оптических параметров в модельных каналах на видеокарте с использованием технологии CUDA. Этот алгоритм является эффективным как по скорости вычислений, так и

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

10 СМ"1 2000 СМ"1

40 СМ"1 50 СМ"1 50 СМ"1 50 СМ"1

Рис. 2. Разбиение участка спектра 10 - 2000 см-1 на интервалы осреднения по 40 - 50 см-1.

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

Алгоритм решения уравнения переноса излучения

Поле собственного и солнечного излучения горизонтально однородной атмосферы описывается 1 -мерными по пространству уравнениями переноса излучения с соответствующими граничными условиями. Для решения этих уравнений используется модификация метода дискретных ординат, детально описанная в работе [5] и имеющая две особенности. Первая из них заключается в том, что расчетная сетка по зенитным углам может быть произвольной. Вторая заключается в том, что для решения возникающей в методе дискретных ординат системы линейных алгебраических уравнений используется метод матричной прогонки. Этот метод является точным и максимально использует структуру матрицы коэффициентов системы для уменьшения объема вычислений. Он является более экономичным и более простым в реализации, чем примененный в программе DISORT метод решения, использующий вычисление собственных чисел и векторов матрицы коэффициентов линейной системы, которая имеет большую размерность. В случае наличия в атмосфере слоев с сильным рассеянием и слабым поглощением (например, слои облаков на Венере и на Земле) итерационные методы могут сходиться медленно и требовать выполнения большого числа итераций для достижения приемлемой точности решения. В этом случае предложенный в данной работе [5] метод имеет преимущество в точности и скорости расчета.

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

" ^и) = -I(*,и) + (1 - ®(г))В(Т(г),у) + S[I](и) , (4)

о (г) аг

а уравнение переноса солнечного излучения имеет вид

и еП>]г (г, и( = -(и, + ^ х(г, , и, р))^о (г) + §[18г ](г, и, р), о (г) аг 4ж

где о(г) и с(г) — коэффициент экстинкции и альбедо 1-кратного рассеяния атмосферного газа на высоте г для излучения с частотой у , £[I](г, и) — перенормированная плотность источника рассеянного излучения с интенсивностью I(г, и) и зенитным углом, имеющим косинус и , %(г, V) — индикатриса рассеяния.

В работе [5] описан метод решения уравнения (4). Уравнение для солнечного излучения решается аналогично. Уравнение переноса излучения вместе с граничными условиями в случае горизонтально однородной атмосферы после дискретизации по оптической толщине и зенитному углу сводится к системе линейных алгебраических уравнений (5) большой размерности с блочно диагональной матрицей коэффициентов [5].

С I - В I = Е

С0 Аг ,0 В0 Аг ,1 Е0

-+ С^г* -В^ГМ1 = Е , при к = 1, ... , М-1 (5)

— АМ^г ,М-1 + CMIr ,М = ЕМ

здесь А, В, С - матрицы размера М х М .

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

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

Результатом расчетов является трехмерный массив значений интенсивности излучения, одно измерение которого - по зенитным углам, другое - по высоте и третье - по модельным каналам.

Оценка быстродействия

Авторами данной работы был создан комплекс программ, в котором реализованы описанные выше алгоритмы. Этот комплекс программ был реализован как на центральном процессоре с использованием технологии распараллеливания Open MP, так и на видеокартах NVIDIA с использованием технологии CUDA.

Таблица 1. Машинное время, затраченное на расчет собственного излучения атмосферы Земли в диапазоне частот от 10 до 2000 см-1 в интервале высот от

поверхности Земли до 76 км

Название устройства Время расчетов, сек

IntelCore i7 (с использованием 8-ми потоков и OpenMP) 840

NVIDIA GTX Titan 45

NVIDIA GTX 1080 20

NVIDIA GTX 1080 Ti 15

Для оценки производительности были проведены расчеты поля собственного и солнечного излучения атмосферы Земли, как на центральном процессоре, так и на графических ускорителях фирмы NVIDIA. Для проведения тестовых расчетов использовался четырех ядерный центральный процессор IntelCore i7, с возможностью запуска восьми параллельных потоков, а также видеокарты NVIDIA GTX Titan, NVIDIA GTX 1080 и NVIDIA GTX 1080 Ti.

В таблице 1 представлено машинное время, затраченное на расчет собственного излучения атмосферы Земли в диапазоне частот от 10 до 2000 см-1 в интервале высот от поверхности Земли до 76 км. Как можно видеть, скорость расчетов на видеокартах NVIDIA GTX 1080 и NVIDIA GTX 1080 Ti существенно выше (более чем в 40 раз) скорости расчетов на центральном процессоре IntelCore i7, использующем 8 потоков и технологии распараллеливания OpenMP. Скорость расчетов на видеокарте NVIDIA GTX Titan, примерно в 2,5 раза ниже скорости расчетов на видеокарте NVIDIA GTX 1080, что объясняется архитектурой этих видеокарт. При этом результаты расчетов на видеокартах и на центральном процессоре сравнимы с относительной погрешностью 10-4 - 10-5, что позволяет делать вывод о том, что по точности видеокарты почти не уступают центральному процессору, а по скорости существенно его опережают.

Заключение

Расчет поля теплового и солнечного излучения в атмосферах Земли и других планет является очень требовательной к вычислительным ресурсам задачей. Авторами данной работы были реализованы алгоритмы решения этой задачи для земной атмосферы с использованием графических ускорителей фирмы

NVIDIA, а также центрального процессора. Анализ результатов расчетов показал, что точность вычислений на графических ускорителях и на центральном процессоре сравнима с относительной погрешностью 10-4 - 10-5, при этом скорость расчетов на видеокартах существенно выше.

Литература

1. Rothman L.S., et al. The HITRAN2012 molecular spectroscopic database // J. Quant. Spectrosc. Rad. Transfer. 2013. Vol. 130, P. 4-50.

2. Mlawer E.J., et al. Development and recent evaluation of the MT CKD model of continuum absorption. Phylosophical Transactions of the Royal Society, 2012, Vol. 370, pp. 2520-2556.

3. M. Kuntz, M. Hopfner. Effcient line-by-line calculation of absorption coefficients // Journal of Quantitative Spectroscopy and Radiative Transfer. 1999. Vol. 63. pp. 97-114.

4. Шильков А.В., Герцев М.Н. Верификация метода лебеговского осреднения // Мат. моделирование. 2015. Т. 27. № 8. C. 13-31.

5. Мингалев И.В. , Федотова Е.А., Игнатьев Н.И., Родин А.В. Новый вариант метода дискретных ординат для расчета собственного излучения в горизонтально однородной атмосфере // ЖВМ и МФ. 2015. т. 55. № 10. с. 109-123.

Сведения об авторах

Орлов Константин Геннадьевич

к.ф.-м.н., ученый секретарь Полярного геофизического института; E-mail: orlov@pgia.ru

Мингалев Игорь Викторович

д.ф.-м.н., ведущий научный сотрудник Полярного геофизического института; Email: mingalev_i@pgia.ru

Федотова Екатерина Алфеевна

младший научный сотрудник Полярного геофизического института; E-mail: godograf87@mail. ru