Моделирование газодинамических течений с использованием технологии cuda Текст научной статьи по специальности «Математика»

ЖАЛНИН Р. В., ЛЕЩАНКИНА Т. М.

МОДЕЛИРОВАНИЕ ГАЗОДИНАМИЧЕСКИХ ТЕЧЕНИЙ С ИСПОЛЬЗОВАНИЕМ ТЕХНОЛОГИИ CUDA Аннотация. В статье описан численный алгоритм решения уравнений газовой динамики с использованием технологии CUDA. Для численного эксперимента были смоделированы условия развития неустойчивости при многократном прохождении ударной волны через контактный разрыв с начальными данными Погги. Результаты численного исследования демонстрируют возможности разработанной параллельной версии программы.

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

ZHALNIN R. V., LESCHANKINA T. M.

MODELING OF GAS DYNAMIC FLOWS USING CUDA TECHNOLOGY Abstract. This article describes a numerical algorithm for solving equations of gas dynamics using CUDA technology. The Poggie's conditions of instability by repeated passage of a shock wave through a contact discontinuity were simulated for the numerical experiment. The experiment results demonstrate the capabilities of the developed parallel version of the program. Keywords: gas dynamics equations, parallel algorithm, CUDA technology.

Введение. В настоящее время большой практический интерес представляют решения задач газовой динамики. Данного вида задачи используются при проектировании самолетов и автомобилей, предсказании погоды, моделировании климата и т. п. Из-за нелинейности уравнений газовой динамики, для их решения необходимо использовать численные методы, которые на данный момент являются универсальными методами решения этих уравнений. В последние годы было проведено множество исследований в сфере разработки эффективных разностных схем высокого разрешения, которые позволяли бы осуществить математическое моделирование вышеуказанных процессов. Реализация данного рода задач на подобных сетках подразумевает использование больших вычислительных ресурсов. И появление новых технологий, ориентированных на использование в своей разработке графических ускорителей в дополнение к уже существующим позволяет использовать новые возможности для проведения интенсивных вычислений. Увеличение производительности видеокарт привело к тому, что количество вычислительных ядер в них достигло нескольких сотен. На сегодняшний день ведущей технологией разработки для работы на GPU (Graphical Processor Unit) является технология CUDA.

CUDA (Compute Unified Device Architecture) - это программная модель, включающая описание вычислительного параллелизма и иерархичной структуры памяти непосредственно в языке программирования [1]. Эта технология позволяет работать с языками программирования C и FORTRAN. С точки зрения программного обеспечения, реализация CUDA представляет собой кроссплатформенную систему компиляции и исполнения программ, части которых работают на CPU и GPU. CUDA предназначена для разработки GPGPU-приложений без привязки к графическим API и поддерживается всеми GPU NVIDIA.

1. Математическая модель и численная схема. Систему уравнений газовой динамики, записанную в переменных Эйлера, можно представить в виде:

др д(ри) d(pv) dt дх ду

д(ри) д(ри2 +р) d(puv)

dt

дх

ду

d(pv) d(puv) d(pv2 + р)

(1)

dt

дх

ду

д(ре) д((ре + р)и) д((ре + р)р) . д1 дх ду

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

показателем адиабаты у: р = (у — 1)ер. В приведенной выше системе уравнений (1)

р = р(х,1) - плотность газа, V = (и, р) - вектор скорости частиц, р = р(х,1) - давление, £ -

удельная внутренняя энергия на единицу массы, е = е + ■

u2+v2

полная энергия, у =

2 С„

показатель адиабаты, ср, су - теплоёмкость при постоянном давлении и постоянном объеме соответственно.

В векторной форме система уравнений, принимает вид:

ди дР(и) дв(и)

dt

+

дх

+

ду

= 0,

(2)

где

U =

1Р ри

pv

\ре.

/

,F(U) =

ри ри2 + р puv \(ре + р)и/

/

,G(U) =

pv puv

pv2 + p \(pe + p)v/

(3)

Для аппроксимации системы уравнений газовой динамики (2) использовалась нелинейная консервативная квазимонотонная дифференциально-разностная схема [3-5]:

Uij - UH + Fi+i/2j - Pj-!/2j | Gjj+j/2 - Gjj-i/2 _ о

т

(4)

где Fi+1/2j, Gij+1/2 - дискретные потоки.

с

р

Дискретные потоки на гранях между ячейками вычислялись как функции двух переменных по интерполированным значениям газодинамических параметров на гранях между ячейками [4; 5]:

F. i = F( Ui,Ur 1

l+2j V i+27 i+27,

G.. i = Y( Ul, 1,Ur. 1

(5)

ij+2

ij+2 ij+2,

Здесь ^¿у+1/2, ^¿у+1/2 - "левые" и " правые" значения вектора и на грани между i и i + 1 ячейками, для которой вычисляется поток ^¿+1/2у, ^¿у+1/2, (я = г, I) "левые" и " правые" значения вектора и на грани между ] и ] + 1 ячейками соответственно. Для того чтобы вычислить значения вектора и, на указанных гранях между ячейками введем новый вектор переменных Y = Y(U). Проведем его интерполяцию на грань между ячейками и, затем, пересчитаем искомое значение вектора и = на данной грани. В расчетах полагалось Y = (р, р, и, V), потоки вычислялись по схеме Лакса-Фридрихса-Русанова:

1

F1 =п l+2J 2

S+Г1

-ai( Ur i -U i

V l+27 l+27;

- a2 ( Ur i - U i

\ 1 ^+2

(6)

где

«1 = max{| u0-| + ciy, | Uj+iy | + q+iy}, a2 = max{| | + ciy, | viy+i | + Qy+i},

(7)

2. Реализация алгоритма на языке CUDA C. Графический ускоритель использовался для проведения расчетов численных потоков на границах ячеек и вычисления ограничителей.

Для проведения расчетов на GPU необходимо выполнить следующие действия: 1. Выделить память для данных и потоков на GPU. Пример: выделяется память под газодинамические параметры:

cudaMalloc((void**)&r_dev, (ntri_new_m[0])*sizeof(float)); cudaMalloc((void**)&p_dev, (ntri_new_m[0])*sizeof(float)); cudaMalloc((void**)&u_dev, (ntri_new_m[0])*sizeof(float)); cudaMalloc((void**)&v_dev, (ntri_new_m[0])*sizeof(float)); где первым параметром идет указатель, в который записывается адрес выделенной памяти, а вторым - размер выделяемой памяти в байтах.

2. Скопировать данные на GPU. Пример: газодинамические параметры (r,p,u,v) копируем на GPU по средствам использования функции:

cudaMemcpy(r_dev, r, (ntri_new_m[0])*sizeof(float), cudaMemcpyHostToDevice); cudaMemcpy(p_dev, p, (ntrinew _m[0])*sizeof(float), cudaMemcpyHostToDevice); cudaMemcpy(u_dev, u, (ntrinew _m[0])*sizeof(float), cudaMemcpyHostToDevice); cudaMemcpy(v_dev, v, (ntri_new_m[0])*sizeof(float), cudaMemcpyHostToDevice); где первым параметром идет указатель, содержащий адрес места-назначения копирования; вторым - указатель, содержащий адрес источника копирования; следом - размер копируемого ресурса в байтах; далее - параметр, указывающий направление копирования.

3. Распределить ядра для вычисления.

Расчет газодинамических параметров состоит из последовательности вызываемых процедур, запускаемых на GPU на каждом шаге по времени:

while (t <= tmax) { boundcond<<<blocks, threads>>> (rdev, pdev, udev, vdev, dldev, d2_dev, d3_dev, d4_dev, d5_dev, d6_dev, otrudev, otruldev, otr_v_dev, otr_v1_dev, alphaldev, alpha2_dev, sootvdev, cm); calcU<<<blocks, threads>>>(r_dev, p dev, u dev, v dev,

rodev, rudev, rvdev, endev, ntridev); ENO_mass<<<blocks, threads>>>(ux_m_dev, uymdev, vxmdev, vymdev,

rxmdev, rymdev, pxmdev, pymdev, ntri dev, trineighldev, cmdev, r dev, u dev, v dev, p dev, ntrinewdev); calcFluxes<<<blocks, threads>>>(r_dev, p dev, u dev, v dev, rtmp_dev, utmp_dev, vtmp_dev, etmp_dev, uxmdev, uy m dev, vx m dev, vy m dev, rx m dev, ry m dev, px m dev, py m dev, xldev, tri neighl dev, tri_v1_dev, cm dev, ntri dev, ntri new dev, nxnewdev); gasDinStep<<<blocks, threads>>>(rn_dev, pndev, undev, vndev, ro dev, ru dev, rv dev, en dev, rtmp_dev, utmp_dev, vtmp dev, etmp dev, s dev, ntri dev); t = t+tau; step++;

}

Процедура boundcond<<<blocks, threads>>>() определяет граничные условия.

4

Процедура calcU<<< blocks, threads>>>() вычисляет консервативные переменные в соответствие с формулой (3).

В процедуре ENO_mass<<<blocks, threads>>>() реализован алгоритм реконструкции газодинамических параметров на границе ячеек.

В процедуре calcFluxes<<<blocks, threads>>>() рассчитываем дискретные потоки по схеме Лакса-Фридрихса в соответствии с формулой (8).

Процедура gasDinStep<<<blocks, threads>>>() рассчитывает значения консервативных, а затем примитивных переменных на каждом шаге по времени.

4. Скопировать полученные данные обратно в CPU. Пример: газодинамические параметры (r,p,u,v) копируем на CPU

cudaMemcpy(r, rndev, (ntri_new_m[0])*sizeof(float), cudaMemcpyDeviceToHost); cudaMemcpy(p, pndev, (ntri_new_m[0])*sizeof(float), cudaMemcpyDeviceToHost); cudaMemcpy(u, undev, (ntri_new_m[0])*sizeof(float), cudaMemcpyDeviceToHost); cudaMemcpy(v, vndev, (ntri_new_m[0])*sizeof(float), cudaMemcpyDeviceToHost);

5. Освободить память GPU. cudaFree(r_dev); cudaFree(pdev); cudaFree(udev); cudaFree(vdev);

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

3. Вычислительный эксперимент. Большой интерес с практической точки зрения представляет задача взаимодействия ударной волны с зоной турбулентного перемешивания. Такие виды задач являются хорошими тестовыми примерами для проверки различных алгоритмов, методов и программ. На практике такие задачи позволяют изучить ход образования вихрей при контактном взаимодействии нескольких газов с разными плотностями и уравнениями состояний [6].

Выполнено численное исследование динамики развития газодинамической неустойчивости в контактной границе элегаз-воздух для условий эксперимента на ударной трубе [7]. Параметры ударной трубы следующие: длина равна 0,4 м, а ширина равна 0,04 м. Труба разделена на три области:

1. область длиной 0,08 м, заполнена элегазом с высоким давлением (SF6_HI);

2. область длиной 0,02 м, заполнена элегазом с низким давлением (SF6_LO);

3. область длиной 0,03 м, заполнена воздухом (Л/ ft).

В качестве расчетных данных была принята следующая физическая постановка задачи: начальная температура 291 К в камере низкого давления давление составляет 10-4

ГПа = 1 бар. Непосредственно за ударной волной в элегазе давление составляет 2,152 бар, плотность - 1,209 * 10-2 г/см3, скорость ударной волны - 195,2 м/с, скорость течения за ударной волной - 97,76 м/с. Начальные плотности элегаза и воздуха в камере низкого давления равны 6,037 * 10 и 1,198 * 10-3 г/смз соответственно. Физические свойства элегаза и воздуха следующие: оба вещества являются невязкими, нетеплопроводными и идеальными газами с показателями адиабаты у = 1,094 (5Р6) и у = 1.4 (воздух), отношение молекулярных масс (5Р6 /воздух) принято равным 5,04.

Необходимо смоделировать развитие неустойчивости при многократном прохождении ударной волны через контактный разрыв.

На рисунке 1 приведен график плотности в начальный момент времени. На рисунках 2 и 3 представлены изменения плотности при первом и повторном прохождении ударной волны через контактную границу.

Рис. 1. Изменение плотности в начальный момент времени.

Рис. 2. Изменение плотности при первом прохождении ударной волны через контактную границу.

Рис. 3. Изменение плотности при втором прохождении ударной волны через контактную границу.

Заключение. Разностная схема и программа расчета были протестированы на задаче с известными экспериментальными данными [7]. Результаты численных расчетов показали соответствие теоретическим предложениям и экспериментальным данным. Таким образом, проведенное тестирование предложенной схемы и программы моделирования

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

ЛИТЕРАТУРА

1. NVIDIA CUDA. Programming Guide [Электронный ресурс]. - Режим доступа: http://developer.nvidia.com/cuda-toolkit-40.

2. Kurganov A., Chi-Tien L. On the reduction of numerical dissipation in central-upwind schemes // Commun. Comput. Phys. February. - 2007. - Vol. 2, No. 1. - pp. 141-163.

3. Жалнин Р. В. О построении параллельного вычислительного алгоритма высокого порядка точности для гиперболических систем уравнений // Журнал Средневолжского математического общества. - 2007. - Т. 9, № 1. - С. 145-153.

4. Жалнин Р. В., Змитренко Н. В., Ладонкина М. Е., Тишкин В. Ф. Численное моделирование развития неустойчивости Рихтмайера-Мешкова с использованием схем высокого порядка точности // Математическое моделирование. - 2007. -Т. 19, № 10. - С. 61-66.

5. Ладонкина М. Е. Численное моделирование турбулентного перемешивания с использованием высокопроизводительных систем: дисс. ... канд. физ.-мат. наук. -М., 2005. - 157 с.

6. Lebo I. G., Nikishin V. V., Rozanov V. B., Tishkin V. F. On the development of instability near the contact boundary between two equal-density gases in passage of a shock wave // Bulletin of the Lebedev Physics Institute (Kratkie Soobsheniya po Fizike). - New York: Allerton Press, Inc, 1997.

7. Poggi F., Thorembey M.-H., Rodriguez G. Velocity measurements in turbulent gaseous mixtures induced by Richtmyer-Meshkov instability // Physics of Fluids. - 1998. -Vol. 10, No. 11. - P. 2698-2700.