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

                 eISSN 2311-2468

EDN:                                 

Национальный исследовательский Мордовский государственный университет,

Уравнение Пуассона относится к числу фундаментальных уравнений математической физики и широко применяется для описания потенциалов, температурных полей, распределения давления и других физических величин. Особое значение оно приобретает при моделировании процессов в задачах с цилиндрической геометрией [1; 2]. При этом использование декартовых сеток для цилиндрических областей приводит к аппроксимации границ «лесенкой», что снижает точность, поэтому переход к цилиндрической системе координат ( r ,ф, z ) является естественным выбором.

Для тестирования разрабатываемых алгоритмов была решена следующая задача

Дирихле:

a ² u

+ —= - 3cos 0, a z²

1 а Г a u ^ 1 a2 u r+ r ar [ar) r2 ap2

< u ( 0 ,0,z ) = u ( 1 , 0 ,z ) = z, u ( r,e, 0 ) = r ( 1 - r ) cos 0, u ( r,0, 1 ) = 1 + r ( 1 - r ) cos 0.

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

С развитием суперкомпьютерных технологий акцент в решении вычислительно емких задач сместился на реализацию параллельных алгоритмов для масштабируемости вычислений1. В частности, в силу доступности и эффективности все чаще для решения таких задач используют графические процессоры (Graphics Processing Unit, GPU)2. Сложность разработки вычислительных алгоритмов для высокопроизводительной реализации с использованием GPU состоит в необходимости адаптации метода к архитектуре GPU и минимизации дивергенции потоков при обработке граничных условий.

Цель исследования – разработать высокопроизводительную реализацию геометрического многосеточного метода для цилиндрических координат на базе технологии единой архитектуры вычислительных устройств ( Compute Unified Device Architecture, CUDA ) и провести анализ ее эффективности.

ОБЗОР ЛИТЕРАТУРЫ

В современной литературе одним из эффективных подходов признан геометрический многосеточный метод ( Geometric Multigrid, GMG ), обеспечивающий оптимальную сложность 0 ( N ). Исследования A. Брандта³ и У. Троттенберга, К. Остерлее и A. Шуллера⁴, заложили теоретический фундамент многосеточных методов. Классические труды, такие как фундаментальное руководство по многосеточным методам⁵ подробно описывают базовые принципы их построения.

Современные исследования направлены на адаптацию многосеточных подходов к сложным численным схемам и параллельным архитектурам. В частности, в работах Р. В. Жалнина, М. С. Нефедова и С. Х. Зининой6 активно исследуется применение многосеточных алгоритмов в сочетании с разрывным методом Галёркина для решения уравнений конвекции-диффузии и нелинейных задач [3; 4]. Авторами показано, что классический h-multigrid демонстрирует более высокую эффективность в подавлении низкочастотных погрешностей по сравнению с p-вариантом, а использование стратегий Ньютона-Крылова на сеточных уровнях позволяет существенно сократить вычислительные затраты при сохранении высокого порядка точности.

Вопросы реализации многосеточных методов на параллельных системах и их интеграции в комплексные алгоритмы рассматриваются в работах исследователей [5; 6]. Продемонстрирована высокая эффективность многосеточных технологий при решении уравнений Навье-Стокса на сотнях процессоров. Параллельно с этим развиваются подходы, использующие мощность графических процессоров (GPU) для ускорения расчетов в областях со сложной конфигурацией границ [7].

Особое место занимают исследования, посвященные моделированию тепломас-сопереноса и химических превращений. В работах Е. В. Мортикова, Е. Е. Писковой, В. Н. Снытникова рассматриваются явные и итерационные методы решения уравнений параболического типа, а также численные модели воздействия лазерного излучения на многокомпонентные реакционные среды [8; 9]. Отмечается, что при моделировании таких процессов, выбор эффективного метода решения СЛАУ является критически важным для достижения приемлемого времени расчета.

МАТЕРИАЛЫ И МЕТОДЫ

В основе многосеточного метода лежит идея об иерархии вложенных сеток, пример вложенной сетки в цилиндрических координатах представлен на рис. 1.

Р и с . 1 . Схематичное представление процесса огрубления сетки: а – подробная сетка Ω_h ; b – грубая сетка Ω_2h F i g . 1 . Schematic representation of the grid coarsening process: a – fine grid Ω_h ; b – coarse grid Ω_2h

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

Source s : hereinafter, figures were prepared by the author based on the results of the study.

• ⁶    Нефедов, М.С., Жалнин Р.В., Зинина С.Х. Исследование многосеточного метода для решения уравнений в частных производных разрывным методом Галёркина. Дифференциальные уравнения и их приложения в математическом моделировании: Сборник материалов XVII Международной научной конференции, Саранск, 29–31 июля 2025 года. – Саранск: Средневолжское математическое общество. 2025:201–205

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

На каждом уровне сетки алгоритм последовательно выполняет несколько шагов. Сначала выполняется итераций сглаживателя на сетке уровня l для устранения высокочастотной ошибки. Затем вычисляется невязка r_l = f_l - a_lu_l . Она переносится на более грубую сетку l - 1 с помощью оператора ограничения I_l ^{l - 1} . На самом грубом уровне система уравнений решается либо прямым алгоритмом, либо большим числом итераций до полной сходимости. После этого найденная поправка e_{l - 1} возвращается обратно на подробную сетку l с помощью интерполяции и решение обновляется с учетом найденной поправки u_l = u_l + e_l . В конце выполняется v ₂ итераций постсглаживания, чтобы убрать новые возмущения, которые могли появиться в процессе интерполяции.

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

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

Во-вторых, W-цикл – двойной рекурсивный вызов метода на грубых уровнях. Достоинства: более надежная сходимость по сравнению с V-циклом, так как ошибка на грубых сетках решается точнее. Недостатки: существенно возрастает количество операций на грубых уровнях. Для GPU это критично, поскольку грубые сетки не обеспечивают достаточной загрузки тысяч ядер видеокарты, что приводит к простою вычислительных мощностей.

В-третьих, F-цикл – гибрид V-цикла и полного многосеточного подхода, где после спуска к грубой сетке и пролонгации коррекции на промежуточных уровнях запускается дополнительный V-цикл для уточнения. Достоинства: обеспечивает лучшую сходимость по сравнению с V-циклом за счет дополнительного уточнения на уровнях; подходит для задач с неравномерными сетками. Недостатки: сложнее в реализации из-за вложенных циклов; на GPU снижается эффективность на очень грубых уровнях из-за низкой загрузки ядер.

Р и с . 2 . Визуализация вариантов обхода сеток: a – V-цикл; b – W-цикл; с – F-цикл F i g . 2 . Visualization of multigrid cycle types: a – V-cycle; b – W-cycle; с – F-cycle

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

U ^{( n )          ( n )             n )}     _,,^{( n )             k )          ( n )}         ,,^{( n )          ( n )}

u i +1 ,j,k ⁺ u i -1 ,j,k , ui+1 ,j,k u i -1 ,j,k ₊ u i,j +1 ,k ⁺ u i,j-1 ,k , u i,j,k +1 ⁺ u i,j,k -1 _ f

u

( n +1 ) i,j,k

h r             2 rtf            ri' %             h i

2	2	2
—	+ _ _	+
h r 2	r i h ,	h z 2

Многосеточные алгоритмы (V- и W-циклы) были реализованы с использованием библиотеки PyCUDA7 для организации параллельных вычислений на графическом ускорителе на базе технологий CUDA8. Использование PyCUDA позволило реализовать сложные рекурсивные структуры GMG, сохраняя при этом гибкость Python.

Все расчеты проводились на вычислительном стенде, оснащенном процессором Intel(R) Core(TM) Ultra 7 155H с тактовой частотой до 3.80 ГГц, 32 ГБ оперативной памяти и графическим ускорителем NVIDIA GeForce RTX 4060 Mobile. Программная среда базировалась на Python 3.12.7, PyCUDA 2025.1 и CUDA Toolkit 12.8.

РЕЗУЛЬТАТЫ ИССЛЕДОВАНИЯ

Для оценки эффективности проводились расчеты на последовательности сгущающихся сеток. В качестве метрики эффективности выбрано время достижения сходимости до относительной невязки £ < 10 5 .

На рисунке 3 приведена зависимость времени выполнения четырех численных методов: метода Якоби и трех вариантов многосеточного метода (V-, W-, F-циклов) от количества расчетных ячеек сетки.

Метод Якоби имеет наибольшее время выполнения, оно резко возрастает с увеличением размерности сетки, что типично для итерационных методов решения систем уравнений и делает его неприменимым для задач с большим числом ячеек. Многосеточные циклы демонстрируют различную эффективность. V- и F-циклы обеспечивают лучшую сходимость и существенно меньшее расчетное время на сетках большой размерности по сравнению с W-циклом. Хотя V- и F-циклы показывают схожие результаты на сетках средней размерности, именно кривая времени выполнения F-цикла имеет наименьший наклон, что указывает на самую низкую вычислительную сложность при увеличении размерности и открывает возможность моделирования реальных задач на детальных сетках

ОБСУЖДЕНИЕ

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

• ⁷    Documentation for the pyCUDA library [Online]: сайт. URL: https://documen.tician.de/pycuda/ (accessed: 20.11.2025)

• ⁸    Информация о технологии CUDA [Электронный ресурс] : сайт. URL: https://developer.nvidia.com/cuda-toolkit/ (дата обращения: 20.11.2025)

Technical sciences 235

Е

CD С

Е о О

ZE Ф 5

IT ex

2 ф

CQ

64 000                 128 000                256 000                512 000 Размер сетки /

Grid size

V-цикл / V-cycle              F-цикл / F-cycle                Якоби / Jacobi                W-цикл / W-cycle

Рис. 3 . Сравнение численных методов решения

• Fig. 3. Comparison of numerical methods

W-цикл, несмотря на теоретически лучшую сходимость, на практике проигрывает V- и F-циклам по времени. Это связано с аппаратными особенностями GPU. W-цикл постоянно спускается на грубые сетки, где количество узлов слишком мало. Из-за этого большинство ядер GPU простаивают, а большая часть времени тратится на накладные расходы API при запуске ядер.

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

ЗАКЛЮЧЕНИЕ

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

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