Solution to the problem of thermal conductivity in a two-dimensional formulation on a graphics processing unit using parallel computing

Full Text

Введение

Программно-аппаратная архитектура параллельных вычислений позволяет существенно повысить скорость вычислений. Технология CUDA применяется в здравоохранении для выявления полезной площади на рентгеновских снимках и сжатия изображения, при этом такой алгоритм работает в 34 раза быстрее, чем последовательный алгоритм, выполняемый на центральном процессоре [1, с. 1144]. Технология параллельного программирования эффективна при решении дифференциальных уравнений в частных производных [2, с. 41]. К сожалению, большинство научных статей про использование технологии параллельных вычислений пишутся без приведения программного кода. Исключением явилась статья [3, с. 9], в которой описывается решение задачи одномерной теплопроводности с использованием явной схемы. Статья [3, с. 6] подтолкнула к изучению технологии параллельного программирования и написанию статей [4, с. 9; 5, с. 24].

В статье [6, с. 346] решается аналогичная задача двухмерной теплопроводности, но с использованием явной схемы при количестве шагов по времени 100. Ускорение составило 48 раз. Явная схема с использованием параллельных вычислений для решения уравнения Навье – Стокса отображена в статье [7, с. 123]. Среднее ускорение составило 32 раза, а максимальное – 100. В статье [8, с. 431] авторы решают задачу теплопроводности явным методом и выделяют его недостаток: для обеспечения устойчивости метода требуется достаточно большое количество узлов. К недостаткам неявной схемы и схемы Кранка – Николсона относится высокая вычислительная сложность на каждой итерации, а преимуществом является устойчивость. Численные результаты решения задачи одномерной теплопроводности [5, с. 37] показали, что среди явного, неявного и метода Кранка – Николсона на графическом процессоре самым быстрым является явный, затем следуют метод Кранка – Николсона и неявный. Так как метод Кранка – Николсона обладает вторым порядком точности, то он был выбран для моделирования задачи двухмерной теплопроводности.

Существуют специальные программные продукты [9, с. 551], позволяющие преобразовать код, написанный на языке C, в код программы CUDA C. За счет таких программ увеличивается скорость выполнения по сравнению с последовательными вычислениями на центральном процессоре (ЦП), но они медленнее, чем если бы код был написан на языке CUDA C. Для написания параллельного кода требуется графический процессор (ГП) фирмы NVIDIA c поддержкой ядер CUDA и пакет программирования CUDA SDK, который позволит взаимодействовать разработчику с видеокартой. Для упрощения использования технологии CUDA в пакете программирования SDK встроено несколько библиотек, которые позволяют использовать возможности CUDA без необходимости писать код самостоятельно [10, с. 90; 11, с. 173–209].

При этом максимальное ускорение программы, выполняющейся параллельно, рассчитывается по формуле Амдала:

$S=\frac{1}{(1-P)+\frac{P}{N}},$ (1)

где P – это часть времени выполнения программы, которая может быть распараллелена на N процессоров.

Практическое использование расчетов на ЦП и ГП показывает, что максимальную производительность можно получить при абсолютно параллельной задаче [4, с. 10], а при решении конечно-разносной схемы Кранка – Николсона требуется синхронизация ячеек (за счет этого тратится общее время выполнения) [5, с. 31].

В статье [5, с. 27] решалась задача одномерной теплопроводности на ГП с использованием технологии параллельного программирования CUDA. Наибольший выигрыш по производительности получился при большем количестве параллельно запущенных нитей с наименьшим количеством итераций.

Результаты и обсуждение

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

• операционная система Windows 10;
• среда программирования Visual Studio 2019 Community Edition;
• средства разработки на языке CUDA (CUDA Toolkit SDK версии 11.5);
• графический драйвер версии 511.79.

Расчеты производились на ЦП AMD Phenom 955 BE и графической карте начального уровня Gigabyte GeForce GT 1030 c частотой ядер 1252 MHz и количеством CUDA ядер 384.

При написании кода на языке CUDA используется как ЦП, так и ГП. Код на ГП пишется на языке CUDA, а на ЦП – на языках С или Python, в данном случае используется язык C.

Сравнение численных методов конечно-разностной аппроксимации показало, что схема Кранка – Николсона является более точной по сравнению с явной и неявной схемами [12, с. 1809], поэтому она использовалась при решении задачи двухмерной теплопроводности несмотря на то, что большая производительность и максимальная нагрузка получаются при реализации явных схем [13, с. 10].

Уравнение Фурье – Кирхгофа описывает нестационарный перенос тепла:

$c\text{ρ}\frac{\partial T}{\partial t}=\frac{\partial }{\partial x}\left(\text{λ}\frac{\partial T}{\partial x}\right)+\frac{\partial }{\partial y}\left(\text{λ}\frac{\partial T}{\partial y}\right)+\frac{\partial }{\partial z}\left(\text{λ}\frac{\partial T}{\partial z}\right)+Q(x,y,z,t,T)$, (2)

где с – удельная теплоемкость, Дж/(кг·К); ρ – плотность, кг/м³; T(x, y, z, τ) – поле температуры в трехмерном пространстве; t – время, с; λ – коэффициент теплопроводности, Вт/(м·К); Q(x, y, z, τ, T) – мощность внутренних источников тепловыделения, Вт.

Уравнение Фурье – Кирхгофа описывает множество вариантов процесса теплопроводности. Рассмотрим решение двухмерной задачи теплопроводности через изолированную пластину с постоянными коэффициентами. На нижней границе пластины, слева и справа поддерживается постоянная температура T_Н = T_Л = T_П, вверху пластины T_В температура границы пластины равна температуре окружающей среды. Начальная температура пластины равна T₀.

Уравнение Фурье – Кирхгофа (2) в двухмерном случае без источников тепла внутри пластины будет выглядеть следующим образом:

$с\text{ρ}\frac{\partial T}{\partial \tau }=\text{λ}\left(\frac{{\partial }^{2}T}{\partial x^2}+\frac{{\partial }^{2}T}{\partial y^2}\right)\mathrm{,0}<x<L\mathrm{,0}<y<L,$ (3)

где L – длина и высота квадратной пластины, м.

Выразив коэффициент теплоотдачи через удельную теплоемкость, плотность и коэффициент температуропроводности $a=\frac{\lambda }{c\text{ρ}}$, получим уравнение (3) в виде

$\frac{\partial T}{\partial \tau }=a\left(\frac{{\partial }^{2}T}{\partial x^2}+\frac{{\partial }^{2}T}{\partial y^2}\right)\mathrm{,0}<x<L\mathrm{,0}<y<L,$ (4)

Начальные условия запишутся следующим образом:

$\tau =0:T=T_0\mathrm{,0}\le x\le L\mathrm{,0}\le y\le L$. (5)

Граничные условия для уравнения (4)

$x\text{}=0:T=T_{Л},\tau >0;$ (6)

$x\text{}=L:T=T_{П},\tau >0;$ (7)

$y\text{}=0:T=T_{В},\tau >0;$ (8)

$y=L:T=T_{Н},\tau >0.$ (9)

Задача решается методом конечных разностей на равномерной квадратной сетке. Для этого пластина разбивается по длине и высоте на N-1 равных промежутков, а для аппроксимации уравнения теплопроводности применяется десятиточечный шаблон конечно-разностной схемы Кранка – Николсона (рис. 1).

Рисунок 1 – Десятиточечный шаблон схемы Кранка – Николсона

При одностороннем нагреве снизу использовались граничные условия первого рода. Для изолированных сторон (левая, верхняя и нижняя) получим:

$\frac{\partial T(x_{\min },y,\text{τ})}{\partial \text{τ}}=\frac{\partial T(x,y_{\max },\text{τ})}{\partial \text{τ}}=\frac{\partial T(x_{\max },y,\text{τ})}{\partial \text{τ}}=0;$ (10)

По всей нижней границе зададим постоянное значение температуры

$\frac{\partial T(x,y_{\min },\text{τ})}{\partial \text{τ}}=100.$ (11)

Заменим дифференциальные операторы в уравнении (4):

$\frac{\partial T}{\partial \text{τ}}=\frac{T_i^{n+1}-T_i^n}{\text{τ}},$ (12)

$\frac{{\partial }^{2}T}{\partial x^2}=\frac{T_{i-1}^{n+1}-2\cdot T_i^{n+1}+T_{i+1}^{n+1}}{\Delta x^2},$ (13)

$\frac{{\partial }^{2}T}{\partial y^2}=\frac{T_{j-1}^{n+1}-2\cdot T_j^{n+1}+T_{j+1}^{n+1}}{\Delta y^2}.$ (14)

Получим

$\frac{T_i^{n+1}-T_i^n}{\text{τ}}=a\left(\frac{T_{i-1}^{n+1}-2\cdot T_i^{n+1}+T_{i+1}^{n+1}}{\Delta x^2}+\frac{T_{j-1}^{n+1}-2\cdot T_j^{n+1}+T_{j+1}^{n+1}}{\Delta y^2}\right)$, (15)

где $T_i^{n+1}$ – значение температуры в ij-м узле разностной сетки в последующий момент времени, $T_i^n$ – значение температуры в ij-м узле разностной сетки в n-й момент времени.

ЦП предназначен для выполнения одного потока инструкций, ГП – для выполнения множества потоков. Для ускорения на ГП используется быстрая общая память, размер которой ограничен. ЦП поддерживает один или два потока на ядро. ГП поддерживает до 1024 потоков на потоковый процессор.

На рисунке 2 показана структура потоков на ГП. Самая маленькая единица в архитектуре CUDA – это поток или нить. Нити объединяются в блоки. Блоки в свою очередь объединяются в сетку.

Рисунок 2 – Структура нитей, блоков и сетки

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

Листинг 1

int idx = blockIdx.x * blockDim.x + threadIdx.x;

где blockIdx.x – индекс блока, blockDim.x – размерность блока, threadIdx.x – индекс нити в блоке.

Реализованная программа на ГП состоит из основного метода mai, который пишется на языке С, в нем происходит: 1) задание начальных параметров; 2) выделение памяти на ЦП и ГП; 3) копирование данных из ЦП в ГП; 4) выполнение расчетов на ГП; 5) копирование данных из ГП в ЦП; 6) вывод результатов и их обработка; 7) освобождение ресурсов на ГП и ЦП.

В статье [5, с. 10] решалась аналогичная задача в одномерной постановке, поэтому приведем отличный от одномерной постановки код.

На ГП используется функция пятиточечного шаблона, и на одной итерации для каждой нити она вызывается дважды для расчета на текущем и следующем шаге. Код приведен в листинге 2, а её графическая интерпретация показана на рисунке 1.

Листинг 2

__device__ double SH5(double L, double R, double T, double D, double M, double dx, double dy)

{

double res;

res = (L – 2 * M + R) / dx / dx + (D – 2 * M + T) / dy / dy;

return res;

}

Как видно из приведенного листинга, данная функция отличается от языка С только наличием специального слова из синтаксиса языка CUDA __device__, которое означает, что функция будет выполняться на ГП.

Основная функция (листинг 3) выполняется на ГП, принимает на вход одномерный массив данных и получает на выходе одномерный массив:

Листинг 3

__global__ void implicit(double* ziro, double* res)

{

//Задание начальных параметров расчета …

…

__shared__ double arr[SIZE]; //n

__shared__ double arrs[SIZE]; //n + 1

__shared__ double arrt[SIZE]; //t

__shared__ double _max; //Переменная, которая доступна всем потокам

_max = 1;

int idx = blockIdx.x * blockDim.x + threadIdx.x;

arr[idx] = ziro[idx];

arrs[idx] = ziro[idx];

arrt[idx] = ziro[idx];

__syncthreads(); //синхронизация потоков

for (int j = 0; j < 1000000; j++) { //количество итераций

_max = 1;

while (_max > eps) { //Пока заданная точность не достигнута

//В центре пластины

if ((idx >= _N) && (idx < SIZE – _N) && (idx % _N != 0) && (idx % _N != _N – 1)) { //Находим значение на следующем шаге

arrs[idx] = arr[idx] + gamma * ((arr[idx] – arrt[idx]) / dt – a * (SH5(arr[idx – 1], arr[idx + 1], arr[idx – _N], arr[idx + _N], arr[idx], _gdx, _gdy) + SH5(arrt[idx – 1], arrt[idx + 1], arrt[idx – _N], arrt[idx + _N], arrt[idx], _gdx, _gdy)) / 2.0);

_max = fabs(arr[_N + 1] – arrs[_N + 1]);

if (_max < fabs(arr[idx] – arrs[idx])) {

_max = fabs(arr[idx] – arrs[idx]);

}

}

__syncthreads();

arr[idx] = arrs[idx];

}

arrt[idx] = arr[idx]; //Получаем значение через dt

__syncthreads();

}

res[idx] = arr[idx]; //Для вывода результат

__syncthreads();

}

Как видно из приведенного кода, __global__ обозначает, что функция будет выполняться на языке CUDA. Основной цикл for определяет количество шагов по времени. А условный цикл while определяет условие достижения заданной точности eps.

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

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

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

В центре пластины используется функция пятиточечного шаблона SH5(), которая приведена в листинге 2. Переход от одномерного массива к квадратному двухмерному представим графически (рис. 3). В данном случае для массива 4 на 4 остаются четыре центральные точки, для которых применяется пятиточечный шаблон. На рисунке 3 приведен шаблон для первой центральной точки с индексом 6.

Рисунок 3 – Переход от одномерного массива к двухмерному

В таблице 1 приведены результаты сравнения скоростей вычислений на ЦП и ГП при размерности двухмерного массива 16 на 16 ячеек и количестве итераций от 10 до 10⁶. Таблицы 2 и 3 показывают скорости выполнения решения задачи двухмерной теплопроводности с количеством ячеек 1024. В таблице 2 используются числа с одинарной точностью (float), а в таблице 3 используются числа с двойной точностью (double). Как видно из таблицы 1, скорость расчета на ГП быстрее, чем на ЦП, от 1,33 до 1,71 раза.

Таблица 1 – Сравнение времени выполнения алгоритма для массива 16 на 16 с использованием чисел одинарной точности

Таблица 2 – Сравнение времени выполнения алгоритма для массива 32 на 32 с использованием чисел одинарной точности

Таблица 3 – Сравнение времени выполнения алгоритма для массива 32 на 32 с использованием чисел двойной точности

Следует отметить, что расчеты с двойной точностью выполняются быстрее, за исключением расчета на ГП при количестве итераций 10⁶. А максимальная производительность ГП по отношению к ЦП получается при меньшем количестве итераций.

Заключение и выводы

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

1. Литературный обзор показал, что быстродействие алгоритмов на графическом процессоре по сравнению с алгоритмами, выполняемыми на центральном процессоре, достигает от нескольких единиц до сотен раз.
2. Решение двухмерной задачи теплопроводности по схеме Кранка –Николсона не является абсолютно параллельной, и максимально возможное ускорение не достигается.
3. При размерности матрицы 16 на 16 и одинарной точности время выполнения на графическом процессоре оказалось в 1,32 – 1,72 быстрее, чем на центральном процессоре. При размерности матрицы 32 на 32 время выполнения на графическом процессоре оказалось в 3,66 – 6,07 быстрее, чем на центральном процессоре.
4. При вычислениях с двойной точностью наблюдается наибольшее ускорение в 71,89 при 10 итерациях расчета, если итераций более 104, то ускорение в расчетах на графическом процессоре с двойной точностью приближается к расчетам с одинарной точностью.

About the authors

Pavel A. Sechenov

Author for correspondence.
Email: pavesa89@mail.ru
ORCID iD: 0000-0003-4551-3822
SPIN-code: 6017-9147
Scopus Author ID: 57192101307
ResearcherId: Y-5965-2019

Candidate of technical sciences, Associate Professor of the Department of Applied Information Technologies and Programming

References

Supplementary files

Supplementary Files

Action

1. JATS XML

Download

2. Figure 1 – A ten– point template of the Krank - Nicholson scheme

Download (72KB)

Indexing metadata

3. Figure 2 – Structure of threads, blocks and grids

Download (36KB)

Indexing metadata

4. Figure 3 – Transition from one-dimensional array to two-dimensional

Download (129KB)

Indexing metadata