Введение. В последнее время на территории Ростовской области отмечается возникновение ряда серьезных экологических проблем. К ним, в частности, относится эвтрофикация вод Азовского моря и Цимлянского водохранилища, которая вызывает рост вредоносных и токсичных видов популяций фитопланктона [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.

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

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

(1)

где

— подмножество характеристик центральных процессоров (CPU) вычислительной системы; — подмножество характеристик видеоадаптеров (GPU) вычислительной системы;

— подмножество характеристик оперативной памяти.

(2)

где

— суммарное число ядер CPU;

— количество потоков, одновременно обрабатываемых одним ядром процессора CPU;

— тактовая частота, МГц;

— частота шины центрального процессора, МГц.

, (3)

где

— множество индексов видеоадаптеров;

— количество видеоадаптеров вычислительной системы;

— индекс видеоадаптера.

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

, (4)

где

— объём видеопамяти видеоадаптера с индексом , Гб;

— количество потоковых мультипроцессоров.

, (5)

где

— суммарный объём оперативной памяти, Гб;

— тактовая частота оперативной памяти, МГц.

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

(6)

где

— подмножество программных потоков, реализующих процесс расчета на CPU;

— подмножество потоковых блоков CUDA, реализующих процесс расчета на потоковых мультипроцессорах GPU;

— количество задействованных программных потоков CPU;

— подмножество потоковых блоков CUDA, реализующих процесс расчета на потоковых мультипроцессорах GPU с индексом ;

— множество индексов GPU;

— количество GPU в вычислительной системе;

— количество задействованных потоковых блоков CUDA, реализующих процесс расчета на потоковых мультипроцессорах GPU с индексом .

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

, (7)

где

— индекс вычислительного устройства в вычислительной системе;

— индекс программного потока CPU или потокового блока GPU.

, (8)

. (9)

Возьмём расчётную область со следующими параметрами:

— характерный размер по оси ;

— по оси ; — по оси .

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

(10)

где

, , — шаги расчётной сетки по соответствующим пространственным направлениям;

, , — количество узлов расчётной сетки по соответствующим пространственным направлениям.

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

(11)

где

— узел расчётной сетки.

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

. (12)

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

, (13)

где

— множество индексов подразделов расчётной сетки ;

— количество подразделов

;

— множество натуральных чисел;

— индекс подраздела .

Так как , тогда:

, (14)

где

— узел подраздела ;

знак обозначает принадлежность к подразделу;

— индекс узла подраздела по координате ;

— количество узлов в подразделе по координате .

(15)

где

— количество узлов по координате -го раздела.

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

, (16)

где

— множество индексов блоков подраздела ;

— количество блоков ;

;

— индекс блока подраздела .

Так как , тогда:

, (17)

где

— узел блока ;

знак обозначает принадлежность к блоку;

— индекс узла блока по координате ;

— количество узлов в блоке по координате .

(18)

где

— количество узлов блока .

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

(19)

где

— множество индексов фрагментов блока подраздела ;

— количество фрагментов ;

;

— индекс фрагмента блока подраздела .

Каждому индексу фрагмента поставим в соответствие кортеж индексов , предназначенный для хранения координат фрагмента в плоскости , где — индекс фрагмента по координате ; — индекс фрагмента по координате .

, (20)

где

— индекс фрагмента по координате ;

– индекс фрагмента по координате .

Количество фрагментов блока вычислим по формуле:

, (21)

где

— количество фрагментов по оси

— количество фрагментов по координате .

Так как , тогда:

, (22)

где

— узел фрагмента;

знак ‿ обозначает принадлежность к фрагменту;

— индексы узла фрагмента по координатам , ;

, — количество узлов расчётной сетки в фрагменте по координатам , ;

, — размеры фрагмента по координатам , .

(23)

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

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

, (24)

где — элемент множества .

Пусть — сопоставление блока программному потоку , тогда:

, (25)

где — программный поток, вычисляющий блок .

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

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

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

, (26)

где

— производительность одного потока CPU;

— число программных потоков, реализующих процесс расчета на CPU;

— производительность GPU с индексом на одном потоковом мультипроцессоре;

— количество потоковых блоков CUDA, реализующих процесс расчета на потоковых мультипроцессорах GPU.

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

. (27)

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

. (28)

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

(29)

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

— количество узлов по координате в последнем блоке расчётной сетки, обрабатываемом потоковыми мультипроцессорами GPU с индексом .

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

(30)

где

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

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

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

. (31)

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

(32)

где

— количество узлов расчётной сетки по координате во всех фрагментах, кроме последнего;

— количество узлов расчётной сетки по координате в последнем фрагменте.

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

(33)

где

— количество узлов расчётной сетки по координате z во всех фрагментах, кроме последнего;

— количество узлов расчётной сетки по координате в последнем фрагменте.

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

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

, (34)

где

— номер шага параллельно-конвейерного метода;

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

. (35)

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

. (36)

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

, (37)

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

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

Вычислительный эксперимент состоял из двух этапов — подготовительного и основного. На подготовительном этапе проверялась корректность декомпозиции расчетной области на подразделы, блоки и фрагменты путем поэлементного сравнения значений в узлах исходной сетки и в фрагментах, полученных в результате декомпозиции. Затем проверялась работа алгоритма управления потоками, в процессе которого фиксировалось время, затрачиваемое на расчет 1, 8, 16 и 32 фрагментов расчетной сетки размерностью 50 узлов по пространственным координатам , и тем же количеством потоков CPU итерационным попеременно-треугольным методом в параллельном режиме. Выполнялось 10 повторов с вычислением среднего арифметического и стандартного отклонения s. По полученным данным вычислялось время , затрачиваемое каждым потоком на обработку одного фрагмента расчетной сетки и ускорение равное отношению времени обработки одного фрагмента потоками к соответствующему времени обработки одним потоком Экспериментальные данные приведены в таблице 1. Проведённый эксперимент показал, что стандартное отклонение имеет наименьшее значение в случае использования 32 параллельных потоков CPU и составляет 0,026 мс, то есть использование 32 параллельных потоков CPU при расчёте 32 фрагментов расчётной сетки даёт более равномерную по времени загруженность программных потоков, что в целом повышает эффективность работы вычислительного узла. При этом среднее значение расчета одного фрагмента составило 4,14 мс. Зависимость ускорения от числа потоков получилась линейная , с коэффициентом детерминации, равным 0,99. Получили, что с увеличением числа потоков ускорение разработанного алгоритма возрастает. Это свидетельствует об эффективном использовании подсистемы при работе с памятью.

На основном этапе вычислительного эксперимента трёхмерная расчетная область, имеющая размеры 1600, 1600, 200 по пространственным координатам , и соответственно, разбивалась на 32 фрагмента по 50 узлов по каждой из координат и . Разбиение на фрагменты по координате приведено в таблице 2. Для каждого варианта декомпозиции с десятикратной повторностью замерялось время обработки всей расчетной сетки предлагаемым параллельно-конвейерным методом и вычислялось его среднее значение . Ускорение вычислялось как отношение ко времени расчета последовательной версией алгоритма, равному 6963 мс. Получено уравнение регрессии с коэффициентом детерминации, равным 0,94. Анализ результатов основного этапа вычислительного эксперимента показал существенное замедление роста при . Поэтому делаем вывод о том, что оптимальным является разбиение на фрагменты по координате на величину, не превышающую 10.

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

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

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