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

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

Содержание

1 Свойства и структура алгоритма
2 Программная реализация алгоритма
3 Литература
4 Свойства и структура алгоритма
5 Программная реализация алгоритма
6 Литература

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

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

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

$\begin{equation} \label{eq:poisson} \sum_{i=1}^{N}\frac{\partial^2 \phi}{\partial x_i^2}=f,~\mathbf{x}\in D, \end{equation}$

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

$\begin{equation} B(\phi)=F,~ \text{на}~ \mathbf{x}\in \Gamma(D), \end{equation}$

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

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

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

$\begin{equation} \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. \end{equation}$

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

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

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

$\begin{equation} \label{eq:poisdiscr} \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, \end{equation}$

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

$\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}$

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

$\begin{equation} \label{eq:fourier} \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)}. \end{equation}$

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

$\begin{equation} \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]}. \end{equation}$

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

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

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

$\begin{equation} \label{eq:fourier2} \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], \end{equation}$

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

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

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

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

Введем для краткости следующие обозначения сеточных функций $\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\}$ и аналогично $f_h$ и $\Phi_h$ . Тогда алгоритм запишется следующим образом:

Вычисление $\text{FFT}_x(f_h)$
Вычисление $\text{FFT}_y(\text{FFT}_x(f_h))$
Вычисление $\text{FFT}_z(\text{FFT}_y(\text{FFT}_x(f_h)))=F_h$
Расчет $\Phi_h$ по $F_h$
Вычисление $\text{IFFT}_x(\Phi_h)$
Вычисление $\text{IFFT}_y(\text{IFFT}_x(\Phi_h))$
Вычисление $\text{IFFT}_z(\text{IFFT}_y(\text{IFFT}_x(\Phi_h)))=\phi_h$

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

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

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

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

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

Зеленые прямоугольники обозначают одномерные сечения входного трехмерного массива. Красным кругом обозначена операция одномерного БПФ, например

$FFT_{x11}$ означает БПФ вдоль оси x при

$j=1,~k=1$ . Желтыми кругами представлены поэлементные преобразования массива для получения коэффициентов Фурье решения (этап 4 алгоритма), например

$L_{2N,1,...,N}$ означает преобразование элементов

$i=2,~j=N,~k=1,...,N$ . Целочисленные индексы для краткости записи перенумерованы от

$0,...,N-1$ до

$1,...,N$

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

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

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

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

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

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

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

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

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
             f_r(i,j,1:k1) = f_r(i,j,1:k1)/(sinx(i)+siny(j)+sinz(1:k1))
             f_i(i,j,1:k1) = f_i(i,j,1:k1)/(sinx(i)+siny(j)+sinz(1:k1))
           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).

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

2.2 Локальность данных и вычислений

В настоящем разделе мы рассмотрим свойства приведенного выше последовательного алгортитма с точки зрения оптимальности взаимодействия с кэш-памятью. С этой точки зрения ключевым фактором является локальность обращений к этой памяти - как пространственная локальность (обращение к близко раположенным адресам), так и временная (частое обращение к одним и тем же адресам). В случае, если алгоритм не обладает достаточной локальностью обращений к памяти, будут происходить частые обращения к более медленной памяти (кэш-памяти более высокого уровня или основной оперативной памяти) - т.н. кэш-промахи. Один кэш-промах на уровне L1 на современных процессорах "стоит" около 10 циклов, а на более высоких уровнях кэша - до 200 циклов.

Ниже мы приводим статистику по кэш-промахам обсуждаемого алгоритма, собранную утилитой cachegrind , входящей в состав пакета valgrind. Расчет произведены на процессоре Intel® Core™ i7-3520M CPU @ 2.90GHz × 4 со следующими уровнями кэш-памяти: L1d 32 КБ, L1i 32 КБ, L2 256 КБ, L3 4096 KБ (буквы i, d означают кэш для инструкций и кэш для данных, соответственно). В расчетах использовалась двойная точность.

Доля кэш-промахов при реализации алгоритма на трехмерных сетках размера

$N\times N\times N$ . Dw - общее число операций записи, Dr - общее число операций чтения, D1mr - промах чтения из кэш 1-го уровня, D1mw - промах записи в кэш 1-го уровня, Dlmr - промах чтения из кэш верхнего (в данном случае, 3-го) уровня, Dlmw - промах записи в кэш верхнего уровня

Из приведенной картинки можно сделать несколько выводов. Во-первых, количество промахов в кэш первого уровня больше, чем в кэш 3-го уровня, и это связано с соотношением размеров этих уровней. Во-вторых, с увеличением размера задачи кэш-промахи учащаются. В связи с этим отметим, что в тестируемом коде используется несколько трехмерных массивов, размер которых даже в задаче $20\times 20 \times 20$ не умещается в кэш 1-го уровня. Поэтому относительно небольшая доля кэш-промахов 1-го уровня свидетельствует об работе аппаратного алгоритма предварительной загрузки (предвыборки) в кэш (prefetching). Однако, при росте размера массивов, этот алгоритм, по-видимому, становится менее эффективным.

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

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

Рассмотрим соотношение вычислений и пересылок данных при MPI-реализации. Вычислительная сложность алгоритма при покоординатном распараллеливании, при запуске на $M$ MPI-процессах составляет $\left(6N^3(\log_2 N) + N^3\right)/M$ , т.е., в лучшем случае, $6N(\log_2 N) + N$ (на $N^2$ процессах). Однако, необходимо принять во внимание накладные расходы, вызванные пересылками данных между процессорами. Для реализации обсуждаемого алгоритма необходимо осуществлять т.н. транспонирования массивов, при которых на каждом из трех этапов преобразования Фурье (этапы 1-3 и 5-7 в п.1.5) обеспечивается нахождение элементов массива вдоль соответствующих направлений в памяти соответствующих процессоров (см. рис.).

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

2.4 Масштабируемость алгоритма и его реализации

2.5 Динамические характеристики и эффективность реализации алгоритма

2.6 Выводы для классов архитектур

2.7 Существующие реализации алгоритма

3 Литература

↑ Марчук Г.И. Методы вычислительной математики, М., 1977, 456с.
↑ Gentleman, W. M.; Sande, G. (1966). "Fast Fourier transforms—for fun and profit". Proc. AFIPS 29: 563–578. doi:10.1145/1464291.1464352

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

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

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