CC BY
71
ДВУХСЕТОЧНЫЕ ПАРАЛЛЕЛЬНЫЕ АЛГОРИТМЫ ДЛЯ РЕШЕНИЯ ДРОБНО-ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ АНОМАЛЬНОЙ ДИФФУЗИИ
С.Ю. Лукащук
Приводятся описание и анализ параллельных алгоритмов решения начальнокраевых задач для уравнений аномальной диффузии, содержащих производные дробного порядка типа Римана-Лиувилля по пространственным и/или временной переменным. Параллельные алгоритмы построены на основе двухсеточного подхода. При этом грубая сетка используется для расчета эффектов пространственного и временного дальнодействия с использованием сплайн-аппроксимации, а мелкая сетка служит для конечно-разностной дискретизации решаемых уравнений. Рассматриваются алгоритмы с декомпозицией как по пространству, так и по времени. Для распараллеливания по времени используется подход, предложенный в известном алгоритме PARAREAL. Приводятся теоретические оценки параллельной эффективности предложенных алгоритмов. Показано, что алгоритмы имеют сверхлинейное ускорение по сравнению с классическим последовательным конечно-разностным алгоритмом и обеспечивают тот же порядок точности вычислений при условии согласованного выбора шагов точной и грубой сеток. Также приводятся некоторые результаты вычислительных экспериментов, подтверждающие эффективность предложенных алгоритмов.
Ключевые слова: двухсеточный параллельный алгоритм, аномальная диффузия, дифференциальное уравнение дробного порядка.
Введение
Во многих неупорядоченных сложных средах, таких как пористые и трещиноватопористые среды, аморфные полупроводники, жидкие полимеры и стекла, перколя-ционные кластеры и самоподобные фрактальные среды, турбулентные потоки жидкости, газа и плазмы, кинетика протекания процессов диффузионного переноса часто отличается от классической, соответствующей нормальному закону распределения [1, 2]. Такие процессы принято называть процессами аномальной диффузии. Аномальный диффузионный перенос обычно обусловлен эффектами памяти среды, пространственной нелокальности и перемежаемости, имеет немарковскую стохастическую природу и поэтому не может быть описан классическими эволюционными уравнениями переноса [3].
Одним из эффективных и широко используемых подходов для описания процессов аномальной диффузии является использование аппарата интегро-дифференцирования дробного порядка [4-6]. В этом случае уравнение процесса является интегро-дифференциальным уравнением, содержащим производные дробного порядка по временной и/или пространственным переменным.
Основной проблемой численного моделирования процессов аномальной диффузии, описываемых уравнениями дробного порядка, является нелокальный характер последних. При использовании конечно-разностной дискретизации обычных диффузионных уравнений значение в текущем узле расчетной сетки зависит только от значений в соседних узлах. При конечно-разностной дискретизации уравнений дробного порядка значение в текущем узле будет зависеть от всех остальных пространственных узлов (при наличии в уравнении дробных производных по пространственным переменным) и всех предыдущих временных слоев (при наличии дробной производной
по времени) [5]. В результате при численных расчетах аномального переноса объем вычислений оказывается существенно (на порядки!) большим по сравнению с объемом вычислений в расчетах классических диффузионных процессов. При наличии дробной производной по времени объем вычислений также растет с увеличением номера временного слоя, поскольку увеличивается диапазон памяти системы, которую нужно учитывать. Существенно возрастает и объем необходимой оперативной памяти для хранения результатов расчета. Использование параллельных алгоритмов, построенных только на принципах пространственной декомпозиции расчетной области (domain decomposition), ситуацию не спасает, поскольку приводит к чрезвычайно большим объемам передачи данных между процессорами.
Указанные проблемы делают актуальной задачу разработки новых последовательных и параллельных численных алгоритмов решения уравнений аномальной диффузии с производными дробного порядка. В данной работе предлагается подход к построению таких алгоритмов, основанный на использовании идеи двухсеточных методов. Грубая сетка используется в этом случае как основа для расчета эффектов пространственного и временного дальнодействия, а мелкая сетка служит для непосредственных вычислений. По узлам грубой сетки может производиться пространственная и, при необходимости, временная декомпозиция расчетной области, что приводит к легко реализуемым параллельным алгоритмам. При распараллеливании по временной оси данный подход приводит к модификации алгоритма PARAREAL для уравнений с производными дробного порядка [7].
1. Постановка задачи
При описании процессов аномальной диффузии наиболее часто используются левосторонняя и правосторонняя производные дробного порядка а типа Римана-Лиувилля [4]:
где п = [а] + 1 и [а] — целая часть числа а.
В дальнейшем для иллюстрации предлагаемого подхода будем рассматривать уравнение аномальной диффузии дробного порядка следующего вида:
где и0(ж), м0(£), и^(£) — заданные функции. Отметим, что постановка таких краевых условий предъявляет определенные дополнительные требования к виду функции f (ж,£), которые однако не являются важными для рассматриваемых алгоритмов и поэтому здесь не приводятся.
(І)
(2)
(0D"u) (x,t) = 7 (0Dl+eu) (x,t) + (1 - y) (xDL+eu) (x,t) + f (x,t) (3)
t Є (0,T], x Є (0,L), a,в,7 Є (0,1).
Для уравнения (3) поставим первую краевую задачу:
u(x, 0) = u0(x), u(0,t) = u0(t), u(L,t)= uL(t), (4)
(4)
При построении численных схем для дробно-дифференциальных уравнений используются различные конечно-разностные аппроксимации дробных производных. Для левосторонней дробной производной (1) в случае равномерной сетки все они могут быть представлены в следующем общем виде:
(здесь к — шаг сетки). Аналогичная формула имеет место и для правосторонней дроб-
классическую аппроксимацию, получаемую из формулы Грюнвальда-Летникова [4], и обеспечивающую первый порядок точности аппроксимации дробной производной для аналитических функций [5]. Повысить порядок точности можно использованием полиномиальной интерполяции искомой функции между узлами расчетной сетки с последующим аналитическим вычислением соответствующих интегралов. Случай линейной межузельной интерполяции, обеспечивающий второй порядок точности, рассмотрен в [8].
2. Декомпозиция задачи на основе двухсеточного подхода
Введем в рассматриваемой пространственно-временной области П = [0,Ь] х [0,Т] равномерную конечно-разностную сетку с шагом АХ по пространственной и АТ по временной переменным:
Сетку (6) в дальнейшем будем называть грубой сеткой. Каждую ячейку, ограниченную смежными узлами грубой сетки, будем рассматривать как отдельную расчетную подобласть:
i=lу=1
Для (£,я) € производные дробного порядка, входящие в уравнение (3), могут быть записаны следующим образом:
к+т
j=о
ной производной (2). При т = 0 и Aj = (—1):?'(а) имеем наиболее часто используемую
шо = {(X = гАХ, Ту = зАТ), г = 0,1,... ,Х,
3 = 0,1,... ,М, АХ = Ь/Х, АТ = Т/М}.
(6)
N М
= [Xi-l,Xi] X [Т/_1,Т], П = у и .
где введены следующие обозначения для интегралов:
Уравнение (3) в области может быть теперь переписано в следующем виде:
— линейный интегральный оператор. Принципиально важным является тот факт, что вычисление Ьи в области о,,- производится по значениям функции и(ж,£) при (х2 ^) ^ 0І,j.
Обозначим через и = {Ц^-, г = 0,1,... , Х, з = 0,1,... , М} сеточную функцию решения задачи на грубой сетке шо В предлагаемом подходе ее значения используются для постановки начальных и граничных условий для уравнения (7), а также вычисления интегралов, входящих в Ьи.
После того, как задача для уравнения (7) замыкается с помощью начальных и граничных условий и указан способ вычисления Ьи по значениям и, численное решение задач в каждой подобласти о,— или соответствующим образом выбранной группе подобластей может выполняться параллельно. Для этого в каждой подобласти вводится точная сетка (для простоты считаем, что разбиение всех подобластей идентично):
С использованием аппроксимаций вида (5) на сетке шр строится конечно-разностная аппроксимация уравнений (7) требуемого уровня точности. Как и в случае классического уравнения диффузии, могут быть построены явные, неявные или полу-явные схемы. Сеточную функцию задачи на точной сетке будем обозначать через и = {ику1, к = 0,1,...,п, I = 0,1,...,т}.
Таким образом, точная сетка шр используется для параллельного решения задач в отдельных подобластях, которые могут выполняться параллельно, а грубая сетка шо используется для постановки краевых условий и вычисления дополнительной функции Ьи.
На основе такого двухсеточного подхода возможно построение алгоритмов двух видов. В алгоритмах первого вида счет ведется последовательно (например, по времени) и в узлы грубой сетки просто копируются значения, полученные при расчете на точной сетке. Затем эти значения рассылаются по всем подобластям, в которых они необходимы для вычисления соответствующих элементов. Второй вид алгоритмов в более полной мере реализует идею двухсеточного подхода и основан на схемах типа ,, предиктор-корректор “. В этом случае каким-либо образом рассчитываются начальные приближения в узлах грубой сетки, а затем проводятся расчеты (параллельные) на точных сетках и уточняются приближения на грубой сетке. К алгоритму второго типа относится, в частности, модификация алгоритма РАИАИЕАЬ для распараллеливания по времени. Далее будут рассмотрены характерные представители алгоритмов обоих видов.
где
Ь = 7 оJxi-1 + (1 - 7) оJxi-1 - о^т-х
шр = {(як = к&т, ^ = Ш), к = 0,1,..., п,
I = 0,1,...,т, = АХ/п, Л = АТ/т}.
3. Описание параллельных алгоритмов
з.1. Параллельный алгоритм с декомпозицией по пространству
Данный алгоритм относится к алгоритмам первого вида. Расчет по времени ведется последовательно, поэтому шаги по времени грубой и точной сеток совпадают: ($£ = АТ, М = т). Характерными особенностями алгоритма, обеспечивающими требуемый уровень точности вычислений, являются: 1) использование кубической сплайн-интерполяции по значениям в узлах грубой сетки на одном временном слое, 2) использование полиномов второй и третьей степени для аппроксимации в окрестности границ подобластей 0,- производных дробного порядка по пространственной переменной.
Сохранение заданной точности вычислений производных дробного порядка по пространству после декомпозиции на подобласти 0,- является главной проблемой при разработке соответствующих параллельных алгоритмов. Использование классических аппроксимаций вида (5) (типа Грюнвальда-Летникова) для смещенных дробных производных, входящих в уравнение (7), может приводить к катастрофической потере точности вычислений в окрестности границ областей 0,- несмотря на сохранение общего порядка точности аппроксимации. Для подавления ошибки предлагается использовать следующие вычислительные приемы.
Во-первых, необходимо выполнить сдвиг вычисляемой функции на величину соответствующего граничного значения:
х4_1 А1+ви = х4_1 А1+в [и(м) - и(Х,-!,^)] + Г-1(-в)(я - Х,-!)-в-1и(Х,-!,^), (9) х^Х+ви = х^Х+в [и(я,г) - и(Х,,£)] + Г-1(-в)(Х, - х)_в-1и(Х,, £). (10)
В результате основная составляющая производной оказывается вычисленной аналитически, а по приближенным формулам рассчитываются уже относительно небольшие поправки.
Во-вторых, при расчете левосторонней производной из правой части (9) в точке я1 и правосторонней производной из правой части (10) в точке жга-1 коэффициенты Аг соответствующих аппроксимаций вида (5) находятся на основе явного вычисления этих производных по трехузловой квадратичной интерполяции. Для следующих точек х2 и жга_2 используется четырехузловая кубическая интерполяция. Для остальных точек может быть использована классическая аппроксимация первого или второго порядка точности.
В-третьих, вычисление интегралов 0^х_ 1 и и х1 JLU производится аналитически для кубической сплайн-интерполяции, построенной по значениям на грубой сетке
и. Известно [9], что оценка погрешности аппроксимации функции f (я) кубическим сплайном з(ж) на равномерной сетке с шагом Н имеет вид |з(ж) - f (я)| < 5АН4/384, где \f1У(х)| < А. Тогда легко получается оценка точности вычисления интеграла
о ^х4-1 и:
(о^х4-1 и (я,^) - з(я,^)( < 5А(£)(АХ)4(£я)_в-1/384,
где |ижжжж(я,^)| < А(£). Аналогичная оценка справедлива для интеграла х1 JLи. Поскольку множитель (£я)_1-в является общим множителем при вычислении дробных производных на точной сетке по аппроксимациям вида (5), то точность вычисления интегралов по сплайн-интерполяции в разностной схеме решения уравнения (7) на
точной сетке имеет порядок О (АХ4). Если для аппроксимации дробных производных используется аппроксимация первого порядка точности по пространству, то для сохранения этого порядка у всей схемы для уравнения (7) необходимо выполнение следующего очевидного неравенства: (АХ)4 < &г. При использовании аппроксимации второго порядка точности неравенство становится более жестким: (АХ)2 < &г.
Далее представлены основные шаги алгоритма при использовании явной схемы для дискретизации уравнений (7).
1. По заданному начальному условию задачи ио(ж) (4) каждый процессор параллельно вычисляет значения сеточной функции и на нулевом временном шаге. Следующие шаги алгоритма являются едиными для временных слоев / = 1, 2,... , т.
2. По значениям и для (/ - 1)-го временного слоя каждый процессор вычисляет коэффициенты соответствующего кубического интерполяционного сплайна.
3. На точной сетке все процессоры параллельно выполняют расчет и^ (к =
1, 2,...,п - 1) на временном слое /. При этом производная по времени аппроксимируется формулой Грюнвальда-Летникова
а\
(о А“и) (ж* ,ъ) - (Л)_“ ^(-1)г(а)им
г=о
а производные по пространству аппроксимируются на (/ - 1)-м временном слое с учетом (9) и (10) по формулам
1 *+1
(х4_1 ^+ви (ж*,^г_1) ~ (^ж)(1+в) X]А (и*+1_г,1_1 - ио,1_1) +
+ ио,1_1
1 П—*+1
Г(-в)(ж* - Х,_1)1+в’
^х4 и^ (ж*,^1_1) ~ (^ж)(1+в) ^ ' Аг (и*_1+г,1_1 - ио,1_1) +
+ ио,1_1
Г(-в)(Х, - ж*)1+в’
Коэффициенты Аг выбираются так, как описано выше. Интегралы о^х-1 и и х1 ^и вычисляются по аналитическим формулам для сплайн-интерполяции, построенной на шаге 2 .
4. Процессор, отвечающий за расчет подобласти 0,+1г (г = 1, 2,..., N - 1), пересылает процессору, отвечающему за расчет соседней подобласти 01;г значение приращения за временной шаг искомой функции в первом узле точной сетки: Ц+1 - и^ (здесь верхний индекс соответствует пространственному номеру области). Процессоры, получившие значения, вычисляют сеточную функцию в узлах грубой сетки на /-м временном слое по формуле
и1,1 = и,,1_1 + О-5 (и1+1 - и1+_1 + иП_1,г - иП_1,г_1) •
5. Выполняется взаимный обмен рассчитанными значениями сеточной функции и таким образом, чтобы полный вектор этих значений находился на каждом процессоре.
Далее переход к следующему временному слою и повтор алгоритма с шага 2.
—г
Отметим, что условие устойчивости использованной явной конечно-разностной схемы имеет вид
(1 + в)(ЭД“ < а(^ж)1+в (11)
(для случая а =1 это условие приведено в [8]). В предельном случае классического уравнения диффузии а = 1, в =1 это условие переходит в хорошо известное условие устойчивости явной разностной схемы для уравнения параболического типа.
3.2. Параллельный алгоритм с декомпозицией по времени
Известно, что параллельные алгоритмы для решения эволюционных задач могут быть построены не только на основе декомпозиции пространственной расчетной области, но и на основе декомпозиции по временной оси. Одним из наиболее широко используемых алгоритмов данного вида является алгоритм РАИАИЕАЬ, впервые предложенный в работе [10]. В настоящее время данный алгоритм хорошо изучен, получены теоретические оценки его точности и сходимости [11-14], предложен целый ряд его модификаций для решения различных классов задач эволюционного типа [11, 15, 16]. В работе [7] была предложена модификация алгоритма РАИАИЕАЬ для решения обыкновенных дифференциальных уравнений с производной дробного порядка а Е (0,1). В данном разделе дается описание новой модификации алгоритма РАИАИЕАЬ, совмещенной с приведенным в предыдущем разделе алгоритмом декомпозиции по пространству, и предназначенной для решения уравнений вида (3).
Алгоритм базируется на пространственно-временной декомпозиции расчетной области 0 на основе грубой сетки шс, представленной в разделе 3, и двумерной сплайн-интерполяции решения по значениям сеточной функции и в узлах грубой сетки.
Пусть с использованием некоторого сеточного метода С, который в дальнейшем будем называть ,,грубым“, построено приближенное решение и краевой задачи (3), (4) в узлах грубой сетки. Рассмотрим теперь задачу для уравнения (7) в области 0,-. Обозначим через и1- функцию решения этого уравнения. С использованием сплайн-интерполяции по значениям и на (3 - 1)-м временном слое можно записать начальное условие для уравнения (7) в виде
и1-(ж,Т_1) = ф(ж, ио-_1, и1-_1,..., Цу-_1). (12)
Пусть интерполяция по временной переменной поля значений и выполняется таким образом, что значения интерполяционных коэффициентов для произвольного (3 - 1)-го слоя зависят только от значений и на этом и предыдущих временных слоях. Эта интерполяция используется для вычисления интегралов о!т,_1. Тогда решение уравнения (7) с начальным условием (12) будет зависеть только от значений и на первых (3 - 1)-м временных слоях:
и1- = и1-(ж, £, и*,г), к = 0,1,..., N / = 0,1,... ,3 - 1, (ж, £) € 0,-. (13)
Отметим, что для получения решений вида (13) необходимо, вообще говоря, на каждом временном слое решать систему уравнений вида (7) для всех подобластей, принадлежащих данному временному слою, с условиями сопряжения решений на границах подобластей. Однако, при численном решении этой задачи на мелкой сетке может быть эффективно использован алгоритм, приведенный в п. 4.1.
Зная решения (13), можно записать уравнения для нахождения значений в узлах
Система (14) является системой нелинейных уравнений относительно сеточной функции и и может быть решена методом Ньютона:
(здесь д — номер итерации метода Ньютона).
Ключевым моментом алгоритма РАИАИЕАЬ является способ вычисления производных, входящих в (15):
где через О обозначено решение задачи на грубой сетке, получаемое с помощью ,,грубого “ метода. Таким образом, для вычисления производных на каждом следующем временном слое используется разность грубых“ решений задачи, полученных по старым и новым, уже вычисленным на предыдущих временных слоях, значениям и. Обозначим через Т(X*, 7}, и) численное решение задачи на точной сетке в области . Полагая и*’-7'(Х*,!}, и) ~ Т(Х*,!}, и), на основе (15), (16) получаем основную расчетную формулу алгоритма РАИЛКЕАЬ:
Данная вычислительная процедура позволяет параллельно по времени итерационно находить приближенные решения задачи в узлах грубой сетки, сопоставимые по точности со значениями, получаемыми при расчете на точной сетке во всей расчетной области.
В результате приходим к следующему параллельному алгоритму (полагаем, что имеется N х М-процессорная решетка {р*,}} и за вычисление каждой области отвечает свой процессор).
1. С использованием классической явной конечно-разностной схемы О, построенной на основе аппроксимаций дробных производных вида (5), каждый процессор находит начальное приближение решения в узлах грубой сетки и0.
2. Организуется глобальная итерационная по д процедура метода Ньютона. На первом шаге д = 0.
3. По значениям Ц?г (к = 0,1,..., N I = 0,1,..., і — 1) каждый процессор строит двумерную кубическую по каждой переменной сплайн-интерполяцию поля грубого решения.
грубой сетки:
и*,7 = и*’7 (Х*,Т} ,им), к = 0,1,..., N I = 0,1,...,і — 1,
і = 1, 2,..., N і = 1, 2,..., М.
(14)
N 7-1
и»;1 = №,7} ,и») + ££
(15)
к=0 1=0
і = 1, 2,..., N і = 1, 2,... ,М
N 7-1
дИ’7 (Х*,Т} ,и»)
ЕЕ
№»+1 — цу « С(Хі,Т,,и«+1) — С(Х„7}, и»), (16)
к=0 1=0
(17)
к = 0,1,..., N I = 0,1,...,і — 1, і = 1, 2,..., N і = 1, 2,..., М.
4. Все процессоры параллельно производят расчет каждый своей подобласти на точной сетке. При этом группы процессоров, отвечающие за расчет одинаковых временных слоев, реализуют алгоритм с декомпозицией по пространству, приведенный в п. 4.1. Интегралы, образующие Ьи в уравнениях (3), рассчитываются по аналитическим формулам, полученным по принятому виду сплайн-интерполяции.
5. Последовательно (по временным слоям) рассчитывается следующее приближение в узлах грубой сетки и9+1 по формуле (17). При этом каждый процессор р,-
Т 7^+1
отвечает за вычисление одного значения и, - и выполняет один дополнительный расчет ,,грубым“ методом С(Х1,Т-, и9+1). После расчета каждого нового временного слоя полученные обновленные значения и9+1 рассылаются каждому процессору.
6. Проверяется условие остановки итерационного процесса: ||и9+1 - и91| < е. Если условие не выполнено, то возврат к шагу 2. При этом временной слой д +1 из расчета исключается (на нем уже получено решение, соответствующее решению на точной сетке, которое не может быть уточнено по (17)).
Описанный алгоритм допускает простую модификацию для случая, когда количество доступных процессоров достаточно для расчета только нескольких временных слоев одновременно. Кроме того, данный алгоритм существенно упрощается в следующих двух случаях. 1) При а = 1, т.е. когда исходное уравнение (3) содержит обыкновенную (не дробную) производную по времени первого порядка — в этом случае отсутствует необходимость интерполяции по времени и алгоритм становится прямой комбинацией алгоритма п. 4.1 и классического алгоритма РАИАИЕАЬ. 2) Если отсутствует пространственная декомпозиция — в этом случае отсутствует необходимость интерполяции по пространственной переменной и алгоритм становится прямой комбинацией модифицированного на уравнения дробного порядка алгоритма РАИАИЕАЬ и классического последовательного алгоритма конечно-разностного решения таких уравнений.
4. Оценка параллельной эффективности алгоритмов
4.1. Оценка эффективности параллельного алгоритма с декомпозицией по пространству
Сначала приведем оценку трудоемкости последовательного алгоритма решения задачи (3), (4), построенного по явной схеме с аппроксимацией дробных производных по формуле Грюнвальда-Летникова. Шаг сетки по пространственной и временной переменным для последовательного алгоритма примем равным соответствующим шагам точной сетки параллельного алгоритма (в этом случае последовательное и параллельное решения будут иметь одинаковый порядок точности). Таким образом, имеем Хп - 1 расчетных узлов по пространству и М расчетных узлов по времени. Расчет /-го временного слоя по последовательному алгоритму требует 2(Хп - 1)(Хп + / + 4) операций. Принимая, как обычно при таких оценках, время выполнения одной операции за единицу, получаем время выполнения последовательного алгоритма равным Т1 = М(Хп - 1)(2Хп + М + 9).
Теперь оценим время выполнения параллельного алгоритма, приведенного в п. 4.1. Будем полагать, что количество используемых для его выполнения процес-
соров равно количеству пространственных подобластей: Р = N. Количество операций, выполняемых каждым процессором на /-м временном слое, складывается из четырех составляющих: 1) операций, необходимых для вычисления дробных производных: 2(п - 1)(п + / + 11), 2) операций, необходимых для вычисления интегралов о 1 и и х ^и по явным аналитическим формулам через сплайн-интерполяцию: 228(п - 1)(Х - 1), 3) операций, необходимых для расчета значений в узлах грубой сетки: 5(Х-1) и 4) операций, необходимых для вычисления коэффициентов интерполяционного сплайна: 22(Х- 1) (для решения системы используется метод прогонки). Таким образом, трудоемкость расчета /-го временного слоя в параллельном алгоритме составляет 2(п - 1)(п + / + 11) + [228(п - 1) + 27](Х - 1) операций, которые выполняются параллельно. Тогда время расчета всех М временных слоев в параллельном алгоритме составит Тр = М[(п - 1)(2п + М + 228Х - 206) + 27(Х - 1)]. Долей последовательного расчета, выполняемого на этапе инициализации алгоритма до начала основной вычислительной процедуры, в данной оценке пренебрегаем.
Тогда для ускорения параллельного алгоритма получаем следующую оценку:
_ (Хп - 1)(2Хп + М + 9) N(2Хп + М)
р Т (п - 1)(2п + М + 228Х - 206) + 27(Х - 1) 2п + М
(последняя оценка становится справедливой при N ^ п). Заметим, что число временных слоев М при использовании явной схемы может оказаться достаточно большим.
При стремлении к нулю шага точной сетки, т.е. при п ^ то, и неизменных остальных параметрах получаем следующую асимптотическую оценку ускорения N2 = Р2. Таким образом, по сравнению с классическим алгоритмом, предлагаемый параллельный алгоритм демонстрирует асимптотически квадратичное ускорение. Причина такого нелинейного ускорения очевидна: в то время как классический алгоритм использует для расчета все пространственные узлы, параллельный алгоритм рассчитывает дальнодействующие нелокальные эффекты с использованием сплайн-аппроксимации, что приводит к существенному уменьшению количества операций в параллельном алгоритме по сравнению с последовательным.
Приведенная оценка не учитывает времени, затрачиваемого на обмен данными между процессорами, и является справедливой для систем с общей памятью. Для систем с распределенной памятью это время необходимо учитывать. Тем не менее, на величину асимптотической оценки это время влияния не оказывает, поскольку пропорционально М ^ - 1).
4.2. Оценка эффективности параллельного алгоритма с декомпозицией по времени
Время выполнения последовательного алгоритма Т1 оценивается аналогично п. 5.1 и составляет Т1 = Мт^п - 1)(2Nn + Мт + 9) (оценка изменилась, поскольку теперь имеем Мт слоев точной сетки по времени).
Оценим время исполнения параллельного алгоритма. Как и ранее, долей последовательного расчета, выполняемого на этапе инициализации алгоритма до начала основной вычислительной процедуры, пренебрегаем. Расчет выполняется на грубой шс и точной шр сетках, определенных в (6) и (8), соответственно.
Наибольшее число операций будут выполнять процессоры, отвечающие за расчет подобластей 01;м, соответствующих последнему М-му временному слою, поэтому
именно они и будут определять время реализации алгоритма. Для каждой итерации д количество операций, выполняемых каждым таким процессором, складывается из четырех составляющих: 1) операций, необходимых для вычисления дробных производных: т(п - 1)(2п + т + 23), 2) операций, необходимых для вычисления интегралов оJхi-1 и, хг^ьи и о/тм-1 по явным аналитическим формулам через сплайн-интерполяцию: 114т(п - 1)[2^ - 1) + (М - 1)], 3) операций, необходимых для расчета на грубой сетке по ,,грубому“ методу и уточнению значений и по формуле (17): М^ - 1)(2N + М + 9) + 2 и 4) операций, необходимых для вычисления коэффициентов интерполяционного сплайна: 22(2MN - М - 1). Сумма этих значений дает трудоемкость одной итерации метода Ньютона в параллельном алгоритме. Если для реализации алгоритма требуется Q итераций, то общее время выполнения параллельного алгоритма равно
Тр = Q[m(n - 1)(2п + т + 128N + 114М - 319) + М^ - 1)^ + М + 53) + 22(М - 1)].
По классической формуле 5Р = Т1 /Тр может быть легко оценено ускорение параллельного алгоритма. Если N ^ п, М ^ т, то Т1 « MNmn(2Nn + Мт), Тр « Qmn(2n + т) и получаем следующую оценку ускорения:
5Р = MN (2^п + Mm)Q-1(2n + т)-1.
Для случая квадратной грубой сетки (М = N) получаем 5Р = N3/Q. Учитывая, что общее количество процессоров, задействованных для выполнения параллельного алгоритма, равно Р = MN = N2, получаем следующую асимптотическую оценку ускорения: 5Р ~ Р3/2^. Максимально возможное количество итераций Q = М = л/Р. Однако на практике количество итераций всегда оказывается существенно меньше, поэтому действительная асимптотическая оценка ускорения получается сверхлинейной.
5. Результаты тестовых расчетов
Точность и эффективность предлагаемых алгоритмов была проверена рядом тестовых расчетов. Ниже представлены некоторые результаты проведенных вычислительных экспериментов.
Решается уравнение (3) в единичном квадрате (Ь = 1, Т =1) при нулевых граничных условиях ио(£) = 0, иь(£) = 0, без функции источника f (ж,£) = 0, с начальным условием ио(ж) = 8т(пж). На рис. 1 показано решение этой задачи в различные моменты времени при а = 0,95, в = 0, 7, 7 = 0, 5.
На рис. 2 показаны ошибки приближенных решений Ди для двух моментов времени £ = 0.1 и £ = 0.3, полученных с помощью классического последовательного алгоритма при т = 105 шагов по времени и количестве узлов сетки по пространству а) п = 200 и б) п = 400. Как и следовало ожидать, измельчение сетки приводит к уменьшению ошибки при неизменном характере ее распределения.
На рис. 3 показаны ошибки приближенных решений Ди для тех же моментов времени, полученных при использовании параллельного алгоритма с декомпозицией по пространству при т = М = 105 и а) количестве процессоров Р = N = 4 и количестве узлов сетки по пространству в каждой подобласти п = 100; б) количестве процессоров
Рис. 1. Решение тестовой задачи для различных моментов времени
0,0012 - \ М>'3 / 0,0007 - \ л
0,0010 0,0006-
0,0008 - / \н>.1N. 0,0005 - / N^<=0.1 N.
Ли 0,0006 - ( \ ди 0,0004 -0,0003 - / \
0,0004 - 0,0002 -
0,0002 - 0,0001 -
0,0000 - о,оооо-
0,0 0,2 0,4 0,6 0,8 1,0
X
а)
б)
0,0 0,2 0,4 0,6 0,8 1,0
X
Рис. 2. Ошибка решения Дм, полученного по классическому последовательному алгоритму: а) п = 200, б) п = 400
а)
б)
Рис. 3. Ошибка решения Ди, полученного с использованием параллельного алгоритма с декомпозицией по пространству: а) Р = N = 4, п = 100, б) Р = N = 8, п = 50
х
X
P = N = 8 и количестве узлов n = 50. Расчеты проводились на одном узле суперкомпьютера Уфимского государственного авиационного технического университета, объединяющего на общей памяти объемом 4 Гб два четырехядерных процессора Intel Xeon 5300. Из приведенных рисунков видно, что параллельный алгоритм показывает тот же порядок точности, что и последовательный.
Заключение
Разработанные двухсеточные параллельные алгоритмы решения уравнений аномальной диффузии дробного порядка показывают сверхлинейное ускорение по сравнению с классическим последовательным конечно-разностным алгоритмом и обеспечивают тот же порядок точности вычислений при условии согласованного выбора шагов точной и грубой сеток. Это свидетельствует о неоптимальности классического алгоритма. Поэтому представляется целесообразным и при последовательных расчетах использовать двухсеточные алгоритмы, получаемые очевидной редукцией предложенных параллельных алгоритмов на однопроцессорный случай.
Статья рекомендована к публикации программным комитетом международной научной конференции «Параллельные вычислительные технологии 2012».
Литература
1. Bouchaud, J.P. Anomalous diffusion in disordered media: statistical mechanisms, models and physical applications / J.P. Bouchaud, A. Georges // Phys. Rep. - 1990.
- Vol. 195. - P. 127-293.
2. Metzler, R. The random walk’s guide to anomalous diffusion: A fractional dynamic approach / R. Metzler, J. Klafter // Phys. Rep. - 2000. - Vol. 339. - P. 1-77.
3. Учайкин, В.В. Автомодельная аномальная диффузия и устойчивые законы /В.В. Учайкин // УФН. - 2003. - Т. 173, № 8. - C. 847-876.
4. Самко, С.Г. Интегралы и производные дробного порядка и некоторые их приложения / С.Г. Самко, А.А. Килбас, О.И. Маричев - Минск: Наука и техника, 1987.
5. Podlubny, I. Fractional differential equations / I. Podlubny - Academic press, San Diego, 1999.
6. Учайкин, В.В. Метод дробных производных / В.В. Учайкин - Ульяновск: Изд-во «Артишок», 2008.
7. Лукащук, С.Ю. Модификация алгоритма PARAREAL для решения дифференциальных уравнений дробного порядка. / С.Ю. Лукащук // Параллельные вычислительные технологии (ПаВТ’2010): Труды международной научной конференции. - Челябинск: Издательский центр ЮУрГУ, 2010. - С. 519-524.
8. Прямые задачи неклассического переноса радионуклидов в геологических формациях / В.М. Головизнин, В.П. Киселев, И.А. Короткин, Ю.И. Юрков // Известия Академии наук, серия «Энергетика». - 2004. - № 4. - C. 121-132.
9. Завьялов, Ю.С. Методы сплайн-функций / Ю.С. Завьялов, Б.И. Квасов, В.Л. Мирошниченко - М.: Наука, 1980.
10. Lions, J.-L. Resolution d’edp par un schema en temps parareal / J.-L. Lions, Y. Maday, G. Turinici // C.R. Acad Sci. Paris. Ser. I Math. - 2001. - Vol. 332. - P. 661-668.
11. Maday, Y. The Parareal in Time Iterative Solver: a Further Direction to Parallel Implementation / Y. Maday, G. Turinici // Lecture Notes in Computational Science and Engineering: Domain Decomposition Methods in Science and Engineering XVIII.
- 2005. - Vol. 40. - P. 441-448.
12. Staff, G.A. Stability of the Parareal Algorithm / G.A. Staff, E.M. Ronquist // Lecture Notes in Computational Science and Engineering: Domain Decomposition Methods in Science and Engineering XVIII. - 2005. - Vol. 40. - P. 449-456.
13. Gander, M.J. On the Superlinear and Linear Convergence of the Parareal Algorithm / M.J. Gander, S. Vandewalle // Lecture Notes in Computational Science and Engineering: Domain Decomposition Methods in Science and Engineering XVI. -2007. - Vol. 55. - P. 291-298.
14. Gander, M.J. Analysis of the parareal time-parallel time-integration method / M.J. Gander, S. Vandewalle // SIAM J. Sci. Comput. - 2007. - Vol. 29, No. 2. - P. 556-578.
15. Fischer, P.F. A Parareal in time semi-implicit approximation of the Navier-Stokes equations / P.F. Fischer, F. Hecht, Y. Maday // Lecture Notes in Computational Science and Engineering: Domain Decomposition Methods in Science and Engineering XVIII. - 2005. - Vol. 40. - P. 433-440.
16. Bal, G. Symplectic Parareal / G. Bal, Q. Wu // Lecture Notes in Computational Science and Engineering: Domain Decomposition Methods in Science and Engineering XVII. - 2008. - Vol. 60. - P. 401-408.
Станислав Юрьевич Лукащук, кандидат физико-математических наук, доцент, кафедра «Высокопроизводительных вычислительных технологий и систем», Уфимский государственный авиационный технический университет, [email protected].
PARALLEL TWO-GRIDS ALGORITHMS FOR SOLUTION OF ANOMALOUS DIFFUSION EQUATIONS OF FRACTIONAL ORDER
S.Yu. Lukashchuk, Ufa State Aviation Technical University (Ufa, Russian Federation)
New parallel algorithms are proposed for solving the initial-boundary value problems for anomalous diffusion equations with the Riemann-Liouville spatial- and/or time-fractional derivatives. A two-grid technique is employed to construct these algorithms. Spline-approximation on a coarse grid is used to compute the spatial and time long-range effects, and a fine grid is used for finite-difference discretization of the fractional diffusion equations. The parallel algorithms with a spatial and a time domain decomposition are discussed separately. The approach originally developed for the Parareal algorithm is used for time domain decomposition. The theoretical estimates of the speed-up and efficiency of the proposed algorithms are given. It has been shown that the algorithms have a superlinear speed-up in comparison with a classical sequential finite-difference algorithm, and have the same accuracy if the size of a fine grid is agreed with the size of a coarse grid. Some computational results are also presented to verify the efficiency of the proposed algorithms.
Keywords: parallel two-grid algorithm, anomalous diffusion, fractional differential equation.
References
1. Bouchaud J.P., Georges A. ANomalous Diffusion in Disordered Media: Statistical Mechanisms, Models and Physical APlications. Phys. Rep., 1990, Vol. 195, P. 127293.
2. Metzler R., Klafter J. The Random Walk’s Guide to ANomalous Diffusion: a
Fractional Dynamic AProach. Phys. Rep., 2000, Vol. 339, P. 1-77.
3. Uchaikin V.V. Self-Similar ANomalous Diffusion and Levy-Stable Laws. Phys. Usp., 2002, Vol. 46, No. 8, P. 821-849.
4. Samko S., Kilbas A., Marichev O. Fractional Integrals and Derivatives. Theory and Applications. Gordon & Breach Sci. Publishers, 1993. 1006 p.
5. Podlubny I. Fractional Differential Equations. Academic press, San Diego, 1999. 340 p.
6. Uchaikin V.V. Metod drobnykh proizvodnykh [The Method of Fractional Derivatives]. Ul’yaNovsk, Artishok Publ., 2008. 512 p.
7. Lukashchuk S.Yu. Modifikatsiya algoritma PARAREAL dlya resheniya differentsial’nikh uravneniy drobNogo poryadka [Modification of the PARAREAL Algorithm for Solution of Fractional Differential Equations]. Trudy MezhdunarodNoy NauchNoy Konferentsii "Parallel’nye Vichislitel’nye TekhNologii (PaVT’2010)" [Proc. Int. Conf. "Parallel Computational TechNologies (PCT’2010)"]. Chelyabinsk, YUrGU Publ., 2010, P. 519-524.
8. Goloviznin V.M., Kisilev V.P., Korotkin I.A., Yurkov Yu.I. Pryamye zadachi neklassicheskogo pereNosa radionuklidov v geologicheskikh formatsiyakh [Direct Problems of Nonclassical Radionuclide Transfer in Geological Formations]. Izv. Ross. Akad. Nauk, Ser. Energetika, 2004, No. 4, P. 121-132.
9. Zav’yalov Yu.S., Kvasov B.I., Miroshnichenko V.L. Metody splain-funktsiy [Methods of Spline-Functions]. Moscow, Nauka, 1980.
10. Lions J.-L., Maday Y., Turinici G. Resolution d’edp par un Schema en Temps Parareal. C.R. Acad Sci. Paris. Ser. I Math., 2001, Vol. 332, P. 661-668.
11. Maday Y., Turinici G. The Parareal in Time Iterative Solver: a Further Direction to Parallel Implementation Lecture Notes in Computational Science and Engineering: Domain Decomposition Methods in Science and Engineering XVIII, 2005, Vol. 40, P. 441-448.
12. Staff G.A., Ronquist E.M. Stability of the Parareal Algorithm Lecture Notes in Computational Science and Engineering: Domain Decomposition Methods in Science and Engineering XVIII, 2005, Vol. 40, P. 449-456.
13. Gander M.J., Vandewalle S. On the Superlinear and Linear Convergence of the Parareal Algorithm Lecture Notes in Computational Science and Engineering: Domain Decomposition Methods in Science and Engineering XVI, 2007, Vol. 55, P. 291-298.
14. Gander M.J., Vandewalle S. Analysis of the Parareal Time-Parallel Time-Integration Method SIAM J. Sci. Comput., 2007, Vol. 29, No 2, P. 556-578.
15. Fischer P.F., Hecht F., Maday Y. A Parareal in Time Semi-implicit AProximation of the Navier-Stokes Equations Lecture Notes in Computational Science and Engineering:
Domain Decomposition Methods in Science and Engineering XVIII, 2005, Vol. 40, P. 433-440.
16. Bal G., Wu Q. Symplectic Parareal Lecture Notes in Computational Science and Engineering: Domain Decomposition Methods in Science and Engineering XVII, 2008, Vol. 60, P. 401-408.
Поступила в редакцию 4 ноября 2012 г.
