ЖУРНАЛ ВЫЧИСЛИТЕЛЬНОЙ МАТЕМАТИКИ И МАТЕМАТИЧЕСКОЙ ФИЗИКИ, 2015, том 55, № 7, с. 1168-1182
УДК 519.6
МНОГОСЕТОЧНЫЙ МЕТОД ДЛЯ ЭЛЛИПТИЧЕСКИХ УРАВНЕНИЙ С АНИЗОТРОПНЫМИ РАЗРЫВНЫМИ КОЭФФИЦИЕНТАМИ1
© 2015 г. В. Т. Жуков, Н. Д. Новикова, О. Б. Феодоритова
(125047Москва, Миусская пл., 4, ИПМРАН) e-mail: zhukov@kiam.ru, nn@kiam.ru, feodor@kiam.ru Поступила в редакцию 03.09.2014 г.
Построен алгоритм многосеточного метода Р. П. Федоренко для разностных эллиптических уравнений. Алгоритм ориентирован на решение трехмерных краевых задач для уравнений с анизотропными разрывными коэффициентами на параллельных компьютерах. Приведены результаты расчетов, подтверждающие работоспособность и параллельную эффективность многосеточного алгоритма, что обеспечивается использованием в качестве многосеточной триады стандартного итерационного метода Чебышёва при решении грубосеточных уравнений, сглаживающих явно-итерационных процедур чебышёвского типа и операторов межсеточных переходов в проблемно-зависимой форме. Библ. 16. Фиг. 3. Табл. 4.
Ключевые слова: трехмерные эллиптические уравнения, анизотропные разрывные коэффициенты, многосеточный метод, итерационный метод Чебышёва, параллельная реализация.
DOI: 10.7868/S0044466915070133
ВВЕДЕНИЕ
В работе изложены принципы и особенности реализации итерационного многосеточного метода Р.П. Федоренко (см. [1]—[3]) решения разностных уравнений, возникающих при аппроксимации трехмерных уравнений эллиптического типа на декартовых сетках. Известно, что для эллиптических уравнений многосеточный метод оптимален, т.е. для нахождения приближенного решения с заданной точностью объем вычислительных затрат прямо пропорционален числу узлов сетки. Это свойство весьма важно, так как многие практические задачи требуют использования сеток с большим числом узлов, до ста миллионов и выше.
Для конкретной задачи число многосеточных итераций и время счета зависит от реализации основных элементов многосеточного метода — операторов межсеточных переходов, алгоритма решения уравнений на самой грубой сетке, а также и в большой степени от сглаживающих процедур. Современные суперкомпьютеры предъявляют к алгоритмическим реализациям дополнительное требование масштабируемости, что в содержательном смысле означает сохранение эффективности выполнения компьютерного кода в широком диапазоне числа используемых процессоров вплоть до экстрамассивного параллелизма. Такому требованию удовлетворяет многосеточный алгоритм при условии его реализации на принципах логической простоты. В [4]—[8] дано детальное численное исследование свойств такого многосеточного алгоритма, в том числе и при решении анизотропных задач с рассмотрением типичных случаев анизотропии по координатным направлениям. Данная статья является продолжением работ [4]—[8] в части рассмотрения анизотропных и разрывных коэффициентов, типичных для задач диффузии и теплопроводности в неоднородных областях.
Построенный многосеточный алгоритм демонстрирует параллельную эффективность при решении краевых задач для эллиптических уравнений. Отметим, что наш опыт реализации многосеточного метода на неструктурированных сетках при решении более сложных уравнений (см. [9], [10]) показывает, что добиться высокого уровня масштабируемости вне ограниченного круга задач и сеток чрезвычайно трудно.
1) Работа выполнена при финансовой поддержке Российского научного фонда (код проекта 14-21-00025).
1. ПОСТАНОВКА ЗАДАЧИ Рассмотрим краевую задачу для эллиптического уравнения
-Шу(^гаёи) + а0 и = /, г е О, -(^гаё и )• п = сти + у, г еГ.
Здесь г = (х, у, I) е О, где О — трехмерная область с границей Г, п — вектор внешней нормали к Г. Функции к(г), /(г), а0(г) > 0, ст(г) > 0, у(г) являются заданными, а функция и(г) — искомой. В общем случае к — тензор диффузии; здесь мы ограничимся случаем, когда к — скалярная положительная функция или диагональная 3 х 3-матрица с положительными диагональными элементами. Считаем, что входные данные обеспечивают разрешимость краевой задачи. В частности, для задачи Неймана (при а0 = 0, ст(г) = 0) выполняется необходимое условие разрешимости
|у( г)йз = /(г) йч, (2)
г п
тогда на подпространстве функций, ортогональных постоянной функции (в скалярном произведении пространства Х2(О)), задача Неймана имеет единственное решение.
Вопросы аппроксимации по пространству задачи (1) здесь не обсуждаются, поэтому для простоты возьмем в качестве области О ^-мерный прямой параллелепипед и в нем неравномерную по каждому координатному направлению декартову сетку Ок = {гп е О, 0 < п < N с сеточной границей Гк. Рассмотрим стандартную схему, построенную методом конечных объемов. Свяжем с каждым узлом сетки Ок дуальную ячейку в виде прямого параллелепипеда, грани которого проходят через середины ребер сетки Ок. Обозначим объемы дуальных ячеек через Кп, 0 < п < N.
Схему записываем на основе закона сохранения, интегрируя исходное уравнение по всем дуальным ячейкам, за исключением ячеек, у которых вершины лежат на границе Гк и в них записано краевое условие Дирихле.
Скалярное произведение в пространстве функций ик, заданных на сетке Ок (и обращающихся в нуль в узлах Гк с условием Дирихле), определяется стандартным образом:
(и, w) = ^UnWnVn, (3)
п
где сумма берется по всем узлам сетки. Сеточная Х2(О)-норма определяется равенством ||и||2 = (и, и). Напомним, что норма линейного оператора Б, действующего из ик в ик, определяется как ||Б|| = = шах||Бг|| при ||г|| = 1, и для самосопряженного оператора норма равна максимальному из модулей его собственных значений.
Дискретизацию краевой задачи в пространство сеточных функций ик будем записывать в операторном виде
Акик = . (4)
Оператор Лк — это разреженная N х ^-матрица, ик, gh — искомая и заданная сеточные функции. Каждая строка матрицыЛк соответствует узлу сетки с номером п и содержит не более семи ненулевых элементов. Правая часть gh схемы (4) в узле с номером п определяется как (як)п = (/к)п/Уп. При таком определении операторЛк является самосопряженным в скалярном произведении (3).
Декартову неравномерную сетку характеризуем параметром к, определяемым равенством к2 = = + Н22 + Н23, где ка — среднеквадратичный шаг неравномерной сетки по каждому из трех координатных направлений, а =1, 2, 3. Точность построенной разностной схемы имеет порядок 0(к2) в норме максимума модуля на последовательности сеток, для которой поверхности разрыва коэффициентов к, а0, / согласованы с сеткой (см. [11]). В практических расчетах мы не требуем такого согласования.
Объясним важный аспект, связанный с вырожденной задачей Неймана. В силу дискретизации, основанной на законе сохранения, аналогом условия разрешимости краевой задачи является соотношение
X УА = ^/пуп, (5)
к п
где суммирование в правой части ведется по всем узлам сетки, а в левой — по граничным узлам, внутренним для граней исходного параллелепипеда О, Бк — площади граней граничных дуальных ячеек; значения Бк берутся со знаком "+" или "—" в зависимости от ориентации внешней нормали к грани относительно декартовой системы координат. Действительно, при суммировании по всем дуальным ячейкам произойдет сокращение потоков на внутренних гранях этих ячеек. В итоге в левой части равенства останутся только суммы потоков на шести гранях области О.
Запишем для вырожденной задачи Неймана при произвольном приближении ий невязку г системы разностных уравнений с компонентами гп = (Д — Айий)п. Тогда условие разрешимости (5) можно записать в виде
г,У, = 0. (6)
В случае вырожденной задачи Неймана точность выполнения сеточного аналога (6) условия (2) на всех сеточных уровнях влияет на эффективность многосеточного алгоритма.
Исследование многосеточного метода опирается на спектральный анализ разностных операторов краевых задач. Здесь отметим, что собственные значения оператора Ай лежат на отрезке [^т1п; ^тах] вещественной оси; > 0, так как оператор Ай является неотрицательно определенным. Для вырожденной задачи Неймана ядро оператора является одномерным подпространством, состоящим из сеточных функций-констант; в этом случае ^т1п = 0, тогда под ^т1п подразумеваем оценку ^¡п > 0 нижней границы оператора на подпространстве, ортогональном ядру оператора. Оценка ^тах достаточно просто получается с помощью теоремы Гершгорина о кругах (см. [12]).
Введем оператор Бк = — Дй, где Дй — разностный оператор Лапласа (при одинаковых с оператором Ай краевых условиях). Считаем, что для операторов Ай и Б известны постоянные 0 < < С2 в неравенстве
С:(Вниь ик) < (Акик, ик) < С2(Окик, ик), (7)
означающем спектральную эквивалентность операторов Ай и Бй.
2. ОБЩАЯ СХЕМА МНОГОСЕТОЧНОГО МЕТОДА
Для решения системы сеточных уравнений (4) возьмем итерационный многосеточный метод в классической форме, используя несколько сеточных уровней (обычно пять—шесть). Каждая многосеточная итерация состоит в переходе с подробной сетки на сетку следующего уровня вплоть до самой грубой сетки и обратно. Удобно объяснить основные элементы метода в двухсе-точном представлении, когда сеток всего две — подробная й-сетка и грубая Н-сетка. Тогда итерирующий оператор многосеточного метода имеет вид
0 = Бр(I- РА^ВА„)Бр. (8)
Здесь Ай, Ан — операторы на подробной и грубой сетках; Р и Я — операторы межсеточных переходов — с грубой сетки на подробную (интерполяции Р) и с подробной на грубую сетку (проектирования Я); Бр — оператор сглаживания с числомр пред- и постсглаживающих шагов. В рассматриваемом алгоритме оператор Ан строится редискретизацией, т.е. на Н-сетке записывается аппроксимация исходного уравнения с однородными краевыми условиями. В качестве правой части системы уравнений Айд = на Н-сетке берется проекция невязки гй= gh — Айий на Н-сетку с помощью оператора Я. Разрешающий оператор многосеточного метода после выполнения
т итераций имеет вид С = (I — 0т)А/ ; с помощью этого оператора можно формально записать
,
полученное приближенное решение как ык = С/ь. Видим, что при 0т —»- 0 разрешающий оператор многосеточного метода приближает обратный кАк: С —» А^ . Подразумевается, что оператор С реализуется с помощью итерационного алгоритма
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.