CC BY
63
УДК 519.642
ПАРАЛЛЕЛЬНЫЙ АЛГОРИТМ РЕШЕНИЯ ДРОБНО-ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ ПЕРЕНОСА НА ОСНОВЕ МОДИФИЦИРОВАННОГО МЕТОДА ШВАРЦА
С.Ю. Лукащук
A PARALLEL ALGORITHM BASED ON AN EXTENDED SCHWARZ DOMAIN DECOMPOSITION METHOD FOR THE SOLUTION OF FRACTIONAL EVOLUTION EQUATIONS
S. Yu. Lukashchuk
Предлагается модификация метода декомпозиции Шварца для дифференциальных уравнений переноса, содержащих производные дробного порядка. Доказывается сходимость метода. Приводятся описание численной схемы и схемы распараллеливания, а также оценка эффективности предлагаемого параллельного алгоритма.
Ключевые слова: параллельный алгоритм, метод декомпозиции Шварца, производная дробного порядка, дробно-дифференциальные уравнения переноса.
An extension of the Schwarz domain decomposition method to evolution equations with fractional time derivatives is considered and its convergence is proved. The numerical scheme and scheme of parallelization for the proposed algorithm are described. An efficiency estimation of the algorithm is given.
Keywords: parallel algorithm, Schwarz domain decomposition method, fractional derivative, fractional evolution equation.
Введение
Методы декомпозиции расчетной области являются основой для построения большинства современных параллельных алгоритмов, нацеленных на решение практических задач в геометрически сложных областях. Обширный класс таких методов образуют методы Шварца, относящиеся к группе итерационных методов подструктур [1]. В методах Шварца решение задачи строится независимо в пересекающихся или непересекающихся подобластях, а для согласования этих решений организуется итерационный процесс. При этом после каждой итерации происходит обмен внутренними граничными условиями между соседними подобластями.
В классическом методе Шварца на внутренних границах выставляются граничные условия первого рода (условия Дирихле). Недостатком такого подхода является низкая скорость сходимости итерационного процесса, а в случае неперекрывающихся областей сходимость может вообще отсутствовать [2]. Эти недостатки преодолены в оптимизированном варианте метода [3]. При этом граничные условия Дирихле заменены на искусственные граничные условия, содержащие линейные дифференциальные операторы специального вида.
Первоначально метод Шварца был разработан для решения уравнений эллиптического типа, однако впоследствии он был распространен и на уравнения гиперболического и параболического типов. В работе [4] показывается, что использование оптимизированного метода Шварца для решения классического уравнения диффузии приводит к тому, что дополнительный оператор, входящий во внутренние граничные условия, становится нелокальным по времени. Для преодоления этого недостатка в [4] предложены различные способы приближения нелокального оператора локальным, сохраняющим сходимость итерационного процесса.
В данной работе показывается, как метод Шварца может быть адаптирован для решения уравнений переноса, содержащих частные производные дробного порядка по времени [5]. В настоящее время наиболее хорошо исследованным видом таких уравнений являются дробно-дифференциальные уравнения аномальной диффузии, позволяющие описывать как аномально медленные (субдиффузия), так и аномально быстрые (супердиффузия) процессы переноса. Такие процессы наиболее часто наблюдаются в неоднородных сложных средах: пористых и трещиновато-пористых средах, перколяционных кластерах и самоподобных фрактальных структурах, турбулентных потоках жидкости, газа и плазмы и многих других [6]. Как правило, аномальный перенос обусловлен эффектами памяти, пространственной нелокальности и перемежаемости.
Следует отметить, что поскольку производная дробного порядка является нелокальным интегро-дифференциальным оператором, численное моделирование по дробнодифференциальным моделям требует весьма значительных вычислительных ресурсов. При этом, если уравнение содержит дробную производную по времени, то объем вычислений будет линейно расти с ростом номера временного слоя. Указанная проблема является одним из факторов, сдерживающих внедрение дробно-дифференциальных моделей в практику современных инженерных расчетов. Поэтому разработка параллельных алгоритмов решения этих задач представляется весьма актуальной.
1. Модификация метода Шварца для уравнения диффузии дробного порядка
Рассмотрим задачу Коши для уравнения диффузии дробного порядка
где О = (—оо, с»), функция у ограничена в О, X? — коэффициент аномальной диффузии,
- левосторонняя частная производная дробного порядка а по і типа Римана - Лиувилля [5].
Разобьем область П на две, возможно неперекрывающиеся, подобласти Пі = (—со, Ь\ и 0,2 = [0, оо), Ь > 0 (см. рисунок). Обозначим приближение к искомой функции у(х,і) в подобласти Лі через и (х, і), а в области О2 _ через V (х,і). Запишем задачу (1) отдельно
(1)
У (х, 0) = д (х), жеП,
а Є (0,1)
для каждой подобласти:
ди (ж, і)
£>-
д2и (ж, і)
дг дга V дх2
и(х,0)=д(х), х Є Гіі,
Эу (х, і) _ да ґд2у (х, і)
дЬ дЬа \ дх2
V (х, 0) = д (х), х Є Г22.
+ і > 0, х Є а Є (0,1);
+ /(х,і), і > 0, х Є 02, а Є (0,1);
(2)
шшшшштшшт,;
о
п1г у(х}0 = и(х,і)
&2, Хх,0 = '{х,д
жшттжжь
Схема декомпозиции области
Решение задачи (1) может быть получено композицией решений «и». Следуя [4], по-
ставим на внутренних границах областей х = 0 и х — Ь следующие граничные условия:
{В + Аі)и(х,і)\х=ь = (В + А і)у{х,і)\х=ь,
(В + Л2) V (ж, і) |я=0 = {В + Л2) и (ж, і) |я=0 ,
(3)
где оператор В имеет вид, соответствующий аномальному диффузионному потоку: В = да ( д \
В-— I —— , а Ах и Л2 - некоторые линейные операторы, действующие по переменной £ и
о1° \дх ) подлежащие определению.
Для системы (2), (3) можно записать следующий итерационный процесс:
*>°- *<*■
х=Ь ]
+ / (х,і), і > 0, х > 0,
(4)
дип (х, *) дЬ
ип (ж, 0) = д (х), х < Ь,
{В + А1)ип(х,г)\х=ь= {В + {х,Щх=1, 1>0;
дьп (х, г) _ да / д2ьп (х, 4) дЬ д1а \ дх2
Vй (х, 0) = д (ж), х > 0,
(£ + Л2)ип(М)|*=о = (В + А2)ип-1(х^)\х=0.
Операторы Ах и Л2 должны выбираться так, чтобы гарантировать сходимость итерационного процесса (4).
Обозначим через <рп (ж, I) = и (х, £) — ип (ж, Ь) и фп (х, 4) = V (ж, £) — Vп (ж, £) соответствующие ошибки на п-й итерации. В силу линейности всех операторов, для определения ошибок из (2) - (4) получаем следующий итерационный процесс:
дфп{х^) да /д2<рп(х,г)'
і > 0, ж < £,
Ы дЬ® ^ дх2
<рп (ж,0) =0, ж < Ь,
(В + Лі) у?” {х,щх=ь = (В + Ах) фп~1 {х,і)\х=ь , г > 0; дфп (ж, і) да ( д2фп (ж, і)
= £>-
(5)
і > 0, ж > 0,
ді діа \ дх2
фп (ж, 0) = 0, ж > 0,
(В + Л2) фп (ж, і) |х=0 = (В + Л2) «р"-1 (ж, і)
!ж=0 '
Сходимость итерационного процесса (5) к нулю независимо от начального приближения гарантирует сходимость исходного процесса (4). Из этого условия и могут быть найдены операторы Лі и Л2.
Обозначим через у (ж, а) преобразование Лапласа от функции у (ж, і) по времени:
Следуя [4] будем полагать, что преобразование Лапласа от А.гу есть Аг (в) у, г = 1,2. Применяя преобразование Лапласа по времени к задаче (5), получаем:
с учетом (7) дает ф2 (х, я) = 0 (х < Ь) и ф2 (х, а) = 0 (ж > 0). Таким образом, условие (9) обеспечивает сходимость итерационного процесса (5) за две итерации независимо от начального приближения и глубины перекрытия Ь.
Применяя к (9) обратное преобразование Лапласа, находим искомые операторы:
Можно показать, что при разбиении исходной области на N подобластей такой способ выбора операторов гарантирует сходимость итерационного процесса метода Шварца за N шагов.
2. Численная схема алгоритма
Теперь рассмотрим некоторые особенности численной реализации описанного в предыдущем разделе метода декомпозиции для случая ограниченной области $7 = [— 1,1].
В качестве модельной будем рассматривать следующую первую начально-краевую задачу для уравнения диффузии дробного порядка:
ОУг К'*-") °) ^° тхх \ і ) > ^ ’
Бзафпх (Ь, а) + Аі_(в) фп (£, а) = ОзафГ1 (Ь, а) + Ах (а) фп~1 (Ь, а); еф71 (ж, в) = Пзаф™х (ж, в), ж > 0,
£>а°^ (0, а) + А2 (а) фп (0, а) - Яа0^"1 (0, а) + А2 (а) фп~1 (0, а).
вф4 (ж, а) = Бзафхх (ж, а), ж < І,
(6)
Разрешая уравнения (6) с учетом ограниченности функций фиф, находим
фп (ж, а) = С?еых, X <Ц фп (ж, а) = С%е~“х, ж > 0; ш = л/а1"4*!)-1.
(7)
Подставляя (7) в граничные условия из (6), приходим к соотношениям
фп+1 (0, з)=р (а) фп~1 (0, а), фп+1 (0, а) = р (а) фп~1 (0, а),
(8)
где
п(я\= (А1 (5) - Рзаш) (А2 (а) + Дааы) 2шЬ
(Ах (а) + £>ааш) (А2 (а) — Из01 ш)
- коэффициент, характеризующий сходимость итерационного процесса. При
(10)
і > 0, ж Є (—1,1), а Є (0,1);
(И)
У (ж, 0) = <? (ж), ж Є [-1,1]; у (-М) = щ (і), у(М) =Мг(<), і > о.
Для простоты будем полагать, что область О разделяется на две равные подобласти без перекрытия: Лі = [—1,0] и = [0,1] • Как и ранее, решение задачи в области Пі обозначим через и (ж, і), а в области Я2 - через V (ж, і). На внутренней границе подобластей ставим стыковочные граничные условия вида (3) с операторами (10).
Каждую подобласть разобьем равномерной конечно-разностной сеткой с шагом Дж и числом узлов М + 1. Будем использовать локальную нумерацию узлов в пределах каждой подобласти: ж і = іАх, і = 0,1,..., М. Шаг по времени Аі будем считать постоянным, временной слой обозначим через tj = ЗАі, j = 0,1,_____Значения сеточных функций обозначим
% (Хі, ^
Производную дробного порядка по времени аппроксимируем на основе формулы Грюн-вальда - Летникова [5]:
Аналогично (4), в каждой подобласти организуем итерационный вычислительный процесс. С учетом конечно-разностной аппроксимации, для ж € О1 имеем:
В качестве начального приближения используем значение с предыдущего временного слоя:
и^- = j = 1,2,___Отметим, что в (12) итерационный процесс организуется только
для временного слоя На всех предыдущих временных слоях процесс считается сошедшимся, поэтому у сеточных функций на этих слоях индекс номера итерации отсутствует.
Аналогично (12), записывается итерационный процесс в области П2 Для функции V (ж, I) .
Легко заметить, что (12) на каждом временном слое tj приводит к системе линейных алгебраических уравнений А1/™ = Щ с трехдиагональной матрицей
^г,0 9іі * 0,1, , М, ^0,^' 3 1,2,...,
(12)
векторами [/" =
А =
\
(1
—а 1
к=1
к=1
(14)
Данная система легко может быть решена методом прогонки.
Аналогичная система записывается в области 02 для вектора неизвестных У™ =
3. Параллельный алгоритм и оценка его эффективности
Теперь приведем формальное описание параллельного алгоритма, порождаемого численной схемой из предыдущего раздела.
Пусть имеется р вычислительных модулей (процессоров или ядер) одинаковой производительности. В этом случае расчетная область разбивается на р неперекрывающихся подобластей одинакового размера, и каждый вычислительный модуль отвечает за расчет своей подобласти. Так как внутреннее граничное условие содержит производную по пространству, каждый вычислитель на каждом временном слое со стороны каждой внутренней границы должен иметь две вспомогательные переменные для хранения значений, передаваемых с соседней подобласти (например, для первой подобласти - это значения г>д у1 и см.(14)). Важно отметить, что в отличие от классических уравнений переноса, при расчете аномальной диффузии необходимо сохранять значения вектора неизвестных на всех временных слоях. Алгоритм состоит из следующих шагов.
1. На основании исходных данных задачи и параметров сетки все вычислители параллельно рассчитывают коэффициенты матрицы А из (13) и неизменные составляющие прогоночных коэффициентов.
2. На каждом вычислителе организуется глобальный цикл по временным шагам. На нулевом шаге вектор искомых значений заполняется в соответствии с начальным усло-
б) На каждом вычислителе организуется внутренний итерационный процесс.
в) Вычислители, отвечающие за расчет соседних подобластей, обмениваются двумя приграничными значениями. При этом на первом итерационном шаге пересылаются значения с предыдущего временного слоя, а на последующих итерационных шагах - значения с предыдущей итерации.
г) На первом итерационном шаге каждый вычислитель по формулам вида (14) параллельно вычисляет вектор правой части системы. На последующих итерациях пересчитываются только первый и последний элементы вектора Н, соответствующие внутренним границам.
д) Каждый вычислитель параллельно решает свою систему методом прогонки.
е) Если номер итерации меньше р, повтор с шага б), иначе переход на новый временной слой и повтор с шага 3.
Проведем теоретическую оценку эффективности этого алгоритма. Пусть, для простоты, в задаче (11) / (х, I) = 0 и граничные условия постоянны. Если число узлов сетки в каждой подобласти равно М + 1, то расчет элементов вектора Н на ,7-ом временном слое на первом итерационном шаге требует для внутренней подобласти 2 + 1) (М + 3) ариф-
метических операций, а на последующих итерациях - только 8 операций, необходимых для пересчета первого и последнего элементов вектора, соответствующих внутренним границам.
вием.
3. Для произвольного у-го временного шага выполняются следующие действия,
а) Вычисляются новые коэффициенты и>^ и
Трудоемкость решения системы методом прогонки с учетом неизменности матрицы А составляет 5 (М + 1). Тогда трудоемкость расчета j-го временного слоя для р вычислителей за р итераций без учета затрат на передачу данных составит
Тр = 2 (2j + 1) (М + 3) + 8 (р - 1) + 5 (М + 1) р.
Трудоемкость реализации алгоритма на одном вычислителе (в этом случае внутренние итерации отсутствуют) составит
Ti = 2 (2j + 1) [(М + 1) р - 2] + 5 [(М + 1) р - 2] + 4.
В результате находим ускорение параллельного алгоритма:
Я =И - Р ~ 5с7^ ~ 2см ~ 6CJCM 1 1 /, гл
р Тр 1 + 5Cjp -f 2см + 8(р — 1)с^сд^ ’ J 2 (2j + 1) ’ М + 1
Как правило, М » р > 1, поэтому см ~ 0 и оценка (15) упрощается. Тогда для больших временных слоев алгоритм имеет асимптотически линейное ускорение: Sp —» р при j —> оо. Эффективность алгоритма = Sp/p в этом случае стремится к единице. Данные асимптотические оценки остаются верными и в случае учета затрат времени на передачу данных.
Предложенный в данной работе алгоритм легко обобщается на случай двумерных и трехмерных пространственных областей.
Статья рекомендована к публикации Программным комитетом Международной научной конференции ^Параллельные вычислительные технологии 2011».
Литература
1. Копысов, С.П. Объектно-ориентированный метод декомпозиции области / С.П. Копы-сов, И.В. Красноперов, В.Н. Рычков // Вычислительные методы и программирование.
- 2003. - Т. 4.- С. 1 - 18.
2. Gander, M.J. A Non-Overlapping Optimized Schwarz Method which Converges with Arbitrarily Weak Dependence on h / M.J. Gander, G.H. Golub // Proc. of the Fourteenth International Conference on Domain Decomposition Methods, Cocoyos, Mexico. - Cocoyos, 2003. - P. 281 - 288.
3. Gander, M. J. Optimized Schwarz methods / M.J. Gander, L. Halpern, F. Nataf // Proc. of the Twelfth International Conference on Domain Decomposition Methods, Chiba, Japan. -Chiba, 2001. - P. 15 - 28.
4. Gander, M. J. Optimal convergence for overlapping and nonoverlapping Schwarz waveform relaxation / M.J. Gander, L. Halpern, F. Nataf // Proc. of the Eleventh international Conference on Domain Decomposition Methods, Greenwich, Great Britain. - Greenwich, 1999. - P. 27 - 36.
5. Podlubny, I. Fractional Differential Equations /I. Podlubny. - San Diego: Academic press, 1999.
6. Metzler, R. The random walk’s guide to anomalous diffusion: A fractional dynamic approach / R. Metzler, J. Klafter // Phys. Rep. - 2000. - V. 339. - P. 1 - 77.
Лукащук Станислав Юрьевич, кандидат физико-математических наук, доцент, кафедра высокопроизводительных вычислительных технологий и систем, ГОУ ВПО «Уфимский государственный авиационный технический университет», [email protected].
Поступила в редакцию 2 марта 2011 г.
