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

Вычислительные технологии

Том 8, № 3, 2003

ПАРАЛЛЕЛЬНЫЙ АЛГОРИТМ ДЛЯ КОМПАКТНЫХ СХЕМ В НЕОДНОРОДНЫХ

ОБЛАСТЯХ*

В. И. ПААСОНЕН

Институт вычислительных технологий СО РАН, Новосибирск, Россия

e-mail: paas@ict.nsc.ru

The method for boundary value problems in inhomogeneous domains is investigated. The compact schemes are applied inside subareas and the multipoint one-dimensional relations are used on the internal and external boundaries. The parallel algorithm is developed for the solution of the problem.

Введение

Известны два способа аппроксимации внешних и внутренних граничных условий. Первый заключается в аппроксимации закона сохранения в балансной ячейке, составленной из однородных частей соседних ячеек сетки [1]. В случае методов высокого порядка точности результирующее соотношение зависит не только от дифференциального уравнения, но и от его следствий и является достаточно сложным [2, 3]. Такой подход всегда приводит к неодномерным граничным условиям, особенно сложным в случае компактных схем, что создает дополнительные проблемы при расщеплении задачи. По этим причинам есть основание применять более универсальный и технологичный алгоритм, способный сочетаться с использованием схем повышенной точности.

Второй способ состоит в аппроксимации именно граничного условия [4, 5] без привлечения самого уравнения. В этом случае граничные условия по построению оказываются универсальными и одномерными, что позволяет единообразно применять их при различных уравнениях и избегать проблем при расщеплении на одномерные задачи. Но с повышением порядка аппроксимации схемы приходится увеличивать число точек для граничного соотношения (длину граничного условия), что приводит к локальным нарушениям диагонального преобладания матриц, подлежащих обращению. Однако численные эксперименты с "длинными" граничными условиями [4, 5] подтверждают применимость алгоритмов такого рода, и это делает актуальным их обоснование.

В настоящей работе предлагается модификация известного алгоритма распараллеливания прогонки [6] применительно к описанному выше сочетанию компактных аппроксимаций уравнения в областях гладкости решения и многоточечных одномерных внутренних и внешних граничных условий. Рассматриваемый параллельный алгоритм осуществляет

* Работа поддержана Министерством образования Российской Федерации по фундаментальным исследованиям, грант Е 02-1.0-25.

© В. И. Паасонен, 2003.

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

1. Постановка одномерной задачи

Аппроксимируем первую производную на равномерной сетке с шагом к, применяя одностороннюю разделенную разность на произвольном множестве точек

Д 1 ^ Т к Т = к ак Т" •

к=О

Здесь — Т, оператор сдвига (Т^и(х) = и(х + к)), а ак — неопределенные коэффициенты разностного оператора. Если они удовлетворяют системе линейных уравнений

^ак ке = 5ц, I = 0

к=О

где 5и — символ Кронекера, то погрешность на достаточно гладких функциях составляет величину 0(кь). Ввиду линейной независимости строк матрицы системы максимально возможный порядок аппроксимации равен в. При этом решение системы вычисляется явно (оно — частный случай более общей системы [9])

(_1)к+1 8

ак =-к-Ск, к = 0, ао = ак •

к= 1

Можно показать, что нашу многоточечную разность "вперед", аппроксимирующую оператор дифференцирования, при Ь = в можно представить в следующих эквивалентных формах:

Д1 = )к+1Ск Тк^ _ Е = ^ (-к)к 1 ( ТН _ Е4 к

к ) С кк ^ к V к

к=1 к=1

Аналогично выглядит разность "назад" порядка в. При тех же выражениях для коэффициентов ак она имеет вид

1Г = Ь "к ^ = _ к ± ак

к=0 к=0

Пусть х0, х1, •••, хг — возрастающая последовательность значений пространственной координаты х, представляющих собой координаты границ однородных слоев х^_1 < х < х^-, (] = 1, •••, г) с постоянными в пределах слоя теплофизическими характеристиками с = с^-, Л = Л^-. Внутри каждого слоя решение и(х,£) удовлетворяет уравнению теплопроводности на интервале времени 0 < £ < Т (или его стационарному варианту)

ди д2и „ .

СЖ = Л"дж2 _ /(х'£)-

На внешних границах х = х0, х = хг поставлены условия первого, второго или третьего рода (объединенные общей формулой)

ди тт , Рз дП + н и = Фз' 3 = 0'г'

где п — внешняя нормаль, а на внутренних границах х = Х1, ...,х = хг-\ заданы условия четвертого рода

ди ди

Хдх) + — [Лдх)_ = 9(х'г)' 3 = 1'...'г— 1

где индексы " — "и " + "означают пределы слева и справа. При нулевой правой части имеет место равенство потоков справа и слева. В противном случае имеется поверхностный источник или сток тепла в зависимости от знака правой части. Это может быть, например, трение между слоями или неидеальный тепловой контакт (щель).

Уравнение теплопроводности аппроксимируем с погрешностью 0(т2 + к4) схемой

ип+1 — ип -2

С л--— = Л Лип — БР+1/2, А = Е + ак2 Л, Б = Е + — Л'

т ' ' 12 '

где Е — тождественный оператор, Л — обычная аппроксимация двойного дифференцирования; параметр осреднения выражается формулой

1 а Лт

~ 12 2' а ск2.

Таким образом, в узле с номером г внутри любого однородного слоя разностное уравнение приводится к трехточечному (временной индекс опущен)

а од— + (1 — 2а) од + аод+1 = Г*.

Здесь Г — результат операций над значениями функции и на нижнем слое и правой части f, определяемый по формуле

Г = 1 (Лип — Sfп

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

В узлах г, совпадающих с границами раздела сред, разностное уравнение запишем с помощью 5-точечных односторонних разностей

аи— + а и+ = 9*.

- з=о + з=о

На внешних границах также используем 5-точечные аналоги первых производных — у аз-з — ^о-о = -фо'^- У азим— — ^-м = —фг.

к1 ъ кг ъ

3=0 3=0

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

Можно рассматривать также уравнение с переменными коэффициентами. В этом случае мы приходим к разностной схеме с переменным коэффициентом а в трехточечном соотношении. Можно допустить также наличие в исходном уравнении членов с первой производной. Это привело бы к более общему асимметричному трехточечному уравнению вида

аг'Шг-1 + ЪгЧ^г + С^+1 =

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

2. Многомерные разностные краевые задачи

В двумерном случае рассматрим краевую задачу для уравнения теплопроводности в области, составленной из подобластей однородности материалов. Каждая из них в свою очередь составлена из прямоугольников со сторонами, параллельными координатным осям. При кусочно-постоянных коэффициентах факторизованное разностное уравнение, аппроксимирующее с погрешностью 0(т2 + к4) уравнение теплопроводности во внутренних узлах подобластей, практически не отличается от классической схемы [10] и имеет вид

иП+1 — ип к2 к2

С А А2-т-= Л Пип - 5/п+1/2, Ак = Е + ак Н2кКк, 5 = Е + ^Л1 + ^Л2,

где

П = (Е + Ц Л2)Л1 + (Е + Ц Л1 )Л2,

Л1, Л2 — аналоги операторов двойного дифференцирования, а параметры осреднения ак имеют вид

1 ак Лт , ак = 12 - у, ак = Скк, к = 1,2

На целых шагах внешних границ поставлены разностные условия первого, второго или третьего рода, а на границах раздела сред — условия четвертого рода. Они имеют порядок в и построены с использованием односторонних разностных операторов, описанных выше. Будем обозначать их Д^к, Д—, где к — номер координатного направления. Для всех трехточечных уравнений расщепленной системы

А1ип+1/2 = Гп, = (ЛПип - Б/п+1/2)1,

ип+1 — ип

А2и,п+1 = ип+1/2, ^п+1 и и

на каждом временном шаге обычным путем формируются одномерные многоточечные граничные условия. Так, обозначив аналог оператора разности потоков справа и слева через

Д-

К\ = А+—1 — А-——1

П1 П1

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

К1ип+1/2 = Е1А2 --—— = А2Е1 --—— = А2 9--.

т т т

Аналогично определяются условия для промежуточной величины ига+1/2 на вертикальных участках внешних границ и для на горизонтальных участках границ. Очевидно, после расщепления задачи и постановки граничных условий на каждом дробном шаге система распадается на одномерные разностные задачи рассмотренного выше типа с почти всюду трехдиагональной структурой.

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

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

^0 - А- + А+

в = Ца' ^ = Г а' ^ = к а •

Первый фрагмент соответствует окрестности левой границы многослойного пакета:

^0 + во в1 в2 вз в4 0 0 ••• —Фо

а1 Ь1 С1 0 0 0 0 ••• Fl

0 0,2 Ь2 С2 0 0 0 ••• F2

0 0 аз Ьз Сз 0 0 ••• Fз

0 0 0 а4 Ь4 С4 0 ••• F4

Второй фрагмент иллюстрирует структуру матрицы в окрестности границы раздела сред:

аг-4 Ъг-4 сг-4 0 0 ^¿-4

0 аг-3 Ъг-3 сг-3 0 0 0 0 0 0 Ъ-3

0 аг-2 Ъг-2 сг-2 0 0 0 0 0 Ъ-2

0 0 аг-1 Ъг-1 Сг-1 0 0 0 0 1

0 7- 7- 7- 7- 7- + 7о+ 7+ 72+ 7+ 74+ 0 дг

0 0 0 0 аг+1 Ъг+1 Сг+1 0 0

0 0 0 0 0 аг+2 Ъг+2 сг+2 0

0 0 0 0 0 0 аг+3 Ъг+3 сг+3 0

0 0 аг+4 Ъг+4 сг+4 ^г+4

Можно считать, что узлов в слое достаточно, чтобы многоточечные односторонние разности не выходили за пределы слоя. Это предположение и необходимо, и естественно, так как в противном случае нарушалась бы аппроксимация ввиду нарушения гладкости на границах раздела сред. Если разбить каждый слой не менее чем на в шагов, то "длинные" условия на соседних границах не будут "переплетаться" между собой, что в принципе позволяет с помощью стандартных операций метода Гаусса понизить их порядок до минимума, вычитая из "длинных" строк комбинации соседних трехточечных строк. В результате таких операций над длинными строками матрица системы превращается в трехдиагональную. Однако мы будем рассматривать иной путь, основанный на параллельных вычислениях в слоях.

3. Система уравнений и параллельный алгоритм

Итак, предположим, что решается система уравнений вида

агиг-1 + Ъгиг + Сгиг+1 = 0 < % < N, % = %к, к = 1, ..., Г — 1, с внешними граничными условиями при % = ¿о = 0 и % = %т = N

^оИо

3=о

фо, ^тПм — аии — = Фг

3=0

и с внутренними граничными условиями

Ак ^^ \

— > Щ-3 + -- N и+ = дг

ТТ Пк+1 —

% = %к к = 1,..., г — 1.

3=о

3=о

Описываемый ниже алгоритм является обобщением и модификацией [6] и отличается от него неоднородностью области и наличием "длинных" граничных условий. В оригинальном алгоритме пространственный промежуток делится на произвольные непересекающиеся отрезки, а в модифицированном требуется, чтобы все границы раздела сред пакета принадлежали множеству точек разбиения области на части. Для простоты разобьем слоистый пакет на части только по границам раздела сред. Следуя [6, 7], предположим, что значения решения в узлах разбиения и на внешних границах уже известны. Обозначим их 'к, к = 0,..., г. Введем в каждом слое Шк = {хк—1 < х < Хк} локальную нумерацию узлов (т = 0,..., М, М = тк), обозначим вектор решения и в слое шк через ук и представим его в виде специальной линейной комбинации

ук = рк Шк-1 + дк Шк + гк,

Л

Л'

т

к

где г к — решение неоднородного уравнения с нулевыми граничными условиями

атгт- 1 + Ьтгт + Стгт+1 ^к +т) г0 гМ °>

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

атРт-1 + ЬтРт + Стрт+1 = Рк = Рм =

&тЯт-1 + ЬтЯт + СтЯт+1 = 0 Як = °> = 1

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

к к

неоднородному уравнению с граничными условиями ук = мк-1 и уМ = мк.

Используя затем граничные условия на внешних и внутренних границах, получим уравнения для независимого вычисления значений Мк. Так, подставляя выражение решения у1 , соответствующее первому слою, в условие на левой границе, получим

^owo - ®j(p]wo + qjjW\ + z1) = фо.

J1

j=0

Очевидно, это равенство представляет собой линейное уравнение для двух неизвестных w0 и w1, коэффициенты которого легко устанавливаются приведением в нем подобных слагаемых.

Из правого граничного условия аналогично устанавливается связь

s

^r Wr - J- aj (Pj Wr-1 + qjj WUr + zj) = фг j=0

для двух неизвестных wr-1 и wr, а из соотношения на k-й границе раздела сред следует равенство

Л s Л s

ТТ У] aj(Pkmk-jwk-i + qtk-jwk + zk -j) + V aj(pk+1wk + qk+1wfc+i + zk+1) = g%k, hk 7=o hk+1 T=o

связывающее значения искомого решения wk-1, wk и wk+1 на трех последовательно расположенных границах слоев.

Таким образом, на втором этапе совокупность значений wk на всех границах может быть определена из системы уравнений с трехдиагональной матрицей

Akwk-1 + Bkwk + Ckwk+1 = Gk, k = 1,..., r - 1,

Bowo + Cow1 = £o, Ar wr-1 + Br wr = Cr, где s s s

Ak = TkX ajpik-j, Bk = TkX ajqk-j + ^ ^ j

k j=o k j=o k+ j=o

С _ ^к+1 а „к+1 г _ п Ак ^ а к Ак+1 ^ а -к+1

к+1 к к+ з=о

8 8 8 Во _ ^о- т° X] ар1' Со _ - т X а 3 Со _ фо + т° X а 3

1 3=о 1 з=о 1 з=о

8 8 8 V- X Л г V- X Л г V- X Л г

Аг _ — I / ртг -з ' Вг _ — I / -3 ' _ + 1 ' — '

I иг ¡иг ¡иг

3=о 3=о з=о

У этой системы уравнений имеется существенное отличие от задач первого этапа, которое заключается в постоянстве ее размерности, не зависящей от степени детальности сетки. Размерность, очевидно, равна общему числу границ (г + 1) в слоистом пакете. В связи с этим нет необходимости исследования системы уравнений второго этапа на устойчивость вычислений при стремлении к нулю шага сетки.

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

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

При решении многомерных задач рассматриваемого класса изложенный алгоритм представляет практический интерес. От распараллеливания можно уверенно ожидать эффективности, если не дробить задачу на одномерные, а остановиться на задачах в подобластях подобно тому, как это сделано в работе [8]. Например, в плоском случае архитектура такого разбиения представляет собой двумерную сеть прямоугольников согласованных размеров. При этом вдоль "линеек" на плоскости (или "брусков" в пространстве), составленных из подобластей, нанизанных одна на другую вдоль координатного направления, возникают задачи рассмотренного выше типа. В отличие от алгоритмов [6, 7] в нашем случае соседние элементы могут быть как из одинаковых, так и из различных материалов. В первом случае на границе возникают трехточечные соотношения, как внутри области однородности, а во втором — многоточечные условия непрерывности потока заданного порядка точности. В первом случае на плоскости контакта подобластей получаются связи, как в приведенных работах, во втором — как описано здесь. Те и другие связи вместе порождают совокупность трехдиагональных систем уравнений малой размерности, не зависящей от Н, для определения решения задачи второго этапа в многомерном случае.

Список литературы

[1] Паасонен В. И. Моделирование тепловых процессов в неоднородных конструкциях с источниками тепла // Изв. АН СССР. Сер. техн. наук. Новосибирск. 1980. № 1. С. 108-112.

[2] Коскин П. И. Схема повышенной точности для уравнения теплопроводности с разрывными коэффициентами // Численные методы механики сплошной среды: Сб. науч. тр. / АН СССР. Сиб. отд-ние. ВЦ; ИТПМ. 1979. Т. 10, № 2. С. 85-96.

[3] Ильин В. П. Балансные аппроксимации повышенной точности для уравнения Пуассона // Сиб. мат. журн. 1996. Т. 37, № 1. 151-169.

[4] ВАлиуллин А.Н., Сафин Р.И., Паасонен В.И. О схеме расщепления с повышенным порядком аппроксимации краевых задач для уравнения Пуассона // Численные методы механики сплошной среды: Сб. науч. тр. / АН СССР. Сиб. отд-ние. ВЦ; ИТПМ. 1972. Т. 3, № 1. С. 17-25.

[5] Ждан С. А., Паасонен В. И. Схема повышенной точности для задачи о непотенциальном течении идеальной жидкости в плоских каналах / / Численные методы механики сплошной среды: Сб. науч. тр. / АН СССР. Сиб. отд-ние. ВЦ; ИТПМ. 1972. Т. 3, № 5. С. 35-39.

[6] Яненко Н.Н., Коновалов А. Н., Бугров А.Н., Шустов Г. В. Об организации параллельных вычислений и распараллеливании прогонки / / Численные методы механики сплошной среды: Сб. науч. тр. / АН СССР. Сиб. отд-ние. ВЦ; ИТПМ. 1978. Т. 9, № 7. С. 136-139.

[7] Вшивков В.А. О распараллеливании вычислительных алгоритмов // Сиб. школа-семинар по параллельным вычислениям. Томск: Изд-во ТГУ, 2002. C. 46-59.

[8] Ильин В. П. Параллельные неявные методы переменных направлений // Журн. вы-числ. математики и мат. физ. 1997. Т. 37, № 8. С. 899-907.

[9] Паасонен В.И. Абсолютно устойчивые схемы повышенной точности для систем гиперболического типа // Численные методы механики сплошной среды: Сб. науч. тр. / АН СССР. Сиб. отд-ние. ВЦ; ИТПМ. 1972. Т. 3, № 3. С. 72-79.

[10] Микеладзе Ш. Е. О численном интегрировании уравнений эллиптического и параболического типов // Изв. АН СССР. Сер. мат. 1941. Т. 5, № 1. С. 85-96.

[11] Самарский А. А. Схемы повышенного порядка точности для многомерного уравнения теплопроводности // Журн. вычисл. математики и мат. физ. 1997. Т. 3, № 5. С. 812-840.

Поступила в редакцию 31 января 2003 г.