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

Barochkin Evgeny Vitalyevich,

Ivanovo State Power Engineering University,

Doctor of Engineering Sciences (Post-Doctoral degree), Professor, Head of the Heat Power Plants Department, address: Ivanovo, Rabfakovskaya St., 34, Building «V» (В), Room 417, telephone (4932) 26-99-13, e-mail: admin@tes.ispu.ru

Ледуховский Григорий Васильевич,

ФГБОУВО «Ивановский государственный энергетический университет имени В.И. Ленина»,

кандидат технических наук, доцент кафедры тепловых электрических станций,

адрес: г. Иваново, ул. Рабфаковская, д. 34, кор. В, ауд. 408,

телефон (4932) 26-99-31,

e-mail: lgv83@yandex.ru

Ledukhovsky Grigory Vasilievich,

Ivanovo State Power Engineering University,

Candidate of Engineering Sciences (PhD), Associate Professor of the Heat Power Plants Department, address: Ivanovo, Rabfakovskays St, 34, Building В (V), Room 408, telephone (4932) 26-99-31, e-mail: lgv83@yandex.ru

УДК 004.942

Методы ускорения расчетов и повышения устойчивости полевой модели гидродинамических турбулентных процессов на CUDA

А.А. Гудухина, И.Ф. Ясинский ФГБОУВО «Ивановский государственный энергетический университет имени В.И.Ленина»,

г. Иваново, Российская Федерация E-mail: igor.yasinskiy@gmail.com

Авторское резюме

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

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

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

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

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

Calculations acceleration and stability improvement methods for hydrodynamic turbulent processes model on CUDA

A.A. Gudukhina, I.F. Yasinskiy Ivanovo State Power Engineering University, Ivanovo, Russian Federation E-mail: igor.yasinskiy@gmail.com

Abstract

Background. High computational complexity of computer calculations remains one of the main problems in the construction of mathematical models of turbulent hydro and gas dynamic environments, since it is necessary to obtain a solution to a non-stationary problem with a grid of such a spatial step, which corresponds to the size of the smallest vortex structures. The authors took into account the recommendations on the use of differential computational schemes for the turbulent dynamics of a liquid, obtained in the work on modeling the bottling of oil products on the water surface. The purpose of this work is to obtain the mathematical model stable calculation of the viscous incompressible fluid dynamics according to the Euler description, taking into account the influence of turbulence, and the acceleration of calculations in a parallel CUDA interface.

Materials and methods. The method of mathematical modeling of physical objects and the method of differential equations numerical integration are used. The method of countercurrent derivatives is used to increase the stability of differential schemes. The Sekundov's method of turbulent viscosity estimation of a continuous medium is used, as well as the analysis of the computational experiments results.

Results. A stable parallel realization has been created for numerical simulation of processes in a continuous medium. It takes into account the influence of turbulizing hydro- and gas-dynamic structures. Stability is achieved by using the counter-current derivative method and replacing traditional difference schemes, when calculating the velocity and pressure fields, with four-point counterparts. The method of uniform distribution of the graphics accelerator computing resources with large sizes of the computational grid is proposed. A stable calculation created of a mathematical model of continuous medium dynamics with the Sekundov's turbulent viscosity estimation at a considerable time interval. Conclusion. Based on the NVidia CUDA interface, this solution allows acceleration of calculations to be achieved from 2 to 8 times, depending on the capabilities of the hardware configuration. The developed system can be a tool for the study of objects in the energy industry, for which acceleration of control decisions is required, compared with standard software, that does not provide for calculations parallelization. Reliability of the results is caused by the compliance of the system conditions with the configuration of two-phase supply systems of power installations.

Key words: continuous media, hydro-gas dynamics, turbulence, parallel computing, CUDA parallel interface

DOI: 10.17588/2072-2672.2018.6.072-080

Введение. Для современной науки представляют интерес природные и технологические процессы механики жидкости и газа, осложненные вихревыми турбулентными течениями. Турбулизация приводит к масштабному воздействию на все характеристики поля течения. Особые свойства таких течений находят применение в энергетическом, теплообмен-ном, технологическом оборудовании ядерной энергетики, аэрокосмической техники, нефтеперерабатывающей промышленности, промышленной теплоэнергетики [1]. Большие перспективы имеет использование вихревого закрученного потока для комплексного решения задач повышения безопасности и эффективности работы ядерно-энергетических установок, конструктивного усовершенствования МГД-генераторов, ускорения протекания химико-технологических процессов, таких как обжиг, плавление, коксование, пиролиз, интенсификации процессов тепло- и массообмена в условиях вынужденной и естественной конвекции [2]. Особый интерес вызывает моделирование осложненной динамики жидкости в контексте проектирования аппаратов для производства высокодисперсных сред [3].

Моделирование вихревой турбулентной динамики является проблемой, которая до сих

пор решена не полностью [4]. При математическом моделировании механики сплошной среды используются следующие способы описания турбулентных режимов:

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

• модели, описывающие турбулентную вязкость в зависимости от поля скоростей, расстояния от твердых поверхностей (модель Секундова);

• параметрические модели турбулентности, в которых выполняется расчет кинетической энергии пульсаций и энергии диссипации (К-Е модель).

Рейнольдсом установлено, что переход от ламинарного течения к турбулентному происходит при достижении безразмерным комплексом Re, названным впоследствии числом Рейнольдса, некоторого предельного значения Reкp. Число Re включает в себя среднерасход-ную скорость среды и, характерный размер d и кинематический коэффициент вязкости потока V [5]. Сначала с увеличением значения числа Re, отражающего склонность сплошной среды к возникновению турбулентных пульсаций в данных условиях, наблюдается смена режима

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

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

Эволюция идей гидродинамики на суперкомпьютерах имеет долгую историю. В качестве заметных итераций можно отметить гибридную параллельную реализацию сеточных уравнений Больцмана в [6]. В скором времени после появления системы CUDA, дающей возможность реализации задач общего назначения на мультипроцессорных видеоускорителях Nvidia, эта система приобрела популярность в области моделирования задач механики сплошной среды. В 2010 году вышла работа Д.Е. Демидова и соавторов [7], в которой обозначены подходы к получению ускорения различных моделей гидродинамических процессов на интерфейсе CUDA.

Известны результаты, изложенные в публикации D. Nishiura и соавторов, 2015 [8], где рассматривается распараллеливание несеточных моделей гидродинамики в парадигме Лагранжа на гетерогенных суперкомпьютерах с общей памятью, включая системы OpenMP и GPGPU. Решение подобных задач уже на массивно-параллельном суперкомпьютере Archer осветили в 2018 году X. Guo и соавторов в [9].

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

Математическое описание задачи о вихревом турбулентном движении вязкой несжимаемой жидкости. Как было отмечено, с ростом числа Рейнольдса размеры вихревых образований в потоке уменьшаются. В этом случае расчетная сетка должна соответствовать размерам самых малых турбулентных структур, что значительно повысит затраты машинного времени на расчет модели. Для оценки возможности ускорения работы была взята область с такими начальными и граничными условиями, которые способствуют образованию турбулентных структур (рис. 1). В ней имеется два боковых впускных отверстия Uвх1, Цвх2 и одно донное выходное ивых. Геометрия источников и стоков хорошо соответствует конфигурации двухфазных систем питания энергетических установок.

Рис. 1. Конфигурация области моделирования динамики сплошной среды

Для описания как ламинарных, так и турбулентных течений применяется система уравнений Рейнольдса [4]:

^ + Zи.^ = УS_ 1 dP.; i 1,2,3; dt t-i 1 дх, t-i dx, p dXj

J=1

J=1

yU = 0;

¿—i Яу

i=1 dXi

Sj = veff

idUi dUjл

dXj dXj

(1)

(2)

(3)

где иI — составляющие вектора скорости по осям х, у, г; — эффективная кинематическая вязкость.

Кинематическая вязкость вычисляется суммой молекулярной вязкости чт0! и турбулентной вязкости УхшЬ. В случае выраженного ламинарного течения с учетом постоянства турбулентной вязкости уравнения (1)—(3) упрощаются до формы Навье-Стокса:

dU,

dU,

d2U 1 dP

— + TU: duL = Vf _1 ; i = 1,2,3. (4)

dt y 1 dXj eff y dx2 -

P dXj

Для расчета турбулентной вязкости в данном случае была взята модель А.Н. Секун-дова [10]:

dv

turb

dt

j=1

dv

turb

d

дХ:

= l— (

1 dXj I

vmol

+ KV

dv

turb)

turb

дХ

(

+vturbf

turb

V 8vmol J

D yLmin (Vmol + PVturb ) V

turb ■

где f() - функция:

f (z ) = 0,2

z2 +1,47 • z + 0,2.

(5)

(6)

г2 -1,47 • г +1 к, р, у - константы, подобранные эмпирическим путем (к = 2; р = 0,06; у = 50); й - величина, описывающая силу деформации среды. Ее получение подразумевает использование формы

(

D =

уЩ

V1 dxj

dU, dU,

w

1/2

dx,

Vj

dXj

(7)

(8)

По модели Кармана длина пути перемешивания £ вычисляется согласно выражению

_К ду/ ду2

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

|г= 0. Для формирования поля турбулентности необходимо установить начальные значения. Данные могут быть получены эмпирическим путем или с использованием модели Прандтля-Кармана.

Модель А.Н. Секундова (5)-(6) показывает достаточно точные результаты, удобна для компьютерного моделирования и применима для описания турбулентных течений в несжимаемых средах. Эту модель удобно сочетать с уравнениями Навье-Стокса, проинтегрировать дифференциальное уравнения (5) можно аналогичным образом (явная схема, метод прогонки и др.).

Распараллеливание задачи. Распараллеливание по точкам системы хорошо ложится на сетку процессов параллельного интерфейса CUDA [11]. Каждый процесс может просчитывать одну точку или группу точек, поставленную ему в соответствие (рис. 2).

Алгоритм работы модели, реализованный в виде главной функции тат(), представлен на рис. 3.

I

ь

"О

tic.: Иск (0,0)

(С. О) (1.0) (2,01

(0,1) (1,4 12,1)

л (0,2) 11,2) (2,2)

blockDim j: Ыа:Шх(0,1)

(О.О) (1.0) (2,0)

(0,1) (1.1) (2.'l)

(0,2) (1,2| «.■2)

blockDImjc/ blockldx(0,2)

(0,0) O.cf) (2.01

(0,1) (1,4 Vt.il)

fa. 2) (1.2) ф)

<0. ф |Т(0| '(2,0)" i ■s о л (0,0) (1.0) (2.0)

Rib V-" (o.i) (1,1) (2,1)

(0.2) (1,2) '(0,2) (1,2) (2,2)

blockDirrji blocklcU (1,1) blockDim jt blockldx (2,1)

(O.U) (i,m ЦО) E S ь (0,0) P,0) (2,0)

(0.1) n.n a-D (0,11 .[1.11 (2.1)

(9,2) (1,2) (2,2) (0,2) (1,2) Pjl

blockDim.K blpckDifiijc

bloctldid.a. _ _ LI ql kid* (2,2)"

(0,0) (1.0) (2,0)

(0,1) 11 Гч (2,1)

(0.2) (1,2) (2,2)

(0,0) ОГО) (2,0)

(0,1) (1,1) (2,1)

(0,2) (1.2) (2.2)

•«-gridDím.x-

Рис. 2. Наложение сетки СиОА на сетку узлов полевой модели

Начало

J

Объявление переменных., выделение глобальной памяти

Заполнение матрицы начальные значений Функцией filljn

Вычисление гола давлений и

поля скоростей & новый момент времени функцией update

Копирование новы* значений в массивы текущих значении

Выполнение шага по времени t f= tau

нет

Вы&Од результатов, освобождение памяти

С

Конец

Рис. 3. Алгоритм главной функции

Блок 1. Начало алгоритма.

Блок 2. Сначала создаются одномерные массивы в памяти CPU для удобства копирования значений в глобальную память. Массивы необходимы для хранения информации об актуальном состоянии моделируемой сплошной среды, а также для хранения значений предыдущих шагов для вычисления новых. Поле давлений сохраняется при помощи единственного массива:

Ux = newfloat[Nx * Ny];

Ux_old = newfloat[Nx * Ny];

Uy = newfloat[Nx * Ny];

Uy_old = newfloat[Nx * Ny];

P = newfloat[Nx * Ny];

Соответствующие им массивы создаются в глобальной памяти видеокарты:

cudaMalloc((void**)&UxDev, size);

cudaMalloc((void**)&Ux_oldDev, size);

cudaMalloc((void**)&UyDev, size);

cudaMalloc((void**)&Uy_oldDev, size);

cudaMalloc((void**)&P_Dev, size);

Блок 3. Массивы рассчитываемых параметров заполняются начальными значениями с помощью функции Ш_тО, которая выполняется на видеокарте и использует параллелизм:

fill_in<<<dim3((Nx / block_size) + 1, (Ny / block_size) + 1), dim3(block_size, block_size) >>> (UxDev, UyDev, P_Dev);

Блоки 4-6. В этом блоке начинается цикл по времени, в рамках которого вызывается функция update(), рассчитывающая обновленные значения полей скорости и давления сплошной среды, которые затем копируются в массив текущих значений и выполняется переход к следующему шагу по времени:

do{

update <<< dim3((Nx / block_size) + 1, (Ny / block_size) + 1), dim3(block_size, block_size) >>> (UxDev, Ux_oldDev, UyDev, Uy_oldDev, P_Dev);

cudaThreadSynchronize(); cudaMemcpy(Ux_oldDev, UxDev, size, cudaMemcpyDeviceToDevice);

cudaMemcpy(Uy_oldDev, UyDev, size, cudaMemcpyDeviceToDevice); t += tau;

} while (t < tmax);

Блоки 7, 8. Цикл заканчивает работу по достижении максимального времени моделирования. После этого выводятся результаты расчетов и освобождается память, выделенная под массивы.

Параллелизм данной реализации заключается в функциях Ш_тО и update().

В общем случае нить с координатами i, j вычисляет параметры точки целевой системы с координатами i и j или если точек в системе больше, чем доступных нитей, то любую точку, которая будет поставлена в соответствие данной нити. В CUDA координаты нити находятся следующим образом [11]:

inti = blockIdx.x*blockDim.x + threadldx.x;

int j = blockIdx.y*blockDim.y + threadldx.y;

где blockldx.x и blockldx.y - координаты блока нитей по осям х и у соответственно; blockDim.x и blockDim.y - размерности блока нитей по осям х и у соответственно; threadldx.x и threadldx.y - координаты нитей в блоке по осям х и у соответственно.

Для равномерного распределения большого количества узлов расчетной сетки среди меньшего количества нитей предлагается определять номер нити для расчета данного узла через операцию остатка от деления номера данного узла на общее число доступных нитей: N = Nc modNf, (9)

где N - номер нити, которая будет обсчитывать данный узел; Nt - общее число доступных нитей; Nc - номер узла, рассчитываемый по формуле

Nc = i ■ Nx + j , (10)

где Nx - размер сетки по оси Х; i, j - координаты ячейки в сетке.

Параллелизм инициализации начальными значениями (функция fill_in()). Данная функция заполняет массивы параметров начальными значениями. Каждой нити ставится в соответствие узел, значения параметров которого она рассчитывает. Если узел попадает на входные или выходные отверстия системы, то начальные значения вычисляются по граничным условиям. В противном случае значения заполняются в соответствии с принятыми начальными услоями (в рассматриваемом случае - нулевыми). Алгоритм выполнения функции представлен на рис. 4.

Начальные скорости и давления в рассматриваемом случае приняты нулевыми:

Ux[j*Nx + i] = 0.0;

Uy[j*Nx + i] = 0.0;

P[j*Nx + i] = p0;

Далее каждая нить вычисляет координаты нижних входных и выходных отверстий. Если эти координаты поставлены в соответствие нити, то они вычисляют скорости в данной точке:

if ((i>= i1) && (i<i2) && (j == Ny - 1)) Uy[j * Nx + i] = F(x);

ихи * Nx + i] = F(x);

Здесь F(x) - функция, по которой вычисляется скорость движения жидкости.

Рис. 4. Блок-схема функции fill-in

Параллелизм расчета новых значений параметров (функция update()). Алгоритм функции представлен на рис. 5. Пересчитыва-ются скорости и давление каждого отдельного узла. Все нити, кроме тех, что находятся на границах, пересчитывают давление по методу слабой сжимаемости [10]:

2 г, дР -с ■ divU = ■

dt

(11)

где с - константа, влияющая на характеристики жидкости в системе.

Эта форма вычисляется следующим образом:

if ((i<Nx - 1) && (j<Ny - 1) && (i> 0) && (j> 0)) {

divU = f1(i, j, Ux, Uy);

P[j * Nx + i] -= c * tau * divU;

}

В целях повышения устойчивости дивергенция скорости D берется по схеме

II k _!_// k _// k _ц k -,k+1 _ +1,j +1 + Uxi+1,j-1 Uxi-1,j+1 Uxi-1,j-1

4h

■ +

Uk

+ UJ

k

Uk

k

_ üj

k

k

yi+1 ,j +1 T Uyi_1 ,j +1 Uyi+1 ,j_1 Uyi_1 ,j_1

(12)

4h

Рис. 5. Алгоритм функции update()

Ее реализует функция f1(), которая выглядит следующим образом:

_device_float f1(int, int, float* J>, float* L )

{

return ((Ux[(j + 1) * Nx + i + 1] + Ux[(j - 1) * Nx + i + 1]) -(Ux[(j - 1) * Nx + i - 1] + Ux[(j + 1) * Nx + i - 1]) +(Uy[(j + 1) * Nx + i - 1] + Uy[(j + 1) * Nx + i + 1]) - (Uy[(j - 1) * Nx + i - 1] + Uy[(j - 1) *

Nx + i + 1])) / 4.0 / h; }

В функцию передаются координаты i и j точки, так что при вызове функции на каждой нити просчитывается отдельно взятый узел системы. Функция будет вызываться в ядре, а выполняться на видеокарте, так как перед ней

стоит спецификатор_device_.

«Граничные» нити находят значения давления копированием значений из соседних клеток.

Для поля давления на всех границах использованы граничные условия II рода, т. е. задается производная давления (плотности) по нормали к границе. В целях упрощения производная давления взята равной нулю на всех

+

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

Проекции вектора скорости среды рассчитываются согласно выражениям:

in k

UkT1 = Uk,

'Xi ,j

'Xij

Uk + \Uki\u k-u k„.

xi.l XiJ\ UXi,j UXi-1, j

uxk,-Uk/ u*-uk.

Xi , 1 I Xi, j| UXi +1 , j UXi, j

2 h

u k + lu k.l + ryi,i\ u k -u k uXi,I uXi, j-1

2 h

u k. -lu k.l \uyi I u k -u k Xi, j+1 Xi, I

+v

u k + u k + u k + u k -4U k

UXi+1,| + UXi-1,| + u Xi, l +1 + u Xi, l-1 4UXi,l h2

k+1 k+1 k+1 k+1 1 Pi+1,|+1 + Pi+1,|-1 Pi-1,|+1 Pi-1,|-1 c k+1 ---+ FXi, j

\

P

uk+1 = uk+t

4h

yi,j

yi,j

uk + \uxk\u k-u *

X', j Xl, j\uyi , j uyi-1, j

uXk/-|uXk/lu k„ -u k

Xi, j I Xi, I uyi +1 , j uyi, j

2 h

u k +u k.l +ryi,i\ [J k -u k uyi,l uyi,l-1

2 h

u k. -lu k.l ryi,j\ k -u k yi, j +1 yi, j

+v-

uyk+M + uyk- и + uyk I+1 + uykj-1 - 4u

yi, j

h

2

. pk+1 pk+1 _ pk+1 _ pk+1

1 Pi+1, |+1 + Pi-1, |+1 Pi+1, 1 Pi-1, j-1 p k+1

'P 4h Fyi,j

\

(13)

(14)

Сначала для всех точек, кроме граничных, вычисляется новое значение скорости по оси X с помощью функции /2():

И ((I <Ых - 1) && 0 <Ыу - 1) && (I > 0) && 0 > 0)) {

I и * Ых + I] = 12(1, 11х_о!с1, 11у_о!с1, Р);

}

Сама функция /2() при этом выполняется на девайсе и имеет вид

Ux[j * Nx + i - 1]) / h + (Ux[j * Nx + i] -abs(Ux[j * Nx + i])) / 2.0 * (Ux[j * Nx + i + 1] -

Ux[j * Nx + i]) / h) - ((Uy[j * Nx + i] + abs(Uy[j * Nx + i])) / 2.0 * (Ux[j * Nx + i] -

Ux[(j - 1) * Nx + i]) / h + (Uy[j * Nx + i] -abs(Uy[j * Nx + i])) / 2.0 * (Ux[(j + 1) * Nx + i] -

Ux[j * Nx + i]) / h) - (P[(j + 1) * Nx + i + 1] + P[(j - 1) * Nx + i + 1] - P[(j + 1) * Nx + i - 1] -

P[(j - 1) * Nx + i - 1]) /4.0 / h / ro + nu * (Ux[j * Nx + i + 1] + Ux[j * Nx + i - 1] +

Ux[(j + 1) * Nx + i] + Ux[(j - 1) * Nx + i] -4 *

Ux[j * Nx + i]) / h / h); }

Аналогично с помощью функции f3() рассчитывается поле скоростей по оси Y, функция вызывается в ядре, а выполняется на видеокарте. В соответствующей функции программы также выполняется обновление поля турбулентной вязкости согласно модели А.Н. Секундова (5).

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

Вычислительные эксперименты по ускорению параллельной реализации математической модели турбулентной динамики сплошной среды. Вычислительные эксперименты проводились на расчетной сетке размерностей 100x100, 200x200, ... , 1000x1000.

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

В таблице представлены результаты по затратам машинного времени для последовательного и параллельного алгоритмов с использованием технологии CUDA.

Последовательная программа запускалась на Intel Core i7 частотой 2,2 Ггц. Параллельные расчеты CUDA проводились на NVIDIA GEFORCE 630M и NVIDIA GTX 1050.

Обсуждение. Сравнительный анализ затрат машинного времени при последовательной и параллельной реализациях программы моделирования механики жидкости иллюстрирует рис. 7.

_device_float f2(int, int, float* Jx, float*

Uy, float* P) {

returnJx[j * Nx + i] + tau * (-((Jx[j * Nx + i] +abs(Jx[j * Nx + i])) / 2.0 * (Jx[j * Nx + i] -

2

h

V

2

h

2

h

\

2

h

k

Рис. 6. Поле скоростей при расчете математической модели движения вязкой несжимаемой жидкости

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

Размерность системы 100 200 300 400 500 600 700 800 900 1000

Последовательная реализация 1,10 4,05 7,50 13,58 21,67 30,58 41,87 54,1 68,22 85,94

Nvidia GeForce GT630M 0,64 1,77 3,85 6,53 10,48 14,46 20,03 25,8 32,96 40,15

Nvidia GeForce GTX 1050 0,31 0,69 1,34 2,25 3,46 4,89 6,52 8,51 11,04 13,44

Рис. 7. Графики затрат машинного времени при последовательной и параллельной реализациях программы моделирования механики жидкости

Параллельный алгоритм с использованием технологии CUDA ускорил вычисления в 8 раз на современной видеокарте NVIDIA GTX 1050 с computer capability 6.1. Менее производительная видеокарта Nvidia Geforce 630M с computer capability 2.1 дала ускорение вычислений в среднем в 2 раза.

Таким образом, построено устойчивое высокопроизводительное решение вычислительно сложной задачи о моделировании вязкой несжимаемой жидкости с включением турбулентных образований на базе интерфейса CUDA графического ускорителя NVidia. Иссле-

дование показало возможности значительного ускорения при расчете описанной математической модели (от 2 до 8 раз) по сравнению с последовательной реализацией в зависимости от производительности аппаратной конфигурации.

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

Список литературы

1. Митрофанова О.В. Гидродинамика и теплообмен закрученных потоков в каналах ядерно-энергетических установок. - М.: ФИЗМАТЛИТ, 2010. -288 с.

2. Лойцянский Л.Г. Механика жидкости и газа: учеб. для вузов. - 7-е изд., испр. - М.: Дрофа, 2003. - 840 с.

3. Кожевников С.О., Калинин Е.Н. К постановке задачи моделирования движения жидкой среды в перемешивающем аппарате // Вестник ЧГУ. -2017. - № 1. - С. 40-47.

4. Моделирование процессов розлива нефтепродуктов по водной поверхности с использованием суперкомпьютеров с графическими ускорителями / Ф.Н. Ясинский, С.Г. Сидоров, И.А. Малый и др. // Вестник ИГЭУ. - 2013. - Вып. 6. - С. 90-94.

5. Ковальногов Н.Н. Прикладная механика жидкости и газа. - Ульяновск: УлГТУ, 2010. - 219 с.

6. Heuveline V., Krause M.J., Latt J. Towards a hybrid parallelization of lattice Boltzmann methods // Computers and Mathematics with Applications. - 2009. -№ 58. - Р. 1071-1080.

7. Демидов Г.Е., Егоров А.Г., Нуриев А.Н. Решение задач вычислительной гидродинамики с применением технологии Nvidia CUDA // Ученые записки КГУ. - 2010. - Т. 152 (1). - С. 142-154.

8. Nishiura D., Furuichi M., Sakaguchi H. Computational performance of a smoothed particle hydrodynamics simulation for shared-memory parallel computing // Computer Physics Communications. -2015. - № 194. - Р. 18-32.

9. New massively parallel scheme for Incompressible Smoothed Particle Hydrodynamics (ISPH) for highly nonlinear and distorted flow / X. Guo, B.D. Rogers, S. Lind, P.K. Stansby // Computer Physics Communications. - 2018. - № 233. - Р. 16-28.

10. Ясинский Ф.Н., Кокорин А.С. Математическое моделирование процессов вентиляции и отопления в больших производственных, культурных и спортивных помещениях // Вестник ИГЭУ. - 2010. -Вып. 3. - С. 90-92.

11. Боресков А.В., Харламов А.А. Основы работы с технологией CUDA. - М.: ДМК Пресс, 2010. - 232 с.

References

1. Mitrofanova, O.V. Gidrodinamika i teploobmen zakruchennykh potokov v kanalakh yaderno-energeticheskikh ustanovok [Hydrodynamics and heat

transfer of swirling flows in channels of nuclear power plants]. Moscow: FIZMATLIT, 2010. 288 p.

2. Loytsyanskiy, L.G. Mekhanika zhidkosti i gaza [Fluid and gas mechanics]. Moscow: Drofa, 2003. 840 p.

3. Kozhevnikov, S.O., Kalinin, E.N. K postanovke zadachi modelirovaniya dvizheniya zhidkoy sredy v peremeshivayushchem apparate [On the formulation of the problem of modeling the motion of a liquid medium in a mixing apparatus]. Vestnik ChGU, 2017, no. 1, pp. 40-47.

4. Yasinsky, F.N., Sidorov, S.G., Malyi, I.A., Po-temkina, O.V., Yasinsky, I.F., Mochalov, A.S., Net-kachev, V.V. Modelirovanie protsessov rozliva neftepro-duktov po vodnoy poverkhnosti s ispol'zovaniem superk-omp'yuterov s graficheskimi uskoritelyami [Modeling of the processes of bottling oil products on the water surface using supercomputers with graphic accelerators]. Vestnik IGEU, 2013, issue 6, pp. 90-94.

5. Koval'nogov, N.N. Prikladnaya mekhanika zhidkosti i gaza [Applied fluid and gas mechanics]. Ulyanovsk: UlGTU, 2010. 219 p.

6. Heuveline, V., Krause, M.J., Latt, J. Towards a hybrid parallelization of lattice Boltzmann methods. Computers and Mathematics with Applications, 2009, no. 58, pp. 1071-1080.

7. Demidov, G.E., Egorov, A.G., Nuriev, A.N. Reshenie zadach vychislitel'noy gidrodinamiki s primeneniem tekhnologii Nvidia CUDA [Solution of the computational hydrodynamics problems by means of Nvidia CUDA technology]. Uchenye zapiski KGU, 2010, vol. 152(1), pp. 142-154.

8. Nishiura, D., Furuichi, M., Sakaguchi, H. Computational performance of a smoothed particle hydrodynamics simulation for shared-memory parallel computing. Computer Physics Communications, 2015, no. 194, pp. 18-32.

9. Guo, X., Rogers, B.D., Lind, S., Stansby, P.K. New massively parallel scheme for Incompressible Smoothed Particle Hydrodynamics (ISPH) for highly nonlinear and distorted flow. Computer Physics Communications, 2018, no. 233, pp. 16-28.

10. Yasinskiy, F.N., Kokorin, A.S. Matemati-cheskoe modelirovanie protsessov ventilyatsii i otopleniya v bol'shikh proizvodstvennykh, kul'turnykh i sportivnykh pomeshcheniyakh [Mathematical modeling of the processes of ventilation and heating in large industrial, cultural and sports facilities]. Vestnik IGEU, 2010, issue 3, pp. 90-92.

11. Boreskov, A.V., Kharlamov, A.A. Osnovy raboty s tekhnoloriey CUDA [Basics of working with CUDA technology]. Moscow: DMK Press, 2010. 232 p.

Ясинский Игорь Федорович,

ФГБОУВО «Ивановский государственный энергетический университет имени В.И. Ленина»,

кандидат технических наук, доцент кафедры высокопроизводительных вычислительных систем,

e-mail: igor2266@yandex.ru.com

Yasinskiy Igor Fedorovich,

Ivanovo State Power Engineering University

Candidate of Engineering Sciences (PhD), Associate Professor of the Department of High-Performance Computer Systems, e-mail: igor2266@yandex.ru

Гудухина Александра Антоновна,

ФГБОУВО «Ивановский государственный энергетический университет имени В.И. Ленина»,

магистрант кафедры высокопроизводительных вычислительных систем,

e-mail: esqulapic@mail.ru

Gudukhina Alexandra Antonovna,

Ivanovo State Power Engineering University,

A Master course student of the Department of High-Performance Computer Systems, e-mail: esqulapic@mail.ru