Реализация метода переменных направлений с использованием технологии CUDA

1. Введение

В настоящей заметке демонстрируется возможность выполнения расчетов на видеокартах (с применением технологии CUDA) при моделировании физических процессов и явлений на примере решения трехмерного уравнения теплопроводности. Приведен сравнительный анализ скорости расчетов на центральном (CPU) и графическом (GPU) процессорах.

2. Описание схемы Дугласа-Рекфорда

При математическом моделировании распространения тепла с учетом фильтрации и фазовых превращений применяется следующее уравнение теплопроводности:

(1)   \begin{equation*}  C_{ev}\frac{\partial T}{\partial t}+\bigtriangledown{(-k\bigtriangledown{T})}+C_w\vec{v}\bigtriangledown{T}-Q=0 \end{equation*}

Физический смысл коэффициентов, участвующих в уравнении (1), приведен в таблице 1.

Таблица 1. Коэффициенты уравнения теплопроводности

Коэффициенты Описание
C_{ev}=C-\rho_{ice}L{\partial w_{ice}}/{\partial T} эффективная теплоемкость грунта
w_{ice}=1-1/(exp(2A_{w}(T_{p}-T))+1) количество замерзшей воды
A_w эмпирический коэффициент
T_p температура фазового перехода
T температура
C=C_{m}w_{ice}+C_{t}(1-w_{ice}) теплоемкость грунта
C_m теплоемкость мерзлого грунта
C_t теплоемкость талого грунта
t время
k=k_{m}w_{ice}+k_{t}(1-w_{ice}) эффективная теплопроводность грунта
k_m теплопроводность мерзлого грунта
k_t теплопроводность талого грунта
C_w теплоемкость воды
\vec{v} вектор скорости фильтрации
Q источник или сток тепла

Уравнение (1) является дифференциальным уравнением с частными производными параболического типа. В силу произвольности расчетной области, в которой заданно уравнение (1), при его решении используют численные методы. Одним из классов численных методов решения параболических дифференциальных уравнений является, так называемые, схемы переменных направлений [1]. Широкое распространение получила схема Дугласа-Рекфорда [1, 2]. Ниже в заметке дано краткое описание упомянутой схемы.

Рассмотрим однородное трехмерное уравнение теплопроводности, которое является упрощенным вариантом уравнения (1):

(2)   \begin{equation*}  \frac{\partial u}{\partial t}=\frac{\partial^{2} u}{\partial x^{2}}+\frac{\partial^{2} u}{\partial y^{2}}+\frac{\partial^{2} u}{\partial z^{2}} \end{equation*}

Пусть расчетная область является параллелепипедом и задана пространственная сетка узлов \omega_{l}=\{(x_{i},y_{j},z_{k}):i=1,...,N_{x};j=1,...,N_{y};k=1,...,N_{z}; \}, а также, временная сетка \omega_{\tau}=\{t_{i}:i=1,...,n\}.

Пусть

H^{x}=\{h^{x}_{i}=x_{i+1}-x_{i}:i=1,...,N_{x}-1\},H^{y}=\{h^{y}_{j}=y_{j+1}-y_{j}:j=1,...,N_{y}-1\}, H^{z}=\{h^{z}_{k}=z_{k+1}-z_{k}:k=1,...,N_{z}-1\}

– шаги расчетной сетки по пространству, а

H^{t}=\{\tau_{i}=t_{i+1}-t_{i}:i=1,...,n-1\}

– шаги по времени.

Предположим, что значение температурного поля в момент времени t_{s}\in{\omega_{\tau}} известно и в узле с координатами (x_{i},y_{j},z_{k}) температура равна y^{s}_{i,j,k}. Для нахождений температурного поля на следующем временном слое, следуя схеме Дугласа-Рекфорда, уравнение (2) аппроксимируется следующим образом:

(3)   \begin{equation*}  \frac{y^{s+1/3}_{ijk}-y^{s}_{ijk}}{\tau}=\Delta_{x}y^{s+1/3}_{ijk}+\Delta_{y}y^{s}_{ijk}+\Delta_{z}y^{s}_{ijk}, j=1,...,N_{x};(j,k) - fix; \end{equation*}

(4)   \begin{equation*}  \frac{y^{s+1}_{ijk}-y^{s+2/3}_{ijk}}{\tau}=\Delta_{z}y^{s+1}_{ijk}-\Delta_{z}y^{s}_{ijk},i=1,...,N_{y};(i,k) - fix; \end{equation*}

(5)   \begin{equation*}  \frac{y^{s+2/3}_{ijk}-y^{s+1/3}_{ijk}}{\tau}=\Delta_{y}y^{s+2/3}_{ijk}-\Delta_{y}y^{s}_{ijk},k=1,...,N_{z};(i,j) - fix; \end{equation*}

где

\Delta_{x}u_{ijk}=\frac{1}{h_{i}^{2}}u_{i+1jk}-\left(\frac{1}{h_{i}^{2}}+\frac{1}{h_{i-1}^{2}}\right)u_{i+1jk}+\frac{1}{h_{i-1}^{2}}u_{i-1jk},
\Delta_{y}u_{ijk}=\frac{1}{h_{j}^{2}}u_{ij+1k}-\left(\frac{1}{h_{j}^{2}}+\frac{1}{h_{j-1}^{2}}\right)u_{ij+1k}+\frac{1}{h_{j-1}^{2}}u_{ij-1k},
\Delta_{z}u_{ijk}=\frac{1}{h_{k}^{2}}u_{ijk+1}-\left(\frac{1}{h_{k}^{2}}+\frac{1}{h_{k-1}^{2}}\right)u_{ijk+1}+\frac{1}{h_{k-1}^{2}}u_{ijk-1}

– конечноразностные производные второго порядка вдоль оси Ox, Oy и Oz, соответственно. При каждой фиксированной паре индексов вдоль осей Oy и Oz (а таких пар ровно N_{y}\times{N_{z}}), N_{x} уравнений (3), собранных вместе, образуют трехдиагональную СЛАУ (систему линейных алгебраических уравнений). Для решения данной СЛАУ методом прогонки требуется O(N_x) операций. Аналогичные утверждения имеют место при решении разностных уравнений (4) и (5).

Так как при каждой фиксированной паре индексов (j,k) ((i,k), (i,j)) уравнения (3) ((4) и (5), соответственно) решаются независимо друг от друга, то говорят, что схема Дугласа - Рекфорда обладает естественным параллелизмом. Данный факт позволил эффективно распараллелить приведенную схему, используя технологию CUDA.

3. Тестовые расчеты

Были проведены расчеты, показывающие прирост производительности для программно-реализованной схемы Дугласа-Рекфорда с использованием технологии CUDA при сравнении с реализацией данной схемы на центральном процессоре компьютера.

Рассмотрим расчетную область P = {(x,y,z):0\leqslant{x}\leqslant{10}, 0\leqslant{y}\leqslant{10}, 0\leqslant{z}\leqslant{10}} (см. рисунок 1), которая имеет форму куба со стороной 10 метров. В области P параллельно оси Oz размещено СОУ (сезонное охлаждающие устройство) S. Нижняя точка A СОУ S имеет координаты (5.0, 5.0, 1.0), а верхняя точка B СОУ S имеет координаты (5.0, 5.0, 9.0).

СОУ S обладает следующими конструктивными параметрами: радиус трубы испарителя равен 16,85 mm, площадь испарителя равна 1.0 m2, площадь конденсатора равна 2.0 m2.

Расчетная область P заполнена материалом M, у которого объемная теплоемкость в талом и мерзлом состояниях равна 2.14\times{10^6}\frac{J}{m^{3}\times{K}}.

Теплопроводность материала M в талом и мерзлом состояниях равна 1.8\frac{W}{m\times{K}}.

Расчетная область

Рисунок 1 – Расчетная область

На границе расчетной области P задано граничное условие второго рода (нулевой поток), на СОУ S задано граничное условие третьего рода, где температура среды равна -20.0 °C, а коэффициент теплообмена равен 10.0\frac{W}{m\times{K}}. Начальная температура материала, которым заполнена расчетная область, равна -1.0 °C.

В приведенном тесте рассчитывается распределение теплового поля через каждые 30 суток на протяжении 360 дней.

Были проведены четыре пары расчетов при разной степени дискретизации расчетной области – первый расчет в каждой паре осуществлялся на одном ядре центрального процессора, а второй – на видеокарте с использованием технологии CUDA.

Для первой пары расчетов был использован шаг по пространству во всех трех направлениях равный 100 mm. Для второй пары расчетов шаг по пространству вдоль осей Ox и Oy равен 100 mm, а вдоль оси Oz – 50 mm. Для третей пары расчетов шаг по пространству вдоль оси Ox равен 100 mm, а вдоль осей Oy и Oz – 50 mm. И для последней пары расчетов был использован шаг по пространству во всех трех направлениях равный 50 mm. Таким образом, в первом случае область P была разбита на 100\times{100}\times{100}=1 000 000 ячеек. Во втором, третьем и четвертом случаях область P разбивалась на 2 000 000, 4 000 000 и 8 000 000 ячеек соответственно.

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

Тепловое распределение около термостабилизатора в цветовой заливке Изолинии вокруг термостабилизатора

Рисунок 2а – тепловое распределение в виде цветовой заливки

Рисунок 2б – тепловое распределение в виде изолиний

Расчеты проводились на процессорах, приведенных в таблице 2.

Таблица 2. Расчетные процессоры

Процессор (CPU) Inlel Core i7 4 ядра
Видеокарта (GPU) GeForce GTX TITAN 2688 ядер

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

Таблица 3. Затраченное время

Количество ячеек Время расчета на Inlel Core i7 (округленно до секунды) Время расчета на GeForce GTX TITAN (округленно до секунды) Ускорение в GPU
(округленно до сотых)
1 000 000 1962 63 31.14
2 000 000 16 296 524 31.10
4 000 000 30 654 1013 30.26
8 000 000 63 858 2079 30.72

Разница скорости вычислений CPU и GPU

Рисунок 3 – Ускорение расчета

4. Вывод

В среднем расчеты на GPU версии ADI решателя с использованием технологии CUDA на видеокарте Nvidia GTX Titan производились в 30 раз быстрее, чем расчеты, осуществляемые на одном ядре центрального процессора Intel Core i7. Отметим, что при проведении тестов решалась относительно несложная тестовая задача (небольшое количество стационарных граничных условий, небольшое количество материалов). При решении более сложных вычислительных задач необходимо ожидать снижение прироста скорости при использовании параллельных вычислений.

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

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

1.Яненко, Н. Н. Метод дробных шагов решения многомерных задач математической физики / Н.Н. Яненко. – Новосибирск: НАУКА, 1967. – 197 с.

2.Thomas, J.W. Numerical partial differential equations: finite difference methods / J.W. Thomas. – USA: Springer 1995. – 437 p.

Реализация метода переменных направлений с использованием технологии CUDA: 2 комментария

  1. Почему сравниваете 1 ядро i7? По факту если взять шестиядерный процессор с гипертрейдингом будет уже 12 потоков, и разница составит только 2-3 раза, а не 30 раз?????

  2. Мы выполняем сравнение с одним потоком/ядром Core i7 для того, чтобы можно было просто определить прирост скорости при сравнении с двумя и более ядрами/потоками. Если изначально сравнивать с шестью, то сложно соотнести прирост при переходе на n-ядерный/потоковый процессор. Ведь в Вашем случае, Вы без труда получили соответствующую калькуляцию. Правда, необходимо ещё ее немного уточнить с поправкой на каналы доступа к памяти. Так, при работе с одним ядром Core i7 и полной разгрузкой процессора для других задач, ядру предоставляется вся шина данных. Это существенно ускоряет процесс расчёта, так как процесс расчета - итерационный, т.е. необходимо постоянно загружать и выгружать данные из оперативной памяти. В случае 6-и ядер с гипертрейдингом процессору не будет предоставлено отдельно 12 полных каналов данных. Поэтому ускорение будет уже не в 12 раз, а хорошо, если в 6-8.

Добавить комментарий

Ваш e-mail не будет опубликован. Обязательные поля помечены *


*

Можно использовать следующие HTML-теги и атрибуты: <a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <strike> <strong>