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

УДК 004.272.22                                                                  Научная статья

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

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

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

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

Финансирование. Работа выполнена при поддержке Российского научного фонда (проект № 21–71–20050).

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

Original article

Model of a Parallel-Pipeline Computational Process for Solving a System of Grid Equations

Vladimir N. Litvinov¹ B , Nelli B. Rudenko² © , Natalya N. Gracheva² ©

• ¹    Don State Technical University, Rostov-on-Don, Russian Federation

• ²    Azov-Black Sea Engineering Institute, Don State Agrarian University, Zernograd, Russian Federation

Introduction. Environmental problems arising in shallow waters and caused by both natural and man-made factors annually do significant damage to aquatic systems and coastal territories. It is possible to identify these problems in a timely manner, as well as ways to eliminate them, using modern computing systems. But earlier studies have shown that the resources of computing systems using only a central processor are not enough to solve large scientific problems, in particular, to predict major environmental accidents, assess the damage caused by them, and determine the possibilities of their elimination. For these purposes, it is proposed to use models of the computing system and decomposition of the computational domain to develop an algorithm for parallel-pipeline calculations. The research objective was to create a model of a parallel-conveyor computational process for solving a system of grid equations by a modified alternating-triangular iterative method using the decomposition of a three-dimensional uniform computational grid that takes into account technical characteristics of the equipment used for calculations.

Materials and Methods. Mathematical models of the computer system and computational grid were developed. The decomposition model of the computational domain was made taking into account the characteristics of a heterogeneous system. A parallel-pipeline method for solving a system of grid equations by a modified alternating-triangular iterative method was proposed.

Results. A program was written in the CUDA C language that implemented a parallel-pipeline method for solving a system of grid equations by a modified alternating-triangular iterative method. The experiments performed showed that with an increase in the number of threads, the computation time decreased, and when decomposing the computational grid, it was rational to split into fragments along coordinate z by a value not exceeding 10. The results of the experiments proved the efficiency of the developed parallel-pipeline method.

Discussion and Conclusion. As a result of the research, a model of a parallel-pipeline computing process was developed using the example of one of the most time-consuming stages of solving a system of grid equations by a modified alternating-triangular iterative method. Its construction was based on decomposition models of a three-dimensional uniform computational grid, which took into account the technical characteristics of the equipment used in the calculations. This program can provide you for the acceleration of the calculation process and even loading of program flows in time. The conducted numerical experiments validated the mathematical model of decomposition of the computational domain.

Acknowledgements: the authors would like to thank the editorial board of the journal and the reviewer for their professional analysis and recommendations for correcting the article.

Funding information. The research was done with the support of the Russian Science Foundation (project No. 21 - 71 - 20050).

Введение. В последнее время на территории Ростовской области отмечается возникновение ряда серьезных экологических проблем. К ним, в частности, относится эвтрофикация вод Азовского моря и Цимлянского водохранилища, которая вызывает рост вредоносных и токсичных видов популяций фитопланктона [1]. Инженерные работы в акваториях рек и морей приводят к загрязнению прилегающих территорий, изменению популяционной структуры биоты и ухудшению условий воспроизводства ценных и промысловых рыб. Изменение климата на юге России привело к увеличению количества случаев затопления некоторых территорий в районе Таганрогского залива и поймы реки Дон, вызванных сгонно-нагонными явлениями. В последнее десятилетие в летний период несколько раз наблюдалось практически полное осушение русла реки Дон, что приводило к полной остановке судоходства. Чтобы спрогнозировать возникновение и развитие подобных случаев, спланировать пути устранения их последствий, оценить нанесенный ими ущерб, требуются современные программные комплексы, построенные с использованием высокоточных математических моделей, численных методов, алгоритмов и структур данных [2].

В основе математических моделей, используемых при прогнозировании природных и техногенных катастроф, лежат системы дифференциальных уравнений в частных производных, например уравнения Пуассона, Навье-Стокса, диффузии-конвекции-реакции, теплопроводности. Численное решение таких систем приводит к необходимости оперативного хранения больших объёмов данных (в массивах различной структуры) и решения систем сеточных уравнений высокой размерности, превышающих 109. Объём оперативной памяти, требуемой для хранения массивов данных при численном решении только одного уравнения Пуассона для трёхмерной области размерностью 103×103×103 попеременно-треугольным итерационным методом, составляет более 64 Гб. В случае численного решения комбинированных задач требуются сотни гигабайт оперативной памяти, которые могут быть доступны лишь при использовании суперкомпьютерных систем.

Проведенное ранее исследование показало, что ресурсов вычислительной системы, использующей только CPU, недостаточно для решения подобных научных задач [3]. Увеличение мощности и видеопамяти GPU позволило использовать для расчетов ресурсы видеоадаптеров [4]. Эффективность использования GPU зависит от применения параллельных алгоритмов для решения вычислительно-трудоемких задач водной экологии [5–7]. Частично решить проблемы нехватки памяти и вычислительной мощности на рабочих станциях можно установкой дополнительных видеоадаптеров в слоты PCI-E X16 непосредственно и в слоты PCI-E X1 с помощью переходников PCI-E X1–PCI-E X16. Таким образом, количество видеоадаптеров, установленных на одной рабочей станции, можно довести до 12 [8–11].

Всё большую популярность в научном сообществе приобретают гетерогенные вычислительные системы, которые позволяют использовать ресурсы CPU и GPU совместно. Применение таких систем дает возможность уменьшить время расчета научных задач [12–14]. Однако эффективное использование гетерогенной вычислительной среды предполагает модернизацию математических моделей, алгоритмов и программ, их численно реализующих. Гетерогенная система позволяет организовать процесс вычислений в параллельном режиме. При этом должны быть учтены принципиальные различия в построении программных систем, использующих совместно CPU и GPU.

Материалы и методы. Опишем предложенные математические модели вычислительной системы, расчетной сетки, а также метод декомпозиции расчетной области.

Пусть D — множество технических характеристик вычислительной системы, тогда:

D = D1U D2 U D3, где D1 — подмножество характеристик центральных процессоров (CPU) вычислительной системы; D2 — подмножество характеристик видеоадаптеров (GPU) вычислительной системы; D3 — подмножество характеристик оперативной памяти.

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

D1 = dd 1-1, d 1,2, d 1-3,d 1,4},(2)

где d ^1,1 — суммарное число ядер CPU; d ^1,2 — количество потоков, одновременно обрабатываемых одним ядром процессора CPU; d ^1,3 — тактовая частота, МГц; d ^1,4 — частота шины центрального процессора, МГц.

D2 = U DkGp„ ={ d 2|3 kGPU € KGPU, d2 € D^},(3)

kGPU € KGPU где KGPU ={1-•••-Ngpu } — множество индексов видеоадаптеров; NGPU — количество видеоадаптеров вычислительной системы; kGPU — индекс видеоадаптера.

Каждый видеоадаптер представим в виде кортежа:

D 2 ²   = ( d^2Д , d ²,2 V                                                    (4)

kGPU       kGPU kGPU где dk2,1   — объём видеопамяти видеоадаптера с индексом kGPU , Гб; dk2,2 — количество потоковых мультипроцессоров.

D ³ = ( d 3,1 , d 3.2),                                                          (5)

где d ^3,1 — суммарный объём оперативной памяти, Гб; d ^3,2 — тактовая частота оперативной памяти, МГц.

Пусть S — множество программных потоков, задействованных в вычислительном процессе, тогда:

S = S 1U S ²,

S1 = {1,..., Ns,},(6)

S2 = U S2 , S2 = (1,..., N., L kkS

GPU        GPU kGPU e KGPU где S1 — подмножество программных потоков, реализующих процесс расчета на CPU; S 2 — подмножество потоковых блоков CUDA, реализующих процесс расчета на потоковых мультипроцессорах GPU; NS1 — количество задействованных программных потоков CPU; Sk2   — подмножество потоковых блоков

CUDA, реализующих процесс расчета на потоковых мультипроцессорах GPU с индексом k_GPU ; K GpU = { 1,-., N_GPU } — множество индексов GPU; N_GPU — количество GPU в вычислительной системе; N_S 2   — количество задействованных потоковых блоков CUDA, реализующих процесс расчета на потоковых

SkGPU мультипроцессорах GPU с индексом kGPU .

Пусть E — множество идентификаторов программных потоков. Тогда для идентификации программных потоков в вычислительной системе каждому элементу множества программных потоков S поставим в соответствие кортеж e из двух элементов:

V s e S 3 e e E : e = ^ n , , n^ ,

где n_d — индекс вычислительного устройства в вычислительной системе; n_t — индекс программного потока CPU или потокового блока GPU.

nd ={'

0,       s e S ¹

.K GPU , s ^{e S 2} ’

I K_S , ,    s e S ¹

ⁿt = ] K,₂    , s e S²

I ^S2GPU           k^GPU

.

Возьмём расчётную область со следующими параметрами: l_x — характерный размер по оси Ox ; l_y — по оси Oy ; l_z — по оси Oz . Сопоставим с указанной областью равномерную расчётную сетку следующего вида:

W = { x i = ih x , y _; = jh_y , z k = kh_z ;

i = 0, n x - 1, j = 0, n y - 1, k = 0, n z - 1;                                               (10)

(nx - 1) hx = lx, (ny-1) hy = ly, (nz - 1) hz = lz}, где hx , hy , hz — шаги расчётной сетки по соответствующим пространственным направлениям; nx , ny , nz количество узлов расчётной сетки по соответствующим пространственным направлениям.

Тогда множество узлов расчётной сетки представим в виде:

^G = { g i , j , k , i = 0, ⁿx ^- 1, ^j = 0, ⁿy ^{- 1, k} = 0, ⁿz ^- 1 } ,

gi,j,k = \xi, yj, zk/, где gi, j,k — узел расчётной сетки.

Число узлов расчётной сетки N_G вычисляется по формуле:

^NG = n x • n y • n z .

Под подразделом расчётной сетки G k с G (далее — подраздел) будем понимать подмножество узлов расчётной сетки G .

G = U G k 1 = { g k 1 | 3 k 1 e K_kx , g^k 1 e G k 1 } , П G k =0 , ^k 1 ^{e K}k1                                                              ^k1 ^{e K}k1

где K_k] ={ 1,..., N_{k ]} } — множество индексов подразделов G^k > расчётной сетки G ; N_ki — количество подразделов G k > ;

K_ki , N_ki с N ; N — множество натуральных чисел; k 1 — индекс подраздела G^k 1 .

Так как G^k 1 с G , тогда:

G k 1 = { g ; ^, _к , i = 0, n x - 1, j = 0, n^k k - 1, k = 0, n z - 1 } ,

где g kj _k — узел подраздела k 1 ; знак ~ обозначает принадлежность к подразделу; j — индекс узла подраздела k ₁ по координате y ; n^k_y ¹ — количество узлов в подразделе k ₁ по координате y .

g kj , k = ( x , y j

Г xi = ihx , yj = 1

, z_k

,

k —1

I ⁿy + j I- ^hy , z k = ^khz ,

где n^b_y ¹ — количество узлов по координате y b ₁ -го раздела.

Под блоком расчётной сетки G^{k 1 , k 2} (далее — блок) будем понимать подмножество узлов расчётной сетки подраздела G^{k 1} .

G ‘ 1 =   и   G k 1 ■ ‘ 2 = { g^ 1 3 ‘ 2 G K_{k k} , g^ 2 G G ‘ 1 > ‘ 2 },   n   G k 1 - ‘ 2 = 0 ,                     (16)

‘ 2 G K ‘ 1, ‘ 2                  ^{1                                1} ’ ²                               ’ ‘ 2 G K ‘1,‘2

где К_к1кг = { 1,..., N_{k k}. ) — множество индексов блоков G^k 1 ■ ^{k 2} подраздела G ‘ 1 ; N^ ^ — количество блоков G ‘ ■ ^k 2 ;

K_{k k z}, N_k, _kz c N ; k ₂ — индекс блока G^k 1 ■ ^k 2 подраздела G ‘ 1 .

Так как G ‘ - ‘ 2 c G ‘ 1 , тогда:

G ^{kp k 2} = { g ^y ‘ ² , i = 0, П х - 1, j = 0, _пм ² - 1, k = 0, n z - 1 } ,                               (17)

где gk1, k2 — узел блока k, k; знак ^ обозначает принадлежность к блоку; ˆj — индекс узла блока k, k по i, j,k координате y ; nyk1,k2 — количество узлов в блоке k1,k2 по координате y .

g

k ₁, k ₂

ˆ

= x^xiy y_P z k

■

I ^k 1 -1 N k 1, ‘ 2 b b ■ ^xi = ih x , У , =I S S n y^-b² + j

^ b =1 b 2 =1

z k = ^k h z ,

где n_y^{b 1 , b 2} — количество узлов блока b ₁, b ₂.

Под фрагментом расчётной сетки G^k 1 ^{, k} 2 ^{, k} 3 (далее — фрагмент) будем понимать подмножество узлов расчётной сетки блока G ^{k 1 , k 2} подраздела G^{k 1} .

G ‘ 1, ‘ 2 = [J    G ‘ 1 , ‘ 2 , ‘ 3 = { g ‘ 1 , ‘ 2 , ‘ 3 | 3 ‘ 3 G K , g ‘ 1 , ‘ 2 , ‘ 3 G G ‘ 1 , ‘ 2 , ‘ 3 } ,

‘ 3 G K_k                                                ‘ 3

Pl   G ‘ 1 , ‘ 2 , ‘ 3 =0 ,

‘ 3 ^G K k ₃

где    ^Kk 1 , k 2 , k 3 = { 1,..., ^Nk 1 , k 2 , k 3 )

множество индексов фрагментов G^k 1 ^{, k} 2 ^{, k} 3 блока G^k 1 ^{, k} 2 подраздела G^k 1 ;

N_{k} kl k}, — количество фрагментов G ‘ ¹ ■ ‘ ² ■ ‘ ³ ; K_k, _{kj ki} , N_k, _{kj kj} c N ; k ₃ — индекс фрагмента G ‘ 1 ■ ‘ ² ■ ‘ ³ блока G^k 1 ■ ‘ ² подраздела G^k 1 .

Каждому индексу k ₃ фрагмента G^{k 1 , k 2 , k 3} поставим в соответствие кортеж индексов k ₄, k ₅, предназначенный для хранения координат фрагмента в плоскости xOz , где k ₄ — индекс фрагмента по координате x ; k ₅ — индекс фрагмента по координате z .

‘3 = ‘4 + Kk4 ■‘5, где k4 — индекс фрагмента по координате x ; k5 – индекс фрагмента по координате z . Количество фрагментов Gk1,k2,k3 блока Gk1,k2 вычислим по формуле:

Kk3 = Kk4 ■ Kk5 , где Kk — количество фрагментов по оси Ox; Kk — количество фрагментов по координате z .

Так как G‘^{1, k}2 ^{, k}3 c G ‘ ^{1, k 2} , тогда:

G‘1,‘2,‘3 =   <•_ . _ i,j,k

_ i = 0, nx -1, j = 0, /7y -

'_______

ˆ

^“       •__.       ------------------

1, k = 0, _ _z -

_____

1 ) , ^_

_ где g_,_ — узел фрагмента; знак „ обозначает принадлежность к фрагменту; i, ‘ — индексы узла фрагмента

_____

по координатам x , z ; n_x , n_z — количество узлов расчётной сетки в фрагменте по координатам х , z ; l_x , l_z — размеры фрагмента по координатам х , z .

_ g _ , • = (x., y, i,j,k i j

, z_k

,

‘ 4 —1

•_____

ˆ

‘ 5 -1

•    •                                   _ xi = I S nb + i I hx, У, = jhy, z‘ = I S nb + ‘ I hz,

где n b — количество узлов b -го фрагмента.

Введем множество сопоставлений блоков расчётной сетки программным потокам M ¹

M ¹ = U I U M 1 1 , k 2

^k 1 ^{G K}k 1 I ‘ 2 ^{g K}k2

,

где M_k ¹ _{, k} — элемент множества M ¹ .

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

Пусть M_k ¹ _{, k} — сопоставление блока G ^k 1 ^{, k} 2 программному потоку s_{k , k} , тогда:

м Ь = Gk^{k 1} , ^k 2 , s V                                         (25)

k 1 , k 2                            k ₁, k ₂

где s^ _k e S — программный поток, вычисляющий блок G k- ^k 2 .

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

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

Предположим, что все вычислительные устройства используются для расчётов. Тогда суммарная производительность вычислительной системы P ^ вычисляется по формуле:

N GPU

P , = P cpu ■ N s 1 + S P ■ N S 2 ,                                       (26)

b =1

где P_CPU — производительность одного потока CPU; N_S 1 — число программных потоков, реализующих процесс расчета на CPU; P^b — производительность GPU с индексом b на одном потоковом мультипроцессоре; N_S^b 2 — количество потоковых блоков CUDA, реализующих процесс расчета на потоковых мультипроцессорах GPU.

Тогда рассчитать количество узлов расчётной сетки n^b_y индексом b можно по формуле:

в подразделе по координате y для каждого GPU с

n by =

Pb

GPU

P l

■ n y .

В процессе вычисления по формуле (27) получим остаток — некоторое количество узлов расчётной сетки. Эти узлы будут располагаться в оперативной памяти. Количество оставшихся узлов n^b_y по координате y рассчитывается по формуле:

N GPU

n yL = n y ^— S n b . b =1

Для вычисления количества узлов по координате y в блоках расчётной сетки, обрабатываемых потоковыми мультипроцессорами GPU, воспользуемся формулами:

^b n yGT

n y ^b

Nb —

S 2

Nb

1 J ^,

—1

b = 1, N^b — 1, S ²

n bGTL = n b

—

^S n bGT , b =1

где n_y^b_GT — количество узлов по координате y в блоках расчётной сетки, обрабатываемых потоковыми мультипроцессорами GPU с индексом b , кроме последнего блока; n_y^b_GTL — количество узлов по координате y в последнем блоке расчётной сетки, обрабатываемом потоковыми мультипроцессорами GPU с индексом b .

Для вычисления количества узлов расчётной сетки по координате y в блоках, обрабатываемых программными потоками, реализующими процесс расчета на CPU, воспользуемся формулами:

nCPU nyCT =

N —

S ¹

CPU nyCTL   ny

¹ J ,

■ nyCT ■ ( N S 1 ^— 1 ) ,

где n_yCT — количество узлов расчётной сетки по координате y , обрабатываемых программными потоками CPU,

кроме последнего потока; n_yCTL — количество узлов расчётной сетки по координате y , обрабатываемых программными потоками CPU, в последнем потоке.

Рассчитаем количество фрагментов расчётной сетки по координате y :

N GPU

N f = N + S N^b . y         s ¹       ь =1      s ²

Пусть задано количество фрагментов N_x^f и N_z^f по координатам x и z соответственно. Тогда количество

узлов расчётной сетки по координате x вычисляется по формулам:

nx

N f - ¹

n f = x

nfL = n - nf ■(Nf -1), x xx x где nx — количество узлов расчётной сетки по координате x во всех фрагментах, кроме последнего; nXL — количество узлов расчётной сетки по координате x в последнем фрагменте.

nz

N f - 1

n zf =

Аналогично вычисляется количество узлов расчётной сетки по координате z :

nfL = nz -nf ■(Nf -1), где nzf — количество узлов расчётной сетки по координате z во всех фрагментах, кроме последнего; nzfL — количество узлов расчётной сетки по координате z в последнем фрагменте.

Опишем модель параллельно-конвейерного метода. Пусть на M ¹ необходимо организовать параллельный процесс вычислений некоторой функции F , причем вычисления в каждом фрагменте G^' k ^,k3 зависят от значений в соседних фрагментах, каждый из которых имеет хотя бы один из индексов по координатам x , y и z на единицу меньший, чем у текущего (рис. 1).

Для организации параллельно-конвейерного метода введём множество кортежей A , задающих соответствия a между идентификаторами программных потоков e , обрабатывающих фрагменты G^ k2, ^k з , номерам шагов параллельно-конвейерного метода r :

V e е E 3 a е A : a = ee , G^k 1 ■ ^k 2 ■ ^k 3 , r ^, (34)

где r = 1, N r — номер шага параллельно-конвейерного метода; N r — число шагов параллельно-конвейерного метода, вычисляемое по формуле:

N = N f N f + N f - 1.                                   (35)

r xz y

Полная загрузка всех вычислителей в предлагаемом параллельно-конвейерном методе начинается с шага r     = N f и заканчивается на шаге r = N f N f . При этом общее количество шагов с полной загрузкой

100 START у                                          100 stop x z вычислителей NrPAR составит:

^NrPAR = r 100 STOP - r 100 START ^{+ 1} = ^{N f N}z - ^Ny ^{+ 11                                  (36)}

Время вычислений некоторой функции F параллельно-конвейерным методом запишем в виде:

N r

T_M , =S max(T),                                       (37)

r =1

где T a — вектор значений затрат времени на обработку фрагментов в параллельном режиме.

		0	1	2	3			N xf -1
e 0 —	* 0	(0,0,0) r =0	(1,0,0) r=1	(2,0,0) r=2	(3,0,0) r =3		N y 0    \ \	( N x -1,0,0) r = N xf -1	x
e 1 —	> 1	(0,1,0) r=1	(1,1,0) r=2	(2,1,0) r=3	(3,1,0)		N y 1	( N x -1,1,0) r = N x f
e 2 —	> 2	(0,2,0) r=2	(1,2,0) r=3	(2,2,0)	(3,2,0)		N y 2	( N x -1,2,0) r = N xf +1
e 3 —	> 3	(0,3,0) r=3	(1,3,0)	(2,3,0)	(3,3,0)	>	N y 3 (  (	( N x -1,3,0) r = N xf +2
	y

Рис. 1. Параллельно-конвейерный вычислительный процесс

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

Результаты исследования . Вычислительные эксперименты проводились на высокопроизводительной вычислительной системе К-60 Института прикладной математики им. М. В. Келдыша Российской академии наук. Использовалась секция с GPU, каждый узел которой оснащён двумя процессорами Intel Xeon Gold 6142 v4, четырьмя видеоадаптерами Nvidia Volta GV100GL и 768 Гб оперативной памяти.

Вычислительный эксперимент состоял из двух этапов — подготовительного и основного. На подготовительном этапе проверялась корректность декомпозиции расчетной области на подразделы, блоки и фрагменты путем поэлементного сравнения значений в узлах исходной сетки и в фрагментах, полученных в результате декомпозиции. Затем проверялась работа алгоритма управления потоками, в процессе которого фиксировалось время, затрачиваемое на расчет 1, 8, 16 и 32 фрагментов расчетной сетки размерностью 50 узлов по пространственным координатам x , y и z тем же количеством потоков CPU N_S 1 итерационным попеременно-треугольным методом в параллельном режиме. Выполнялось 10 повторов с вычислением среднего арифметического T a и стандартного отклонения а . По полученным данным вычислялось время T a ¹ = T a I N S 1 , затрачиваемое каждым потоком на обработку одного фрагмента расчетной сетки и ускорение E = T a ¹ ( N S 1 ) / T a ¹ (1), равное отношению времени T_a ¹ ( N_S 1 ) обработки одного фрагмента N_S 1 потоками к соответствующему времени обработки одним потоком T_a ¹ (1). Экспериментальные данные приведены в таблице 1. Проведённый эксперимент показал, что стандартное отклонение имеет наименьшее значение в случае использования 32 параллельных потоков CPU и составляет 0,026 мс, то есть использование 32 параллельных потоков CPU при расчёте 32 фрагментов расчётной сетки даёт более равномерную по времени загруженность программных потоков, что в целом повышает эффективность работы вычислительного узла. При этом среднее значение расчета одного фрагмента составило 4,14 мс. Зависимость ускорения E от числа потоков получилась линейная E = 0,603 + 0,804 N S 1 , с коэффициентом детерминации, равным 0,99. Получили, что с увеличением числа потоков ускорение разработанного алгоритма возрастает. Это свидетельствует об эффективном использовании подсистемы при работе с памятью.

Таблица 1

Результаты подготовительного этапа вычислительного эксперимента

NS1 max (Ta ) , мс а, мс Ta, мс E 1 3,38 0,141 3,38 1,00 8 3,66 0,042 0,46 7,39 16 3,94 0,028 0,25 13,73 32 4,14 0,026 0,13 26,13

На основном этапе вычислительного эксперимента трёхмерная расчетная область, имеющая размеры 1600, 1600, 200 по пространственным координатам x , y и z соответственно, разбивалась на 32 фрагмента по 50 узлов по каждой из координат x и y . Разбиение на фрагменты по координате z приведено в таблице 2. Для каждого варианта декомпозиции с десятикратной повторностью замерялось время обработки всей расчетной сетки предлагаемым параллельно-конвейерным методом и вычислялось его среднее значение T_pm . Ускорение E_pm вычислялось как отношение T_pm ко времени T_sm расчета последовательной версией алгоритма, равному 6963 мс. Получено уравнение регрессии E pm = 7,35 + 1,97 ln( N Z ) с коэффициентом детерминации, равным 0,94. Анализ результатов основного этапа вычислительного эксперимента показал существенное замедление роста E_pm при N / >  10 . Поэтому делаем вывод о том, что оптимальным является разбиение на фрагменты по координате z на величину, не превышающую 10.

Таблица 2

Результаты основного этапа вычислительного эксперимента

N z f	nzf	Tpm , мс	E pm
1	200	1033,20	6,74
2	100	779,00	8,94
4	50	651,90	10,68
8	25	588,35	11,84
20	10	550,22	12,66

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

Результаты, полученные в ходе вычислительных экспериментов, подтверждают эффективность разработанного метода. Подтверждена и корректность декомпозиции расчетной области на подразделы, блоки и фрагменты. Проверена работа алгоритма управления потоками, при этом выявлено, что стандартное отклонение имеет наименьшее значение в случае использования 32 параллельных потоков CPU и составляет 0,026 мс, то есть использование 32 параллельных потоков CPU при расчёте 32 фрагментов расчётной сетки даёт более равномерную по времени загруженность программных потоков. При этом среднее значение расчета одного фрагмента составило 4,14 мс.

Результаты обработки замеров времени расчетов предлагаемым параллельно-конвейерным методом показали существенное замедление роста ускорения при разбиении на фрагменты по координате z при N_z^f >  10 . Получили, что оптимальным является разбиение на фрагменты по координате z на величину, не превышающую 10.