Уравнение Пуассона, решение дискретным преобразованием Фурье

Основные авторы описания: В.М.Степаненко, Е.В.Мортиков.

1 Свойства и структура алгоритма

1.1 Общее описание алгоритма

Уравнение Пуассона для многомерного пространства имеет вид:

[math] \sum_{i=1}^{N}\frac{\partial^2 \phi}{\partial x_i^2}=f,~\mathbf{x}\in D, [/math]

где [math]D \in \mathbb{R}^N[/math] - область определения решения [math]\phi(\mathbf{x})[/math], [math]\mathbf{x}=(x_1,...,x_N)^T[/math] - вектор независимых переменных. Уравнение Пуассона дополняется граничными условиями:

[math] B(\phi)=F, \mathbf{x} \in \Gamma(D), [/math]

где [math]\Gamma(D)[/math] - граница области [math]D[/math], a [math]B(\phi)[/math] - оператор, определяющий граничные условия. [math]B(\phi)=\phi[/math] соответствует граничным условиям Дирихле, [math]B(\phi)=\partial\phi/\partial n[/math] ([math]\mathbf{n}[/math] - внешняя нормаль к границе [math]\Gamma(D)[/math]) - условиям Неймана. Иногда задают также смешанное граничное условие: [math]B(\phi)=C\phi+\partial\phi/\partial n[/math] ([math]C[/math] - константа). Встречаются также так называемые "периодические граничные условия", при которых задача ставится для бесконечной области, но предполагается периодичность решения по подмножеству переменных из [math]\mathbf{x}[/math].

Уравнение Пуассона возникает во многих задачах математической физики, например, в электростатике (в этом случае [math]\phi[/math] - потенциал электрической силы) и гидродинамике ([math]\phi[/math] - давление жидкости или газа); при этом [math]N=2,3[/math] для плоской и трехмерной задач, соответственно.

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

[math] \sum_{i=1}^{N}\frac{\phi_{k_1,...,k_i+1,...,k_N}-2\phi_{k_1,...,k_i,...,k_N}+\phi_{k_1,...,k_i-1,...,k_N}}{\Delta x_i^2}=f_{k_1,...,k_N},~(k_1,...,k_N) \in D_N. [/math]

Здесь вторые производные заменены конечно-разностной аппроксимацией второго порядка точности (образуя т.н. схему "крест" для плоской задачи), а решение ищется на дискретном множестве точек [math]N[/math]-мерного пространства, [math]D_N[/math]. Конечными разностями аппроксимируются при этом также граничные условия.

1.2 Математическое описание алгоритма

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

[math] \frac{\phi_{i+1,j,k}-2\phi_{i,j,k}+\phi_{i-1,j,k}}{ \Delta x^2}\,+\, \frac{\phi_{i,j+1,k}-2\phi_{i,j,k}+\phi_{i,j-1,k}}{ \Delta y^2}\,+\, \frac{\phi_{i,j,k+1}-2\phi_{i,j,k}+\phi_{i,j,k-1}}{ \Delta z^2} = f_{i,j,k},~(i,j,k) \in D_h, [/math]

где [math]D_h=\{0:N_x-1\}\times\{0:N_y-1\}\times\{0:N_z-1\} [/math] - параллелепипед в сеточной области. В качестве граничных условий для простоты примем так называемые троякопериодические граничные условия

[math] \begin{align} \phi_{0,j,k} &= \phi_{N_x,j,k},\\ \phi_{i,0,k} &= \phi_{i,N_y,k},\\ \phi_{i,j,0} &= \phi_{i,j,N_z}. \end{align} [/math]

Периодические граничные условия удовлетворяются "автоматически", если представить решение в виде стандартного дискретного обратного преобразования Фурье:

[math] \phi_{i,j,k}=\frac{1}{N_x N_y N_z}\sum_{l=0}^{N_x-1}\sum_{m=0}^{N_y-1}\sum_{n=0}^{N_z-1}\Phi_{l,m,n}e^{2\pi \overline{i} \left(\frac{il}{N_x}+\frac{jm}{N_y}+\frac{kn}{N_z}\right)}. [/math]

Здесь [math]\overline{i}=\sqrt{-1}[/math]. Аналогичное разложение применяется также к правой части уравнения, [math]f_{i,j,k}[/math]. Решение дискретного уравнения Пуассона преобразованием Фурье удобно тем, что базисные функции разложения Фурье являются собственными функциями дискретного оператора Лапласа. А именно, подставляя разложение Фурье неизвестной функции [math]\phi_{i,j,k}[/math] и правой части в исходное уравнение, получаем:

[math] \Phi_{l,m,n}=-\frac{F_{l,m,n}}{4\left[\sin^2\left(\frac{\pi l}{N_x}\right) + \sin^2\left(\frac{\pi m}{N_y}\right) + \sin^2\left(\frac{\pi n}{N_z}\right) \right]}, [/math]

где [math]F_{l,m,n}[/math] - преобразование Фурье правой части. Отсюда очевиден алгоритм решения уравнения: сначала правая часть уравнения разлагается в ряд Фурье, затем по приведенной выше формуле вычисляются коэффициенты Фурье решения и, наконец, решение восстанавливается обратным преобразованием Фурье.

1.3 Вычислительное ядро алгоритма

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

[math] \phi_{i,j,k}=\frac{1}{N_x} \sum_{l=0}^{N_x-1} \left[ e^{2\pi \overline{i}\left(\frac{il}{N_x}\right)} \frac{1}{N_y} \sum_{m=0}^{N_y-1} \left[ e^{2\pi \overline{i}\left(\frac{jm}{N_y}\right)} \frac{1}{N_z} \sum_{n=0}^{N_z-1} \Phi_{l,m,n}e^{2\pi \overline{i} \left(\frac{kn}{N_z}\right)}\right]\right], [/math]

откуда видно, что трехмерное преобразование Фурье сводится в последовательности трех одномерных. Для вычисления одномерного преобразования широко используется эффективный алгоритм быстрого преобразования Фурье (БПФ, FFT, Fast Fourier Transform).

1.4 Макроструктура алгоритма

Из приведенного выше ясно, что основной макрооперацией алгоритма решения уравнения Пуассона является одномерное быстрое преобразование Фурье. В дальнейшем будем обозначать ее [math]\text{FFT}_i,~i=x,y,z[/math] по каждому из трех направлений, соответственно, а обратное преобразование Фурье - [math]\text{IFFT}_i,~i=x,y,z[/math].

1.5 Схема реализации последовательного алгоритма

Введем для краткости следующие обозначения сеточных функций: [math]\phi_h=\{\phi_{i,j,k},~(i,j,k) \in D_h\},~F_h=\{F_{l,m,n},~l=0,...,N_x-1,~m=0,...,N_y-1,~n=0,...,N_z-1\}[/math] и аналогично [math]f_h[/math] и [math]\Phi_h[/math]. Тогда алгоритм запишется следующим образом:

1. Вычисление [math]\text{FFT}_x(f_h)[/math]
2. Вычисление [math]\text{FFT}_y(\text{FFT}_x(f_h))[/math]
3. Вычисление [math]\text{FFT}_z(\text{FFT}_y(\text{FFT}_x(f_h)))=F_h[/math]
4. Расчет [math]\Phi_h[/math] по [math]F_h[/math]
5. Вычисление [math]\text{IFFT}_x(\Phi_h)[/math]
6. Вычисление [math]\text{IFFT}_y(\text{IFFT}_x(\Phi_h))[/math]
7. Вычисление [math]\text{IFFT}_z(\text{IFFT}_y(\text{IFFT}_x(\Phi_h)))=\phi_h[/math]

Примечательно, что весь описанный алгоритм реализуется с помощью одного трехмерного массива [math]N_x \times N_y \times N_z[/math], поскольку результат каждого одномерного преобразования Фурье можно записать во входной массив и затем использовать его на входе для следующего преобразования, а операция (4) - это поэлементное изменение массива.

1.6 Последовательная сложность алгоритма

Рассмотрим случай куба, [math]N_x=N_y=N_z=N=2^k[/math]. Тогда сложность одномерного БПФ по каждому из направлений составляет (по порядку) [math]N(\text{log}_2 N)[/math] операций. Поскольку одномерное БПФ на каждом из этапов алгоритма, кроме четвертого, осуществляется [math]N^2[/math] раз, а количество операций на этапе 4 составляет (снова по порядку) [math]N^3[/math], то общее количество арифметических действий составляет [math]N^3(6\text{log}_2 N+1)[/math].

1.7 Информационный граф

Представим информационный граф алгоритма для этапов 1-4, поскольку этапы 5-7 реализуются аналогично этапам 1-3.

Зеленые прямоугольники обозначают одномерные сечения входного трехмерного массива. Красным кругом обозначена операция одномерного БПФ, например [math]FFT_{x11}[/math] означает БПФ вдоль оси x при [math]j=1,~k=1[/math]. Желтыми кругами представлены поэлементные преобразования массива для получения коэффициентов Фурье решения (этап 4 алгоритма), например [math]L_{2N,1,...,N}[/math] означает преобразование элементов [math]i=2,~j=N,~k=1,...,N[/math]. Для удобства целочисленные индексы нумеруются числами [math]1,...,N[/math] вместо [math]0,...,N-1[/math]

Информационная зависимость между ярусами графа заключается в том, что одномерное БПФ зависит от результата произведенного ранее БПФ по перпендикулярному направлению (т.е. с предыдущего яруса), только если обрабатываемые этими БПФ одномерные сечения массива пересекаются. Следовательно, каждое из БПФ по второй координате, [math]FFT_{i,y,k},~i=1,...,N[/math], зависит от результатов каждого из БПФ по первой координате [math]FFT_{x,j,k},~j=1,...,N[/math], и только от них, если эти преобразования производятся в одной плоскости [math]k=K[/math]. Соответственно, каждое БПФ по третьей координате [math]FFT_{i,j,z},~j=1,...,N[/math], зависит от каждого БПФ по второй координате [math]FFT_{i,y,k},~k=1,...,N[/math], в одной плоскости [math]i=I[/math].

1.8 Ресурс параллелизма алгоритма

В описанном выше алгоритме можно выделить, по меньшей мере, два уровня параллелизма. Во-первых, вычисление каждого одномерного БПФ может быть распределено между вычислительными ядрами. И, во-вторых, на каждом этапе алгоритма (ярусе графа) одномерные БПФ независимы друг от друга и также могут быть выполнены паралельно (координатный параллелизм). Для сопоставления параллельной (при координатном распараллеливании) и последовательной сложности алгоритма определим в качестве основных операций одномерное БПФ и поэлементный перерасчет массива (этап 4 алгоритма). Тогда последовательный алгоритм будет состоять из последовательно выполняемых [math]6N^2[/math] одномерных БПФ и [math]N^3[/math] операций над отдельными элементами массива, параллельный же алгоритм будет выполнен за 6 шагов одномерного БПФ и 1 операцию поэлементного преобразования. Если же выразить введенные макрооперации через элементарные арифметические операции, то получим, что сложность последовательного алгоритма - [math]6N^3(\text{log}_2N)+N^3[/math], а сложность параллельного - [math]6N(\text{log}_2N)+1[/math]. Тогда при больших [math]N[/math] отношение последовательной сложности к параллельной становится равным [math]\propto N^2(6\text{log}_2N \text{ln}2+\text{ln}2+5)[/math].

1.9 Входные и выходные данные алгоритма

В качестве входных данных выступает трехмерный массив [math]N \times N \times N[/math] правой части уравнения. Выходной массив - трехмерный массив решения той же размерности.

1.10 Свойства алгоритма

Как было показано выше, при неограниченном количестве вычислительных устройств и больших объемах входных данных вычислительная сложность параллельного алгоритма убывает быстрее, чем [math]N^2[/math].

Поскольку объем входных данных составляет [math]N^3[/math], вычислительная мощность последовательного алгоритма становится равной [math]6\text{log}_2N+1[/math], а параллельного алгоритма - [math](6N\text{log}_2N+1)/N^3 \propto N^{-2}[/math] (при больших [math]N[/math]).

2 Программная реализация алгоритма

2.1 Особенности реализации последовательного алгоритма

Приведем пример реализации алгоритма на языке Фортран-90 с применением библиотеки FFTW. В этом случае трехмерное прямое преобразование Фурье правой части уравнения реализуется как последовательность трех одномерных (этапы 1-3 в п.1.5). Полученные коэффициенты Фурье делятся на собственное число оператора Лапласа (этап 4 в п.1.5) (см. комментарии в коде):

         ! x - transform
         do k = 1, k1
           do j = 1, j1
             workfor_x(1:i1) = cmplx(f_r(1:i1,j,k),f_i(1:i1,j,k))
             call FFTW_EXECUTE_DFT(planfor_x,workfor_x,workfor_x)
             f_r(1:i1,j,k) = real(real (workfor_x(1:i1)) )
             f_i(1:i1,j,k) = real(aimag(workfor_x(1:i1)) )
           enddo
         enddo
         ! y - transform
         do k = 1, k1
           do i = 1, i1
             workfor_y(1:j1) = cmplx(f_r(i,1:j1,k),f_i(i,1:j1,k))
             call FFTW_EXECUTE_DFT(planfor_y,workfor_y,workfor_y)
             f_r(i,1:j1,k) = real(real (workfor_y(1:j1)) )
             f_i(i,1:j1,k) = real(aimag(workfor_y(1:j1)) )
           enddo
         enddo
         ! z - transform and dividing by eigenvalue
         do j = 1, j1
           do i = 1, i1
             workfor_z(1:k1) = cmplx(f_r(i,j,1:k1),f_i(i,j,1:k1))
             call FFTW_EXECUTE_DFT(planfor_z,workfor_z,workfor_z)
             f_r(i,j,1:k1) = real(real (workfor_z(1:k1)) )
             f_i(i,j,1:k1) = real(aimag(workfor_z(1:k1)) )
             ! Dividing by eigenvalue of Laplace operator
             sinx = sin((i-1) * Pi / i1)^2 * 4
             siny = sin((j-1) * Pi / j1)^2 * 4
             do k = 1, k1
               sinz = sin((k-1) * Pi / k1)^2 * 4
               f_r(i,j,1:k1) = f_r(i,j,1:k1)/(sinx + siny + sinz)
               f_i(i,j,1:k1) = f_i(i,j,1:k1)/(sinx + siny + sinz)
             enddo
           enddo
         enddo

Как только коэффициенты Фурье решения найдены, аналогично производится обратное преобразование Фурье:

         ! inverse x - transform
         do k = 1, k1
           do j = 1, j1
             workback_x(1:i1) = cmplx(f_r(1:i1,j,k),f_i(1:i1,j,k))
             call FFTW_EXECUTE_DFT(planback_x,workback_x,workback_x)
             f_r(1:i1,j,k) = real(real (workback_x(1:i1)) )/real(i1)
             f_i(1:i1,j,k) = real(aimag(workback_x(1:i1)) )/real(i1)
           enddo
         enddo
         ! inverse y - transform
         do k = 1, k1
           do i = 1, i1
             workback_y(1:j1) = cmplx(f_r(i,1:j1,k),f_i(i,1:j1,k))
             call FFTW_EXECUTE_DFT(planback_y,workback_y,workback_y)
             f_r(i,1:j1,k) = real(real (workback_y(1:j1)) )/real(j1)
             f_i(i,1:j1,k) = real(aimag(workback_y(1:j1)) )/real(j1)
           enddo
         enddo
         ! inverse z - transform
         do j = 1, j1
           do i = 1, i1
             workback_z(1:k1) = cmplx(f_r(i,j,1:k1),f_i(i,j,1:k1))
             call FFTW_EXECUTE_DFT(planback_z,workback_z,workback_z)
             f_r(i,j,1:k1) = real(real (workback_z(1:k1)) )
             f_i(i,j,1:k1) = real(aimag(workback_z(1:k1)) )
           enddo
         enddo

Как видим, вложенность циклов соответствует порядку расположения элементов массивов в языке Фортран - по столбцам - это позволяет существенно сократить время обращения к памяти. Многократно используются одни и те же вспомогательные массивы, что при небольших размерах сетки должно обеспечить быстроту работы с кэш-памятью. Однако, данная реализация не является оптимальной. Так, вместо многократного вызова подпрограммы FFTW_EXECUTE_DFT в циклах, можно было бы обращаться к одной подпрограмме, выполняющей заданное множество одномерных БПФ.

Из приведенных выше листингов программы становится очевиден также и ресурс параллелизма алгоритма: каждое одномерное преобразование Фурье выполняется независимо. Самым простым способом использования этого обстоятельства будет вставка соответствующих директив OpenMP перед внутренними циклами - в случае использования многоядерной архтитектуры компьютера на общей памяти. В случае же параллельного выполнения на кластерах с распределенной памятью, необходимо будет позаботиться о распределении входного массива (правой части уравнения) между процессорами, а также "транспозиции" массивов при переходе от БПФ по одному из направлений к БПФ по следующему направлению.

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

Важным обстоятельством является то, что большинство вариантов БПФ, включая алгоритм Кули-Тьюки, в арифметике с плавающей запятой имеют очень низкую погрешность. Верхняя граница относительной ошибки для алгоритма Кули-Тьюки составляет [math]O(\epsilon \log N)[/math] (для сравнения - относительная ошибка преобразования Фурье непосредственно по формуле дискретного преобразования Фурье - [math]O(\epsilon N^{3/2})[/math]) . Поэтому для ускорения алгоритма можно использовать одинарную точность представления переменных без существенного увеличения ошибки.

2.2 Возможные способы и особенности параллельной реализации алгоритма

Как было показано выше, существует, по меньшей мере, два уровня параллелизма в описываемом алгоритме. Во-первых, это параллельная реализация одномерных БПФ и, во-вторых, координатный параллелизм - параллельное исполнение независимых одномерных БПФ . Эти два уровня могут быть скомбинированы, например, путем распределения разных БПФ между MPI-процессами, а каждого отдельного БПФ - между нитями на общей памяти (POSIX, OpenMP, и т.д.). Ниже мы приведем детали базового алгоритма распараллеливания с использованием координатного параллелизма, рассмотрим его недостатки и перечислим возможные способы повышения его масштабируемости.

Координатный параллелизм в контексте данной задачи реализуется в результате разбиения трехмерной области на подобласти, в каждой из которых можно независимо от расчетов в других подобластях, выполнить одномерные БПФ. К настоящему времени предложены две реализации такого подхода: одномерное разбиение области (slab decomposition) и двумерное разбиение (pencil decomposition) . При одномерном разбиении требуется меньшее количество MPI-пересылок, однако, максимальное количество процессоров/ядер ограничено максимальным по трем измерениям размером области. Поэтому с точки зрения массивного параллелизма двумерная декомпозиция предпочтительнее, и ниже мы остановимся на этом подходе.

При MPI-реализации необходимо принять во внимание накладные расходы, вызванные пересылками данных между процессорами. Для реализации обсуждаемого алгоритма необходимо осуществлять т.н. транспонирования массивов, при которых на каждом из трех этапов преобразования Фурье (этапы 1-3 и 5-7 в п.1.5) обеспечивается нахождение элементов массива вдоль соответствующих направлений в памяти соответствующих процессоров:

Транспонирование трехмерного массива при параллельной реализации трехмерного БПФ

Рассмотрим соотношение вычислений и пересылок данных при MPI-реализации. Вычислительная сложность алгоритма при покоординатном распараллеливании, при запуске на [math]M[/math] MPI-процессах составляет [math] \left(6N^3(\log_2 N) + N^3\right)/M[/math], т.е., в лучшем случае, [math] 6N(\log_2 N) + N[/math] (на [math]N^2[/math] процессах). При этом, суммарное количество отсылаемых и получаемых MPI-пакетов составляет [math]8N_{px}+4N_{py}[/math], а суммарное количество отсылаемых и получаемых элементов массива на каждом процессе - [math]12N^3/(N_{px}N_{py})[/math] (здесь - [math]N_{px},~N_{py}[/math] - количества процессов в двумерной MPI-декомпозиции области по осям [math]x[/math] и [math]y[/math], соответственно - см. рис. выше). Учитывая, что [math]M=N_{px}N_{py}[/math], получаем, что отношение количества пересылаемых элементов массива к количеству арифметических операций составляет [math]\alpha=12/(6\log_2 N + 1)[/math]. Таким образом, отношение объема пересылаемых элементов к количеству операций не зависит от количества MPI-процессов и уменьшается с ростом размера задачи. При больших [math]N[/math], однако, эта зависимость становится слабой, т.е. "насыщается". В то же время, при фиксированном размере задачи, рост количества используемых процессов, не изменяя параметра [math]\alpha[/math], вызывает рост числа пакетов ([math]\propto \sqrt{M}[/math], в силу [math]N_{px} \propto N_{py} \propto \sqrt{M}[/math] и формул, приведенных выше). В силу латентности сети, это приводит к росту времени пересылки единицы данных и, следовательно, замедлению ускорения алгоритма.

Приведенный здесь анализ может служить только для ориентировочных оценок, поскольку он не дает времен вычислений и MPI-обменов. Более детальные расчеты для случаев одномерной и двумерной декомпозиции области обширно представлены в литературе .

В профилируемом в следующих разделах программном коде параллельной реализации алгоритма вызовы одномерных БПФ по трем направлениям, показанные в листинге п.2.1, перемежаются вызовами специальных процедур, осуществляющих транспонирование массива. Отдельно транспонируются действительная и мнимая части массива (что уже не оптимально - лучше пересылать комплексные числа) с помощью коммуникаций "точка-точка". При этом используются асинхронные пересылки с помощью MPI_ISEND и MPI_IRECV.