Параллельные алгоритмы разделения области для параболических задач Текст научной статьи по специальности «Математика»

УДК 519.63

ПАРАЛЛЕЛЬНЫЕ АЛГОРИТМЫ РАЗДЕЛЕНИЯ ОБЛАСТИ ДЛЯ ПАРАБОЛИЧЕСКИХ ЗАДАЧ*)

П. Е, Захаров

1. Введение

Теория и практика методов декомпозиции области (domain decomposition, DD) для стационарных краевых задач подробно рассмотрены в [1-4]. Различают методы декомпозиции области с налеганием и без налегания подобластей. Приближенное решение нестационарной задачи может быть получено после стандартной неявной аппроксимации по времени и решения соответствующих дискретных уравнений на каждом временном слое одним из вариантов методов DD для стационарной задачи. Учитывая переходящее свойство нестационарной задачи (см., например, реализацию, основанную на методе Шварца [5,6]), можно построить оптимальный итеративный метод DD, где количество итераций не зависит от размера сетки времени и пространства.

Такая особенность нестационарных задач полностью учитывается в безытерационных схемах DD. В некоторых случаях возможно [7,8] использование только одной итерации альтернирующего метода Шварца для параболического уравнения второго порядка без потери точности приближенного решения. Безытерационные схемы DD связаны с особыми вариантами аддитивных схем (схем расщепления) — регионально-аддитивными схемами [9].

Схемы декомпозиции области для нестационарных задач могут быть выделены методом декомпозиции области, выбором оператора де-

Работа выполнена при поддержке Российского фонда фундаментальных исследований (проект № 13—01-00719А).

© 2013 Захаров П. Е.

композиции (обменом граничных условий) и схемой расщепления. Для дифференциальной задачи используют методы декомпозиции области

с налеганием подобластей (Пав = Па П Q в ф 0) или без налегания (Пав = 0) [2)4]. Методы без налегания подобластей связаны с явной формулировкой граничных условий в интерфейсных границах. Эти методы широко используются в задачах, где для каждой отдельной подобласти вводится своя собственная вычислительная сетка (триангуляция). В общем, также используются схемы DD с налеганием подобластей. Методы декомпозиции области с минимальным налеганием, где ширина налегания равна размеру сетки (Qa^ = #(h)), часто могут быть интерпретированы как методы без налегания, дополненные соответствующими граничными условиями на интерфейсах.

В данной работе рассматривается двукомпонентная схема декомпозиции области без налегания для приближенного решения краевой задачи для параболического уравнения.

2. Модельная параболическая задача

Рассмотрим модельную начально-краевую задачу для параболического уравнения второго порядка. В ограниченной области Q неизвестная функция w(x, t) удовлетворяет уравнению

где &(х) > к > 0, х € О. Уравнение (1) дополняется однородным граничным условием Дирихле

p

О = Qa, Qa = Qa U dQa, a = 1,2,... ,p,

u(x,t) = 0, xe 6П, 0<t<T.

(2)

Также дополнительно вводим начальное условие

u(x,0) = u0(x), xeû.

(3)

Таким образом, рассматривается нестационарная задача диффузии (1)-(3) для множества функций u(x, t), удовлетворяющих граничным условиям (2). Вместо (1) и (2) будем использовать дифференциально-операторную запись

du

— + = f(t), 0 < t < Т, (4)

где

—¿¿НУ-

a= 1 4 у

Рассматривается задача Коши для эволюционного уравнения (4) с начальным условием

u(0) = u0. (5)

Для множества функций (2) определим гильбертово пространство Ж = Jz?2(0) со скалярным произведением и нормой

M = /U(xMx)dX' М = (U'U)1 ^. n

В Ж оператор диффузионного переноса л/ самосопряжен и положительно определен:

«# = «#* > кМ, S = ¿(0) > О, (6)

где S — единичный оператор в Ж.

Для решения задачи (4)-(6) определим простейшую априорную оценку, которая будет отправной точкой для рассматриваемых сеточных задач. Самосопряженный положительно определенный оператор 'З может быть связан с гильбертовым пространством Ж® со скалярным произведением и нормой

1

(u,v)® = {@u,v), ||u||®= (u,u)^"\

u

\jt\W\^+\Wuf = {f,^u). (7)

Принимая во внимание неравенство

из (7) получаем

а

1

Используя лемму Гропуолла, приходим к желаемой оценке

+ ^ I \\т\\*м,

(8)

о

которая выражает устойчивость решения задачи (4)-(6) относительно начального значения и правой части.

Детальное изучение аппроксимаций по пространству и времени будем проводить на примере начально-краевой задачи в прямоугольной

Приближенное решение задается в узлах равномерной прямоугольной сетки области Л:

ш = {х | X = (Ж1,Ж2), Ха = гаЬа, ¿а = 0Д,...,Жа, МаЬ,у = 1а, а=1,2},

а через ш и дш обозначим множество внутренних и граничных узлов соответственно (ш = ш и дш). Для сеточных функций у(х) = 0, х € дш, определим гильбертово пространство Н = Ь2(ш) со скалярным произведением и нормой

3. Разностная аппроксимация

области

П = {х | х = (ж1, ж2), 0 < жа < 1а, а = 1, 2}.

Предполагая, что коэффициент &(х) в Л достаточно гладкий, используем дискретный оператор диффузии

Ау = --ргк(х 1 + /11/2, х2)(у(х 1 + Ъ,\, х2) - у(х 1,х2))

Ч

+ -т^Нхх - ^1/2, x2)(y(xl,x2) - у(х\ - Ь^,х2))

ч

— -Г^к(х1, х2 + И2/2)(у(х1,х2 + 1г2) - у(х 1, ж2)) щ

+ - h2/2)(y(xl, х2) - у(х\,х2 - h2)). (9)

щ

В пространстве И оператор А самосопряжен и положительно определен:

4 пп

А = А*^к(61 + 62)Е, 6а = — шп2 —а= 1,2. (10)

После аппроксимации по пространству из уравнений (1), (2) переходим к дифференциально-разностному уравнению

= хеш, о <КТ. (11)

аЬ

С учетом (3) дополним (11) начальным условием

у(х,0) = м°(х), хе ш. (12)

Для решения дифференциально-разностной задачи Коши (11), (12) выполняется следующая априорная оценка (см. (8)):

г

ЫГа^ИГаЦ / \\т\\2м. (13)

о

Теперь уделим внимание аппроксимации по времени. При построении схемы декомпозиции области для задачи (11), (12) как отправная точка используется стандартная двухслойная схема. Пусть т — шаг по времени и уп = у(Ьп), Ьп = ит, п = 0,1,..., N, Мт = Т. Уравнение

(11) аппроксимируем двухслойной схемой с весами уп уп

--— + А(ауп+1 + (1 -а)уп) = <рп, п = 0,1,..., N - 1, (14)

где = f(tn + ат), и дополним начальным условием

У° = и°. (15)

На достаточно гладких решениях задачи разностная схема (14), (15) имеет ошибку аппроксимации #(т2 + (а — 1/2)т + h2), где h2 = (h\ + hl )/2.

Теорема 3.1. Разностная схема (14), (15) безусловно устойчива при а ^ О.5, и для численного решения выполняется оценка

IIyn+1 IIb < IIЛ1Ь + ¿М2, n = о, i,...,лг-1, (16)

где

D = А + (а — 1/2)тА2. Доказательство. Перепишем разностную схему (14):

yn+l _ yn yn+l I yn

(E + (а — 1/2)т А)—-g— + А--= vn,

т

скалярно умножим на тА(уп+1 + уп). Поскольку а > 1/2, а значит, D ^ А, имеем

1|уп+1||Ь - IIЛ1Ь + l\\A{yn+1 + у")||2 = r((fin,A(yn+1 + у«)). Учитывая, что

+ у«)) < \\\А{у^ + уп)\\2 +

получаем требуемую оценку (16).

Априорная оценка (16) для решения задачи (14), (15) является дискретным аналогом априорной оценки (13) для решения дпфферен-

D А т

4. Разделение области

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

_

I*

Л

Рис. 1. Декомпозиция сетки: и = £ и 7.

Сетку и делим та множество прямоугольных подобластей которые соединены между собой множеством внутренних граничных узлов 7. Подобласти £ не связаны между собой и могут быть рассмотрены параллельно. Фрагмент сетки показан на рис. 1. Решение задачи при такой декомпозиции можно разделить на 3 этапа.

Этап 1. Вычисление условий на границах подобластей явной схемой

Решение исходной задачи с использованием только схем (17), (18) является условно устойчивым, поэтому дополнительно вводится этап коррекции в 7.

Этап 3. Коррекция условий на границе подобластей неявной схемой

(17)

(18)

Уп+1 =уп+1/2, хе С,

уп уп

»-у— + Ауп+1 = ч>п, х е 7- (19)

т

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

5. Исследование устойчивости и сходимости

Введем функцию принадлежности подобласти 7:

х(х> = {); хе ? <эд

10, х е С.

Схему (17) и (18) умножаем на х и (1 - х) соответственно и суммируем: у.п+1/2 _ у п

У--У- + (1 - х)Ауп+1/2 + ХАуп = <Рп. (21)

т

Аналогично сложим схемы (18) и (19):

уп уп

-У- + ХАуп+1 + (1 - х)Ауп+1/2 = (22)

т

Схема (21), (22) представляет собой регионально-аддитивную схему, которую можно интерпретировать как разностную схему Дугласа — Рекфорда

уп / уп

-У— + А1Уп+1/2 + А2уп = уп, (23)

т

уп уп /

----+ А2 (уп+1 - у") = 0, (24)

т

где

^^ Х1А, XI = X, = Х2А, хъ = 1 - X.

Схему Дугласа — Рекфорда (23), (24) записываем в более общем

виде как факторизованную схему:

уп уп

£>1 £>2 -— + Ауп = <рп, (25)

т

где

Ба = Е + атха А, (26)

с правой частью, определенной в виде ф" = ¡{<гЬп+1 + ( 1 — я)Ьп). Прямая подстановка подтверждает, что схема (25), (26) совпадает со схемой (23), (24) при а=1.

Сформулируем основную теорему устойчивости для факторизо-ванной схемы (25), (26) [10].

Теорема 5.1. Факторизованная регионально-аддитивная схема (25), (26) безусловно устойчива для а ^ 1/2, н для дискретного решения выполняется оценка

\\В2уп+1Щ < ||В2упЩ + т\\В-фп||А, п = 0,1,...,Ж — 1. (27)

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

Точность анализируется стандартным способом — через рассмотрение соответствующей задачи для погрешности

гп(х) = уп(х) — ип(х), х е и,

где ип(х) =и(я.,1п) — точное решение дифференциальной задачи (1)-(3). Из (25), (26) приходим к задаче для погрешности

гП+1 _ г"

ВгВ2-+ Агп = фп, (28)

т

= 0. (29)

Учитывая (27) из задачи (28), (29), получаем

п

\\В2гп+1 \\А ^тфк\\А, п = 0,1,...,М — 1. (30) к=0

Для погрешности аппроксимации имеем

ип ип

фп = - ВгВ2--Аип. (31)

т

Учитывая (26), из (31) выводим

Фп = Ф" + Ф",

ф? = срп - (Е+ (а - 1/2)тА)--А-

т

ип+1 _ ип

Ф2 = -ег2 т2 XI Ах2 А

т

Первый член погрешности стандартный для схем с весами, тогда как второй член получается из разделения на подобласти. Для достаточно гладких решений для задачи (1), (3) получаем

фп = 0(Н2 + т2 + (а - 1/2)т).

Рассмотрим член фп более подробно. Учитывая, что ||Ба|| > 1 при а > 0, имеем

Б-Фп\\А < ^т2

Х1АХ2А-

< а2т2М||(1- Х2)АХ21| А

А

т

< аЧ*МцАх*||а +**т*М||х2ах2||а = °*т*М||*2||а +°*т*М||*2||а-

|Б-1 Фп\\А = ^ (**т2 |Ь |Ы .

Таким образом,

\ Б А

Данные аргументы позволяют нам сформулировать следующее утверждение.

Теорема 5.2. Для погрешности факторизованной регионально-аддитивной разностной схемы (25), (26) при а ^ 1/2 для задачи (1)-(3) получаем оценку

||Б2гп+1 ||а < М(Ь2+т2 + (а - 1/2)т + ^^|Ь|Ы. (32)

Для схем декомпозиции области без налегания, рассматриваемых здесь, с дискретными эллиптическими операторами второго порядка (9) оценка (32) дает неравенство

||Б2гп+1 ||а < М(Н2 + т2 + (а - 1/2)т + ^т2Н-/2Н-/2). (33)

а / т

сматрнваемая факторизованная схема при Н = 1 соответствует схеме с весами, и погрешность аппроксимации равна #(Н2 + т2 + (а - 1/2)т).

6. Численные эксперименты

Численные эксперименты проводились для параболического уравнения (1), где

Цх) = 1, /(х,£)=0, хеП, 0<£ ^ Т = 0.05. (34)

Задача (1)-(3), (34) рассматривается в единичном квадрате = 12 = 1 с начальным условием

и0(х) = х е П. (35)

Решение задачи (1)-(3), (34), (35) записывается в виде

и(х.,Ь) = ехр(—2

■ ■Я

- яи ни^^^н

■ V А I НИ

..........

...........

..................Иг'Ж......

ШШЬШ......Й,: ...

.........

... ■■

Рис. 2. Разница точного и приближенного решения при разделении области на 9 частей.

Погрешность приближенного решения вычисляем по формуле

е(Г) = Цуп(х) — ип(х )Ц

в норме Н = Ь2(и). На рис. 2 изображена норма Цуп(х) — ип(х)Ц решения факторизованной регионально-аддитивной схемой, значение увеличивается от светлого до темного.

0.1

0.01

0.001

0.0001

1е-05

1е-06 1е

01

Рис. 3. Зависимость от шага по времени при и = 1. 0.01

0.001

«О 0.0001

1е-05

1е-06

0.0001 0.001 0.01 0.1

к

Рис. 4. Зависимость от шага по пространству при и = 1.

Численные эксперименты проводились на вычислительном кластере Северо-восточного федерального университета им. М. К. Аммосова. Для выявления асимптотической зависимости погрешности от параметров задачи вычисления проводились в большом диапазоне параметров. Все графики для погрешности е = е(Т) заданы в логарифмическом масштабе.

К= 1.0с • 00 Л к = 2.5е-01 Ь = 6.2е-02 Л =1.бс-02

¡-07

1е-06

1е-05 0.0001

0.001

К= 1.0е+00 Ь = 2.5е-01 и = б.2е-02 Л = 1.6е-02

Рис. 5. Зависимость от размера подобластей при a = 1.

0.1 0.01 0.001 ю 0.0001 1е-05 1е-06

1е"°1е1о7 1^06 1е-05 0.0001 0.001 0.01

Рис. 6. Зависимость от шага по времени при a = 1/2.

Для факторизованной схемы первого порядка (a = 1) при увеличении шага по времени сперва преобладает член погрешности ^(т), потом начинает преобладать ^(т2) для h > 1 (рис. 3). При уменьшении шага по пространству наблюдается увеличение погрешности на конечный момент времени (рис. 4). Погрешность с определенного момента увеличивается при уменьшении размера подобласти (рис. 5).

Рис. 7. Зависимость от шага по пространству при a = 1/2.

Рис. 8. Зависимость от размера подобластей при a = 1/2.

Для факторизованной схемы второго порядка (a = 1/2) при уве-тт теоретической оценке (рис. 6). На рис. 7 легко увидеть член погрешно-h

ретическую оценку (33).

Для исследования зависимости погрешности от временного шага

т варьировали от 2_9 (26 временных слоев) до 2-20 ( 52 4 29 временных слоев) при h\ = h2 = 2-8 (65536 узлов в сетке) (см. рис. 3 и 6). Исследование зависимости погрешности от шага по пространству проводилось для hi и h2, равных от 2-4 (256 узлов в сетке) до 2-11 (4194304 узлов в сетке) при т = 2-13 (410 временных слоев) (см. рис. 4 и 7). Для исследования зависимости погрешности от размера подобласти использовалось до 1024 вычислительных ядер при hi = h2 = 2-8 (см. рис. 5 и 8).

ЛИТЕРАТУРА

1. Matbew Т. Domain decomposition methods for the numerical solution of parabolic and elliptic differential equations. Berlin: Springer-Verl, 2008.

2. Quarteroni A., VaHi A. Domain decomposition methods for partial differential equations. Oxford: Clarendon Press, f999.

3. Самарский А. А., Гулин А. В. Устойчивость разностных схем. M.: Наука, f973.

4. Toselli A., Widlund О. Domain decomposition methods — algorithm and theory. Berlin: Springer-Verl., 2005.

5. Cai X.-C. Additive Schwarz algorithms for parabolic convection-diffusion equations 11 Numer. Math. 1991. V. 60, N 1. P. 41-61.

6. Cai X.-C. Multiplicative Schwarz methods for parabolic problems // SIAM J. Sci. Comput. 1994. V. 15, N 3. P. 587-603.

7. Kuznetsov Yu. A. New algoritm for approximate realization of implicit difference schemes 11 Sov. J. Number. Methods (Eng.) 2009. V. 25, N 7. P. 810-826.

8. Kuznetsov Yu. A. Overlapping domain decomposition methods for fe-problems with elliptic singular perturbed operators // Proc. 4th Int. Symp. Domain decomposition methods for partial differential equations. Moscow, 1990. P. 223241.

9. Samarskii A. A., Mat us P. P., Vabishchevich P. N. Difference schemes with operator factors. Mathematics and its Applications. Dordrecht: Kluwer Acad. Publ., 2002.

10. Вабищевич П. H. Аддитивные оиераторно-разностные схемы (схемы расщепления). M.: URSS, 2013.

г. Якутск

6 августа 2013 г.