Уравнение Лапласа/Пуассона — одно из фундаментальных уравнений математической физики, описывающее стационарные состояния(не зависящие от времени) в различных физических процессах, таких как теплопроводность, электростатика, гидродинамика и теория упругости.
Численные методы его решения подразделяются на итерационные и прямые. В итерационных методах решаются системы линейных алгебраических уравнений(СЛАУ) за некоторое количество итераций до достижения заданной точности. Прямые методы находят решения за одну итерацию, но они не универсальны, в том смысле, что неприменимы для непрямоугольных областей или в случае переменных коэффициентов.

Программа визуализирует решение и промежуточные результаты задачи Дирихле для уравнения Лапласа/Пуассона в прямоугольнике на CPU и GPU.

Уравнение Пуассона(при f=0 уравнение Лапласа) :

\begin{gather}\frac{\partial ^{2}}{\partial x}u + \frac{\partial ^{2}}{\partial y}u = -f\end{gather}

Задача Дирихле(краевые условия 1-го рода) - заданы значения u на границе области(в данном случае в прямоугольнике).

Результаты визуализируются при помощи тепловой карты(HeatMap). Используются цвета и оттенки синего, сине-зелёного(циан), зеленого, желтого и красного цветов(от меньшего к большему значению). Нажав на кнопку '2D -> 3D' можно получить 3D изображение вычисленного решения. В 3D режиме перемещения мыши при нажатой левой кнопке приводят к изменению проекции 3D изображения. А если ещё и зажать Shift, то будет смещаться всё изображение.

Программа написана на C#(WPF). В OpenCL используется Cloo, в CUDA - ManagedCUDA. Точность всех методов O(h^{2}).

Итерационные методы

Простой итерции(Якоби).
Скользящей итерации(Гаусса-Зейделя).
Верхней релаксации.
Расщепления.
Переменных направлений(ADI).
Попеременно-треугольный.
Скорейшего спуска.
Минимальных невязок.
Сопряжённых градиентов
Бисопряжённых градиентов
3-х слойный Чебышевский
Многосеточный (MultiGrid)

Прямые методы

Циклической(полной) редукции(CR)(2 варианта)
Разделения переменных(FA)(4 варианта, 2 основных - FFT и прогонка, 2 дополнительных - редукция и маршевый скалярный)
Неполной редукции(FACR)
Матричной прогонки
Маршевый(GMA)(параметрический)

Особенности методов

  • Методы простой и скользящей итерции простые, но медленные, сами по себе не используются. Зато они используются в качестве сглаживателей в многосеточном методе. Без их использования многосеточный метод не сходится.
    Для простой итерации есть возможность использовать Чебышевский набор параметров - ускорение в 10 раз.

  • Методы расщепления и переменных направлений - счёт на установление с оптимальным итерационным параметром.
    Для переменных направлений можно использовать ускорение Жордана(придумал Вашпресс(Wachspress)).

  • Для попеременно-треугольного метода также можно использовать Чебышевский набор параметров.

  • Многосеточный метод интересен тем, что уже 1-я итерация визуально даёт картинку, близкую к окончательному решению.

Быстрые методы:
Маршевый, Разделения переменных, Неполной редукции, Циклической редукции, Многосеточный, Переменых направлений(с Жордановым ускорением).

Медленные методы:
Простой итерации, Скользящей итерации, Скорейшего спуска, Минимальных невязок.

Длинная арифметика

Числа битностей 32(float), 64(double), 128(DD128), 256(QD256) поддерживаются благодаря интерфейсу INumber в C#. Для битностей 128 и 256 за основу был взят NetQD, который в свою очередь основан на https://github.com/aoki-t/QD, http://web.mit.edu/tabbott/Public/quaddouble-debian/qd-2.3.4-old/docs/qd.pdf. Добавил туда поддержку интерфейса INumber и сопутствующие ему(ITrigonometricFunctions, ILogarithmicFunctions и т.д.). Сейчас DD128, QD256 лежат на ГитХабе отдельным проектом https://github.com/NickUtenkov/DD128_QD256_INumber

Опции программы

  1. Для итерационных методов можно выбрать один из 4-х способов интерполяции начального приближения :
    Среднее (min+max)/2
    Средне-арифметическая
    Кусочно-линейная(средне-арифметическая)
    Кусочно-линейная(весовая)

  2. Количество итерациий, можно также задавать в файле входных данных.

  3. Визуализация - позволяет отображать промежуточные результаты решения. Если опция выключена, то отображаются только начальное приближение и результат.

  4. Шаг визуализации - можно менять количество промежуточных изображений

  5. Сохранять анимированный GIF - в OutputData папку сохраняются анимированные GIF'ы(визуализация лучше чтобы была включена).

  6. Многопоточность - отключение опции позволит узнать насколько медленнее работает метод без распараллеливания.

Расход памяти

Память измеряется в КИПах(коэффициент использования памяти). Рабочий массив размера MxN(или NxN) имеет КИП равный 1. Причём массив может быть двумерным, одномерным(адресуемый как двумерный) или ступенчатым. Размеры массивов, которые намного меньше рабочего - не учитаваем.

Значения КИП для уравнения Лапласа(для Пуассона, как правило, на единицу больше) :
Итерационные методы
Простой итерции - 2(текущая и следующая итерация)
Скользящей итерации и верхней релаксации - 1(используется красно-чёрное чередование)
Расщепления - 3
Переменных направлений - 3
Попеременно-треугольный - 4
Скорейшего спуска - 3
Минимальных невязок - 3
Сопряжённых градиентов - 4
Бисопряжённых градиентов - 6(для стабилизированного - 8)
3-х слойный Чебышевский - 4
Многосеточный примерно(из-за нюансов) равен (1 + 1/4 + 1/16 + 1/32 + ...) * 3

Прямые методы
Циклической редукции - 1-1.5
Разделения переменных - 1
Неполной редукции - 1
Матричной прогонки - 1
Маршевый - 1.75

Прямые методы и быстрее и памяти потребляют меньше.

Описание папок и формата файлов

<папка пользователя>/VLP2D_Data(Перед 1-м запуском программы перенесите папку VLP2D_Data в <Каталог пользователя>)
Подкаталоги :
InputData - содержит файлы данных. Можно создавать свои файлы данных.
Ini - хранятся параметры вычислений той или иной задачи.
CUDA, OCL - содержат откомпиллированные ядра GPU. Их можно удалять, вносить изменения в исходный текст ядер и они создадутся заново при следующем запуске программы.
OutputData - в эту папку сохраняются анимированные GIF'ы, при включенной опции.

Описание параметров файлов данных

xMin - левая граница по кординате X
yMin - нижяя граница по кординате Y
xMax - правая граница по кординате X
yMax - верхняя граница по кординате Y
segmentsX - число отрезков разбиения по X(число узлов по X равно (segmentsX + 1))
segmentsY - число отрезков разбиения по Y(число узлов по Y равно (segmentsY + 1))
funcRHS - выражение для функции правой части(для уравнения Лапласа равна нулю)(RHS - Right hand side)
funcLeft - выражение для функции на левой границе
funcRight - выражение для функции на правой границе
funcTop - выражение для функции на верхней границе
funcBottom - выражение для функции на нижней границе
funcBoundary - выражение для функции на всей границе
funcAnalitic - выражение для функции решения уравнения(если известна), можно использовать для сравнения с результатом программы.
epsilon - точность решения(для итерационных методов)
epsilon можно задавать общий для всех разрядностей или индивидуально epsilonSingle, epsilonDouble, epsilonDD128, epsilonQD256.
Для функции на границе надо задавать либо funcBoundary, либо funcLeft, funcRight, funcBottom, funcTop.
Если задана funcAnalitic, то дополнительно отображается HeatMap отклонений аналитического решения от вычисленного.

xMin, yMin равны нулю, если не указаны
xMax, yMax равны 1, если не указаны
segmentsX, segmentsY равны 100, если не указаны
xMin, yMin, xMax, yMax можно задавать в виде функции от константы.
Например, в выражении xMax=2*asin(1)(или xMax=pi(1), о pi(x) см. ниже) будет вычислено число Пи(с текущей разрядностью чисел).

Следующие параметры задают число итераций для 4-х вариантов вычисления начального приближения, битности и GPU.
SimpleIteration - для простой итерации
SlidingIteration - для скользящей итерации
SOR - для верхней релаксации
Splitting - для расщепления
VarDir - для переменых направлений
PTM - для попеременно-треугольного
GradientDescent - для скорейшего спуска
MinimumResidual - для минимальных невязок
ConjugateGradient - для сопряженных градиентов
BiconjugateGradient - для бисопряженных градиентов
Chebishev3Layers - для 3-слойной Чебышевской
MultiGrid - для многосеточного

Примеры :
SimpleIterationQD256=1,2,3,4 - CPU, QD256 и 4 значения количества итераций
SimpleIterationOCLDD128 - OpenCL, DD128
SimpleIterationCUDADouble - CUDA, Double

Число итераций нужно в том числе для точного отображения ProgressBar. Если точный ProgressBar не важен, то можно задавать большие числа, например 1,000,000 . Для Чебышевского набора параметров и Жорданова ускорения число итераций известно заранее. Лучше задавать количество итераций на 1 больше, тогда если окончательное количество на 1 меньше, то итерации сошлись.

Для использования выражений для функций(и значений xMin, yMin,xMax, yMax) используется глубоко модернизированная библиотека ELW(на Хабре была статья о ней).

Примеры файлы данных в папке VLP2D_DataInputData
x^3 - y^2 + 10 x y + 10 x - 10 y + 1 - функция 2-х переменных, x^3 это x в степени 3
x^0.5 нужно использовать вместо sqrt(x)

Помимо арифмтических операторов(+, -, *, /, унарный плюс, унарный минус) поддерживаются ещё два :
^ - возведение в степень
Тернарный оператор ?:

Поддерживаемые математические функции

sin - синус
cos - косинус
tan - тангенс
asin - арксинус
acos - арккосинус
atan - арктангенс
sh - синус гиперболический
ch - косинус гиперболический
th - тангенс гиперболический
arsh - арксинус гиперболический
arch - арккосинус гиперболический
arth - арктангенс гиперболический
exp - экспонента
log(x,y) - логарифм x по основанию y
ln - логарифм натуральный
lg - логарифм десятичный
pi(x) - число Пи умноженное на x(хотя можно использовать 2*asin(1)*x )

Свои функции можно добавлять в Utils.addCustomFunctions(см. реализацию функции pi(x)).

Сходимость альфы

Для решения трёхдиагональных СЛАУ используется метод прогонки. Во всех методах, кроме маршевого, матрица системы задается одним числом(диагональ, d). Над и поддиагонали равны 1(или -1). В этих случаях имеет место сходимость коэффициента альфа при |d| > 2. В некоторых случаях альфа сходится уже на 11-м номере(для double), ~22(для DD128), ~42(для QD256).

FFT

В прямых методах FA, FACR, GMA для вычисления сумм синусов используется БПФ посредством библиотеки FFTW(а для битностей 128 и 256 используется Lomont на C#).
Можно использовать и интеловский MKL FFT, но его размер в 100 раз больше(6 DLL общим размером в 266 Мб против 2.6 Мб у FFTW).
В CUDA используется встроенный FFT, в OpenCL - АМДшный.

Реализация длинной арифметики на GPU

Поскольку в CUDA есть компилятор C++, использовать длинную арифметику в ней проще, чем в OpenCL. Для этого надо добавить операторы +, -, *, / и необходимые функции. В OpenCL вместо использования C++(есть там всякие SPIR, SYCL, cl_ext_cxx_for_opencl) было решено использовать трюк с макросом. Выражение a = b + c(все имеют тип, скажем, DD128) записываются как a = HP(b + c). HP(High Precision) - это макрос, после обработки которого арифметические операторы заменяются на вызовы соответствующих функций. a = b + c становится a = add_HH(b, c), функция add_HH - сложение чисел типа DD128 или QD256.
Интеловские процессоры(даже не сильно старые, типа i5-13500H) поддерживают только тип float в OpenCL, CUDA поддерживает double(хотя и медленнее, бывает, что медленнее в 64 раза).

Для отладки в OpenCL помогает OCLgrind Simulator. И можно использовать такую вещь как https://github.com/intel/opencl-intercept-layer. Если поставить древний AMD-APP-SDKInstaller, то процессор(и интеловский тоже) будет отображатся как устройство OpenCL.

Внешние зависимости(не идут вместе с проектом)

OpenCL - OpenCL.dll(обычно находится в C:\Windows\SysWOW64 и/или C:\Windows\System32)
CUDA - названия файлов в SDK 12.5 cufft64_11.dll, nvrtc-builtins64_125.dll, nvrtc64_120_0.dll, находятся в C:\Program Files\NVIDIA GPU Computing Toolkit\CUDA\v12.5\bin, версия SDK может быть другой.

Ссылки (в исходных текстах программы) на книги и статьи

[SNR] Samarskii-Nikolaev russian - А. А. Самарский, Е. С. Николаев Методы решения сеточных уравнений.
[SNE] Samarskii-Nikolaev english - Aleksandr A. Samarskii Evgenii S. Nikolaev Numerical Methods for Grid Equations Volume I Direct Methods.
[S_VVCM] Самарский А. А. Введение в численные методы.
[IL] Ильин В.П., Кузнецов Ю.И. Трехдиагональные матрицы и их приложения.

Заключение

Файлы проекта находятся по адресу. Готовый исполняемый файл находится в каталоге VisualLaplacePoisson2D\bin\x64\Release\net9.0-windows. Некоторые анимированые Gif находятся в каталоге ResultPictures.