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

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

Том 7, № 2, 2002

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

МАТЕМАТИКИ*

Ю. М. ЛАЕВСКИй Институт вычислительной математики и математической геофизики СО РАИ, Новосибирск, Россия e-mail: laev@labchem.sscc.ru

The paper deals with discussing some features of numerical methods for solution large problems of mathematical physics. The modern point of view to approximation problems, iterative methods, efficiency preconditions design, solution of nonstationary problems a presented.

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

Под большими задачами математической физики в настоящее время понимаются задачи, дискретные модели которых имеют порядка 106 и более степеней свободы. При этом речь идет о моделях, описываемых системами дифференциальных уравнений в частных производных. Алгоритмы решения интегральных уравнений здесь не рассматриваются, хотя, говоря о методах декомпозиции области, мы коснемся и их. Работа с задачами таких размерностей потребовала принципиально новых подходов, что, в конечном итоге, позволяет говорить о появлении алгоритмов нового поколения (термин в достаточной степени условный). Следует отметить, что для больших задач чрезвычайно важными становятся вопросы усвоения исходных данных и обработки результатов, организации информационных потоков в процессе вычислений, способы хранения информации и т. д. — все то, что принято называть вычислительной технологией. В данной работе указанные вопросы затрагиваться не будут. При их решении самую существенную роль играет архитектура ЭВМ, на которой проводится моделирование, но существуют и общие технологические принципы, которых следует придерживаться при разработке вычислительных программ.

* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований, грант №01-01-00819, и программы "Университеты России", грант ЗН 336-99.

© Ю.М. Лаевский, 2002.

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

1. Дискретизация и сопутствующие постановки задач

Основным инструментом пространственной дискретизации для уравнений математической физики сейчас следует считать метод конечных элементов (МКЭ). Связано это, с одной стороны, с переходом к решению больших реальных задач в областях, имеющих как правило, очень сложную геометрию и требующих использования существенно неравномерных неструктурированных сеток, и, с другой стороны, с возможностью достижения практически полной универсализации и автоматизации построения сеточных систем. При этом современный взгляд на МКЭ не исключает различные отклонения от строгого метода Галеркина: диагонализацию матрицы масс, использование аналогов направленных разностей при аппроксимации конвективных слагаемых, комбинацию с методом конечных объемов (МКО) и т. д.

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

(и,у)ы = / (и).

Указанная постановка не всегда очевидна. Например, для задачи Стокса с естественными краевыми условиями путем к такой постановке является выписывание уравнения движения в виде

— а(и,р) = ¥, где а(и,р) — симметрический тензор напряжений

а(и,р) = n(gгad и + (grad и)т) — р1;

I — единичный тензор; и — вектор скорости; р — давление. В случае традиционной записи уравнения движения в виде

—пДи + grad р = ¥

успех достигается только для главных краевых условий и = 0. Дело в том, что в традиционной записи мы "забыли" о слагаемом grad(divи), равном нулю в силу условия несжимаемости div и = 0. Отметим, что в линейной теории упругости такого затруднения не возникает именно из-за отсутствия условия несжимаемости, и уравнение движения имеет вид

—^Ди — (Л + ^)grad ^у и) = ¥.

Применение МКЭ, по существу, приводит к задаче о представлении линейного ограниченного функционала в конечномерном подпространстве У^ С Н с локальным финитным базисом (как правило, кусочно-полиномиальным). С точки зрения анализа погрешности МКЭ — это поиск ортогональной проекции решения исходной задачи на подпространство.

В рамках проекционных постановок большое внимание в последние годы уделяется смешанному МКЭ. Собственно смешанный МКЭ подразумевает новую исходную постановку, позволяющую снизить априорную гладкость искомого решения, а следовательно, и уменьшить размерность аппроксимирующего подпространства. Это хорошо видно на примере уравнений движения вязкой жидкости. Если, например, компоненты скорости в задаче Стокса с краевыми условиями и = 0 искать в пространстве соленоидальных функций (условие несжимаемости ввести в определение пространства разрешимости)

Н = {и е (я01(П))т, а1у и = 0},

то, с одной стороны, получим обобщенную задачу для вектора скорости, в которой давление исключено,

а(и, V) = п J grad и : grad = J f • vdП, V е Н, п п

а с другой — для давления получим уравнение Пуассона, т. е. давление ищется как элемент пространства Н 1(П)/еоп81. Аппроксимация по МКЭ для такой постановки приводит к весьма сложной конструкции, поскольку необходимо, по крайней мере приближенно, удовлетворить условию У^ С Н. Но задачу Стокса можно рассмотреть и в другой постановке: найти функции и е (Н1(П))т и р е ¿2(П)/сопв1 такие, что

а(и, V) + Ь(р, V) = У f • vdП, V е (Н1(П))т, п

Ь(д, и) = 0, д е Ь2(П)/сопв1,

где

Ь(д, и) = — J д div vdП.

п

В этой постановке соленоидальность поля скоростей уже не априорное требование, а результат решения задачи, а давление — суммируемая в квадрате функция. Аппроксимация по смешанному МКЭ осуществляется на паре пространств У^ С (Н0(П))т, С Ь2(П)/сопв1 и приводит к алгебраической задаче о поиске седловой точки

А В и ] = ( Р вт о ) V р ) V о

с симметричной, но не положительно определенной матрицей. При этом корректность приведенной системы является следствием условия Ладыженской — Бабушки-Бреци на конечно-элементные пространства У^ и Б^: существует не зависящее от Н число 7 > 0 такое, что

• г Ь(д, v) ^

ш! йиР ........ ^ 7-

дея По NV Н1

К совершенно аналогичной алгебраической системе приводит смешанный МКЭ при применении к уравнению Пуассона, бигармоническому уравнению и т. д. Единственное отличие от задачи Стокса состоит в необходимости предварительного получения смешанной постановки задачи введением векторной переменной.

И, наконец, весьма перспективным, на мой взгляд, для решения векторных задач (в том числе полученных из скалярных уравнений смешанной постановкой) выглядит векторный МКЭ. Суть метода заключается в использовании векторных базисных функций и привязке степеней свободы не к вершинам разбиения, а к ребрам или граням. При этом нетрудно построить базисы, для которых на каждом элементе либо div v = 0, либо rot v = 0. В первом случае оказываются непрерывными тангенциальные, а во втором — нормальные компоненты базисных векторов при переходе от одного элемента к другому. Проблематика, связанная с применением векторного МКЭ к различным классам задач, в настоящее время интенсивно развивается и содержит много открытых вопросов.

2. Методы решения СЛАУ

Рассмотрим только итерационные процедуры, предложенные и активно исследовавшиеся в последние 10-15 лет. Что же касается прямых методов, то в настоящее время на Западе сложилась точка зрения, что применение задач с более чем 105 степеней свободы слишком дорого, даже для современных суперЭВМ. Но для задач, полученных применением МКЭ, с нерегулярной структурой матрицы до порядков 104 — 5 • 104 (в зависимости от используемой ЭВМ) предпочтение следует отдать прямым методам.

Среди итерационных методов долгое время наиболее употребимыми были методы релаксационного типа: Якоби, Гаусса — Зейделя, SOR, SSOR и их блочные модификации. Теория сходимости этих методов основана на понятии регулярного расщепления матриц и хорошо развита. Одновременно с релаксационными появились методы, основанные на проектировании на последовательность подпространств Крылова, остававшиеся какое-то время на втором плане. Под подпространством Крылова степени n понимается линейная оболочка векторов из RN:

Kn(A; r) = {r, Ar,..., An-1r}

(A — квадратная матрица размерности N x N). Интересно отметить, что в современных обзорах методы, опирающиеся на информацию о спектре матрицы, вообще не обсуждаются. Для итерационного решения систем линейных алгебраических уравнений (СЛАУ)

Au = f

с использованием подпространств Крылова Kn(A; ro), где r0 = f — Auo, а u0 — начальное приближение, имеются два основных подхода:

1. Подход Ритца — Галеркина, когда очередное приближение uk ищется на основе ортогонализации невязки текущему подпространству Крылова:

f — Auk ± Kk(A; ro).

2. Минимизация невязки, когда очередное приближение uk ищется на основе минимизации эвклидовой нормы |f — AukЦ2 на текущем подпространстве Крылова Kk(A; r0).

Для симметричных матриц в рамках подхода 1 наиболее предпочтительными долгое время оставались метод Ланцоша и в случае положительно определенной матрицы его двучленный вариант — метод сопряженных градиентов (CG). К этой же группе относятся известные методы FOM или GENCG. Для решения несимметричных задач был предложен метод Bi-CG, использующий дополнительно пространства Крылова Kk(AT; s0). Однако этот метод весьма не надежен. Наиболее известные методы второй группы — это

MINR.ES и GMRES. Алгоритм MINRES пригоден для решения симметричных, но зна-конеопределенных систем и использует тот же прием построения базиса, что и СО. Метод GMRES использует метод вращений Гивенса при приведении матрицы Хессенберга к верхнетреугольному виду и является в настоящее время основным инструментом решения несимметричных систем. Широко применяется алгоритм GMRES(ш), осуществляющий перезапуск GMRES через каждые ш шагов. Существуют его различные модификации: FGMRES, GMRESR и т.д. И, наконец, широко известен применяемый к несимметричным системам метод Bi-CGSTAB, являющийся комбинацией методов Bi-CG и GMRES(1). По мнению Голуба и Ван дер Ворста, процесс формирования нового поколения алгоритмов линейной алгебры, основанных на пространствах Крылова, практически завершен, хотя остается ряд открытых вопросов, таких, например, как проблема останова итерационного процесса, разработка параллельных модификаций.

3. Проблемы переобусловливания

Хорошо известно, что скорость сходимости итерационных процессов зависит от спектральных свойств матрицы линейной системы. Дискретные аналоги уравнений в частных производных дают плохую обусловленность, но, к счастью, как правило, известную с точки зрения порядка СЛАУ. Вопрос о переобусловливании системы замечательно сформулировал Ланцош (1952): "Построение обратной матрицы эквивалентно линейному преобразованию, переводящему данную матрицу в единичную. Единичную матрицу можно вообразить как хорошо обусловленную матрицу со всеми собственными числами, равными единице". Собственно цель переобусловливания состоит в построении матрицы В (пе-реобусловливателя), в некотором смысле близкой к исходной матрице А (матрица АВ-1 близка к единичной), а операция В-1г реализуется экономично (либо прямым методом, либо некоторым внутренним итерационным процессом, скорость сходимости которого из-за специальной структуры В много выше, чем для матрицы А).

Среди методов алгебраического переобусловливания, когда матрица В определяется только элементами матрицы А и никак не использует информацию об исходной дифференциальной задаче, наиболее популярной является неполная Ш факторизация, осуществляющая приближенное Ш-разложение с занулением внедиагональных малых элементов верхней и нижней треугольных матриц. В случае симметричной положительно определенной матрицы А неполное разложение Холесского, используемое в качестве пе-реобусловливателя в методе CG, дает хорошие результаты для довольно больших задач (ICCG). Существуют различные варианты этого метода, например, метод с диагональной компенсацией.

В последние годы появились работы, в которых матрица М = В-1 (собственно именно эта матрица и нужна в процессе вычислений) ищется как разреженная матрица, делающая малой ¡I — АМ| для некоторой подходящей нормы.

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

о современной вычислительной математике

79

тальным понятием при построении переобусловливателя является спектральная эквивалентность матриц: матрицы А и В спектрально эквивалентны, если для любого вектора и выполняются неравенства

(здесь положительные числа а и в не зависят от размерности матриц). В этом случае собственные числа матрицы АВ-1 принадлежат отрезку [а, в] и задача становится хорошо обусловленной.

В настоящее время интенсивно изучаются вопросы построения эффективных переобус-ловливателей для задач с седловой точкой, в том числе несимметричных. Ответы на них имеют первостепенное значение для задач механики жидкости.

4. Методы декомпозиции и фиктивных областей

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

В своей центральной идее методы декомпозиции области (DDM) являются развитием и обобщением классического метода Шварца альтернирования по подобластям. Аддитивный вариант метода Шварца (ASM), являющийся основой построения и анализа DDM с налеганием подобластей для решения задачи о представлении линейного ограниченного функционала в гильбертовом пространстве H = Hi + • • • + Hm, можно сформулировать в виде следующего итерационного процесса:

При этом сходимость процесса со скоростью, не зависящей от размерности пространства Н, определяется соответствующей независимостью констант из следующих условий: для любого элемента V е Н существует представление

a (Bu, u) < (Au, u) < в (Bu, u)

u

k+i=uk - t (zk+•••+zm),

где Zi £ Hi и для любого v £ Hi

(Zi, v)Hi = (uk, v)h - f (v).

V = vi +----+ V;

m

такое, что

(vi, vi)h +-----+ (vm, vm)H < Y (v, v)h ,

и для любого w £ Hi имеют место неравенства

и матрицы А и В спектрально эквивалентны.

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

Н = Но + Нг, Но = Но,1 ф • • • ф Но,т

(Нг — пространство продолжений сеточных функций).

В векторно-матричных терминах СЛАУ сводится к уравнению для дополнения Шура:

Биз = (А33 - АТзА-11А1з - АТзА-1 А2з)из = дз-

Непрерывный аналог этого равенства является интегральным уравнением, а матрица Б — сеточным аналогом интегрального оператора Пуанкаре — Стеклова. Имеется многочисленная литература по построению переобусловливателей для этой системы.

Современная трактовка метода фиктивных областей наиболее полно отражена в подходе, который принято называть методом фиктивного пространства (ГБИ):

Но и Н — гильбертовы пространства со скалярными произведениями (ио,ио)о и (и, и),

А : Но ^ Но, В : Н ^ Н

— линейные самосопряженные положительно определенные операторы,

Я : Н ^ Но, Т : Но ^ Н

— линейные операторы такие, что оператор ЯТ : Но ^ Но тождественный и для любых ио е Но и и е Н

(АЯи, Яи)о < сД(Ви,и), ст(ВТио,Тио) < (Аио,ио)о, где сд и ст — положительные константы. Тогда

ст(А-1ио,ио)о < (ЯВ-1Я*ио,ио)о < сд(А-1ио,ио)о, где Я* : Но ^ Н — оператор, сопряженный к Я.

5. Многосеточные и многоуровневые методы

Многосеточные методы (МО) являются одними из самых эффективных средств решения уравнений в частных производных. В основе этих методов лежат вычисления на последовательности вложенных сеток

(Н0 С • • • С ,

порождающих последовательность вложенных конечномерных пространств

УНо с ••• с УНз

и на последовательности операторных уравнений

Аз из = ¡з, из е , з = 0,...,'1.

Классические многосеточные методы начали развиваться с начала 60-х годов. Суть их заключается в следующей рекуррентной процедуре:

— на j-м уровне проводится m итераций для j-й системы, задающих приближенное решение Vj и невязку Tj;

— невязка Tj £ Vhj проектируется на пространство Vhj_1: Tj-i = RjTj;

— на основе описываемой процедуры находится приближенное решение Wj— системы Aj-iwj-i = Tj-i;

— погрешность wj-1 интерполируется на j-й уровень: Wj = Ijwj-1;

— на j-м уровне полагается Vj = Vj — Wj, и, начиная с Vj как с начального приближения, проводятся l итераций для j-й системы (подавление ошибки интерполяции).

Выше был описан так называемый V-цикл. Широко распространена его модификация — W-цикл. Наиболее простая версия классических MG — каскадный метод, состоящий только в процедуре интерполяции с итерационным сглаживанием, получил теоретическое обоснование совсем недавно. В основе каскадного метода лежит "правильное" количество итераций на каждом уровне.

В конце 80-х годов начал формироваться современный взгляд на многосеточные методы — многоуровневые методы. Их основу составляет МКЭ на иерархических базисах (HB). Вводя стандартные в МКЭ интерполирующие операторы

П : Vhj ^ Vhj, j = 0,...,J,

используем представление

nj = Ro + ■ ■ ■ + Rj ,

где

Ro = По, Rj = П — П-i, j = 1, ••• ,J.

С учетом того, что в пространстве VJ оператор nj является тождественным, имеет место разложение

uj = vo +-----+ Vj,

где Vj = Rj uj, и применим ASM. При этом

Bo+ = RTA-iRo, B+ = RTD-iRj, j = 1, ■ ■ ■ , J,

где Dj = diag(Rj ARj) и переобусловливатель для решения системы на верхнем уровне есть

B-i = B+ + ••• + B+.

Для приведенного метода число обусловленности cond (B-i/2AB-i/2) зависит от hj. Этого недостатка лишен другой HB-метод — BPX, скорость сходимости которого не зависит не только от hj, но и от J. Для BPX операторы Rj определяются как

Ro = Qo, Rj = Qj — Qj-i, j = 1, ■ ■ ■ , J,

где Qj — Ь2-проекторы:

Qj : Vhj ^ Vhj, (Qju, v)L2 = (u, v)La Vv £ Vhj.

При этом

B+ = h2-dRT Rj, j = 1, ••• ,J,

где hj = 2-jh0; d — размерность исходной дифференциальной задачи. Отметим, что многоуровневые методы, основанные на ASM, являются аддитивными вариантами классических MG-методов и, следовательно, прекрасно распараллеливаются.

По мнению Саада, в настоящее время пользователи итерационной техники разделены на два лагеря — на сторонников исключительно MG-алгоритмов и сторонников чисто алгебраических методов в пространствах Крылова. Имеются компромиссные варианты (AMLI), но они используются достаточно редко. Методы в пространствах Крылова обладают значительной универсальностью. Но если речь идет о дискретном аналоге уравнения в частных производных, для которого имеется MG-техника (гладкое решение, регулярные сетки), то ее использование много эффективнее.

6. Методы решения нестационарных задач

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

В последнее десятилетие проводятся исследования, связанные с безытерационными DDM. Были построены и изучены так называемые регионально аддитивные методы, суть которых состоит в использовании схем расщепления не по координатным направлениям, а по подобластям. При этом были рассмотрены схемы как с налеганием, так и без налегания подобластей. И, наконец, совсем недавно был предложен безытерационный вариант ASM.

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

7. Задачи с особенностями

Коснемся задач, возникающих вследствие анизотропии, малых включений и особенностей, при наличии композитных материалов, пограничных слоев, многомасштабных пространственно-временных процессов. Численное решение таких задач основано на адаптивных стратегиях, включающих анализаторы апостериорной информации о погрешности, генерацию и адаптацию расчетных сеток, что позволяет использовать сетки, существенно неравномерные и неструктурированные. При этом плохая обусловленность возникающих СЛАУ связана не только с неравномерностью сетки, но и с собственно зависимостью от параметров, близких к своим предельным значениям. Для этого, в общем случае, требуются параметрические методы переобусловливания нового типа и использование эффективных

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

Выводы

1. Продолжаются интенсивные исследования по МКЭ с целью поиска конечномерных пространств, наиболее адекватно отвечающих исходной дифференциальной задаче. Упор делается на описание векторных полей и использование неструктурированных сеток.

2. В целом для решения СЛАУ имеется достаточно богатый набор алгоритмов, и основная задача состоит в разработке качественных программ с учетом распараллеливания, эффективных критериев останова и организации широкомасштабного тестирования на больших реальных задачах.

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

— декомпозиции области;

— фиктивного пространства;

— многоуровневой декомпозиции;

— полуаналитическим (с явным выделением особенностей, погранслоев, внутренних фронтов и т. д.).

4. При моделировании сложных нестационарных процессов следует ориентироваться на новые возможности явных схем и использовать различные явно-неявные методики.

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